/* 
    This file was generated by "mtxrun --script "mtx-wtoc.lua" from the metapost cweb files but
    now maintained as C file.
*/

/* 
    todo: check the random code, we might as well use the lua randomizer for all number models. 

    done: no longer intercept zero because a test is worse for branch prediction than a redundant action
*/

# include "mpmathdouble.h"

# define  PI                      3.1415926535897932384626433832795028841971

# define  fraction_multiplier     4096.0
# define  angle_multiplier        16.0
# define  coef_bound              ((7.0 / 3.0) * fraction_multiplier)
# define  fraction_threshold      0.04096
# define  half_fraction_threshold (fraction_threshold / 2)
# define  scaled_threshold        0.000122
# define  half_scaled_threshold   (scaled_threshold / 2)
# define  near_zero_angle         (0.0256*angle_multiplier)
# define  p_over_v_threshold      0x80000
# define  equation_threshold      0.001

# define  warning_limit           pow(2.0,52.0)
# define  epsilon                 pow(2.0,-52.0)

# define  unity                   1.0
# define  two                     2.0
# define  three                   3.0
# define  half_unit               0.5
# define  three_quarter_unit      0.75

# define  EL_GORDO                (DBL_MAX / 2.0 - 1.0)
# define  negative_EL_GORDO       (-EL_GORDO)
# define  one_third_EL_GORDO      (EL_GORDO / 3.0)

# define  fraction_half           (0.5 * fraction_multiplier)
# define  fraction_one            (1.0 * fraction_multiplier)
# define  fraction_two            (2.0 * fraction_multiplier)
# define  fraction_three          (3.0 * fraction_multiplier)
# define  fraction_four           (4.0 * fraction_multiplier)

# define  no_crossing             (fraction_one + 1)
# define  one_crossing            fraction_one
# define  zero_crossing           0

# define  one_eighty_deg          ( 180.0 * angle_multiplier)
# define  negative_one_eighty_deg (-180.0 * angle_multiplier)
# define  three_sixty_deg         ( 360.0 * angle_multiplier)

# define  two_to_the(A)           (1 << (unsigned)(A))
# define  set_cur_cmd(A)          mp->cur_mod_->command = (A)
# define  set_cur_mod(A)          mp->cur_mod_->data.n.data.dval = (A)

/*tex 

Apart from it looking weird in results a |-0.0| can also give wrong results, for instance in a 
|tan2|, see there for a comment. This is why we now do some checking on zero which in some cases
also might be a bit more efficient as it avoids a multiplication or division, so likely we break 
even. We can use |signbit| but equally well just set or add zero.

*/

static inline double mp_double_make_fraction(double p, double q)
{
    return (p / q) * fraction_multiplier;
}

static inline double mp_double_take_fraction(double p, double q)
{
    return (p * q) / fraction_multiplier;
}

static inline double mp_double_make_scaled(double p, double q)
{
    return p / q;
}

static inline void mp_double_allocate_number(MP mp, mp_number *n, mp_number_type t)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = 0.0,
        .type      = t,
    };
}

static inline void mp_double_allocate_clone(MP mp, mp_number *n, mp_number_type t, const mp_number *v)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = v->data.dval,
        .type      = t,
    };
}

static inline void mp_double_allocate_abs(MP mp, mp_number *n, mp_number_type t, const mp_number *v)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = fabs(v->data.dval),
        .type      = t,
    };
}

static inline void mp_double_allocate_div(MP mp, mp_number *n, mp_number_type t, const mp_number *a, const mp_number *b)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = a->data.dval == 0.0 ? 0.0 : a->data.dval / b->data.dval,
        .type      = t,
    };
}

static inline void mp_double_allocate_mul(MP mp, mp_number *n, mp_number_type t, const mp_number *a, const mp_number *b)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = a->data.dval == 0.0 || b->data.dval == 0.0 ? 0.0 : a->data.dval * b->data.dval,
        .type      = t,
    };
}

static inline void mp_double_allocate_add(MP mp, mp_number *n, mp_number_type t, const mp_number *a, const mp_number *b)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = a->data.dval + b->data.dval,
        .type      = t,
    };
}

static inline void mp_double_allocate_sub(MP mp, mp_number *n, mp_number_type t, const mp_number *a, const mp_number *b)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = a->data.dval - b->data.dval,
        .type      = t,
    };
}

static inline void mp_double_allocate_double(MP mp, mp_number *n, double v)
{
    (void) mp;
    *n = (mp_number) {
        .data.dval = v,
        .type      = mp_scaled_type,
    };
}

static inline void mp_double_free_number(MP mp, mp_number *n)
{
    (void) mp;
    n->type = mp_nan_type;
}

static inline void mp_double_set_from_int(mp_number *A, mp_scaled_t B)
{
    A->data.dval = B;
}

static inline void mp_double_set_from_boolean(mp_number *A, mp_scaled_t B)
{
    A->data.dval = B;
}

static inline void mp_double_set_from_scaled(mp_number *A, mp_scaled_t B)
{
    static const double INV_65536 = 1.0 / 65536.0;
    A->data.dval = B * INV_65536;
}

static inline void mp_double_set_from_double(mp_number *A, double B)
{
    A->data.dval = B;
}

static inline void mp_double_set_from_addition(mp_number *A, const mp_number *B, const mp_number *C)
{
    A->data.dval = B->data.dval + C->data.dval;
}

static inline void mp_double_set_half_from_addition(mp_number *A, const mp_number *B, const mp_number *C)
{
    A->data.dval = (B->data.dval + C->data.dval) / 2.0;
}

