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

# include "mpmathdecimal.h"

# define DECNUMDIGITS 1000

# include "decNumber.h"

# define E_STRING                "2.7182818284590452353602874713526624977572470936999595749669676277240766303535"
# define PI_STRING               "3.1415926535897932384626433832795028841971693993751058209749445923078164062862"

# define fraction_multiplier     4096
# define angle_multiplier        16
# define unity                   1
# define two                     2
# define three                   3
# define four                    4

# define half_unit               0.5
# define three_quarter_unit      0.75

# 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 epsilon                 pow(2.0,-173.0)
# define epsilonf                pow(2.0,-52.0)

# define EL_GORDO                "1E1000000"
# define negative_EL_GORDO       "-1E1000000"
# define warning_limit           "1E1000000"

# define DECPRECISION_DEFAULT    34
# define FACTORIALS_CACHESIZE    50

# define too_precise(a)          (a == (DEC_Inexact + DEC_Rounded))
# define too_large(a)            (a & DEC_Overflow)

# define fraction_half           (fraction_multiplier / 2)
# define fraction_one            (1 * fraction_multiplier)
# define fraction_two            (2 * fraction_multiplier)
# define fraction_three          (3 * fraction_multiplier)
# define fraction_four           (4 * fraction_multiplier)

# define no_crossing             mp_decimal_data.fraction_one_plus_decNumber
# define one_crossing            mp_decimal_data.fraction_one_decNumber
# define zero_crossing           mp_decimal_data.zero

# define odd(A)                  (abs(A) % 2 == 1)
# define set_cur_cmd(A)          mp->cur_mod_->type = (A)
# define set_cur_mod(A)          decNumberCopy((decNumber *) (mp->cur_mod_->data.n.data.num), A)

/* This one saves some typing and also looks better: */

# define  decNumberIsPositive(A)  (! (decNumberIsZero(A) || decNumberIsNegative(A)))

typedef struct mp_decimal_info {
    decContext   set;
    decContext   limitedset;
    decNumber    zero;
    decNumber    one;
    decNumber    minusone;
    decNumber    two_decNumber;
    decNumber    three_decNumber;
    decNumber    four_decNumber;
    decNumber    fraction_multiplier_decNumber;
    decNumber    angle_multiplier_decNumber;
    decNumber    fraction_one_decNumber;
    decNumber    fraction_one_plus_decNumber;
    decNumber    PI_decNumber;
    decNumber    epsilon_decNumber;
    decNumber    EL_GORDO_decNumber;
    decNumber    negative_EL_GORDO_decNumber;
    decNumber  **factorials;
    int          last_cached_factorial;
    int          initialized;
} mp_decimal_info;

static mp_decimal_info mp_decimal_data = {
    .factorials            = NULL,
    .last_cached_factorial = 0,
    .initialized           = 0,
};

static void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q);
static void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q);

static void mp_decnumber_check(MP mp, decNumber *dec, decContext *context)
{
    int test = 0;
    (void) mp;
    if (context->status & DEC_Overflow) {
        test = 1;
        context->status &= ~DEC_Overflow;
    }
    if (context->status & DEC_Underflow) {
        test = 1;
        context->status &= ~DEC_Underflow;
    }
    if (context->status & DEC_Errors) {
        test = 1;
        decNumberZero(dec);
    }
    context->status = 0;
    if (decNumberIsSpecial(dec)) {
        test = 1;
        if (decNumberIsInfinite(dec)) {
            if (decNumberIsNegative(dec)) {
                decNumberCopyNegate(dec, &mp_decimal_data.EL_GORDO_decNumber);
            } else {
                decNumberCopy(dec, &mp_decimal_data.EL_GORDO_decNumber);
            }
        } else {
            /* Nan */
            decNumberZero(dec);
        }
    }
    if (decNumberIsZero(dec) && decNumberIsNegative(dec)) {
        decNumberZero(dec);
    }
    mp->arithmic_error = test;
}

/* See mpmathdouble for documentation. */

static void checkZero(decNumber *ret)
{
    if (decNumberIsZero(ret) && decNumberIsNegative(ret)) {
        decNumberZero(ret);
    }
}

static int decNumberLess(decNumber *a, decNumber *b)
{
    decNumber comp;
    decNumberCompare(&comp, a, b, &mp_decimal_data.set);
    return decNumberIsNegative(&comp);
}

static int decNumberGreater(decNumber *a, decNumber *b)
{
    decNumber comp;
    decNumberCompare(&comp, a, b, &mp_decimal_data.set);
    return decNumberIsPositive(&comp);
}

static void decNumberFromDouble(decNumber *A, double B)
{
    char buffer[1000];
    char *c = buffer;
    snprintf(buffer, 1000, "%-650.325lf", B);
    while (*c++) {
        if (*c == ' ') {
            *c = '\0';
            break;
        }
    }
    decNumberFromString(A, buffer, &mp_decimal_data.set);
}

static double decNumberToDouble(decNumber *A)
{
    char *buffer = mp_memory_allocate(A->digits + 14);
    double res = 0.0;
    decNumberToString(A, buffer);
    if (sscanf(buffer, "%lf", &res)) {
        mp_memory_free(buffer);
        return res;
    } else {
        mp_memory_free(buffer);
        /* |mp->arithmic_error = 1;| */
        return 0.0;
    }
}

/*tex

    \startformula
    \arctan(x) = x - \frac {x^3}{3} + \frac {x^5{5} - \frac {x^7}{7} + \ldots
    \stopformula

    This power series works well, if $x$ is close to zero ($|x|<0.5$). If x is larger, the series 
    converges too slowly, so in order to get a smaller x, we apply the identity

    \startformula
    \arctan(x) = 2 \arctan \left (\frac {\sqrt{1 + x^2}-1} {x} \right)
    \stopformula

    twice. The first application gives us a new $x$ with $x < 1$. The second application gives us 
    a new x with $x < 0.4142136$. For that $x$, we use the power series and multiply the result by 
    four.

*/

static void decNumberAtan(decNumber *result, decNumber *x_orig, decContext *localset)
{
    decNumber x;
    decNumberCopy(&x, x_orig);
    if (decNumberIsZero(&x)) {
        decNumberCopy(result, &x);
    } else {
        decNumber f, g, mx2, term;
        for (int i = 0; i < 2; i++) {
            decNumber y;
            decNumberMultiply(&y, &x, &x, localset);                   /* $ y = x^2     $ */
            decNumberAdd(&y, &y, &mp_decimal_data.one, localset);      /* $ y = y + 1   $ */
            decNumberSquareRoot(&y, &y, localset);                     /* $ y = sqrt(y) $ */
            decNumberSubtract(&y, &y, &mp_decimal_data.one, localset); /* $ y = y - 1   $ */
            decNumberDivide(&x, &y, &x, localset);                     /* $ x = y / x   $ */
            if (decNumberIsZero(&x)) {
                decNumberCopy(result, &x);
                return;
            }
        }
        decNumberCopy(&f, &x);                     /* $ f(0) =  x   $ */
        decNumberCopy(&g, &mp_decimal_data.one);   /* $ g(0) =  1   $ */
        decNumberCopy(&term, &x);                  /* $ term =  x   $ */
        decNumberCopy(result, &x);                 /* $ sum  =  x   $ */
        decNumberMultiply(&mx2, &x, &x, localset); /* $ mx2  =  x^2 $ */
        decNumberMinus (&mx2, &mx2, localset);     /* $ mx2  = -x^2 $ */
        for (int i = 0; i < 2 * localset->digits; i++) {
            decNumberMultiply(&f, &f, &mx2, localset);
            decNumberAdd(&g, &g, &mp_decimal_data.two_decNumber, localset);
            decNumberDivide(&term, &f, &g, localset);
            decNumberAdd(result, result, &term, localset);
        }
        decNumberAdd(result, result, result, localset);
        decNumberAdd(result, result, result, localset);
    }
}

