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

/*
    We cannot predefine (share) as much as we like because some numbers (constants) depend
    on the precission. Of course sharing doesn't really save much, it is more about
    abstraction and the (future) posibility to fetch some of these properties from the
    engine (as we do with \TEX).

*/

# include "mpmathdecimal.h"

# define DECNUMDIGITS 1000

# include "decNumber.h"

# define E_STRING                 "2.7182818284590452353602874713526624977572470936999595749669676277240766303535"
# define PI_STRING                "3.1415926535897932384626433832795028841971693993751058209749445923078164062862"
# define EL_GORDO_STRING          "1E1000000"
# define NEGATIVE_EL_GORDO_STRING "-1E1000000"
# define WARNING_LIMIT_STRING     "1E1000000"

# define mp_fraction_multiplier   4096
# define mp_angle_multiplier      16

# define mp_epsilon               pow(2.0,-173.0)
# define mp_epsilonf              pow(2.0,-52.0)

/* */

# 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 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)))

/*tex
    We have some more \quote {constants} that can go here so that we can use copies
    instead. I will do that when we have some performance critical example. In the end
    this only makes sense for number systems that stored numbers in allocated blobs.
*/

typedef struct mp_decimal_info {
    decContext   set;
    decContext   limitedset;
    decNumber    zero;
    decNumber    one;
    decNumber    minusone;
    decNumber    two;
    decNumber    three;
    decNumber    four;
    decNumber    five;
    decNumber    sixteen;
    decNumber    fraction_multiplier;
    decNumber    angle_multiplier;
    decNumber    fraction_one;
    decNumber    fraction_two;
    decNumber    fraction_three;
    decNumber    fraction_four;
    decNumber    fraction_one_plus;
    decNumber    fraction_half;
 // decNumber    dp90;
 // decNumber    dp180;
 // decNumber    dp270;
 // decNumber    dp360;
 // decNumber    dm90;
 // decNumber    dm180;
 // decNumber    dm270;
 // decNumber    dm360;
 // decNumber    d16;
 // decNumber    d64;
 // decNumber    d256;
 // decNumber    d4096;
 // decNumber    d65536;
    decNumber    pi;
    decNumber    epsilon;
    decNumber    EL_GORDO;
    decNumber    negative_EL_GORDO;
    decNumber    no_crossing;
    decNumber    one_crossing;
    decNumber    zero_crossing;
    /* */
    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, const decNumber *p, const decNumber *q);
static void mp_decimal_take_fraction (MP mp, decNumber *ret, const decNumber *p, const decNumber *q);

static inline 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);
            } else {
                decNumberCopy(dec, &mp_decimal_data.EL_GORDO);
            }
        } else {
            /* Nan */
            decNumberZero(dec);
        }
    }
    if (decNumberIsZero(dec) && decNumberIsNegative(dec)) {
        decNumberZero(dec);
    }
    mp->arithmic_error = test;
}

/* See mpmathdouble for documentation. */

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

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

static inline int decNumberGreater(const decNumber *a, const 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(const 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 = mp_error_code(mp, 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, const 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, 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, const decNumber *y, const 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, localset);
            } else {
                decNumberAdd(result, result, &mp_decimal_data.pi, 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, &mp_decimal_data.four, 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, &mp_decimal_data.two, localset);
        } else {
            /* Otherwise it is +0 if x is positive, +pi if x is neg */
            if (decNumberIsNegative(x)) {
                decNumberCopy(result, &mp_decimal_data.pi);
            } 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, const 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, const 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, const mp_number *a, const 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, const mp_number *a, const 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, const mp_number *a, const 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, const mp_number *a, const 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, (int32_t) B);
}

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