static inline void mp_double_set_from_subtraction(mp_number *A, const mp_number *B, const mp_number *C)
{
    A->data.dval = B->data.dval - C->data.dval;
}

static inline void mp_double_set_half_from_subtraction(mp_number *A, const mp_number *B, const mp_number *C)
{
    A->data.dval = (B->data.dval - C->data.dval) / 2.0;
}

static inline void mp_double_set_from_div(mp_number *A, const mp_number *B, const mp_number *C)
{
    A->data.dval = B->data.dval / C->data.dval;
}

static inline void mp_double_set_from_mul(mp_number *A, const mp_number *B, const mp_number *C)
{
    A->data.dval = B->data.dval * C->data.dval;
}

static inline void mp_double_set_from_int_div(mp_number *A, const mp_number *B, mp_scaled_t C)
{
    A->data.dval = B->data.dval == 0.0 ? 0.0 : B->data.dval / C;
}

static inline void mp_double_set_from_int_mul(mp_number *A, const mp_number *B, mp_scaled_t C)
{
    A->data.dval = B->data.dval * (double) C;
}

static inline void mp_double_set_from_of_the_way(MP mp, mp_number *A, const mp_number *t, const mp_number *B, const mp_number *C)
{
    (void) mp;
    A->data.dval = B->data.dval - mp_double_take_fraction(B->data.dval - C->data.dval, t->data.dval);
}

static inline void mp_double_negate(mp_number *A)
{
    /* the addition gets rid of the sign bit when we have zero: no -ffast-math */
    A->data.dval = -A->data.dval + 0.0;
}

static inline void mp_double_add(mp_number *A, const mp_number *B)
{
    A->data.dval += B->data.dval;
}

static inline void mp_double_subtract(mp_number *A, const mp_number *B)
{
    A->data.dval -= B->data.dval;
}

static inline void mp_double_half(mp_number *A)
{
    A->data.dval /= 2.0;
}

static inline void mp_double_double(mp_number *A)
{
    A->data.dval *= 2.0;
}

static inline void mp_double_add_scaled(mp_number *A, mp_scaled_t B)
{
    /* also for negative B */
    static const double INV_65536 = 1.0 / 65536.0;
    A->data.dval += (B * INV_65536);
}

static inline void mp_double_multiply_int(mp_number *A, mp_scaled_t B)
{
    A->data.dval *= (double) B;
}

static inline void mp_double_divide_int(mp_number *A, mp_scaled_t B)
{
    A->data.dval /= (double) B;
}

static inline void mp_double_abs(mp_number *A)
{
    A->data.dval = fabs(A->data.dval);
}

static inline void mp_double_clone(mp_number *A, const mp_number *B)
{
    *A = *B;
}

static inline void mp_double_negated_clone(mp_number *A, const mp_number *B)
{
    /* the addition gets rid of the sign bit when we have zero: no -ffast-math */
    A->data.dval = -B->data.dval + 0.0;
}

static inline void mp_double_abs_clone(mp_number *A, const mp_number *B)
{
    A->data.dval = fabs(B->data.dval);
}

static inline void mp_double_swap(mp_number *A, mp_number *B)
{
    mp_number tmp = *A;
    *A = *B;
    *B = tmp;
}

static inline void mp_double_fraction_to_scaled(mp_number *A)
{
    A->type = mp_scaled_type;
    A->data.dval /= fraction_multiplier;
}

static inline void mp_double_angle_to_scaled(mp_number *A)
{
    A->type = mp_scaled_type;
    A->data.dval /= angle_multiplier;
}

static inline void mp_double_scaled_to_fraction(mp_number *A)
{
    A->type = mp_fraction_type;
    A->data.dval *= fraction_multiplier;
}

static inline void mp_double_scaled_to_angle(mp_number *A)
{
    A->type = mp_angle_type;
    A->data.dval *= angle_multiplier;
}

static mp_scaled_t mp_double_to_scaled(const mp_number *A)
{
    return (mp_scaled_t) mpscaledround(A->data.dval * 65536.0);
}

static mp_scaled_t mp_double_to_int(const mp_number *A)
{
    return (mp_scaled_t) (A->data.dval);
}

static mp_scaled_t mp_double_to_boolean(const mp_number *A)
{
    return (mp_scaled_t) (A->data.dval);
}

static double mp_double_to_double(const mp_number *A)
{
    return A->data.dval;
}

static int mp_double_odd(const mp_number *A)
{
    return odd_long(lround(A->data.dval));
}

static int mp_double_equal(const mp_number *A, const mp_number *B)
{
    return A->data.dval == B->data.dval;
}

static int mp_double_greater(const mp_number *A, const mp_number *B)
{
    return A->data.dval > B->data.dval;
}

static int mp_double_less(const mp_number *A, const mp_number *B)
{
    return A->data.dval < B->data.dval;
}

static int mp_double_non_equal_abs(const mp_number *A, const mp_number *B)
{
    return fabs(A->data.dval) != fabs(B->data.dval);
}

static char *mp_double_number_tostring(MP mp, const mp_number *n)
{
    static char set[64];
    int l = 0;
    char *ret = mp_memory_allocate(64);
    (void) mp;
    snprintf(set, 64, mp->less_digits ? "%.3g" : "%.17g", n->data.dval);
    while (set[l] == ' ') {
        l++;
    }
    strcpy(ret, set+l);
    return ret;
}

static inline void mp_double_print_number(MP mp, const mp_number *n)
{
    char *str = mp_double_number_tostring(mp, n);
    mp_print_e_str(mp, str);
    mp_memory_free(str);
}