static void decNumberAtan2(decNumber *result, decNumber *y, decNumber *x, decContext *localset)
{
    if (! decNumberIsInfinite(x) && ! decNumberIsZero(y) && ! decNumberIsInfinite(y) && ! decNumberIsZero(x)) {
        decNumber temp;
        decNumberDivide(&temp, y, x, localset);
        decNumberAtan(result, &temp, localset);
        /*
            decNumberAtan doesn't quite return the values in the ranges we
            want for x < 0. So we need to do some correction
        */
        if (decNumberIsNegative(x)) {
            if (decNumberIsNegative(y)) {
                decNumberSubtract(result, result, &mp_decimal_data.PI_decNumber, localset);
            } else {
                decNumberAdd(result, result, &mp_decimal_data.PI_decNumber, localset);
            }
        }
    } else {
        if (decNumberIsInfinite(y) && decNumberIsInfinite(x)) {
            /* If x and y are both inf, the result depends on the sign of x */
            decNumberDivide(result, &mp_decimal_data.PI_decNumber, &mp_decimal_data.four_decNumber, localset);
            if (decNumberIsNegative(x) ) {
                decNumber a;
                decNumberFromDouble(&a, 3.0);
                decNumberMultiply(result, result, &a, localset);
            }
        } else if (! decNumberIsZero(y) && ! decNumberIsInfinite(x) ) {
            /* If y is non-zero and x is non-inf, the result is +-pi/2 */
            decNumberDivide(result, &mp_decimal_data.PI_decNumber, &mp_decimal_data.two_decNumber, localset);
        } else {
            /* Otherwise it is +0 if x is positive, +pi if x is neg */
            if (decNumberIsNegative(x)) {
                decNumberCopy(result, &mp_decimal_data.PI_decNumber);
            } else {
                decNumberZero(result);
            }
        }
        /* Atan2 will be negative if y < 0 */
        if (decNumberIsNegative(y)) {
            decNumberMinus(result, result, localset);
        }
    }
}

static void mp_decimal_allocate_number(MP mp, mp_number *n, mp_number_type t)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = t;
    decNumberZero(n->data.num);
}

static void mp_decimal_allocate_clone(MP mp, mp_number *n, mp_number_type t, mp_number *v)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = t;
    decNumberZero(n->data.num);
    decNumberCopy(n->data.num, v->data.num);
}

static void mp_decimal_allocate_abs(MP mp, mp_number *n, mp_number_type t, mp_number *v)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = t;
    decNumberZero(n->data.num); /* not needed */
    decNumberAbs(n->data.num, v->data.num, &mp_decimal_data.set);
}

static void mp_decimal_allocate_div(MP mp, mp_number *n, mp_number_type t, mp_number *a, mp_number *b)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = t;
    decNumberZero(n->data.num); /* not needed */
    decNumberDivide(n->data.num, a->data.num, b->data.num, &mp_decimal_data.set);
}

static void mp_decimal_allocate_mul(MP mp, mp_number *n, mp_number_type t, mp_number *a, mp_number *b)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = t;
    decNumberZero(n->data.num); /* not needed */
    decNumberMultiply(n->data.num, a->data.num, b->data.num, &mp_decimal_data.set);
}

static void mp_decimal_allocate_add(MP mp, mp_number *n, mp_number_type t, mp_number *a, mp_number *b)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = t;
    decNumberZero(n->data.num); /* not needed */
    decNumberAdd(n->data.num, a->data.num, b->data.num, &mp_decimal_data.set);
}

static void mp_decimal_allocate_sub(MP mp, mp_number *n, mp_number_type t, mp_number *a, mp_number *b)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = t;
    decNumberZero(n->data.num); /* not needed */
    decNumberSubtract(n->data.num, a->data.num, b->data.num, &mp_decimal_data.set);
}

static void mp_decimal_allocate_double(MP mp, mp_number *n, double v)
{
    (void) mp;
    n->data.num = mp_memory_allocate(sizeof(decNumber));
    n->type = mp_scaled_type;
    decNumberZero(n->data.num); /* not needed */
    decNumberFromDouble(n->data.num, v);
}

static void mp_decimal_free_number(MP mp, mp_number *n)
{
    (void) mp;
    if (n->data.num) {
        mp_memory_free(n->data.num);
        n->data.num = NULL;
        n->type = mp_nan_type;
    }
}

static void mp_decimal_set_from_int(mp_number *A, mp_scaled_t B)
{
    decNumberFromScaled(A->data.num, B);
}

static void mp_decimal_set_from_boolean(mp_number *A, mp_scaled_t B)
{
    decNumberFromScaled(A->data.num, B);
}

