/* 
    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. 
*/

# 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  odd(A)                  (abs(A) % 2 == 1)
# 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. 

*/

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

// printf("%f %f %f %f\n",p,q,p*q,(p * q) / fraction_multiplier);
// -4096.000000 0.000000 -0.000000 -0.000000

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

static inline double mp_double_make_scaled(double p, double q) { 
 // return p / q; 
    return p == 0.0 ? 0.0 : p / q; 
}

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

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

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

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

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

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

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

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

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

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

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

static void mp_double_set_from_scaled(mp_number *A, int B)
{
    A->data.dval = B == 0 ? 0.0 : B / 65536.0;
}

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

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

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

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

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

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

static void mp_double_set_from_mul(mp_number *A, mp_number *B, mp_number *C)
{
    A->data.dval = B->data.dval == 0.0 || C->data.dval == 0.0 ? 0.0 : B->data.dval * C->data.dval;
}

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

static void mp_double_set_from_int_mul(mp_number *A, mp_number *B, int C)
{
    A->data.dval = B->data.dval == 0.0 || C == 0 ? 0.0 : B->data.dval * C;
}

static void mp_double_set_from_of_the_way(MP mp, mp_number *A, mp_number *t, mp_number *B, 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 void mp_double_negate(mp_number *A)
{
    if (A->data.dval == -0.0) { /* already checked */
        A->data.dval = 0.0;
    } else if (A->data.dval != 0.0) {
        A->data.dval = -A->data.dval;
    }
}

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

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

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

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

static void mp_double_add_scaled(mp_number *A, int B)
{
    /* also for negative B */
    A->data.dval += (B / 65536.0);
}

static void mp_double_multiply_int(mp_number *A, int B)
{
    if (A->data.dval != 0.0) {
        A->data.dval *= (double) B;
    }
}

static void mp_double_divide_int(mp_number *A, int B)
{
    if (A->data.dval != 0.0) {
        A->data.dval /= (double) B;
    }
}

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

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

static void mp_double_negated_clone(mp_number *A, mp_number *B)
{
    A->data.dval = -B->data.dval;
    if (A->data.dval == -0.0) {
        A->data.dval = 0.0;
    }
}

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

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

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

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

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

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

static int mp_double_to_scaled(mp_number *A)
{
    return (int) lround(A->data.dval * 65536.0);
}

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

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

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

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

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

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

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

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

static char *mp_double_number_tostring(MP mp, 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 void mp_double_print_number(MP mp, 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, mp_number *x_orig, 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 = 1;
            ret->data.dval = EL_GORDO;
        }
    } else if (-y <= EL_GORDO + x) {
        ret->data.dval = x + y;
    } else {
        mp->arithmic_error = 1;
        ret->data.dval = negative_EL_GORDO;
    }
}

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

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

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

static void mp_double_number_make_scaled(MP mp, mp_number *ret, mp_number *p_orig, 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, int 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, int 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++;
    }
    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, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, 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 int mp_double_ab_vs_cd(mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, 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, mp_number *aa, mp_number *bb, 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 int mp_double_unscaled(mp_number *x_orig)
{
    return (int) lround(x_orig->data.dval);
}

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

static 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, 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, mp_number *a_orig, mp_number *b_orig)
{
    double a = fabs(a_orig->data.dval);
    double b = fabs(b_orig->data.dval);
    errno = 0;
    ret->data.dval = sqrt(a*a + b*b);
    if (errno) {
        mp->arithmic_error = 1;
        ret->data.dval = EL_GORDO;
    }
}

static void mp_double_pyth_sub(MP mp, mp_number *ret, mp_number *a_orig, 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, mp_number *a_orig, mp_number *b_orig)
{
    errno = 0;
    ret->data.dval = pow(a_orig->data.dval, b_orig->data.dval);
    if (errno) {
        mp->arithmic_error = 1;
        ret->data.dval = EL_GORDO;
    }
}

static void mp_double_m_log(MP mp, mp_number *ret, 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, 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 = 1;
            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, mp_number *x_orig, mp_number *y_orig)
{
    /*        
        if (x_orig->data.dval == -0.0) {
            x_orig->data.dval = 0.0;
        }
        if (y_orig->data.dval == -0.0) {
            y_orig->data.dval = 0.0;
        }
    */
    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;
        if (ret->data.dval == -0.0) {
            ret->data.dval = 0.0;
        }
// printf("D x=%f y=%f atan=%f\n",y_orig->data.dval,x_orig->data.dval,ret->data.dval);
    }
}