static void mp_double_slow_add(MP mp, mp_number *ret, const mp_number *x_orig, const mp_number *y_orig)
{
    double x = x_orig->data.dval;
    double y = y_orig->data.dval;
    if (x >= 0.0) {
        if (y <= EL_GORDO - x) {
            ret->data.dval = x + y;
        } else {
            mp->arithmic_error = mp_error_code(mp, 1);
            ret->data.dval = EL_GORDO;
        }
    } else if (-y <= EL_GORDO + x) {
        ret->data.dval = x + y;
    } else {
        mp->arithmic_error = mp_error_code(mp, 2);
        ret->data.dval = negative_EL_GORDO;
    }
}

static void mp_double_slow_sub(MP mp, mp_number *ret, const mp_number *x_orig, const mp_number *y_orig)
{
    double x = x_orig->data.dval;
    double y = y_orig->data.dval;
    if (x >= 0.0) {
        if (-y <= EL_GORDO - x) {
            ret->data.dval = x - y;
        } else {
            mp->arithmic_error = mp_error_code(mp, 1);
            ret->data.dval = EL_GORDO;
        }
    } else if (y <= EL_GORDO + x) {
        ret->data.dval = x - y;
    } else {
        mp->arithmic_error = mp_error_code(mp, 2);
        ret->data.dval = negative_EL_GORDO;
    }
}

static inline void mp_double_number_make_fraction(MP mp, mp_number *ret, const mp_number *p, const mp_number *q)
{
    (void) mp;
    ret->data.dval = mp_double_make_fraction(p->data.dval, q->data.dval);
}

static inline void mp_double_number_take_fraction(MP mp, mp_number *ret, const mp_number *p, const mp_number *q)
{
   (void) mp;
   ret->data.dval = mp_double_take_fraction(p->data.dval, q->data.dval);
}

static inline void mp_double_number_take_scaled(MP mp, mp_number *ret, const mp_number *p_orig, const mp_number *q_orig)
{
    (void) mp;
    ret->data.dval = p_orig->data.dval * q_orig->data.dval;
}

static inline void mp_double_number_make_scaled(MP mp, mp_number *ret, const mp_number *p_orig, const mp_number *q_orig)
{
    (void) mp;
    ret->data.dval = p_orig->data.dval / q_orig->data.dval;
}

static void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop)
{
    double result;
    char *end = (char *) stop;
    errno = 0;
    result = strtod((char *) start, &end);
    if (errno == 0) {
        set_cur_mod(result);
        if (result >= warning_limit) {
            if (internal_value(mp_warning_check_internal).data.dval > 0 && (mp->scanner_status != mp_tex_flushing_state)) {
                char msg[256];
                snprintf(msg, 256, "Number is too large (%g)", result);
                mp_error(
                    mp,
                    msg,
                    "Continue and I'll try to cope with that big value; but it might be dangerous."
                    "(Set warningcheck := 0 to suppress this message.)"
                );
            }
        }
    } else if (mp->scanner_status != mp_tex_flushing_state) {
        mp_error(
            mp,
            "Enormous number has been reduced.",
            "I could not handle this number specification probably because it is out of"
            "range."
        );
        set_cur_mod(EL_GORDO);
    }
    set_cur_cmd(mp_numeric_command);
}

static void mp_double_aux_find_exponent(MP mp)
{
    if (mp->buffer[mp->cur_input.loc_field] == 'e' || mp->buffer[mp->cur_input.loc_field] == 'E') {
        mp->cur_input.loc_field++;
        if (!(mp->buffer[mp->cur_input.loc_field] == '+'
           || mp->buffer[mp->cur_input.loc_field] == '-'
           || mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class)) {
            mp->cur_input.loc_field--;
            return;
        }
        if (mp->buffer[mp->cur_input.loc_field] == '+'
         || mp->buffer[mp->cur_input.loc_field] == '-') {
            mp->cur_input.loc_field++;
        }
        while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
            mp->cur_input.loc_field++;
        }
    }
}

static void mp_double_scan_fractional_token(MP mp, mp_scaled_t n) /* n is scaled */
{
    unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
    unsigned char *stop;
    (void) n;
    while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
        mp->cur_input.loc_field++;
    }
    mp_double_aux_find_exponent(mp);
    stop = &mp->buffer[mp->cur_input.loc_field-1];
    mp_wrapup_numeric_token(mp, start, stop);
}

/*tex
    The input format is the same as for the C language, so we just collect valid bytes in the 
    buffer, then call |strtod()|. It looks like we have no buffer overflow check here. 
*/

static void mp_double_scan_numeric_token(MP mp, mp_scaled_t n)
{
    unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
    unsigned char *stop;
    (void) n;
    while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
        mp->cur_input.loc_field++;
    }
    if (mp->buffer[mp->cur_input.loc_field] == '.' && mp->buffer[mp->cur_input.loc_field+1] != '.') {
        mp->cur_input.loc_field++;
        while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
            mp->cur_input.loc_field++;
        }
    }
    mp_double_aux_find_exponent(mp);
    stop = &mp->buffer[mp->cur_input.loc_field-1];
    mp_wrapup_numeric_token(mp, start, stop);
}

static void mp_double_velocity(MP mp, mp_number *ret, const mp_number *st, const mp_number *ct, const mp_number *sf, const mp_number *cf, const mp_number *t)
{
    double acc, num, denom; /* registers for intermediate calculations */
    (void) mp;
    acc = mp_double_take_fraction(st->data.dval - (sf->data.dval / 16.0), sf->data.dval - (st->data.dval / 16.0));
    acc = mp_double_take_fraction(acc, ct->data.dval - cf->data.dval);
    num = fraction_two + mp_double_take_fraction(acc, sqrt(2)*fraction_one);
    denom = fraction_three
          + mp_double_take_fraction(ct->data.dval, 3*fraction_half*(sqrt(5.0)-1.0))
          + mp_double_take_fraction(cf->data.dval, 3*fraction_half*(3.0-sqrt(5.0)));
    if (t->data.dval != unity) {
        num = mp_double_make_scaled(num, t->data.dval);
    }
    if (num / 4 >= denom) {
        ret->data.dval = fraction_four;
    } else {
        ret->data.dval = mp_double_make_fraction(num, denom);
    }
}