static void mp_decimal_set_from_scaled(mp_number *A, mp_scaled_t B)
{
    decNumber c;
    decNumberFromScaled(&c, 65536);
    decNumberFromScaled(A->data.num, B);
    decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_double(mp_number *A, double B)
{
    decNumberFromDouble(A->data.num, B);
}

static void mp_decimal_set_from_addition(mp_number *A, mp_number *B, mp_number *C)
{
    decNumberAdd(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
}

static void mp_decimal_set_half_from_addition(mp_number *A, mp_number *B, mp_number *C)
{
    decNumber c;
    decNumberAdd(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
    decNumberFromScaled(&c, 2);
    decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
{
    decNumberSubtract(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
}

static void mp_decimal_set_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
{
    decNumber c;
    decNumberSubtract(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
    decNumberFromScaled(&c, 2);
    decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_div(mp_number *A, mp_number *B, mp_number *C)
{
    decNumberDivide(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
}

static void mp_decimal_set_from_mul(mp_number *A, mp_number *B, mp_number *C)
{
    decNumberMultiply(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
}

static void mp_decimal_set_from_int_div(mp_number *A, mp_number *B, mp_scaled_t C)
{
    decNumber c;
    decNumberFromScaled(&c, C);
    decNumberDivide(A->data.num, B->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_int_mul(mp_number *A, mp_number *B, mp_scaled_t C)
{
    decNumber c;
    decNumberFromScaled(&c, C);
    decNumberMultiply(A->data.num, B->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_of_the_way(MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C)
{
    decNumber c;
    decNumber r1;
    decNumberSubtract(&c, B->data.num, C->data.num, &mp_decimal_data.set);
    mp_decimal_take_fraction(mp, &r1, &c, t->data.num);
    decNumberSubtract(A->data.num, B->data.num, &r1, &mp_decimal_data.set);
    mp_decnumber_check(mp, A->data.num, &mp_decimal_data.set);
}

static void mp_decimal_negate(mp_number *A)
{
    decNumberCopyNegate(A->data.num, A->data.num);
    checkZero(A->data.num);
}

static void mp_decimal_add(mp_number *A, mp_number *B)
{
    decNumberAdd(A->data.num, A->data.num, B->data.num, &mp_decimal_data.set);
}

static void mp_decimal_subtract(mp_number *A, mp_number *B)
{
    decNumberSubtract(A->data.num, A->data.num, B->data.num, &mp_decimal_data.set);
}

static void mp_decimal_half(mp_number *A)
{
    decNumber c;
    decNumberFromScaled(&c, 2);
    decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_double(mp_number *A)
{
    decNumber c;
    decNumberFromScaled(&c, 2);
    decNumberMultiply(A->data.num, A->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_add_scaled(mp_number *A, mp_scaled_t B)
{
    decNumber b, c;
    decNumberFromScaled(&c, 65536);
    decNumberFromScaled(&b, B);
    decNumberDivide(&b, &b, &c, &mp_decimal_data.set);
    decNumberAdd(A->data.num, A->data.num, &b, &mp_decimal_data.set);
}

static void mp_decimal_multiply_int(mp_number *A, mp_scaled_t B)
{
    decNumber b;
    decNumberFromScaled(&b, B);
    decNumberMultiply(A->data.num, A->data.num, &b, &mp_decimal_data.set);
}

static void mp_decimal_divide_int(mp_number *A, mp_scaled_t B)
{
    decNumber b;
    decNumberFromScaled(&b, B);
    decNumberDivide(A->data.num, A->data.num, &b, &mp_decimal_data.set);
}

static void mp_decimal_abs(mp_number *A)
{
    decNumberAbs(A->data.num, A->data.num, &mp_decimal_data.set);
}

static void mp_decimal_clone(mp_number *A, mp_number *B)
{
    decNumberCopy(A->data.num, B->data.num);
}

static void mp_decimal_negated_clone(mp_number *A, mp_number *B)
{
    decNumberCopyNegate(A->data.num, B->data.num);
    checkZero(A->data.num);
}

static void mp_decimal_abs_clone(mp_number *A, mp_number *B)
{
    decNumberAbs(A->data.num, B->data.num, &mp_decimal_data.set);
}

static void mp_decimal_swap(mp_number *A, mp_number *B)
{
    decNumber swap_tmp;
    decNumberCopy(&swap_tmp,   A->data.num);
    decNumberCopy(A->data.num, B->data.num);
    decNumberCopy(B->data.num, &swap_tmp);
}

static void mp_decimal_fraction_to_scaled(mp_number *A)
{
    A->type = mp_scaled_type;
    decNumberDivide(A->data.num, A->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set);
}

static void mp_decimal_angle_to_scaled(mp_number *A)
{
    A->type = mp_scaled_type;
    decNumberDivide(A->data.num, A->data.num, &mp_decimal_data.angle_multiplier_decNumber, &mp_decimal_data.set);
}

static void mp_decimal_scaled_to_fraction(mp_number *A)
{
    A->type = mp_fraction_type;
    decNumberMultiply(A->data.num, A->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set);
}

static void mp_decimal_scaled_to_angle(mp_number *A)
{
    A->type = mp_angle_type;
    decNumberMultiply(A->data.num, A->data.num, &mp_decimal_data.angle_multiplier_decNumber, &mp_decimal_data.set);
}

static mp_scaled_t mp_decimal_to_scaled(mp_number *A)
{
    mp_scaled_t result;
    decNumber corrected;
    decNumberFromScaled(&corrected, 65536);
    decNumberMultiply(&corrected, &corrected, A->data.num, &mp_decimal_data.set);
    decNumberReduce(&corrected, &corrected, &mp_decimal_data.set);
    result = mpscaledround(decNumberToDouble(&corrected));
    return result;
}

static mp_scaled_t mp_decimal_to_int(mp_number *A)
{
    mp_scaled_t result;
    mp_decimal_data.set.status = 0;
    result = decNumberToScaled(A->data.num, &mp_decimal_data.set);
    if (mp_decimal_data.set.status == DEC_Invalid_operation) {
        mp_decimal_data.set.status = 0;
        /* |mp->arithmic_error = 1;| */
        return 0;
    } else {
        return result;
    }
}

static mp_scaled_t mp_decimal_to_boolean(mp_number *A)
{
    mp_scaled_t result;
    mp_decimal_data.set.status = 0;
    result = decNumberToScaled(A->data.num, &mp_decimal_data.set);
    if (mp_decimal_data.set.status == DEC_Invalid_operation) {
        mp_decimal_data.set.status = 0;
        /* |mp->arithmic_error = 1;| */
        return 0;
    } else {
        return result;
    }
}

static double mp_decimal_to_double(mp_number *A)
{
    char *buffer = mp_memory_allocate((size_t) ((decNumber *) A->data.num)->digits + 14);
    double res = 0.0;
    decNumberToString(A->data.num, buffer);
    if (sscanf(buffer, "%lf", &res)) {
        mp_memory_free(buffer);
        return res;
    } else {
        mp_memory_free(buffer);
        /* |mp->arithmic_error = 1;| */
        return 0.0;
    }
}

static int mp_decimal_odd(mp_number *A)
{
    decNumber r1, r2;
    decNumberAbs(&r1, A->data.num, &mp_decimal_data.set);
    decNumberRemainder(&r2, &r1, &mp_decimal_data.two_decNumber, &mp_decimal_data.set);
    decNumberCompare(&r1, &r2, &mp_decimal_data.one, &mp_decimal_data.set);
    return decNumberIsZero(&r1);
}

static int mp_decimal_equal(mp_number *A, mp_number *B)
{
    decNumber res;
    decNumberCompare(&res, A->data.num, B->data.num, &mp_decimal_data.set);
    return decNumberIsZero(&res);
}

static int mp_decimal_greater(mp_number *A, mp_number *B)
{
    decNumber res;
    decNumberCompare(&res, A->data.num, B->data.num, &mp_decimal_data.set);
    return decNumberIsPositive(&res);
}

static int mp_decimal_less(mp_number *A, mp_number *B)
{
    decNumber res;
    decNumberCompare(&res, A->data.num, B->data.num, &mp_decimal_data.set);
    return decNumberIsNegative(&res);
}

static int mp_decimal_non_equal_abs(mp_number *A, mp_number *B)
{
    decNumber res, a, b;
    decNumberCopyAbs(&a, A->data.num);
    decNumberCopyAbs(&b, B->data.num);
    decNumberCompare(&res, &a, &b, &mp_decimal_data.set);
    return ! decNumberIsZero(&res);
}

static char *mp_decnumber_tostring(decNumber *n)
{
    decNumber corrected;
    char *buffer = mp_memory_allocate((size_t) ((decNumber *) n)->digits + 14);
    decNumberCopy(&corrected, n);
    decNumberTrim(&corrected);
    decNumberToString(&corrected, buffer);
    return buffer;
}

static char *mp_decimal_number_tostring(MP mp, mp_number *n)
{
    (void) mp;
    return mp_decnumber_tostring(n->data.num);
}

static void mp_decimal_print_number(MP mp, mp_number *n)
{
    char *str = mp_decnumber_tostring(n->data.num);
    mp_print_e_str(mp, str);
    mp_memory_free(str);
}

static void mp_decimal_slow_add(MP mp, mp_number *ret, mp_number *A, mp_number *B)
{
    (void) mp;
    decNumberAdd(ret->data.num, A->data.num, B->data.num, &mp_decimal_data.set);
}

static void mp_decimal_make_fraction(MP mp, decNumber *ret, decNumber *p, decNumber *q)
{
    decNumberDivide(ret, p, q, &mp_decimal_data.set);
    mp_decnumber_check(mp, ret, &mp_decimal_data.set);
    decNumberMultiply(ret, ret, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set);
}

static void mp_decimal_number_make_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q)
{
    mp_decimal_make_fraction(mp, ret->data.num, p->data.num, q->data.num);
}

static void mp_decimal_take_fraction(MP mp, decNumber *ret, decNumber *p, decNumber *q)
{
    (void) mp;
    decNumberMultiply(ret, p, q, &mp_decimal_data.set);
    decNumberDivide(ret, ret, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set);
}

static void mp_decimal_number_take_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q)
{
    mp_decimal_take_fraction(mp, ret->data.num, p->data.num, q->data.num);
}

static void mp_decimal_number_take_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
{
    (void) mp;
    decNumberMultiply(ret->data.num, p_orig->data.num, q_orig->data.num, &mp_decimal_data.set);
}

static void mp_decimal_number_make_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
{
    decNumberDivide(ret->data.num, p_orig->data.num, q_orig->data.num, &mp_decimal_data.set);
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_decimal_aux_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop)
{
    decNumber result;
    size_t l = stop-start+1;
    char *buf = mp_memory_allocate(l + 1);
    buf[l] = '\0';
    (void) strncpy(buf, (const char *) start, l);
    mp_decimal_data.set.status = 0;
    decNumberFromString(&result,buf, &mp_decimal_data.set);
    mp_memory_free(buf);
    if (mp_decimal_data.set.status == 0) {
        set_cur_mod(&result);
    } else if (mp->scanner_status != mp_tex_flushing_state) {
        if (too_large(mp_decimal_data.set.status)) {
            mp_decnumber_check(mp, &result, &mp_decimal_data.set);
            set_cur_mod(&result);
            mp_error(
                mp,
                "Enormous number has been reduced",
                "I could not handle this number specification because it is out of range."
            );
        } else if (too_precise(mp_decimal_data.set.status)) {
            set_cur_mod(&result);
            if (decNumberIsPositive((decNumber *) internal_value(mp_warning_check_internal).data.num) && (mp->scanner_status != mp_tex_flushing_state)) {
                char msg[256];
                snprintf (msg, 256, "Number is too precise (numberprecision = %d)", mp_decimal_data.set.digits);
                mp_error(
                    mp,
                    msg,
                    "Continue and I'll round the value until it fits the current numberprecision\n"
                    "(Set warningcheck:=0 to suppress this message.)"
                );
            }
        } else {
            /* this also captures underflow */
            mp_error(
                mp,
                "Erroneous number specification changed to zero",
                "I could not handle this number specification"
            );
            decNumberZero(&result);
            set_cur_mod(&result);
        }
    }
    set_cur_cmd((mp_variable_type) mp_numeric_command);
}

static void mp_decimal_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_decimal_scan_fractional_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++;
    }
    mp_decimal_aux_find_exponent(mp);
    stop = &mp->buffer[mp->cur_input.loc_field-1];
    mp_decimal_aux_wrapup_numeric_token(mp, start, stop);
}

static void mp_decimal_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_decimal_aux_find_exponent(mp);
    stop = &mp->buffer[mp->cur_input.loc_field-1];
    mp_decimal_aux_wrapup_numeric_token(mp, start, stop);
}

static void mp_decimal_velocity(MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t)
{
    decNumber acc, num, denom; /* registers for intermediate calculations */
    decNumber r1, r2;
    decNumber arg1, arg2;
    decNumber i16, fone, fhalf, ftwo, sqrtfive;
    decNumberFromScaled(&i16, 16);
    decNumberFromScaled(&fone, fraction_one);
    decNumberFromScaled(&fhalf, fraction_half);
    decNumberFromScaled(&ftwo, fraction_two);
    decNumberFromScaled(&sqrtfive, 5);
    decNumberSquareRoot(&sqrtfive, &sqrtfive, &mp_decimal_data.set);

    decNumberDivide(&arg1, sf->data.num, &i16, &mp_decimal_data.set);   /* arg1 = sf / 16*/
    decNumberSubtract(&arg1, st->data.num,&arg1, &mp_decimal_data.set); /* arg1 = st - arg1*/
    decNumberDivide(&arg2, st->data.num, &i16, &mp_decimal_data.set);   /* arg2 = st / 16*/
    decNumberSubtract(&arg2, sf->data.num,&arg2, &mp_decimal_data.set); /* arg2 = sf - arg2*/
    mp_decimal_take_fraction(mp, &acc, &arg1, &arg2);                   /* acc = (arg1 * arg2) / fmul*/

    decNumberCopy(&arg1, &acc);
    decNumberSubtract(&arg2, ct->data.num, cf->data.num, &mp_decimal_data.set); /* arg2 = ct - cf*/
    mp_decimal_take_fraction(mp, &acc, &arg1, &arg2);                           /* acc = (arg1 * arg2 ) / fmul*/

    decNumberSquareRoot(&arg1, &mp_decimal_data.two_decNumber, &mp_decimal_data.set); /* arg1 = $\sqrt{2}$*/
    decNumberMultiply(&arg1, &arg1, &fone, &mp_decimal_data.set);                     /* arg1 = arg1 * fmul*/
    mp_decimal_take_fraction(mp, &r1, &acc, &arg1);                                   /* r1 = (acc * arg1) / fmul*/
    decNumberAdd(&num, &ftwo, &r1, &mp_decimal_data.set);                             /* num = ftwo + r1*/

    decNumberSubtract(&arg1,&sqrtfive, &mp_decimal_data.one, &mp_decimal_data.set);        /* arg1 = $\sqrt{5}$ - 1*/
    decNumberMultiply(&arg1,&arg1,&fhalf, &mp_decimal_data.set);                           /* arg1 = arg1 * fmul/2*/
    decNumberMultiply(&arg1,&arg1,&mp_decimal_data.three_decNumber, &mp_decimal_data.set); /* arg1 = arg1 * 3*/

    decNumberSubtract(&arg2,&mp_decimal_data.three_decNumber, &sqrtfive, &mp_decimal_data.set); /* arg2 = 3 - $\sqrt{5}$*/
    decNumberMultiply(&arg2,&arg2, &fhalf, &mp_decimal_data.set);                               /* arg2 = arg2 * fmul/2*/
    decNumberMultiply(&arg2,&arg2, &mp_decimal_data.three_decNumber, &mp_decimal_data.set);     /* arg2 = arg2 * 3*/
    mp_decimal_take_fraction(mp, &r1, ct->data.num, &arg1) ;                                    /* r1 = (ct * arg1) / fmul*/
    mp_decimal_take_fraction(mp, &r2, cf->data.num, &arg2);                                     /* r2 = (cf * arg2) / fmul*/

    decNumberFromScaled(&denom, fraction_three);              /* denom = 3fmul*/
    decNumberAdd(&denom, &denom, &r1, &mp_decimal_data.set); /* denom = denom + r1*/
    decNumberAdd(&denom, &denom, &r2, &mp_decimal_data.set); /* denom = denom + r1*/

    decNumberCompare(&arg1, t->data.num, &mp_decimal_data.one, &mp_decimal_data.set);
    if (! decNumberIsZero(&arg1)) { /* t != r1*/
        decNumberDivide(&num, &num, t->data.num, &mp_decimal_data.set); /* num = num / t */
    }
    decNumberCopy(&r2, &num); /* r2 = num / 4*/
    decNumberDivide(&r2, &r2, &mp_decimal_data.four_decNumber, &mp_decimal_data.set);
    if (decNumberLess(&denom, &r2)) {
        decNumberFromScaled(ret->data.num, fraction_four); /* num/4 >= denom => denom < num/4*/
    } else {
        mp_decimal_make_fraction(mp, ret->data.num, &num, &denom);
    }
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static int mp_decimal_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig)
{
    decNumber a, b, c, d;
    decNumber ab, cd;
    decNumberCopy(&a, (decNumber *) a_orig->data.num);
    decNumberCopy(&b, (decNumber *) b_orig->data.num);
    decNumberCopy(&c, (decNumber *) c_orig->data.num);
    decNumberCopy(&d, (decNumber *) d_orig->data.num);
    decNumberMultiply(&ab, (decNumber *) a_orig->data.num, (decNumber *)b_orig->data.num, &mp_decimal_data.set);
    decNumberMultiply(&cd, (decNumber *) c_orig->data.num, (decNumber *)d_orig->data.num, &mp_decimal_data.set);
    if (decNumberLess(&ab, &cd)) {
        return -1;
    } else if (decNumberGreater(&ab, &cd)) {
        return 1;
    } else {
        return 0;
    }
}

static void mp_decimal_crossing_point(MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc)
{
    decNumber a, b, c;
    double d; /* recursive counter */
    decNumber x, xx, x0, x1, x2; /* temporary registers for bisection */
    decNumber scratch, scratch2;
    decNumberCopy(&a, (decNumber *) aa->data.num);
    decNumberCopy(&b, (decNumber *) bb->data.num);
    decNumberCopy(&c, (decNumber *) cc->data.num);
    if (decNumberIsNegative(&a)) {
        decNumberCopy(ret->data.num, &zero_crossing);
        goto RETURN;
    }
    if (! decNumberIsNegative(&c)) {
        if (! decNumberIsNegative(&b)) {
            if (decNumberIsPositive(&c)) {
                decNumberCopy(ret->data.num, &no_crossing);
            } else if (decNumberIsZero(&a) && decNumberIsZero(&b)) {
                decNumberCopy(ret->data.num, &no_crossing);
            } else {
                decNumberCopy(ret->data.num, &one_crossing);
            }
            goto RETURN;
        }
        if (decNumberIsZero(&a)) {
            decNumberCopy(ret->data.num, &zero_crossing);
            goto RETURN;
        }
    } else if (decNumberIsZero(&a) && ! decNumberIsPositive(&b)) {
        decNumberCopy(ret->data.num, &zero_crossing);
        goto RETURN;
    }
    /* Use bisection to find the crossing point... */
    d = epsilonf;
    decNumberCopy(&x0, &a);
    decNumberSubtract(&x1, &a, &b, &mp_decimal_data.set);
    decNumberSubtract(&x2, &b, &c, &mp_decimal_data.set);
    /* not sure why the error correction has to be >= 1E-12 */
    decNumberFromDouble(&scratch2, 1E-12);
    do {
        decNumberAdd(&x, &x1, &x2, &mp_decimal_data.set);
        decNumberDivide(&x, &x, &mp_decimal_data.two_decNumber, &mp_decimal_data.set);
        decNumberAdd(&x, &x, &scratch2, &mp_decimal_data.set);
        decNumberSubtract(&scratch, &x1, &x0, &mp_decimal_data.set);
        if (decNumberGreater(&scratch, &x0)) {
            decNumberCopy(&x2, &x);
            decNumberAdd(&x0, &x0, &x0, &mp_decimal_data.set);
            d += d;
        } else {
            decNumberAdd(&xx, &scratch, &x, &mp_decimal_data.set);
            if (decNumberGreater(&xx,&x0)) {
                decNumberCopy(&x2,&x);
                decNumberAdd(&x0, &x0, &x0, &mp_decimal_data.set);
                d += d;
            } else {
                decNumberSubtract(&x0, &x0, &xx, &mp_decimal_data.set);
                if (! decNumberGreater(&x,&x0)) {
                    decNumberAdd(&scratch, &x, &x2, &mp_decimal_data.set);
                    if (! decNumberGreater(&scratch, &x0)) {
                        decNumberCopy(ret->data.num, &no_crossing);
                        goto RETURN;
                    }
                }
                decNumberCopy(&x1,&x);
                d = d + d + epsilonf;
            }
        }
    } while (d < fraction_one);
    decNumberFromDouble(&scratch, d);
    decNumberSubtract(ret->data.num,&scratch, &mp_decimal_data.fraction_one_decNumber, &mp_decimal_data.set);
  RETURN:
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static mp_scaled_t mp_decimal_round_unscaled(mp_number *x_orig)
{
    return mpscaledround(mp_decimal_to_double(x_orig));
}

static void mp_decimal_floor(mp_number *i)
{
    int round = mp_decimal_data.set.round;
    mp_decimal_data.set.round = DEC_ROUND_FLOOR;
    decNumberToIntegralValue(i->data.num, i->data.num, &mp_decimal_data.set);
    mp_decimal_data.set.round = round;
}

static void mp_decimal_fraction_to_round_scaled(mp_number *x_orig)
{
    x_orig->type = mp_scaled_type;
    decNumberDivide(x_orig->data.num, x_orig->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set);
}

static void mp_decimal_square_rt(MP mp, mp_number *ret, mp_number *x_orig)
{
    decNumber x;
    decNumberCopy(&x, x_orig->data.num);
    if (! decNumberIsPositive(&x)) {
        if (decNumberIsNegative(&x)) {
            char msg[256];
            char *xstr = mp_decimal_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."
            );
        }
        decNumberZero(ret->data.num);
    } else {
        decNumberSquareRoot(ret->data.num, &x, &mp_decimal_data.set);
    }
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_decimal_pyth_add(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
{
    decNumber a, b;
    decNumber asq, bsq;
    decNumberCopyAbs(&a, a_orig->data.num);
    decNumberCopyAbs(&b, b_orig->data.num);
    decNumberMultiply(&asq, &a, &a, &mp_decimal_data.set);
    decNumberMultiply(&bsq, &b, &b, &mp_decimal_data.set);
    decNumberAdd(&a, &asq, &bsq, &mp_decimal_data.set);
    decNumberSquareRoot(ret->data.num, &a, &mp_decimal_data.set);
    /*
    if (set.status != 0) {
      mp->arithmic_error = 1;
      decNumberCopy(ret->data.num, &mp_decimal_data.EL_GORDO_decNumber);
    }
    */
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_decimal_pyth_sub(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
{
    decNumber a, b;
    decNumberCopyAbs(&a, a_orig->data.num);
    decNumberCopyAbs(&b, b_orig->data.num);
    if (! decNumberGreater(&a, &b)) {
        if (decNumberLess(&a, &b)) {
            char msg[256];
            char *astr = mp_decimal_number_tostring(mp, a_orig);
            char *bstr = mp_decimal_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."
            );
        }
        decNumberZero(&a);
    } else {
        decNumber asq, bsq;
        decNumberMultiply(&asq, &a, &a, &mp_decimal_data.set);
        decNumberMultiply(&bsq, &b, &b, &mp_decimal_data.set);
        decNumberSubtract(&a, &asq, &bsq, &mp_decimal_data.set);
        decNumberSquareRoot(&a, &a, &mp_decimal_data.set);
    }
    decNumberCopy(ret->data.num, &a);
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_decimal_power_of(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
{
    decNumberPower(ret->data.num, a_orig->data.num, b_orig->data.num, &mp_decimal_data.set);
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_decimal_m_log(MP mp, mp_number *ret, mp_number *x_orig)
{
    if (! decNumberIsPositive((decNumber *) x_orig->data.num)) {
        char msg[256];
        char *xstr = mp_decimal_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."
        );
        decNumberZero(ret->data.num);
    } else {
        decNumber twofivesix;
        decNumberFromScaled(&twofivesix, 256);
        decNumberLn(ret->data.num, x_orig->data.num, &mp_decimal_data.limitedset);
        mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.limitedset);
        decNumberMultiply(ret->data.num, ret->data.num, &twofivesix, &mp_decimal_data.set);
    }
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_decimal_m_exp(MP mp, mp_number *ret, mp_number *x_orig)
{
    decNumber temp, twofivesix;
    decNumberFromScaled(&twofivesix, 256);
    decNumberDivide(&temp, x_orig->data.num, &twofivesix, &mp_decimal_data.set);
    mp_decimal_data.limitedset.status = 0;
    decNumberExp(ret->data.num, &temp, &mp_decimal_data.limitedset);
    if (mp_decimal_data.limitedset.status & DEC_Clamped) {
        if (decNumberIsPositive((decNumber *) x_orig->data.num)) {
            mp->arithmic_error = 1;
            decNumberCopy(ret->data.num, &mp_decimal_data.EL_GORDO_decNumber);
        } else {
            decNumberZero(ret->data.num);
        }
    }
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.limitedset);
    mp_decimal_data.limitedset.status = 0;
}

static void mp_decimal_n_arg(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
{
    if (decNumberIsZero((decNumber *) x_orig->data.num) && decNumberIsZero((decNumber *) y_orig->data.num)) {
        if (decNumberIsNegative((decNumber *) internal_value(mp_default_zero_angle_internal).data.num)) {
            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."
            );
            decNumberZero(ret->data.num);
        } else { 
            decNumberCopy(ret->data.num, internal_value(mp_default_zero_angle_internal).data.num);
        }
    } else {
        decNumber atan2val, oneeighty_angle;
        ret->type = mp_angle_type;
        decNumberFromScaled(&oneeighty_angle, 180 * angle_multiplier);
        decNumberDivide(&oneeighty_angle, &oneeighty_angle, &mp_decimal_data.PI_decNumber, &mp_decimal_data.set);
        checkZero(y_orig->data.num);
        checkZero(x_orig->data.num);
        decNumberAtan2(&atan2val, y_orig->data.num, x_orig->data.num, &mp_decimal_data.set);
        decNumberMultiply(ret->data.num,&atan2val, &oneeighty_angle, &mp_decimal_data.set);
        checkZero(ret->data.num);
    }
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void sinecosine(decNumber *theangle, decNumber *c, decNumber *s)
{
    int prec = mp_decimal_data.set.digits/2;
    decNumber p, pxa, fac, cc;
    decNumber n1, n2, p1;
    decNumberZero(c);
    decNumberZero(s);
    if (prec < DECPRECISION_DEFAULT) {
        prec = DECPRECISION_DEFAULT;
    }
    for (int n = 0; n < prec; n++) {
        decNumberFromScaled(&p1, n);
        decNumberFromScaled(&n1, 2*n);
        decNumberPower(&p,  &mp_decimal_data.minusone, &p1, &mp_decimal_data.limitedset);
        if (n == 0) {
            decNumberCopy(&pxa, &mp_decimal_data.one);
        } else {
            decNumberPower(&pxa, theangle, &n1, &mp_decimal_data.limitedset);
        }
        if (2*n < mp_decimal_data.last_cached_factorial) {
            decNumberCopy(&fac,mp_decimal_data.factorials[2*n]);
        } else {
            decNumberCopy(&fac,mp_decimal_data.factorials[mp_decimal_data.last_cached_factorial]);
            for (int i = mp_decimal_data.last_cached_factorial+1; i <= 2*n; i++) {
                decNumberFromScaled(&cc, i);
                decNumberMultiply (&fac, &fac, &cc, &mp_decimal_data.set);
                if (i < FACTORIALS_CACHESIZE) {
                    mp_decimal_data.factorials[i] = mp_memory_allocate(sizeof(decNumber));
                    decNumberCopy(mp_decimal_data.factorials[i], &fac);
                    mp_decimal_data.last_cached_factorial = i;
                }
            }
        }
        decNumberDivide(&pxa, &pxa, &fac, &mp_decimal_data.set);
        decNumberMultiply(&pxa, &pxa, &p, &mp_decimal_data.set);
        decNumberAdd(s, s, &pxa, &mp_decimal_data.set);
        decNumberFromScaled(&n2, 2*n+1);
        decNumberMultiply(&fac, &fac, &n2, &mp_decimal_data.set); /* fac = fac * (2*n+1)*/
        decNumberPower(&pxa, theangle, &n2, &mp_decimal_data.limitedset);
        decNumberDivide(&pxa, &pxa, &fac, &mp_decimal_data.set);
        decNumberMultiply(&pxa, &pxa, &p, &mp_decimal_data.set);
        decNumberAdd(c, c, &pxa, &mp_decimal_data.set);
    }
}

static void mp_decimal_sin_cos(MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin)
{
    decNumber rad;
    decNumber one_eighty;
    double tmp = mp_decimal_to_double(z_orig)/16.0;
    if ((tmp == 90.0)||(tmp == -270)){
        decNumberZero(n_cos->data.num);
        decNumberCopy(n_sin->data.num, &mp_decimal_data.fraction_multiplier_decNumber);
    } else if ((tmp == -90.0)||(tmp == 270.0)) {
        decNumberZero(n_cos->data.num);
        decNumberCopyNegate(n_sin->data.num, &mp_decimal_data.fraction_multiplier_decNumber);
    } else if ((tmp == 180.0) || (tmp == -180.0)) {
        decNumberCopyNegate(n_cos->data.num, &mp_decimal_data.fraction_multiplier_decNumber);
        decNumberZero(n_sin->data.num);
    } else {
        decNumberFromScaled(&one_eighty, 180 * 16);
        decNumberMultiply(&rad, z_orig->data.num, &mp_decimal_data.PI_decNumber, &mp_decimal_data.set);
        decNumberDivide(&rad, &rad, &one_eighty, &mp_decimal_data.set);
        sinecosine(&rad, n_sin->data.num, n_cos->data.num);
        decNumberMultiply(n_cos->data.num, n_cos->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set);
        decNumberMultiply(n_sin->data.num, n_sin->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set);
    }
    mp_decnumber_check(mp, n_cos->data.num, &mp_decimal_data.set);
    mp_decnumber_check(mp, n_sin->data.num, &mp_decimal_data.set);
}

# 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 QUALITY       1009               /* recommended quality level for high-res use */
# define TT            70                 /* guaranteed separation between streams */
# define is_odd(x)     ((x)&1)            /* units bit of x */

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

static mp_decimal_random_info mp_decimal_random_data = {
    .dummy   = -1,
    .started = -1,
    .ptr     = &mp_decimal_random_data.dummy
};

static void ran_array(long aa[],int n)
{
    int i, j;
    for (j = 0; j < KK;j++) {
        aa[j] = mp_decimal_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_decimal_random_data.x[i] = mod_diff(aa[j - KK], aa[j - LL]);
    }
    for (;i < KK; i++, j++) {
        mp_decimal_random_data.x[i] = mod_diff(aa[j - KK], mp_decimal_random_data.x[i - LL]);
    }
}

static void 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;
        ss <<= 1;
        if (ss >= MM) {
            /* cyclic shift 29 bits */
            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_decimal_random_data.x[j + KK -LL] = x[j];
    }
    for (; j < KK; j++) {
        mp_decimal_random_data.x[j - LL] = x[j];
    }
    for (j = 0; j < 10; j++) {
        /* warm things up */
        ran_array(x, KK + KK - 1);
    }
    mp_decimal_random_data.ptr = &mp_decimal_random_data.started;
}

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

static void mp_decimal_init_randoms(MP mp, int seed)
{
    int k = 1; /* more or less random integers */
    int j = abs(seed);
    while (j >= fraction_one) {
        j = j/2;
    }
    for (int i = 0; i <= 54; i++) {
        int jj = k;
        k = j - k;
        j = jj;
        if (k < 0) {
            k += fraction_one;
        }
        decNumberFromScaled(mp->randoms[(i * 21) % 55].data.num, j);
    }
    /* \quote {warm up} the array */
    mp_new_randoms(mp);
    mp_new_randoms(mp);
    mp_new_randoms(mp);
    ran_start((unsigned long) seed);
}

static void mp_decimal_number_modulo(mp_number *a, mp_number *b)
{
    decNumberRemainder(a->data.num, a->data.num, b->data.num, &mp_decimal_data.set);
}

static void mp_next_unif_random(MP mp, mp_number *ret)
{
    decNumber a;
    decNumber b;
    unsigned long int op = (unsigned) (*mp_decimal_random_data.ptr>=0? *mp_decimal_random_data.ptr++: ran_arr_cycle());
    (void) mp;
    decNumberFromScaled(&a, op);
    decNumberFromScaled(&b, MM);
    decNumberDivide(&a, &a, &b, &mp_decimal_data.set); /* a = a/b */
    decNumberCopy(ret->data.num, &a);
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_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_decimal_clone(ret, &(mp->randoms[mp->j_random]));
}

static void mp_decimal_m_unif_rand(MP mp, mp_number *ret, mp_number *x_orig)
{
    mp_number x, abs_x, u, y; /* |y| is trial value */
    mp_decimal_allocate_number(mp, &y, mp_fraction_type);
    mp_decimal_allocate_clone(mp, &x, mp_scaled_type, x_orig);
    mp_decimal_allocate_abs(mp, &abs_x, mp_scaled_type, &x);
    mp_decimal_allocate_number(mp, &u, mp_scaled_type);
    mp_next_unif_random(mp, &u);
    decNumberMultiply(y.data.num, abs_x.data.num, u.data.num, &mp_decimal_data.set);
    if (mp_decimal_equal(&y, &abs_x)) {
        mp_decimal_clone(ret, &((math_data *)mp->math)->md_zero_t);
    } else if (mp_decimal_greater(&x, &((math_data *)mp->math)->md_zero_t)) {
        mp_decimal_clone(ret, &y);
    } else {
        mp_decimal_negated_clone(ret, &y);
    }
    mp_decimal_free_number(mp, &x);
    mp_decimal_free_number(mp, &abs_x);
    mp_decimal_free_number(mp, &y);
    mp_decimal_free_number(mp, &u);
}

static void mp_decimal_m_norm_rand(MP mp, mp_number *ret)
{
    mp_number abs_x, u, r, la, xa;
    mp_decimal_allocate_number(mp, &la, mp_scaled_type);
    mp_decimal_allocate_number(mp, &xa, mp_scaled_type);
    mp_decimal_allocate_number(mp, &abs_x, mp_scaled_type);
    mp_decimal_allocate_number(mp, &u, mp_scaled_type);
    mp_decimal_allocate_number(mp, &r, mp_scaled_type);
    do {
        do {
            mp_number v; /* maybe move outside loop */
            mp_decimal_allocate_number(mp, &v, mp_scaled_type);
            mp_next_random(mp, &v);
            mp_decimal_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t);
            mp_decimal_number_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v);
            mp_decimal_free_number(mp, &v);
            mp_next_random(mp, &u);
            mp_decimal_clone(&abs_x, &xa);
            mp_decimal_abs(&abs_x);
        } while (! mp_decimal_less(&abs_x, &u));
        mp_decimal_number_make_fraction(mp, &r, &xa, &u);
        mp_decimal_clone(&xa, &r);
        mp_decimal_m_log(mp, &la, &u);
        mp_decimal_set_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la);
    } while (mp_decimal_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0);
    mp_decimal_clone(ret, &xa);
    mp_decimal_free_number(mp, &r);
    mp_decimal_free_number(mp, &abs_x);
    mp_decimal_free_number(mp, &la);
    mp_decimal_free_number(mp, &xa);
    mp_decimal_free_number(mp, &u);
}

static void mp_decimal_set_precision(MP mp)
{
    int i = decNumberToInt32((decNumber *) internal_value(mp_number_precision_internal).data.num, &mp_decimal_data.set);
    mp_decimal_data.set.digits = i;
    mp_decimal_data.limitedset.digits = i;
}

static void mp_decimal_free_math(MP mp)
{
    mp_decimal_free_number(mp, &(mp->math->md_three_sixty_deg_t));
    mp_decimal_free_number(mp, &(mp->math->md_one_eighty_deg_t));
    mp_decimal_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t));
    mp_decimal_free_number(mp, &(mp->math->md_fraction_one_t));
    mp_decimal_free_number(mp, &(mp->math->md_zero_t));
    mp_decimal_free_number(mp, &(mp->math->md_half_unit_t));
    mp_decimal_free_number(mp, &(mp->math->md_three_quarter_unit_t));
    mp_decimal_free_number(mp, &(mp->math->md_unity_t));
    mp_decimal_free_number(mp, &(mp->math->md_two_t));
    mp_decimal_free_number(mp, &(mp->math->md_three_t));
    mp_decimal_free_number(mp, &(mp->math->md_one_third_inf_t));
    mp_decimal_free_number(mp, &(mp->math->md_inf_t));
    mp_decimal_free_number(mp, &(mp->math->md_negative_inf_t));
    mp_decimal_free_number(mp, &(mp->math->md_warning_limit_t));
    mp_decimal_free_number(mp, &(mp->math->md_one_k));
    mp_decimal_free_number(mp, &(mp->math->md_sqrt_8_e_k));
    mp_decimal_free_number(mp, &(mp->math->md_twelve_ln_2_k));
    mp_decimal_free_number(mp, &(mp->math->md_coef_bound_k));
    mp_decimal_free_number(mp, &(mp->math->md_coef_bound_minus_1));
    mp_decimal_free_number(mp, &(mp->math->md_fraction_threshold_t));
    mp_decimal_free_number(mp, &(mp->math->md_half_fraction_threshold_t));
    mp_decimal_free_number(mp, &(mp->math->md_scaled_threshold_t));
    mp_decimal_free_number(mp, &(mp->math->md_half_scaled_threshold_t));
    mp_decimal_free_number(mp, &(mp->math->md_near_zero_angle_t));
    mp_decimal_free_number(mp, &(mp->math->md_p_over_v_threshold_t));
    mp_decimal_free_number(mp, &(mp->math->md_equation_threshold_t));
    /*
        For sake of speed, we accept this memory leak:

        for (i = 0; i <= mp_decimal_data.last_cached_factorial; i++) {
            mp_memory_free(mp_decimal_data.factorials[i]);
        }
        mp_memory_free(mp_decimal_data.factorials);
    */
    mp_memory_free(mp->math);
}

math_data *mp_initialize_decimal_math(MP mp)
{
    math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data));
    decContextDefault(&mp_decimal_data.set, DEC_INIT_BASE);        /* initialize */
    mp_decimal_data.set.traps = 0;                                 /* no traps, thank you */
    decContextDefault(&mp_decimal_data.limitedset, DEC_INIT_BASE); /* initialize */
    mp_decimal_data.limitedset.traps = 0;                          /* no traps, thank you */
    mp_decimal_data.limitedset.emax = 999999;
    mp_decimal_data.limitedset.emin = -999999;
    mp_decimal_data.set.digits = DECPRECISION_DEFAULT;
    mp_decimal_data.limitedset.digits = DECPRECISION_DEFAULT;
    if (! mp_decimal_data.initialized) {
        mp_decimal_data.initialized = 1 ;
        decNumberFromScaled(&mp_decimal_data.one, 1);
        decNumberFromScaled(&mp_decimal_data.minusone, -1);
        decNumberFromScaled(&mp_decimal_data.zero, 0);
        decNumberFromScaled(&mp_decimal_data.two_decNumber, two);
        decNumberFromScaled(&mp_decimal_data.three_decNumber, three);
        decNumberFromScaled(&mp_decimal_data.four_decNumber, four);
        decNumberFromScaled(&mp_decimal_data.fraction_multiplier_decNumber, fraction_multiplier);
        decNumberFromScaled(&mp_decimal_data.fraction_one_decNumber, fraction_one);
        decNumberFromScaled(&mp_decimal_data.fraction_one_plus_decNumber, (fraction_one+1));
        decNumberFromScaled(&mp_decimal_data.angle_multiplier_decNumber, angle_multiplier);
        decNumberFromString(&mp_decimal_data.PI_decNumber, PI_STRING, &mp_decimal_data.set);
        decNumberFromDouble(&mp_decimal_data.epsilon_decNumber, epsilon);
        decNumberFromString(&mp_decimal_data.EL_GORDO_decNumber, EL_GORDO, &mp_decimal_data.set);
        decNumberFromString(&mp_decimal_data.negative_EL_GORDO_decNumber, negative_EL_GORDO, &mp_decimal_data.set);
        mp_decimal_data.factorials = (decNumber **) mp_memory_allocate(FACTORIALS_CACHESIZE * sizeof(decNumber *));
        mp_decimal_data.factorials[0] = (decNumber *) mp_memory_allocate(sizeof(decNumber));
        decNumberCopy(mp_decimal_data.factorials[0], &mp_decimal_data.one);
    }
    math->md_allocate        = mp_decimal_allocate_number;
    math->md_free            = mp_decimal_free_number;
    math->md_allocate_clone  = mp_decimal_allocate_clone;
    math->md_allocate_abs    = mp_decimal_allocate_abs;
    math->md_allocate_div    = mp_decimal_allocate_div;
    math->md_allocate_mul    = mp_decimal_allocate_mul;
    math->md_allocate_add    = mp_decimal_allocate_add;
    math->md_allocate_sub    = mp_decimal_allocate_sub;
    math->md_allocate_double = mp_decimal_allocate_double;
    mp_decimal_allocate_number(mp, &math->md_precision_default, mp_scaled_type);
    decNumberFromScaled(math->md_precision_default.data.num, DECPRECISION_DEFAULT);
    mp_decimal_allocate_number(mp, &math->md_precision_max, mp_scaled_type);
    decNumberFromScaled(math->md_precision_max.data.num, DECNUMDIGITS);
    mp_decimal_allocate_number(mp, &math->md_precision_min, mp_scaled_type);
    decNumberFromScaled(math->md_precision_min.data.num, 2);
    /* here are the constants for scaled objects */
    mp_decimal_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type);
    decNumberCopy(math->md_epsilon_t.data.num, &mp_decimal_data.epsilon_decNumber);
    mp_decimal_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
    decNumberCopy(math->md_inf_t.data.num, &mp_decimal_data.EL_GORDO_decNumber);
    mp_decimal_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type);
    decNumberCopy(math->md_negative_inf_t.data.num, &mp_decimal_data.negative_EL_GORDO_decNumber);
    mp_decimal_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
    decNumberFromString(math->md_warning_limit_t.data.num, warning_limit, &mp_decimal_data.set);
    mp_decimal_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type);
    decNumberDivide(math->md_one_third_inf_t.data.num, math->md_inf_t.data.num, &mp_decimal_data.three_decNumber, &mp_decimal_data.set);
    mp_decimal_allocate_number(mp, &math->md_unity_t, mp_scaled_type);
    decNumberCopy(math->md_unity_t.data.num, &mp_decimal_data.one);
    mp_decimal_allocate_number(mp, &math->md_two_t, mp_scaled_type);
    decNumberFromScaled(math->md_two_t.data.num, two);
    mp_decimal_allocate_number(mp, &math->md_three_t, mp_scaled_type);
    decNumberFromScaled(math->md_three_t.data.num, three);
    mp_decimal_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type);
    decNumberFromString(math->md_half_unit_t.data.num, "0.5", &mp_decimal_data.set);
    mp_decimal_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type);
    decNumberFromString(math->md_three_quarter_unit_t.data.num, "0.75", &mp_decimal_data.set);
    mp_decimal_allocate_number(mp, &math->md_zero_t, mp_scaled_type);
    decNumberZero(math->md_zero_t.data.num);
    /* fractions */
    {
        decNumber fourzeroninesix;
        decNumberFromScaled(&fourzeroninesix, 4096);
        mp_decimal_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type);
        decNumberDivide(math->md_arc_tol_k.data.num, &mp_decimal_data.one, &fourzeroninesix, &mp_decimal_data.set);         /* quit when change in arc length estimate reaches this */
    }
    mp_decimal_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type);
    decNumberFromScaled(math->md_fraction_one_t.data.num, fraction_one);
    mp_decimal_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
    decNumberFromScaled(math->md_fraction_half_t.data.num, fraction_half);
    mp_decimal_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
    decNumberFromScaled(math->md_fraction_three_t.data.num, fraction_three);
    mp_decimal_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
    decNumberFromScaled(math->md_fraction_four_t.data.num, fraction_four);
    /* angles */
    mp_decimal_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type);
    decNumberFromScaled(math->md_three_sixty_deg_t.data.num, 360  * angle_multiplier);
    mp_decimal_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type);
    decNumberFromScaled(math->md_one_eighty_deg_t.data.num, 180 * angle_multiplier);
    mp_decimal_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type);
    decNumberFromScaled(math->md_negative_one_eighty_deg_t.data.num, -180 * angle_multiplier);
    /* various approximations */
    mp_decimal_allocate_number(mp, &math->md_one_k, mp_scaled_type);
    decNumberFromDouble(math->md_one_k.data.num, 1.0/64);
    mp_decimal_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type);
    decNumberFromDouble(math->md_sqrt_8_e_k.data.num, 112428.82793 / 65536.0);               /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
    mp_decimal_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type);
    decNumberFromDouble(math->md_twelve_ln_2_k.data.num, 139548959.6165 / 65536.0);          /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
    mp_decimal_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type);
    decNumberFromDouble(math->md_coef_bound_k.data.num,coef_bound);
    mp_decimal_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type);
    decNumberFromDouble(math->md_coef_bound_minus_1.data.num,coef_bound - 1 / 65536.0);
    mp_decimal_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type);
    decNumberFromDouble(math->md_twelvebits_3.data.num, 1365 / 65536.0);                     /* $1365\approx 2^{12}/3$ */
    mp_decimal_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type);
    decNumberFromDouble(math->md_twentysixbits_sqrt2_t.data.num, 94906265.62 / 65536.0);     /* $2^{26}\sqrt2\approx94906265.62$ */
    mp_decimal_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type);
    decNumberFromDouble(math->md_twentyeightbits_d_t.data.num, 35596754.69 / 65536.0);       /* $2^{28}d\approx35596754.69$ */
    mp_decimal_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type);
    decNumberFromDouble(math->md_twentysevenbits_sqrt2_d_t.data.num, 25170706.63 / 65536.0); /* $2^{27}\sqrt2\,d\approx25170706.63$ */
    /* thresholds */
    mp_decimal_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type);
    decNumberFromDouble(math->md_fraction_threshold_t.data.num, fraction_threshold);
    mp_decimal_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
    decNumberFromDouble(math->md_half_fraction_threshold_t.data.num, half_fraction_threshold);
    mp_decimal_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
    decNumberFromDouble(math->md_scaled_threshold_t.data.num, scaled_threshold);
    mp_decimal_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
    decNumberFromDouble(math->md_half_scaled_threshold_t.data.num, half_scaled_threshold);
    mp_decimal_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
    decNumberFromDouble(math->md_near_zero_angle_t.data.num, near_zero_angle);
    mp_decimal_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type);
    decNumberFromDouble(math->md_p_over_v_threshold_t.data.num, p_over_v_threshold);
    mp_decimal_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
    decNumberFromDouble(math->md_equation_threshold_t.data.num, equation_threshold);
    /* functions */
    math->md_from_int                 = mp_decimal_set_from_int;
    math->md_from_boolean             = mp_decimal_set_from_boolean;
    math->md_from_scaled              = mp_decimal_set_from_scaled;
    math->md_from_double              = mp_decimal_set_from_double;
    math->md_from_addition            = mp_decimal_set_from_addition;
    math->md_half_from_addition       = mp_decimal_set_half_from_addition;
    math->md_from_subtraction         = mp_decimal_set_from_subtraction;
    math->md_half_from_subtraction    = mp_decimal_set_half_from_subtraction;
    math->md_from_of_the_way          = mp_decimal_set_from_of_the_way;
    math->md_from_div                 = mp_decimal_set_from_div;
    math->md_from_mul                 = mp_decimal_set_from_mul;
    math->md_from_int_div             = mp_decimal_set_from_int_div;
    math->md_from_int_mul             = mp_decimal_set_from_int_mul;
    math->md_negate                   = mp_decimal_negate;
    math->md_add                      = mp_decimal_add;
    math->md_subtract                 = mp_decimal_subtract;
    math->md_half                     = mp_decimal_half;
    math->md_do_double                = mp_decimal_double;
    math->md_abs                      = mp_decimal_abs;
    math->md_clone                    = mp_decimal_clone;
    math->md_negated_clone            = mp_decimal_negated_clone;
    math->md_abs_clone                = mp_decimal_abs_clone;
    math->md_swap                     = mp_decimal_swap;
    math->md_add_scaled               = mp_decimal_add_scaled;
    math->md_multiply_int             = mp_decimal_multiply_int;
    math->md_divide_int               = mp_decimal_divide_int;
    math->md_to_int                   = mp_decimal_to_int;
    math->md_to_boolean               = mp_decimal_to_boolean;
    math->md_to_scaled                = mp_decimal_to_scaled;
    math->md_to_double                = mp_decimal_to_double;
    math->md_odd                      = mp_decimal_odd;
    math->md_equal                    = mp_decimal_equal;
    math->md_less                     = mp_decimal_less;
    math->md_greater                  = mp_decimal_greater;
    math->md_non_equal_abs            = mp_decimal_non_equal_abs;
    math->md_round_unscaled           = mp_decimal_round_unscaled;
    math->md_floor_scaled             = mp_decimal_floor;
    math->md_fraction_to_round_scaled = mp_decimal_fraction_to_round_scaled;
    math->md_make_scaled              = mp_decimal_number_make_scaled;
    math->md_make_fraction            = mp_decimal_number_make_fraction;
    math->md_take_fraction            = mp_decimal_number_take_fraction;
    math->md_take_scaled              = mp_decimal_number_take_scaled;
    math->md_velocity                 = mp_decimal_velocity;
    math->md_n_arg                    = mp_decimal_n_arg;
    math->md_m_log                    = mp_decimal_m_log;
    math->md_m_exp                    = mp_decimal_m_exp;
    math->md_m_unif_rand              = mp_decimal_m_unif_rand;
    math->md_m_norm_rand              = mp_decimal_m_norm_rand;
    math->md_pyth_add                 = mp_decimal_pyth_add;
    math->md_pyth_sub                 = mp_decimal_pyth_sub;
    math->md_power_of                 = mp_decimal_power_of;
    math->md_fraction_to_scaled       = mp_decimal_fraction_to_scaled;
    math->md_scaled_to_fraction       = mp_decimal_scaled_to_fraction;
    math->md_scaled_to_angle          = mp_decimal_scaled_to_angle;
    math->md_angle_to_scaled          = mp_decimal_angle_to_scaled;
    math->md_init_randoms             = mp_decimal_init_randoms;
    math->md_sin_cos                  = mp_decimal_sin_cos;
    math->md_slow_add                 = mp_decimal_slow_add;
    math->md_sqrt                     = mp_decimal_square_rt;
    math->md_print                    = mp_decimal_print_number;
    math->md_tostring                 = mp_decimal_number_tostring;
    math->md_modulo                   = mp_decimal_number_modulo;
    math->md_ab_vs_cd                 = mp_decimal_ab_vs_cd;
    math->md_crossing_point           = mp_decimal_crossing_point;
    math->md_scan_numeric             = mp_decimal_scan_numeric_token;
    math->md_scan_fractional          = mp_decimal_scan_fractional_token;
    /* housekeeping */
    math->md_free_math                = mp_decimal_free_math;
    math->md_set_precision            = mp_decimal_set_precision;
    return math;
}