static void mp_decimal_set_from_scaled(mp_number *A, mp_scaled_t B)
{
    decNumber c;
    decNumberFromScaled(&c, 65536); /* can be copy */
    decNumberFromScaled(A->data.num, (int32_t) 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, const mp_number *B, const 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, const mp_number *B, const mp_number *C)
{
    decNumber c;
    decNumberAdd(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
    decNumberFromScaled(&c, 2); /* can be copy */
    decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_subtraction(mp_number *A, const mp_number *B, const 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, const mp_number *B, const mp_number *C)
{
    decNumber c;
    decNumberSubtract(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set);
    decNumberFromScaled(&c, 2); /* can be copy */
    decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_div(mp_number *A, const mp_number *B, const 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, const mp_number *B, const 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, const mp_number *B, mp_scaled_t C)
{
    decNumber c;
    decNumberFromScaled(&c, (int32_t) C);
    decNumberDivide(A->data.num, B->data.num, &c, &mp_decimal_data.set);
}

static void mp_decimal_set_from_int_mul(mp_number *A, const mp_number *B, mp_scaled_t C)
{
    decNumber c;
    decNumberFromScaled(&c, (int32_t) 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, const mp_number *t, const mp_number *B, const 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, const 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, const 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); /* can be copy */
    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); /* can be copy */
    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); /* can be copy */
    decNumberFromScaled(&b, (int32_t) 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, (int32_t) 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, (int32_t) 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 inline void mp_decimal_clone(mp_number *A, const mp_number *B)
{
    decNumberCopy(A->data.num, B->data.num);
}

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

static inline void mp_decimal_abs_clone(mp_number *A, const 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, &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, &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, &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, &mp_decimal_data.set);
}

static mp_scaled_t mp_decimal_to_scaled(const mp_number *A)
{
    mp_scaled_t result;
    decNumber corrected;
    decNumberFromScaled(&corrected, 65536); /* can be copy */
    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(const 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 = mp_error_code(mp, 2);| */
        return 0;
    } else {
        return result;
    }
}

static mp_scaled_t mp_decimal_to_boolean(const 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 = mp_error_code(mp, 3);| */
        return 0;
    } else {
        return result;
    }
}

static double mp_decimal_to_double(const 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 = mp_error_code(mp, 4);| */
        return 0.0;
    }
}

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

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

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

static inline int mp_decimal_less(const mp_number *A, const 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(const mp_number *A, const 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(const decNumber *n)
{
    decNumber corrected;
    char *buffer = mp_memory_allocate((size_t) n->digits + 14);
    decNumberCopy(&corrected, n);
    decNumberTrim(&corrected);
    decNumberToString(&corrected, buffer);
    return buffer;
}

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

static void mp_decimal_print_number(MP mp, const 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, const mp_number *A, const mp_number *B)
{
    (void) mp;
    decNumberAdd(ret->data.num, A->data.num, B->data.num, &mp_decimal_data.set);
}

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

static void mp_decimal_make_fraction(MP mp, decNumber *ret, const decNumber *p, const 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, &mp_decimal_data.set);
}

static void mp_decimal_number_make_fraction(MP mp, mp_number *ret, const mp_number *p, const 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, const decNumber *p, const decNumber *q)
{
    (void) mp;
    decNumberMultiply(ret, p, q, &mp_decimal_data.set);
    decNumberDivide(ret, ret, &mp_decimal_data.fraction_multiplier, &mp_decimal_data.set);
}

static void mp_decimal_number_take_fraction(MP mp, mp_number *ret, const mp_number *p, const 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, const mp_number *p_orig, const 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, const mp_number *p_orig, const 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, const unsigned char *start, const unsigned char *stop)
{
    decNumber result;
    size_t l = (size_t)(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)
{
    const unsigned char *start = &mp->buffer[mp->cur_input.loc_field - 1];
    const 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)
{
    const unsigned char *start = &mp->buffer[mp->cur_input.loc_field - 1];
    const 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, const mp_number *st, const mp_number *ct, const mp_number *sf, const mp_number *cf, const mp_number *t)
{
    decNumber acc, num, denom; /* registers for intermediate calculations */
    decNumber r1, r2;
    decNumber arg1, arg2;
    decNumber i16, fone, fhalf, ftwo, sqrtfive;
    decNumberCopy(&i16, &mp_decimal_data.sixteen);
    decNumberCopy(&fone, &mp_decimal_data.fraction_one);
    decNumberCopy(&fhalf, &mp_decimal_data.fraction_half);
    decNumberCopy(&ftwo, &mp_decimal_data.fraction_two);
    decNumberCopy(&sqrtfive, &mp_decimal_data.sixteen);
    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, &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, &mp_decimal_data.set);      /* arg1 = arg1 * 3*/

    decNumberSubtract(&arg2,&mp_decimal_data.three, &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, &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*/

    decNumberCopy(&denom, &mp_decimal_data.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, &mp_decimal_data.set);
    if (decNumberLess(&denom, &r2)) {
        decNumberCopy(ret->data.num, &mp_decimal_data.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 (const mp_number *a_orig, const mp_number *b_orig, const mp_number *c_orig, const mp_number *d_orig)
{
    decNumber a, b, c, d;
    decNumber ab, cd;
    decNumberCopy(&a, (const decNumber *) a_orig->data.num);
    decNumberCopy(&b, (const decNumber *) b_orig->data.num);
    decNumberCopy(&c, (const decNumber *) c_orig->data.num);
    decNumberCopy(&d, (const decNumber *) d_orig->data.num);
    decNumberMultiply(&ab, (const decNumber *) a_orig->data.num, (const decNumber *) b_orig->data.num, &mp_decimal_data.set);
    decNumberMultiply(&cd, (const decNumber *) c_orig->data.num, (const 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, const mp_number *aa, const mp_number *bb, const 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, (const decNumber *) aa->data.num);
    decNumberCopy(&b, (const decNumber *) bb->data.num);
    decNumberCopy(&c, (const decNumber *) cc->data.num);
    if (decNumberIsNegative(&a)) {
        decNumberCopy(ret->data.num, &mp_decimal_data.zero_crossing);
        goto RETURN;
    }
    if (! decNumberIsNegative(&c)) {
        if (! decNumberIsNegative(&b)) {
            if (decNumberIsPositive(&c)) {
                decNumberCopy(ret->data.num, &mp_decimal_data.no_crossing);
            } else if (decNumberIsZero(&a) && decNumberIsZero(&b)) {
                decNumberCopy(ret->data.num, &mp_decimal_data.no_crossing);
            } else {
                decNumberCopy(ret->data.num, &mp_decimal_data.one_crossing);
            }
            goto RETURN;
        }
        if (decNumberIsZero(&a)) {
            decNumberCopy(ret->data.num, &mp_decimal_data.zero_crossing);
            goto RETURN;
        }
    } else if (decNumberIsZero(&a) && ! decNumberIsPositive(&b)) {
        decNumberCopy(ret->data.num, &mp_decimal_data.zero_crossing);
        goto RETURN;
    }
    /* Use bisection to find the crossing point... */
    d = mp_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, &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, &mp_decimal_data.no_crossing);
                        goto RETURN;
                    }
                }
                decNumberCopy(&x1,&x);
                d = d + d + mp_epsilonf;
            }
        }
    } while (d < mp_fraction_multiplier);
    decNumberFromDouble(&scratch, d);
    decNumberSubtract(ret->data.num,&scratch, &mp_decimal_data.fraction_one, &mp_decimal_data.set);
  RETURN:
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static mp_scaled_t mp_decimal_round_unscaled(const 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, &mp_decimal_data.set);
}

static void mp_decimal_square_rt(MP mp, mp_number *ret, const 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, const mp_number *a_orig, const 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 = mp_error_code(mp, 5);
      decNumberCopy(ret->data.num, &mp_decimal_data.EL_GORDO);
    }
    */
    mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set);
}