static inline int mp_double_ab_vs_cd(const mp_number *a_orig, const mp_number *b_orig, const mp_number *c_orig, const mp_number *d_orig)
{
    double ab = a_orig->data.dval * b_orig->data.dval;
    double cd = c_orig->data.dval * d_orig->data.dval;
    if (ab > cd) {
        return 1;
    } else if (ab < cd) {
        return -1;
    } else {
        return 0;
    }
}

static void mp_double_crossing_point(MP mp, mp_number *ret, const mp_number *aa, const mp_number *bb, const mp_number *cc)
{
    double a = aa->data.dval;
    double b = bb->data.dval;
    double c = cc->data.dval;
    (void) mp;
    if (a < 0.0) {
        ret->data.dval = zero_crossing;
        return;
    } else if (c >= 0.0) {
        if (b >= 0.0) {
            if (c > 0.0) {
                ret->data.dval = no_crossing;
            } else if ((a == 0.0) && (b == 0.0)) {
                ret->data.dval = no_crossing;
            } else {
                ret->data.dval = one_crossing;
            }
            return;
        } else if (a == 0.0) {
            ret->data.dval = zero_crossing;
            return;
        }
    } else if ((a == 0.0) && (b <= 0.0)) {
        ret->data.dval = zero_crossing;
        return;
    }
    /* Use bisection to find the crossing point... */
    {
        double d = epsilon;
        double x0 = a;
        double x1 = a - b;
        double x2 = b - c;
        do {
            /* not sure why the error correction has to be >= 1E-12 */
            double x = (x1 + x2) / 2 + 1E-12;
            if (x1 - x0 > x0) {
                x2 = x;
                x0 += x0;
                d += d;
            } else {
                double xx = x1 + x - x0;
                if (xx > x0) {
                    x2 = x;
                    x0 += x0;
                    d += d;
                } else {
                    x0 = x0 - xx;
                    if ((x <= x0) && (x + x2 <= x0)) {
                        ret->data.dval = no_crossing;
                        return;
                    }
                    x1 = x;
                    d = d + d + epsilon;
                }
            }
        } while (d < fraction_one);
        ret->data.dval = (d - fraction_one);
    }
}

/*tex 

    This is the same one as in |mpmath.c|:

*/

static inline mp_scaled_t mp_double_unscaled(const mp_number *x_orig)
{
    return (mp_scaled_t) mpscaledround(x_orig->data.dval);
}

static inline void mp_double_floor(mp_number *i)
{
    i->data.dval = floor(i->data.dval);
}

static inline void mp_double_fraction_to_round_scaled(mp_number *x_orig)
{
    double x = x_orig->data.dval;
    x_orig->type = mp_scaled_type;
    x_orig->data.dval = x / fraction_multiplier;
}

static void mp_double_square_rt(MP mp, mp_number *ret, const mp_number *x_orig) /* return, x: scaled */
{
    double x = x_orig->data.dval;
    if (x > 0) {
        ret->data.dval = sqrt(x);
    } else {
        if (x < 0) {
            char msg[256];
            char *xstr = mp_double_number_tostring(mp, x_orig);
            snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr);
            mp_memory_free(xstr);
            mp_error(
                mp,
                msg,
                "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
                "Proceed, with fingers crossed."
            );
        }
        ret->data.dval = 0;
    }
}

static void mp_double_pyth_add(MP mp, mp_number *ret, const mp_number *a_orig, const mp_number *b_orig)
{
    double a = fabs(a_orig->data.dval);
    double b = fabs(b_orig->data.dval);
    errno = 0;
    ret->data.dval = hypot(a, b);
    if (errno) {
        mp->arithmic_error = mp_error_code(mp, 3);
        ret->data.dval = EL_GORDO;
    }
}

static void mp_double_pyth_add3(MP mp, mp_number *ret, const mp_number *a_orig, const mp_number *b_orig, const mp_number *c_orig)
{
    double a = fabs(a_orig->data.dval);
    double b = fabs(b_orig->data.dval);
    double c = fabs(c_orig->data.dval);
    errno = 0;
    ret->data.dval = sqrt(a*a + b*b + c*c);
    if (errno) {
        mp->arithmic_error = mp_error_code(mp, 3);
        ret->data.dval = EL_GORDO;
    }
}

static void mp_double_pyth_sub(MP mp, mp_number *ret, mp_number *a_orig, const mp_number *b_orig)
{
    double a = fabs(a_orig->data.dval);
    double b = fabs(b_orig->data.dval);
    if (a > b) {
        a = sqrt(a*a - b*b);
    } else {
        if (a < b) {
            char msg[256];
            char *astr = mp_double_number_tostring(mp, a_orig);
            char *bstr = mp_double_number_tostring(mp, b_orig);
            snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, bstr);
            mp_memory_free(astr);
            mp_memory_free(bstr);
            mp_error(
                mp,
                msg,
                "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
                "Proceed, with fingers crossed."
            );
        }
        a = 0;
    }
    ret->data.dval = a;
}