static void mp_double_sin_cos(MP mp, 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.
*/

# define KK            100                /* the long lag  */
# define LL            37                 /* the short lag */
# define MM            (1L<<30)           /* the modulus   */
# define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
# define TT            70                 /* guaranteed separation between streams */
# define is_odd(x)     ((x)&1)            /* units bit of x */
# define QUALITY       1009               /* recommended quality level for high-res use */

/*tex 
    The destination array length (must be at least KK).
*/

typedef struct mp_double_random_info {
    long  x[KK];
    long  buf[QUALITY];
    long  dummy;
    long  started;
    long *ptr;
} mp_double_random_info;

static mp_double_random_info mp_double_random_data = {
    .dummy   = -1,
    .started = -1,
    .ptr     = &mp_double_random_data.dummy
};

static void mp_double_aux_ran_array(long aa[], int n)
{
    int i, j;
    for (j = 0; j < KK; j++) {
        aa[j] = mp_double_random_data.x[j];
    }
    for (; j < n; j++) {
        aa[j] = mod_diff(aa[j - KK], aa[j - LL]);
    }
    for (i = 0; i < LL; i++, j++) {
        mp_double_random_data.x[i] = mod_diff(aa[j - KK], aa[j - LL]);
    }
    for (; i < KK; i++, j++) {
        mp_double_random_data.x[i] = mod_diff(aa[j - KK], mp_double_random_data.x[i - LL]);
    }
}

static void mp_double_aux_ran_start(long seed)
{
    int t, j;
    long x[KK + KK - 1]; /* the preparation buffer */
    long ss = (seed+2) & (MM - 2);
    for (j = 0; j < KK; j++) {
        /* bootstrap the buffer */
        x[j] = ss;
        /* cyclic shift 29 bits */
        ss <<= 1;
        if (ss >= MM) {
            ss -= MM - 2;
        }
    }
    /* make x[1] (and only x[1]) odd */
    x[1]++;
    for (ss = seed & (MM - 1), t = TT - 1; t;) {
        for (j = KK - 1; j > 0; j--) {
            /* "square" */
            x[j + j] = x[j];
            x[j + j - 1] = 0;
        }
        for (j = KK + KK - 2; j >= KK; j--) {
            x[j - (KK -LL)] = mod_diff(x[j - (KK - LL)], x[j]);
            x[j - KK] = mod_diff(x[j - KK], x[j]);
        }
        if (is_odd(ss)) {
            /* "multiply by z" */
            for (j = KK; j>0; j--) {
                x[j] = x[j-1];
            }
            x[0] = x[KK];
            /* shift the buffer cyclically */
            x[LL] = mod_diff(x[LL], x[KK]);
        }
        if (ss) {
            ss >>= 1;
        } else {
            t--;
        }
    }
    for (j = 0; j < LL; j++) {
        mp_double_random_data.x[j + KK - LL] = x[j];
    }
    for (;j < KK; j++) {
        mp_double_random_data.x[j - LL] = x[j];
    }
    for (j = 0; j < 10; j++) {
        /* warm things up */
        mp_double_aux_ran_array(x, KK + KK - 1);
    }
    mp_double_random_data.ptr = &mp_double_random_data.started;
}

static long mp_double_aux_ran_arr_cycle(void)
{
    if (mp_double_random_data.ptr == &mp_double_random_data.dummy) {
        /* the user forgot to initialize */
        mp_double_aux_ran_start(314159L);
    }
    mp_double_aux_ran_array(mp_double_random_data.buf, QUALITY);
    mp_double_random_data.buf[KK] = -1;
    mp_double_random_data.ptr = mp_double_random_data.buf + 1;
    return mp_double_random_data.buf[0];
}

static void mp_double_init_randoms(MP mp, int seed)
{
    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((unsigned long) seed);
}

static void mp_double_modulo(mp_number *a, 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_double_random_data.ptr>=0? *mp_double_random_data.ptr++: mp_double_aux_ran_arr_cycle());
    double a = op / (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, 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_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_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;
}