static void mp_decimal_pyth_add3(MP mp, mp_number *ret, const mp_number *a_orig, const mp_number *b_orig, const mp_number *c_orig)
{
    decNumber a, b, c;
    decNumber asq, bsq, csq;
    decNumberCopyAbs(&a, a_orig->data.num);
    decNumberCopyAbs(&b, b_orig->data.num);
    decNumberCopyAbs(&c, c_orig->data.num);
    decNumberMultiply(&asq, &a, &a, &mp_decimal_data.set);
    decNumberMultiply(&bsq, &b, &b, &mp_decimal_data.set);
    decNumberMultiply(&csq, &c, &c, &mp_decimal_data.set);
    decNumberAdd(&a, &asq, &bsq, &mp_decimal_data.set);
    decNumberAdd(&b, &a, &csq, &mp_decimal_data.set);
    decNumberSquareRoot(ret->data.num, &b, &mp_decimal_data.set);
    /*
    if (set.status != 0) {
      mp->arithmic_error = mp_error_code(mp, 5);
      decNumberCopy(ret->data.num, &mp_decimal_data.EL_GORDO);
    }
    */
    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, const 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, const mp_number *a_orig, const 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, const mp_number *x_orig)
{
    if (! decNumberIsPositive((const 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); /* can be copy */
        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, const mp_number *x_orig)
{
    decNumber temp, twofivesix;
    decNumberFromScaled(&twofivesix, 256); /* can be copy */
    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((const decNumber *) x_orig->data.num)) {
            mp->arithmic_error = mp_error_code(mp, 6);
            decNumberCopy(ret->data.num, &mp_decimal_data.EL_GORDO);
        } 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, const mp_number *x_orig, const mp_number *y_orig)
{
    if (decNumberIsZero((const decNumber *) x_orig->data.num) && decNumberIsZero((const decNumber *) y_orig->data.num)) {
        if (decNumberIsNegative((const 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 * mp_angle_multiplier); /* can be copy */
        decNumberDivide(&oneeighty_angle, &oneeighty_angle, &mp_decimal_data.pi, &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(const 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, const 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);
    } else if ((tmp == -90.0)||(tmp == 270.0)) {
        decNumberZero(n_cos->data.num);
        decNumberCopyNegate(n_sin->data.num, &mp_decimal_data.fraction_multiplier);
    } else if ((tmp == 180.0) || (tmp == -180.0)) {
        decNumberCopyNegate(n_cos->data.num, &mp_decimal_data.fraction_multiplier);
        decNumberZero(n_sin->data.num);
    } else {
        decNumberFromScaled(&one_eighty, 180 * 16); /* can be copy */
        decNumberMultiply(&rad, z_orig->data.num, &mp_decimal_data.pi, &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, &mp_decimal_data.set);
        decNumberMultiply(n_sin->data.num, n_sin->data.num, &mp_decimal_data.fraction_multiplier, &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);
}

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

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

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

static void mp_decimal_init_randoms(MP mp, int seed)
{
    mp->random_data = (mp_random_info) {
        .dummy   = -1,
        .started = -1,
        .ptr     = &mp->random_data.dummy
    };
    int k = 1; /* more or less random integers */
    int j = abs(seed);
    while (j >= mp_fraction_multiplier) {
        j = j/2;
    }
    for (int i = 0; i <= 54; i++) {
        int jj = k;
        k = j - k;
        j = jj;
        if (k < 0) {
            k += mp_fraction_multiplier;
        }
        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(mp, (unsigned long) seed);
}

static void mp_decimal_number_modulo(mp_number *a, const 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->random_data.ptr >= 0 ? *mp->random_data.ptr++ : ran_arr_cycle(mp));
    (void) mp;
    decNumberFromScaled(&a, op);
    decNumberFromScaled(&b, mp_random_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, const 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((const 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) {
        /* these are shared */
        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, 2);
        decNumberFromScaled(&mp_decimal_data.three, 3);
        decNumberFromScaled(&mp_decimal_data.four, 4);
        decNumberFromScaled(&mp_decimal_data.five, 5);
        decNumberFromScaled(&mp_decimal_data.sixteen, 16);
        decNumberFromScaled(&mp_decimal_data.fraction_multiplier, mp_fraction_multiplier);
        decNumberFromScaled(&mp_decimal_data.fraction_one, mp_fraction_multiplier);
        decNumberFromScaled(&mp_decimal_data.fraction_two, 2 * mp_fraction_multiplier);
        decNumberFromScaled(&mp_decimal_data.fraction_three, 3 * mp_fraction_multiplier);
        decNumberFromScaled(&mp_decimal_data.fraction_four, 4 * mp_fraction_multiplier);
        decNumberFromScaled(&mp_decimal_data.fraction_half, mp_fraction_multiplier / 2);
        decNumberFromScaled(&mp_decimal_data.fraction_one_plus, mp_fraction_multiplier + 1);
        decNumberFromScaled(&mp_decimal_data.angle_multiplier, mp_angle_multiplier);
        decNumberFromString(&mp_decimal_data.pi, PI_STRING, &mp_decimal_data.set);
        decNumberFromDouble(&mp_decimal_data.epsilon, mp_epsilon);
        decNumberFromString(&mp_decimal_data.EL_GORDO, EL_GORDO_STRING, &mp_decimal_data.set);
        decNumberFromString(&mp_decimal_data.negative_EL_GORDO, NEGATIVE_EL_GORDO_STRING, &mp_decimal_data.set);
        decNumberFromScaled(&mp_decimal_data.no_crossing, mp_fraction_multiplier + 1);
        decNumberFromScaled(&mp_decimal_data.one_crossing, mp_fraction_multiplier);
        decNumberFromScaled(&mp_decimal_data.zero_crossing, 0);
        /* */
        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);
    decNumberCopy(math->md_precision_min.data.num, &mp_decimal_data.two);
    /* 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);
    mp_decimal_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
    decNumberCopy(math->md_inf_t.data.num, &mp_decimal_data.EL_GORDO);
    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);
    mp_decimal_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
    decNumberFromString(math->md_warning_limit_t.data.num, WARNING_LIMIT_STRING, &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, &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);
    decNumberCopy(math->md_two_t.data.num, &mp_decimal_data.two);
    mp_decimal_allocate_number(mp, &math->md_three_t, mp_scaled_type);
    decNumberCopy(math->md_three_t.data.num, &mp_decimal_data.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);
    decNumberCopy(math->md_fraction_one_t.data.num, &mp_decimal_data.fraction_one);
    mp_decimal_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
    decNumberCopy(math->md_fraction_half_t.data.num, &mp_decimal_data.fraction_half);
    mp_decimal_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
    decNumberCopy(math->md_fraction_three_t.data.num, &mp_decimal_data.fraction_three);
    mp_decimal_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
    decNumberCopy(math->md_fraction_four_t.data.num, &mp_decimal_data.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  * mp_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 * mp_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 * mp_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);
    # define mp_coef_bound ((7.0 / 3.0) * mp_fraction_multiplier)
    decNumberFromDouble(math->md_coef_bound_k.data.num, mp_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, mp_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, 0.04096);
    mp_decimal_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
    decNumberFromDouble(math->md_half_fraction_threshold_t.data.num, 0.04096 / 2);
    mp_decimal_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
    decNumberFromDouble(math->md_scaled_threshold_t.data.num, 0.000122);
    mp_decimal_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
    decNumberFromDouble(math->md_half_scaled_threshold_t.data.num, 0.000122 / 2);
    mp_decimal_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
    decNumberFromDouble(math->md_near_zero_angle_t.data.num, 0.0256 * mp_angle_multiplier);
    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, 0x80000);
    mp_decimal_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
    decNumberFromDouble(math->md_equation_threshold_t.data.num, 0.001);
    /* 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_add3                = mp_decimal_pyth_add3;
    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_slow_sub                 = mp_decimal_slow_sub;
    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;
}