static void mp_double_power_of(MP mp, mp_number *ret, const mp_number *a_orig, const mp_number *b_orig)
{
    errno = 0;
    ret->data.dval = pow(a_orig->data.dval, b_orig->data.dval);
    if (errno) {
        mp->arithmic_error = mp_error_code(mp, 4);
        ret->data.dval = EL_GORDO;
    }
}

static void mp_double_m_log(MP mp, mp_number *ret, const mp_number *x_orig)
{
    if (x_orig->data.dval > 0) {
        ret->data.dval = log(x_orig->data.dval)*256.0;
    } else {
        char msg[256];
        char *xstr = mp_double_number_tostring(mp, x_orig);
        snprintf(msg, 256, "Logarithm of %s has been replaced by 0", xstr);
        mp_memory_free(xstr);
        mp_error(
            mp,
            msg,
            "Since I don't take logs of non-positive numbers, I'm zeroing this one.\n"
            "Proceed, with fingers crossed."
        );
        ret->data.dval = 0;
    }
}

static void mp_double_m_exp(MP mp, mp_number *ret, const mp_number *x_orig)
{
    errno = 0;
    ret->data.dval = exp(x_orig->data.dval / 256.0);
    if (errno) {
        if (x_orig->data.dval > 0) {
            mp->arithmic_error = mp_error_code(mp, 5);
            ret->data.dval = EL_GORDO;
        } else {
            ret->data.dval = 0;
        }
    }
}

/*tex

    We (MS/HH) ran into an issue @ 2024-03-09 that showed up because a file from early 2023 gave 
    different results: |draw (50,0) {dir 120} .. (-50,0)  ..{dir 60} (50,0);| and in the end the 
    issue was due to a |-0| not making |atan2| happy. It was ok in scaled (what was the default
    in \CONTEXT\ before 2024), posit and decimal, but failed in double (which became default in 
    mid 2023) and in binary (both also tested in \LUATEX). 

    We now try to catch the |-0.0| upstream but we keep the patch commented. However, it was 
    decided that in mplib (in \LUAMETATEX) we patch the |n_arg| function, so there the commented 
    code below is used. 

*/

static void mp_double_n_arg(MP mp, mp_number *ret, const mp_number *x_orig, const mp_number *y_orig)
{
    if (x_orig->data.dval == 0.0 && y_orig->data.dval == 0.0) {
        if (internal_value(mp_default_zero_angle_internal).data.dval < 0) {
            mp_error(
                mp,
                "angle(0,0) is taken as zero",
                "The 'angle' between two identical points is undefined. I'm zeroing this one.\n"
                "Proceed, with fingers crossed."
            );
            ret->data.dval = internal_value(mp_default_zero_angle_internal).data.dval;
        } else { 
            ret->data.dval = 0;
        }
    } else {
        ret->type = mp_angle_type;
        ret->data.dval = atan2(y_orig->data.dval, x_orig->data.dval) * (180.0 / PI) * angle_multiplier;
        /* we need to get rid of the sign bit */
        if (ret->data.dval == 0.0) {
            ret->data.dval = 0.0;
        }
    }
}

static void mp_double_sin_cos(MP mp, const mp_number *z_orig, mp_number *n_cos, mp_number *n_sin)
{
    double rad = (z_orig->data.dval / angle_multiplier); /* still degrees */
    (void) mp;
    if ((rad == 90.0) || (rad == -270)){
        n_cos->data.dval = 0.0;
        n_sin->data.dval = fraction_multiplier;
    } else if ((rad == -90.0) || (rad == 270.0)) {
        n_cos->data.dval = 0.0;
        n_sin->data.dval = -fraction_multiplier;
    } else if ((rad == 180.0) || (rad == -180.0)) {
        n_cos->data.dval = -fraction_multiplier;
        n_sin->data.dval = 0.0;
    } else {
        rad = rad * PI / 180.0;
        n_cos->data.dval = cos(rad) * fraction_multiplier;
        n_sin->data.dval = sin(rad) * fraction_multiplier;
    }
}

/*tex
    This is the |http://www-cs-faculty.stanford.edu/~uno/programs/rng.c| with small cosmetic 
    modifications. The code only documented here as the other non scaled number models use the 
    same method.
*/

static void mp_double_aux_ran_array(MP mp, long aa[], int n)
{
    int i, j;
    for (j = 0; j < mp_random_KK; j++) {
        aa[j] = mp->random_data.x[j];
    }
    for (; j < n; j++) {
        aa[j] = mp_random_mod_diff(aa[j - mp_random_KK], aa[j - mp_random_LL]);
    }
    for (i = 0; i < mp_random_LL; i++, j++) {
        mp->random_data.x[i] = mp_random_mod_diff(aa[j - mp_random_KK], aa[j - mp_random_LL]);
    }
    for (; i < mp_random_KK; i++, j++) {
        mp->random_data.x[i] = mp_random_mod_diff(aa[j - mp_random_KK], mp->random_data.x[i - mp_random_LL]);
    }
}

static void mp_double_aux_ran_start(MP mp, long seed)
{
    int t, j;
    long x[mp_random_KK + mp_random_KK - 1]; /* the preparation buffer */
    long ss = (seed + 2) & (mp_random_MM - 2);
    for (j = 0; j < mp_random_KK; j++) {
        /* bootstrap the buffer */
        x[j] = ss;
        /* cyclic shift 29 bits */
        ss <<= 1;
        if (ss >= mp_random_MM) {
            ss -= mp_random_MM - 2;
        }
    }
    /* make x[1] (and only x[1]) odd */
    x[1]++;
    for (ss = seed & (mp_random_MM - 1), t = mp_random_TT - 1; t;) {
        for (j = mp_random_KK - 1; j > 0; j--) {
            /* "square" */
            x[j + j] = x[j];
            x[j + j - 1] = 0;
        }
        for (j = mp_random_KK + mp_random_KK - 2; j >= mp_random_KK; j--) {
            x[j - (mp_random_KK - mp_random_LL)] = mp_random_mod_diff(x[j - (mp_random_KK - mp_random_LL)], x[j]);
            x[j - mp_random_KK] = mp_random_mod_diff(x[j - mp_random_KK], x[j]);
        }
        if (odd_long(ss)) {
            /* "multiply by z" */
            for (j = mp_random_KK; j > 0; j--) {
                x[j] = x[j-1];
            }
            x[0] = x[mp_random_KK];
            /* shift the buffer cyclically */
            x[mp_random_LL] = mp_random_mod_diff(x[mp_random_LL], x[mp_random_KK]);
        }
        if (ss) {
            ss >>= 1;
        } else {
            t--;
        }
    }
    for (j = 0; j < mp_random_LL; j++) {
        mp->random_data.x[j + mp_random_KK - mp_random_LL] = x[j];
    }
    for (;j < mp_random_KK; j++) {
        mp->random_data.x[j - mp_random_LL] = x[j];
    }
    for (j = 0; j < 10; j++) {
        /* warm things up */
        mp_double_aux_ran_array(mp, x, mp_random_KK + mp_random_KK - 1);
    }
    mp->random_data.ptr = &mp->random_data.started;
}

static long mp_double_aux_ran_arr_cycle(MP mp)
{
    if (mp->random_data.ptr == &mp->random_data.dummy) {
        /* the user forgot to initialize */
        mp_double_aux_ran_start(mp, 314159L);
    }
    mp_double_aux_ran_array(mp, mp->random_data.buf, mp_random_QUALITY);
    mp->random_data.buf[mp_random_KK] = -1;
    mp->random_data.ptr = mp->random_data.buf + 1;
    return mp->random_data.buf[0];
}

static void mp_double_init_randoms(MP mp, int seed)
{
    mp->random_data = (mp_random_info) {
        .dummy   = -1,
        .started = -1,
        .ptr     = &mp->random_data.dummy
    };
    int k = 1;
    int j = abs(seed);
    int f = (int) fraction_one; /* avoid warnings */
    while (j >= f) {
        j = j / 2;
    }
    for (int i = 0; i <= 54; i++) {
        int jj = k;
        k = j - k;
        j = jj;
        if (k < 0) {
            k += f;
        }
        mp->randoms[(i * 21) % 55].data.dval = j;
    }
    mp_new_randoms(mp);
    mp_new_randoms(mp);
    mp_new_randoms(mp);
    /* warm up the array */
    mp_double_aux_ran_start(mp, (long) seed);
}

static inline void mp_double_modulo(mp_number *a, const mp_number *b)
{
    double tmp;
    a->data.dval = modf((double) a->data.dval / (double) b->data.dval, &tmp) * (double) b->data.dval;
}

static void mp_double_aux_next_unif_random(MP mp, mp_number *ret)
{
    unsigned long int op = (unsigned) (*mp->random_data.ptr >= 0 ? *mp->random_data.ptr++ : mp_double_aux_ran_arr_cycle(mp));
    double a = op / (mp_random_MM * 1.0);
    (void) mp;
    ret->data.dval = a;
}

static void mp_double_aux_next_random(MP mp, mp_number *ret)
{
    if (mp->j_random == 0) {
        mp_new_randoms(mp);
    } else {
        mp->j_random = mp->j_random - 1;
    }
    mp_double_clone(ret, &(mp->randoms[mp->j_random]));
}

static void mp_double_m_unif_rand(MP mp, mp_number *ret, const mp_number *x_orig)
{
    mp_number x, abs_x, u, y; /* |y| is trial value */
    mp_double_allocate_number(mp, &y, mp_fraction_type);
    mp_double_allocate_clone(mp, &x, mp_scaled_type, x_orig);
    mp_double_allocate_abs(mp, &abs_x, mp_scaled_type, &x);
    mp_double_allocate_number(mp, &u, mp_scaled_type);
    mp_double_aux_next_unif_random(mp, &u);
    y.data.dval = abs_x.data.dval * u.data.dval;
    mp_double_free_number(mp, &u);
    if (mp_double_equal(&y, &abs_x)) {
        mp_double_clone(ret, &((math_data *)mp->math)->md_zero_t);
    } else if (mp_double_greater(&x, &((math_data *)mp->math)->md_zero_t)) {
        mp_double_clone(ret, &y);
    } else {
        mp_double_negated_clone(ret, &y);
    }
    mp_double_free_number(mp, &abs_x);
    mp_double_free_number(mp, &x);
    mp_double_free_number(mp, &y);
}

static void mp_double_m_norm_rand(MP mp, mp_number *ret)
{
    mp_number abs_x, u, r, la, xa;
    mp_double_allocate_number(mp, &la, mp_scaled_type);
    mp_double_allocate_number(mp, &xa, mp_scaled_type);
    mp_double_allocate_number(mp, &abs_x, mp_scaled_type);
    mp_double_allocate_number(mp, &u, mp_scaled_type);
    mp_double_allocate_number(mp, &r, mp_scaled_type);
    do {
        do {
            mp_number v;
            mp_double_allocate_number(mp, &v, mp_scaled_type);
            mp_double_aux_next_random(mp, &v);
            mp_double_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t);
            mp_double_number_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v);
            mp_double_free_number(mp, &v);
            mp_double_aux_next_random(mp, &u);
            mp_double_clone(&abs_x, &xa);
            mp_double_abs(&abs_x);
        } while (! mp_double_less(&abs_x, &u));
        mp_double_number_make_fraction(mp, &r, &xa, &u);
        mp_double_clone(&xa, &r);
        mp_double_m_log(mp, &la, &u);
        mp_double_set_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la);
    } while (mp_double_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0);
    mp_double_clone(ret, &xa);
    mp_double_free_number(mp, &r);
    mp_double_free_number(mp, &abs_x);
    mp_double_free_number(mp, &la);
    mp_double_free_number(mp, &xa);
    mp_double_free_number(mp, &u);
}

static void mp_double_set_precision(MP mp)
{
    (void) mp;
}

static void mp_double_free_math(MP mp)
{
    mp_double_free_number(mp, &(mp->math->md_three_sixty_deg_t));
    mp_double_free_number(mp, &(mp->math->md_one_eighty_deg_t));
    mp_double_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t));
    mp_double_free_number(mp, &(mp->math->md_fraction_one_t));
    mp_double_free_number(mp, &(mp->math->md_zero_t));
    mp_double_free_number(mp, &(mp->math->md_half_unit_t));
    mp_double_free_number(mp, &(mp->math->md_three_quarter_unit_t));
    mp_double_free_number(mp, &(mp->math->md_unity_t));
    mp_double_free_number(mp, &(mp->math->md_two_t));
    mp_double_free_number(mp, &(mp->math->md_three_t));
    mp_double_free_number(mp, &(mp->math->md_one_third_inf_t));
    mp_double_free_number(mp, &(mp->math->md_inf_t));
    mp_double_free_number(mp, &(mp->math->md_negative_inf_t));
    mp_double_free_number(mp, &(mp->math->md_warning_limit_t));
    mp_double_free_number(mp, &(mp->math->md_one_k));
    mp_double_free_number(mp, &(mp->math->md_sqrt_8_e_k));
    mp_double_free_number(mp, &(mp->math->md_twelve_ln_2_k));
    mp_double_free_number(mp, &(mp->math->md_coef_bound_k));
    mp_double_free_number(mp, &(mp->math->md_coef_bound_minus_1));
    mp_double_free_number(mp, &(mp->math->md_fraction_threshold_t));
    mp_double_free_number(mp, &(mp->math->md_half_fraction_threshold_t));
    mp_double_free_number(mp, &(mp->math->md_scaled_threshold_t));
    mp_double_free_number(mp, &(mp->math->md_half_scaled_threshold_t));
    mp_double_free_number(mp, &(mp->math->md_near_zero_angle_t));
    mp_double_free_number(mp, &(mp->math->md_p_over_v_threshold_t));
    mp_double_free_number(mp, &(mp->math->md_equation_threshold_t));
    mp_memory_free(mp->math);
}

math_data *mp_initialize_double_math(MP mp)
{
    math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data));
    /* alloc */
    math->md_allocate        = mp_double_allocate_number;
    math->md_free            = mp_double_free_number;
    math->md_allocate_clone  = mp_double_allocate_clone;
    math->md_allocate_abs    = mp_double_allocate_abs;
    math->md_allocate_div    = mp_double_allocate_div;
    math->md_allocate_mul    = mp_double_allocate_mul;
    math->md_allocate_add    = mp_double_allocate_add;
    math->md_allocate_sub    = mp_double_allocate_sub;
    math->md_allocate_double = mp_double_allocate_double;
    /* precission */
    mp_double_allocate_number(mp, &math->md_precision_default, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_precision_max, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_precision_min, mp_scaled_type);
    /* here are the constants for |scaled| objects */
    mp_double_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_unity_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_two_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_three_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_zero_t, mp_scaled_type);
    /* |fractions| */
    mp_double_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
    /* |angles| */
    mp_double_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type);
    mp_double_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type);
    mp_double_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type);
    /* various approximations */
    mp_double_allocate_number(mp, &math->md_one_k, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type);
    /* thresholds */
    mp_double_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
    mp_double_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
    mp_double_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type);
    mp_double_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
    /* initializations */
    math->md_precision_default.data.dval         = 16 * unity; 
    math->md_precision_max.data.dval             = 16 * unity;
    math->md_precision_min.data.dval             = 16 * unity;
    math->md_epsilon_t.data.dval                 = epsilon;
    math->md_inf_t.data.dval                     = EL_GORDO;
    math->md_negative_inf_t.data.dval            = negative_EL_GORDO;
    math->md_one_third_inf_t.data.dval           = one_third_EL_GORDO;
    math->md_warning_limit_t.data.dval           = warning_limit;
    math->md_unity_t.data.dval                   = unity;
    math->md_two_t.data.dval                     = two;
    math->md_three_t.data.dval                   = three;
    math->md_half_unit_t.data.dval               = half_unit;
    math->md_three_quarter_unit_t.data.dval      = three_quarter_unit;
    math->md_arc_tol_k.data.dval                 = (unity/4096);                /* quit when change in arc length estimate reaches this */
    math->md_fraction_one_t.data.dval            = fraction_one;
    math->md_fraction_half_t.data.dval           = fraction_half;
    math->md_fraction_three_t.data.dval          = fraction_three;
    math->md_fraction_four_t.data.dval           = fraction_four;
    math->md_three_sixty_deg_t.data.dval         = three_sixty_deg;
    math->md_one_eighty_deg_t.data.dval          = one_eighty_deg;
    math->md_negative_one_eighty_deg_t.data.dval = negative_one_eighty_deg;
    math->md_one_k.data.dval                     = 1.0/64 ;
    math->md_sqrt_8_e_k.data.dval                = 1.71552776992141359295;       /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
    math->md_twelve_ln_2_k.data.dval             = 8.31776616671934371292 * 256; /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
    math->md_coef_bound_k.data.dval              = coef_bound;
    math->md_coef_bound_minus_1.data.dval        = coef_bound - 1/65536.0;
    math->md_twelvebits_3.data.dval              =     1365 / 65536.0;           /* $1365\approx 2^{12}/3$ */
    math->md_twentysixbits_sqrt2_t.data.dval     = 94906266 / 65536.0;           /* $2^{26}\sqrt2\approx94906265.62$ */
    math->md_twentyeightbits_d_t.data.dval       = 35596755 / 65536.0;           /* $2^{28}d\approx35596754.69$ */
    math->md_twentysevenbits_sqrt2_d_t.data.dval = 25170707 / 65536.0;           /* $2^{27}\sqrt2\,d\approx25170706.63$ */
    math->md_fraction_threshold_t.data.dval      = fraction_threshold;
    math->md_half_fraction_threshold_t.data.dval = half_fraction_threshold;
    math->md_scaled_threshold_t.data.dval        = scaled_threshold;
    math->md_half_scaled_threshold_t.data.dval   = half_scaled_threshold;
    math->md_near_zero_angle_t.data.dval         = near_zero_angle;
    math->md_p_over_v_threshold_t.data.dval      = p_over_v_threshold;
    math->md_equation_threshold_t.data.dval      = equation_threshold;
    /* functions */
    math->md_from_int                 = mp_double_set_from_int;
    math->md_from_boolean             = mp_double_set_from_boolean;
    math->md_from_scaled              = mp_double_set_from_scaled;
    math->md_from_double              = mp_double_set_from_double;
    math->md_from_addition            = mp_double_set_from_addition;
    math->md_half_from_addition       = mp_double_set_half_from_addition;
    math->md_from_subtraction         = mp_double_set_from_subtraction;
    math->md_half_from_subtraction    = mp_double_set_half_from_subtraction;
    math->md_from_of_the_way          = mp_double_set_from_of_the_way;
    math->md_from_div                 = mp_double_set_from_div;
    math->md_from_mul                 = mp_double_set_from_mul;
    math->md_from_int_div             = mp_double_set_from_int_div;
    math->md_from_int_mul             = mp_double_set_from_int_mul;
    math->md_negate                   = mp_double_negate;
    math->md_add                      = mp_double_add;
    math->md_subtract                 = mp_double_subtract;
    math->md_half                     = mp_double_half;
    math->md_do_double                = mp_double_double;
    math->md_abs                      = mp_double_abs;
    math->md_clone                    = mp_double_clone;
    math->md_negated_clone            = mp_double_negated_clone;
    math->md_abs_clone                = mp_double_abs_clone;
    math->md_swap                     = mp_double_swap;
    math->md_add_scaled               = mp_double_add_scaled;
    math->md_multiply_int             = mp_double_multiply_int;
    math->md_divide_int               = mp_double_divide_int;
    math->md_to_int                   = mp_double_to_int;
    math->md_to_boolean               = mp_double_to_boolean;
    math->md_to_scaled                = mp_double_to_scaled;
    math->md_to_double                = mp_double_to_double;
    math->md_odd                      = mp_double_odd;
    math->md_equal                    = mp_double_equal;
    math->md_less                     = mp_double_less;
    math->md_greater                  = mp_double_greater;
    math->md_non_equal_abs            = mp_double_non_equal_abs;
    math->md_round_unscaled           = mp_double_unscaled;
    math->md_floor_scaled             = mp_double_floor;
    math->md_fraction_to_round_scaled = mp_double_fraction_to_round_scaled;
    math->md_make_scaled              = mp_double_number_make_scaled;
    math->md_make_fraction            = mp_double_number_make_fraction;
    math->md_take_fraction            = mp_double_number_take_fraction;
    math->md_take_scaled              = mp_double_number_take_scaled;
    math->md_velocity                 = mp_double_velocity;
    math->md_n_arg                    = mp_double_n_arg;
    math->md_m_log                    = mp_double_m_log;
    math->md_m_exp                    = mp_double_m_exp;
    math->md_m_unif_rand              = mp_double_m_unif_rand;
    math->md_m_norm_rand              = mp_double_m_norm_rand;
    math->md_pyth_add                 = mp_double_pyth_add;
    math->md_pyth_add3                = mp_double_pyth_add3;
    math->md_pyth_sub                 = mp_double_pyth_sub;
    math->md_power_of                 = mp_double_power_of;
    math->md_fraction_to_scaled       = mp_double_fraction_to_scaled;
    math->md_scaled_to_fraction       = mp_double_scaled_to_fraction;
    math->md_scaled_to_angle          = mp_double_scaled_to_angle;
    math->md_angle_to_scaled          = mp_double_angle_to_scaled;
    math->md_init_randoms             = mp_double_init_randoms;
    math->md_sin_cos                  = mp_double_sin_cos;
    math->md_slow_add                 = mp_double_slow_add;
    math->md_slow_sub                 = mp_double_slow_sub;
    math->md_sqrt                     = mp_double_square_rt;
    math->md_print                    = mp_double_print_number;
    math->md_tostring                 = mp_double_number_tostring;
    math->md_modulo                   = mp_double_modulo;
    math->md_ab_vs_cd                 = mp_double_ab_vs_cd;
    math->md_crossing_point           = mp_double_crossing_point;
    math->md_scan_numeric             = mp_double_scan_numeric_token;
    math->md_scan_fractional          = mp_double_scan_fractional_token;
    /* housekeeping */
    math->md_free_math                = mp_double_free_math;
    math->md_set_precision            = mp_double_set_precision;
    return math;
}