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

# include "mpmathscaled.h"
# include "mpstrings.h"

# define  coef_bound               0x25555555
# define  fraction_threshold       2685
# define  half_fraction_threshold  1342
# define  scaled_threshold         8
# define  half_scaled_threshold    4
# define  near_zero_angle          26844
# define  p_over_v_threshold       0x80000
# define  equation_threshold       64

/*tex
    Fixed-point arithmetic is done on {\sl scaled integers} that are multiples of $2^{-16}$. In 
    other words, a binary point is assumed to be sixteen bit positions from the right end of a 
    binary computer word.
*/

# define  unity                    0x10000
# define  two                      (2*unity)
# define  three                    (3*unity)
# define  half_unit                (unity/2)
# define  three_quarter_unit       (3*(unity/4))
# define  EL_GORDO                 0x7FFFFFFF
# define  negative_EL_GORDO        (-EL_GORDO)
# define  one_third_EL_GORDO       0x2AAAAAAA

/*tex We need these preprocessor values */

# define  TWEXP31                  2147483648.0
# define  TWEXP28                  268435456.0
# define  TWEXP16                  65536.0
# define  TWEXP_16                 (1.0/65536.0)
# define  TWEXP_28                 (1.0/268435456.0)
# define  no_crossing              (fraction_one + 1)
# define  one_crossing             fraction_one
# define  zero_crossing            0

/*tex
    The |scaled| quantities in \MP\ programs are generally supposed to be less than $2^{12}$ in 
    absolute value, so \MP\ does much of its internal arithmetic with 28~significant bits of 
    precision. A |fraction| denotes a scaled integer whose binary point is assumed to be 28 bit 
    positions from the right.
*/

# define  fraction_half            0x08000000
# define  fraction_one             0x10000000
# define  fraction_two             0x20000000
# define  fraction_three           0x30000000
# define  fraction_four            0x40000000

/*tex
    Octants are represented in a \quote {Gray code,} since that turns out to be computationally 
    simplest.
*/

# define  negate_x                 1
# define  negate_y                 2
# define  switch_x_and_y           4
# define  first_octant             1
# define  second_octant            (first_octant + switch_x_and_y)
# define  third_octant             (first_octant + switch_x_and_y + negate_x)
# define  fourth_octant            (first_octant                  + negate_x)
# define  fifth_octant             (first_octant                  + negate_x + negate_y)
# define  sixth_octant             (first_octant + switch_x_and_y + negate_x + negate_y)
# define  seventh_octant           (first_octant + switch_x_and_y + negate_y)
# define  eighth_octant            (first_octant                  + negate_y)
# define  forty_five_deg           0x02D00000
# define  ninety_deg               0x05A00000
# define  one_eighty_deg           0x0B400000
# define  negative_one_eighty_deg -0x0B400000
# define  three_sixty_deg          0x16800000
# define  odd(A)                   (abs(A)%2==1)
# define  two_to_the(A)            (1<<(unsigned)(A))
# define  set_cur_cmd(A)           mp->cur_mod_->command = (A)
# define  set_cur_mod(A)           mp->cur_mod_->data.n.data.val = (A)

static const int mp_m_spec_log[29] = {
    0, 93032640, 38612034, 17922280, 8662214, 4261238, 2113709, 1052693, 525315,
    262400, 131136, 65552, 32772, 16385, 8192, 4096, 2048, 1024, 512, 256, 128,
    64, 32, 16, 8, 4, 2, 1, 1
};

static const int mp_m_spec_atan[27] = {
    0, 27855475, 14718068, 7471121, 3750058, 1876857, 938658, 469357, 234682,
    117342, 58671, 29335, 14668, 7334, 3667, 1833, 917, 458, 229, 115, 57, 29,
    14, 7, 4, 2, 1
};

static int mp_scaled_aux_take_fraction (MP mp, int p, int q);
static int mp_scaled_aux_make_fraction (MP mp, int p, int q);
static int mp_scaled_round_unscaled    (mp_number *x_orig);

static void mp_scaled_allocate_number(MP mp, mp_number *n, mp_number_type t)
{
    (void) mp;
    n->data.val = 0;
    n->type = t;
}

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

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

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

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

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

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

static void mp_scaled_allocate_double(MP mp, mp_number *n, double v)
{
    (void) mp;
    n->type = mp_scaled_type;
    n->data.val = (int) (v * 65536.0);
}

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

static void mp_scaled_set_from_int(mp_number *A, int B)
{
    A->data.val = B * 65536;
}

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

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

static void mp_scaled_set_from_double(mp_number *A, double B)
{
    A->data.val = (int) (B * 65536.0);
}

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

static void mp_scaled_set_half_from_addition(mp_number *A, mp_number *B, mp_number *C)
{
    A->data.val = (B->data.val + C->data.val) / 2;
}

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

static void mp_scaled_set_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
{
    A->data.val = (B->data.val - C->data.val) / 2;
}

static void mp_scaled_set_from_div(mp_number *A, mp_number *B, mp_number *C)
{
    A->data.val = B->data.val / C->data.val;
}

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

static void mp_scaled_set_from_int_div(mp_number *A, mp_number *B, int C)
{
    A->data.val = B->data.val / C;
}

static void mp_scaled_set_from_int_mul(mp_number *A, mp_number *B, int C)
{
    A->data.val = B->data.val * C;
}

static void mp_scaled_set_from_of_the_way(MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C)
{
    (void) mp;
    A->data.val = B->data.val - mp_scaled_aux_take_fraction(mp, (B->data.val - C->data.val), t->data.val);
}

static void mp_scaled_negate(mp_number *A)
{
    A->data.val = -A->data.val;
}

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

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

static void mp_scaled_half(mp_number *A)
{
    A->data.val = A->data.val / 2;
}

static void mp_scaled_double(mp_number *A)
{
    A->data.val = A->data.val + A->data.val;
}

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

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

static void mp_scaled_divide_int(mp_number *A, int B)
{
    A->data.val = A->data.val / B;
}

static void mp_scaled_abs(mp_number *A)
{
    A->data.val = abs(A->data.val);
}

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

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

static void mp_scaled_abs_clone(mp_number *A, mp_number *B)
{
    A->data.val = abs(B->data.val);
}

static void mp_scaled_swap(mp_number *A, mp_number *B)
{
    int swap_tmp = A->data.val;
    A->data.val = B->data.val;
    B->data.val = swap_tmp;
}

static void mp_scaled_fraction_to_scaled(mp_number *A)
{
    A->type = mp_scaled_type;
    A->data.val = A->data.val / 4096;
}

static void mp_scaled_angle_to_scaled(mp_number *A)
{
    A->type = mp_scaled_type;
    if (A->data.val >= 0) {
        A->data.val =   ( A->data.val + 8) / 16;
    } else {
        A->data.val = -((-A->data.val + 8) / 16);
    }
}

static void mp_scaled_scaled_to_fraction(mp_number *A)
{
    A->type = mp_fraction_type;
    A->data.val = A->data.val * 4096;
}

static void mp_scaled_scaled_to_angle(mp_number *A)
{
    A->type = mp_angle_type;
    A->data.val = A->data.val * 16;
}

static int mp_scaled_to_int(mp_number *A)
{
    return A->data.val;
}

static int mp_scaled_to_scaled(mp_number *A)
{
    return A->data.val;
}

static int mp_scaled_to_boolean(mp_number *A)
{
    return A->data.val;
}

static double mp_scaled_to_double(mp_number *A)
{
    return A->data.val / 65536.0;
}

static int mp_scaled_odd(mp_number *A)
{
    return odd(mp_scaled_round_unscaled(A));
}

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

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

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

static int mp_scaled_non_equal_abs(mp_number *A, mp_number *B)
{
    return abs(A->data.val) != abs(B->data.val);
}

/*tex 

    One of \MP's most common operations is the calculation of $\lfloor {a + b \over 2}\rfloor$, 
    the midpoint of two given integers |a| and~|b|. The most decent way to do this is to write 
    |(a + b) / 2|; but on many machines it is more efficient to calculate |(a + b) >> 1|.

    Therefore the midpoint operation will always be denoted by |half(a + b)| in this program. If 
    \MP\ is being implemented with languages that permit binary shifting, the |half| macro should 
    be changed to make this operation as efficient as possible. Since some systems have shift 
    operators that can only be trusted to work on positive numbers, there is also a macro |halfp| 
    that is used only when the quantity being halved is known to be positive or zero.

    Here is a procedure analogous to |print_int|. If the output of this procedure is subsequently 
    read by \MP\ and converted by the |round_decimals| routine above, it turns out that the 
    original value will be reproduced exactly. A decimal point is printed only if the value is not 
    an integer. If there is more than one way to print the result with the optimum number of digits
    following the decimal point, the closest possible value is given.

    The invariant relation in the |repeat| loop is that a sequence of decimal digits yet to be 
    printed will yield the original number if and only if they form a fraction~$f$ in the range 
    $s - \delta \L 10 \cdot 2^{16} f < s$. We can stop if and only if $f = 0$ satisfies this 
    condition; the loop will terminate before $s$ can possibly become zero. We round to five 
    digits.

*/

/*tex Todo: check if we can replace this: */

static char *mp_string_scaled(MP mp, int s)
{
    static char scaled_string[32];
    int i = 0;
    (void) mp;
    if (s < 0) {
        scaled_string[i++] = '-';
        s = -s;
    }
    /*tex The integer part: */
    snprintf((scaled_string+i), 12, "%d", (int) (s / unity));
    while (*(scaled_string+i)) {
        i++;
    }
    s = 10 * (s % unity) + 5;
    if (s != 5) {
        /*tex The amount of allowable inaccuracy, scaled: */
        int delta = 10;
        scaled_string[i++] = '.';
        do {
            /*tex Round the final digit: */
            if (delta > unity) {
                s = s + 0x8000 - (delta / 2);
            }
            scaled_string[i++] = (char) ((s / unity) + '0');
            s = 10 * (s % unity);
            delta = delta * 10;
        } while (s > delta);
    }
    scaled_string[i] = '\0';
    return scaled_string;
}

/*tex
    Addition is not always checked to make sure that it doesn't overflow, but in places where 
    overflow isn't too unlikely the |slow_add| routine is used.
*/

static void mp_scaled_slow_add(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
{
    int x = x_orig->data.val;
    int y = y_orig->data.val;
    if (x >= 0) {
        if (y <= EL_GORDO - x) {
            ret->data.val = x + y;
        } else {
            mp->arithmic_error = 1;
            ret->data.val =  EL_GORDO;
        }
    } else if (-y <= EL_GORDO + x) {
        ret->data.val = x + y;
    } else {
        mp->arithmic_error = 1;
        ret->data.val = negative_EL_GORDO;
    }
}

/*tex

    The |make_fraction| routine produces the |fraction| equivalent of |p/q|, given integers |p| 
    and~|q|; it computes the integer $f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are
    positive. If |p| and |q| are both of the same scaled type |t|, the \quote {type relation}
    |make_fraction(t,t) = fraction| is valid; and it's also possible to use the subroutine \quote 
    {backwards,} using the relation |make_fraction(t,fraction) = t| between scaled types.

    If the result would have magnitude $2^{31}$ or more, |make_fraction| sets |arithmic_error := 
    true|. Most of \MP's internal computations have been designed to avoid this sort of error.

    If this subroutine were programmed in assembly language on a typical machine, we could simply
    compute |(*p)div q|, since a double-precision product can often be input to a fixed-point 
    division instruction. But when we are restricted to integer arithmetic it is necessary either
    to resort to multiple-precision maneuvering or to use a simple but slow iteration. The 
    multiple-precision technique would be about three times faster than the code adopted here, but 
    it would be comparatively long and tricky, involving about sixteen additional multiplications 
    and divisions.

    This operation is part of \MP's \quote {inner loop}; indeed, it will consume nearly 10\pct! of 
    the running time (exclusive of input and output) if the code below is left unchanged. A 
    machine-dependent recoding will therefore make \MP\ run faster. The present implementation is 
    highly portable, but slow; it avoids multiplication and division except in the initial stage. 
    System wizards should be careful to replace it with a routine that is guaranteed to produce 
    identical results in all cases. 

    As noted below, a few more routines should also be replaced by machine-dependent code, for 
    efficiency. But when a procedure is not part of the \quote {inner loop,} such changes aren't 
    advisable; simplicity and robustness are preferable to trickery, unless the cost is too high. 

*/

/* 
    Todo: check if we can replace this (see mpmathdouble): 
*/

static int mp_scaled_aux_make_fraction(MP mp, int p, int q)
{
    if (q == 0) {
        mp_confusion (mp, "division by zero");
        return 0;
    } else {
        double d = TWEXP28 * (double) p / (double) q;
        if ((p ^ q) >= 0) {
            d += 0.5;
            if (d >= TWEXP31) {
                mp->arithmic_error = 1;
                return EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? -q : q) & 0x7FFF) * (((i & 0x3FFF) << 1) - 1) & 0x800) != 0) {
                    --i;
                }
                return i;
            }
        } else {
            d -= 0.5;
            if (d <= -TWEXP31) {
                mp->arithmic_error = 1;
                return -negative_EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? q : -q) & 0x7FFF) * (((i & 0x3FFF) << 1) + 1) & 0x800) != 0) {
                    ++i;
                }
                return i;
            }
        }
    }
}

static void mp_scaled_make_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q)
{
    ret->data.val = mp_scaled_aux_make_fraction (mp, p->data.val, q->data.val);
}

/*tex
    The dual of |make_fraction| is |take_fraction|, which multiplies a given integer~|q| by a 
    fraction~|f|. When the operands are positive, it computes $p = \lfloor qf/2^{28} + { 1\over 2}
    \rfloor$, a symmetric function of |q| and~|f|.

    This routine is even more \quote {inner loopy} than |make_fraction|; the present implementation 
    consumes almost 20\pct! of \MP's computation time during typical jobs, so a machine-language 
    substitute is advisable.
*/

static int mp_scaled_aux_take_fraction(MP mp, int p, int q)
{
    double d = (double) p *(double) q *TWEXP_28;
    if ((p ^ q) >= 0) {
        d += 0.5;
        if (d >= TWEXP31) {
            if (d != TWEXP31 || (((p & 0x7FFF) * (q & 0x7FFF)) & 040000) == 0) {
                mp->arithmic_error = 1;
            }
            return EL_GORDO;
        } else {
            int i = (int) d;
            if (d == (double) i && (((p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
                --i;
            }
            return i;
        }
    } else {
        d -= 0.5;
        if (d <= -TWEXP31) {
            if (d != -TWEXP31 || ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) == 0) {
                mp->arithmic_error = 1;
            }
            return -negative_EL_GORDO;
        } else {
            int i = (int) d;
            if (d == (double) i && ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
                ++i;
            }
            return i;
        }
    }
}

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

/*tex

    When we want to multiply something by a |scaled| quantity, we use a scheme analogous to 
    |take_fraction| but with a different scaling. Given positive operands, |take_scaled| computes 
    the quantity $p=\lfloor qf/2^{16}+{1\over2}\rfloor$.

    Once again it is a good idea to use a machine-language replacement if possible; otherwise 
    |take_scaled| will use more than 2\pct! of the running time when the Computer Modern fonts are
    being generated. 

*/

static int mp_take_scaled(MP mp, int p, int q)
{ 
    /* q = scaled */
    double d = (double) p * (double) q * TWEXP_16;
    if ((p ^ q) >= 0) {
        d += 0.5;
        if (d >= TWEXP31) {
            if (d != TWEXP31 || (((p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) == 0) {
                mp->arithmic_error = 1;
            }
            return EL_GORDO;
        } else {
            int i = (int) d;
            if (d == (double) i && (((p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
                --i;
            }
            return i;
        }
    } else {
        d -= 0.5;
        if (d <= -TWEXP31) {
            if (d != -TWEXP31 || ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) == 0) {
                mp->arithmic_error = 1;
            }
            return -negative_EL_GORDO;
        } else {
            int i = (int) d;
            if (d == (double) i && ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
                ++i;
            }
            return i;
        }
    }
}

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

/*tex 

    For completeness, there's also |make_scaled|, which computes a quotient as a |scaled| number 
    instead of as a |fraction|. In other words, the result is $\lfloor2^{16}p/q+{1\over2}\rfloor$, 
    if the operands are positive. \ (This procedure is not used especially often, so it is not 
    part of \MP's inner loop.)

*/

static int mp_make_scaled(MP mp, int p, int q)
{
    if (q == 0) {
        mp_confusion(mp, "division by zero");
        return 0;
    } else {
        double d = TWEXP16 * (double) p / (double) q;
        if ((p ^ q) >= 0) {
            d += 0.5;
            if (d >= TWEXP31) {
                mp->arithmic_error = 1;
                return EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? -q : q) & 0x7FFF) * (((i & 0x3FFF) << 1) - 1) & 0x800) != 0) {
                    --i;
                }
                return i;
            }
        } else {
            d -= 0.5;
            if (d <= -TWEXP31) {
                mp->arithmic_error = 1;
                return -negative_EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? q : -q) & 0x7FFF) * (((i & 0x3FFF) << 1) + 1) & 0x800) != 0) {
                    ++i;
                }
                return i;
            }
        }
    }
}

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

/*tex
    The following function is used to create a scaled integer from a given decimal fraction 
    $(.d_0d_1\ldots d_{k-1})$, where |0<=k<=17|. This converts a decimal fraction.
*/

static int mp_round_decimals(MP mp, unsigned char *b, int k)
{
    unsigned a = 0;
    (void) mp; 
    for (int l = k - 1; l >= 0; l--) {
        if (l < 16) {
            /*tex digits for |k >= 17| cannot affect the result */
            a = (a + (unsigned) (*(b + l) - '0') * two) / 10;
        }
    }
    return (int) (a + 1)/2; /* (a + 1)  >> 1 */
}

/*tex 

    We no longer have these character mapping so we can just as well do the next without class 
    checking. Also because signs are hard checked.
*/

static void mp_scaled_wrapup_numeric_token(MP mp, int n, int f)
{
    if (n < 32768) {
        int mod = (n * unity + f); /* scaled */
        set_cur_mod(mod);
        if (mod >= fraction_one) {
            if (internal_value(mp_warning_check_internal).data.val > 0 && (mp->scanner_status != mp_tex_flushing_state)) {
                char msg[256];
                snprintf(msg, 256, "Number is too large (%s)", mp_string_scaled(mp,mod));
                mp_error(
                    mp,
                    msg,
                    "It is at least 4096. Continue and I'll try to cope with that big value;\n"
                    "but it might be dangerous. (Set warningcheck:=0 to suppress this message.)"
                );
            }
        }
    } else if (mp->scanner_status != mp_tex_flushing_state) {
        mp_error(
            mp,
            "Enormous number has been reduced",
            "I can\'t handle numbers bigger than 32767.99998; so I've changed your constant\n"
            "to that maximum amount."
        );
        set_cur_mod(EL_GORDO);
    }
    set_cur_cmd(mp_numeric_command);
}

static void mp_scaled_scan_fractional_token(MP mp, int n)
{ 
    /* n: scaled */
    int f; /* scaled */
    int k = 0;
    do {
        k++;
        mp->cur_input.loc_field++;
    } while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class);
    f = mp_round_decimals(mp, (unsigned char *)(mp->buffer + mp->cur_input.loc_field - k), (int) k);
    if (f == unity) {
        n++;
        f = 0;
    }
    mp_scaled_wrapup_numeric_token(mp, n, f);
}

static void mp_scaled_scan_numeric_token(MP mp, int n)
{
    while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
        if (n < 32768) {
            n = 10 * n + mp->buffer[mp->cur_input.loc_field] - '0';
        }
        mp->cur_input.loc_field++;
    }
    if (! (mp->buffer[mp->cur_input.loc_field] == '.' && mp->char_class[mp->buffer[mp->cur_input.loc_field + 1]] == mp_digit_class)) {
        mp_scaled_wrapup_numeric_token(mp, n, 0);
    } else {
        mp->cur_input.loc_field++;
        mp_scaled_scan_fractional_token(mp, n);
    }
}

/*tex

    Here is a typical example of how the routines above can be used. It computes the function 

    $$
        {1 \over 3 \tau}f(\theta,\phi)= {\tau^{-1} \bigl( 2 + \sqrt2\, (\sin \theta -{1 \over 16}
        \sin \phi) (\sin \phi - {1 \over 16} \sin \theta) (\cos \theta - \cos \phi) \bigr) \over
        3\, \bigl( 1 + {1 \over 2} (\sqrt5 - 1) \cos \theta + {1 \over 2} (3 - \sqrt 5 \,) \cos 
        \phi \bigr)},
    $$

    where $\tau$ is a |scaled| \quote {tension} parameter. This is \MP's magic fudge factor for 
    placing the first control point of a curve that starts at an angle $\theta$ and ends at an 
    angle $\phi$ from the straight path. (Actually, if the stated quantity exceeds 4, \MP\ reduces 
    it to~4.)

    The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$. (It's a sum of 
    eight terms whose absolute values can be bounded using relations such as $\sin \theta \cos 
    \theta |1 \over 2|$.) Thus the numerator is positive; and since the tension $\tau$ is 
    constrained to be at least $3\over4$, the numerator is less than $16 \over 3$. The denominator 
    is nonnegative and at most~6. Hence the fixed-point calculations below are guaranteed to stay 
    within the bounds of a 32-bit computer word.

    The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction| arguments |st|, 
    |ct|, |sf|, and |cf|, representing $\sin \theta$, $\cos \theta$, $\sin \phi$, and $\cos \phi$, 
    respectively.

*/

static void mp_scaled_velocity(MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t)
{
    int acc, num, denom; 
    acc = mp_scaled_aux_take_fraction(mp, st->data.val - (sf->data.val / 16), sf->data.val - (st->data.val / 16));
    acc = mp_scaled_aux_take_fraction(mp, acc, ct->data.val - cf->data.val);
    num = fraction_two + mp_scaled_aux_take_fraction(mp, acc, 379625062);
    /*tex
        $2^{28}\sqrt2\approx379625062.497$
    */
    denom = fraction_three + mp_scaled_aux_take_fraction(mp, ct->data.val, 497706707) + mp_scaled_aux_take_fraction(mp, cf->data.val, 307599661);
    /*tex
        $3\cdot2^{27}\cdot(\sqrt5-1)\approx497706706.78$ and
        $3\cdot2^{27}\cdot(3-\sqrt5\,)\approx307599661.22$
    */
    if (t->data.val != unity) {
        /*tex 
            |make_scaled(fraction,scaled) = fraction| 
        */
        num = mp_make_scaled(mp, num, t->data.val);
    }
    if (num / 4 >= denom) {
        ret->data.val = fraction_four;
    } else {
        ret->data.val = mp_scaled_aux_make_fraction(mp, num, denom);
    }
}

/*tex
    The following somewhat different subroutine tests rigorously if $ab$ is greater than, equal to,
    or less than~$cd$, given integers $(a,b,c,d)$. In most cases a quick decision is reached. The 
    result is $+1$, 0, or~$-1$ in the three respective cases.
*/

static int mp_scaled_ab_vs_cd(mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig)
{
    int a = a_orig->data.val;
    int b = b_orig->data.val;
    int c = c_orig->data.val;
    int d = d_orig->data.val;
    if (a < 0) {
        a = -a;
        b = -b;
    }
    if (c < 0) {
        c = -c;
        d = -d;
    }
    if (d <= 0) {
        if (b >= 0) {
            if ((a == 0 || b == 0) && (c == 0 || d == 0)) {
                return 0;
            } else {
                return 1;
            }
        } else if (d == 0) {
            return a == 0 ? 0 : -1;
        } else {
            int q = a;
            a = c;
            c = q;
            q = -b;
            b = -d;
            d = q;
        }
    } else if (b <= 0) {
        if (b < 0 && a > 0) {
            return -1;
        } else {
            return c == 0 ? 0 : -1;
        }
    }
    while (1) {
        int q = a / d;
        int r = c / b;
        if (q != r) {
            return q > r ? 1 : -1;
        } else {
            q = a % d;
            r = c % b;
            if (r == 0) {
                return q ? 1 : 0;
            } else if (q == 0) {
                return -1;
            } else {
                a = b;
                b = q;
                c = d;
                d = r;
            }
        }
    }
    /*tex Now |a > d > 0| and |c > b > 0|. */
}

/* 
    Now here's a subroutine that's handy for all sorts of path computations: Given a quadratic
    polynomial $B(a,b,c;t)$, the |crossing_point| function returns the unique |fraction| value |t|
    between 0 and~1 at which $B(a,b,c;t)$ changes from positive to negative, or returns |t = 
    fraction_one + 1| if no such value exists. If |a < 0| (so that $B(a,b,c;t)$ is already negative
    at |t = 0|), |crossing_point| returns the value zero.

    The general bisection method is quite simple when $n=2$, hence |crossing_point| does not take 
    much time. At each stage in the recursion we have a subinterval defined by |l| and~|j| such 
    that $B(a,b,c;2^{-l}(j+t))=B(x_0,x_1,x_2;t)$, and we want to \quote {zero in} on the 
    subinterval where $x_0\G0$ and $\min(x_1,x_2)<0$.

    It is convenient for purposes of calculation to combine the values of |l| and~|j| in a single 
    variable $d=2^l+j$, because the operation of bisection then corresponds simply to doubling 
    $d$ and possibly adding~1. Furthermore it proves to be convenient to modify our previous 
    conventions for bisection slightly, maintaining the variables $X_0 = 2^lx_0$, $X_1 = 2^l(x_0 -
    x_1)$, and $X_2 = 2^l(x_1 - x_2)$. With these variables the conditions $x_0 \ge 0$ and
    $\min(x_1,x_2) < 0$ are equivalent to $\max(X_1,X_1+X_2) > X_0 \ge 0$.

    The following code maintains the invariant relations $0\L|x0| < \max(|x1|,|x1| + |x2|)$, 
    $\vert|x1|\vert < 2^{30}$, $\vert|x2|\vert < 2^{30}$; it has been constructed in such a way 
    that no arithmetic overflow will occur if the inputs satisfy $a < 2^{30}$, $\vert a - b\vert <
    2^{30}$, and $\vert b - c\vert < 2^{30}$.
*/

static void mp_scaled_crossing_point(MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc)
{
    int x, xx, x0, x1, x2; /* temporary registers for bisection */
    int a = aa->data.val;
    int b = bb->data.val;
    int c = cc->data.val;
    int d;                 /* recursive counter */
    (void) mp;
    if (a < 0) {
        ret->data.val = zero_crossing;
        return;
    } else if (c >= 0) {
        if (b >= 0) {
            if (c > 0) {
                ret->data.val = no_crossing;
            } else if ((a == 0) && (b == 0)) {
                ret->data.val = no_crossing;
            } else {
                ret->data.val = one_crossing;
            }
            return;
        } else if (a == 0) {
            ret->data.val = zero_crossing;
            return;
        }
    } else if (a == 0) {
        if (b <= 0) {
            ret->data.val = zero_crossing;
            return;
        }
    }
    /* Use bisection to find the crossing point... */
    d = 1;
    x0 = a;
    x1 = a - b;
    x2 = b - c;
    do {
        x = (x1 + x2) / 2;
        if (x1 - x0 > x0) {
            x2 = x;
            x0 += x0;
            d += d;
        } else {
            xx = x1 + x - x0;
            if (xx > x0) {
                x2 = x;
                x0 += x0;
                d += d;
            } else {
                x0 = x0 - xx;
                if ((x <= x0) && (x + x2 <= x0)) {
                    ret->data.val = no_crossing;
                    return;
                } else {
                    x1 = x;
                    d = d + d + 1;
                }
            }
        }
    } while (d < fraction_one);
    ret->data.val = d - fraction_one;
}

/*tex
    We conclude this set of elementary routines with some simple rounding and truncation 
    operations; |round_unscaled| rounds a |scaled| and converts it to |int|
*/

static int mp_scaled_round_unscaled(mp_number *x_orig)
{
    int x = x_orig->data.val;
    if (x >= 32768) {
        return 1 + ((x-32768) / 65536);
    } else if (x >= -32768) {
        return 0;
    } else {
        return -(1+((-(x+1)-32768) / 65536));
    }
}
static void mp_scaled_floor(mp_number *i)
{
    i->data.val = i->data.val & -65536;
}

static void mp_scaled_fraction_to_round_scaled(mp_number *x_orig)
{
    int x = x_orig->data.val;
    x_orig->type = mp_scaled_type;
    x_orig->data.val = (x>=2048 ? 1+((x-2048) / 4096)  : ( x>=-2048 ? 0 : -(1+((-(x+1)-2048) / 4096))));
}

/*tex
    \MP\ computes all of the necessary special functions from scratch, without relying on |real| 
    arithmetic or system subroutines for sines, cosines, etc.
    
    To get the square root of a |scaled| number |x|, we want to calculate $s = \lfloor 2^8\! \sqrt 
    x + {1\over2} \rfloor$. If $x > 0$, this is the unique integer such that $2^{16}x - s\L s^2 < 
    2^{16}x + s$. The following subroutine determines $s$ by an iterative method that maintains the
    invariant relations $x=2^{46 - 2k}x_0 \bmod 2^{30}$, $0 < y = \lfloor 2^{16 - 2k}x_0\rfloor -
    s^2 + s\L q = 2s$, where $x_0$ is the initial value of $x$. The value of~$y$ might, however, be 
    zero at the start of the first iteration.
*/

static void mp_scaled_sqrt(MP mp, mp_number *ret, mp_number *x_orig)
{
    int x = x_orig->data.val;
    if (x <= 0) {
        if (x < 0) {
            char msg[256];
            snprintf(msg, 256, "Square root of %s has been replaced by 0", mp_string_scaled(mp, x));
            mp_error(
                mp,
                msg,
                "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
                "Proceed, with fingers crossed."
            );
        }
        ret->data.val = 0;
    } else {
        int k = 23; /* iteration control counter */
        int y;
        int q = 2;
        while (x < fraction_two) {  /* i.e., |while x<|\unskip */
            k--;
            x = x + x + x + x;
        }
        if (x < fraction_four)
            y = 0;
        else {
            x = x - fraction_four;
            y = 1;
        }
        do {
            /*tex Decrease |k| by 1, maintaining the invariant relations between |x|, |y|, and~|q|: */
            x += x;
            y += y;
            if (x >= fraction_four) {
                /* note that |fraction_four=| */
                x = x - fraction_four;
                y++;
            };
            x += x;
            y = y + y - q;
            q += q;
            if (x >= fraction_four) {
                x = x - fraction_four;
                y++;
            };
            if (y > (int) q) {
                y -= q;
                q += 2;
            } else if (y <= 0) {
                q -= 2;
                y += q;
            };
            k--;
        } while (k != 0);
        ret->data.val = (int) (q/2);
    }
}

/*tex
    Pythagorean addition $\psqrt{a^2 + b^2}$ is implemented by an elegant iterative scheme due to 
    Cleve Moler and Donald Morrison [{\sl IBM Journal   of Research and Development \bf27} (1983),
    577--581]. It modifies |a| and~|b| in such a way that their Pythagorean sum remains invariant, 
    while the smaller argument decreases.
*/

static void mp_scaled_pyth_add(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
{
    int a = abs(a_orig->data.val);
    int b = abs(b_orig->data.val);
    if (a < b) {
        int r = b;
        b = a;
        a = r;
    }
    /* now |0<=b<=a| */
    if (b > 0) {
        int big;  
        /*tex Is the result dangerously near $2^{31}$? */
        if (a < fraction_two) {
            big = 0;
        } else {
            a = a / 4;
            b = b / 4;
            big = 1;
        }
        /*tex 
            We reduced the precision to avoid arithmetic overflow. Replace |a| by an approximation 
            to $\psqrt{a^2+b^2}:$ 
        */
        while (1) {
            int r = mp_scaled_aux_make_fraction(mp, b, a);
            r = mp_scaled_aux_take_fraction(mp, r, r);
            /*tex Now $r\approx b^2/a^2$: */
            if (r == 0) {
                break;
            } else {
                r = mp_scaled_aux_make_fraction(mp, r, fraction_four + r);
                a = a + mp_scaled_aux_take_fraction(mp, a + a, r);
                b = mp_scaled_aux_take_fraction(mp, b, r);
            }
        }
        if (big) {
            if (a < fraction_two) {
                a = a + a + a + a;
            } else {
                mp->arithmic_error = 1;
                a = EL_GORDO;
            }
        }
    }
    ret->data.val = a;
}

/*tex
    The key idea here is to reflect the vector $(a,b)$ about the line through $(a,b/2)$. Here is a
    similar algorithm for $\psqrt{a^2-b^2}$. It converges slowly when $b$ is near $a$, but 
    otherwise it works fine.
*/

static void mp_scaled_pyth_sub(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
{
    int a = abs(a_orig->data.val);
    int b = abs(b_orig->data.val);
    if (a <= b) {
        /*tex Handle erroneous |pyth_sub| and set |a:=0|: */
        if (a < b) {
            char msg[256];
            char *astr = mp_strdup(mp_string_scaled(mp, a));
            snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, mp_string_scaled(mp, b));
            mp_memory_free(astr);
            mp_error(
                mp,
                msg,
                "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
                "Proceed, with fingers crossed."
            );
        }
        a = 0;
    } else {
        int big;  
        /*tex Is the result dangerously near $2^{31}$? */
        if (a < fraction_four) {
            big = 0;
        } else {
            a = (int) a/2;
            b = (int) b/2;
            big = 1;
        }
        /*tex Replace |a| by an approximation to $\psqrt{a^2-b^2}$ */
        while (1) {
            int r = mp_scaled_aux_make_fraction(mp, b, a);
            r = mp_scaled_aux_take_fraction(mp, r, r);
            /*tex Now $r\approx b^2/a^2$. */
            if (r == 0) {
                break;
            } else {
                r = mp_scaled_aux_make_fraction(mp, r, fraction_four - r);
                a = a - mp_scaled_aux_take_fraction(mp, a + a, r);
                b = mp_scaled_aux_take_fraction(mp, b, r);
            }
        }
        if (big) {
            a *= 2;
        }
    }
    ret->data.val = a;
}

/*tex 
    We just abuse doubles here.
*/

static void mp_scaled_power_of(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
{
    double p = pow(mp_scaled_to_double(a_orig), mp_scaled_to_double(b_orig));
    long r = lround(p * 65536.0);
    if (r > 0) {
        if (r >= EL_GORDO) {
            mp->arithmic_error = 1;
            r = EL_GORDO;
        }
    } else if (r < 0) {
        if (r <= - EL_GORDO) {
            mp->arithmic_error = 1;
            r = - EL_GORDO;
        }
    }
    ret->data.val = r;
}

/*tex

    The subroutines for logarithm and exponential involve two tables. The first is simple: 
    |two_to_the[k]| equals $2^k$. The second involves a bit more calculation, which the author 
    claims to have done correctly: |mp_m_spec_log[k]| is $2^{27}$ times $\ln\bigl(1/(1-2^{-k})
    \bigr)= 2^{-k}+{1\over2}2^{-2k}+{1\over3}2^{-3k}+\cdots\,$, rounded to the nearest integer.
    Here is the routine that calculates $2^8$ times the natural logarithm of a |scaled| quantity; 
    it is an integer approximation to $2^{24}\ln(x/2^{16})$, when |x| is a given positive integer.

    The method is based on exercise 1.2.2--25 in {\sl The Art of Computer Programming}: During the
    main iteration we have $1\L 2^{-30}x<1/(1-2^{1-k})$, and the logarithm of $2^{30}x$ remains to 
    be added to an accumulator register called~$y$. Three auxiliary bits of accuracy are retained 
    in~$y$ during the calculation, and sixteen auxiliary bits to extend |y| are kept in~|z| during 
    the initial argument reduction. (We add $100\cdot2^{16}=6553600$ to~|z| and subtract 100 
    from~|y| so that |z| will not become negative; also, the actual amount subtracted from~|y| 
    is~96, not~100, because we want to add~4 for rounding before the final division by~8.)

*/

static void mp_scaled_m_log(MP mp, mp_number *ret, mp_number *x_orig)
{
    int x = x_orig->data.val;
    if (x <= 0) {
        /* Handle non-positive logarithm: */
        char msg[256];
        snprintf(msg, 256, "Logarithm of %s has been replaced by 0", mp_string_scaled(mp, x));
        mp_error(
            mp,
            msg,
            "Since I don't take logs of non-positive numbers, I'm zeroing this one.\n"
            "Proceed, with fingers crossed."
        );
        ret->data.val = 0;
    } else {
        int k = 2;                    /* iteration counter, starts at 2 */
        int y = 1302456956 + 4 - 100; /* $14\times2^{27}\ln2\approx1302456956.421063$ */
        int z = 27595 + 6553600;      /* and $2^{16}\times .421063\approx 27595$ */
        /* $2^{27}\ln2\approx 93032639.74436163$ and $2^{16}\times.74436163\approx 48782$ */
        while (x < fraction_four) {
            x = 2*x;
            y -= 93032639;
            z -= 48782;
        }
        y = y + (z / unity);
        while (x > fraction_four + 4) {
            /*tex
                Increase |k| until |x| can be multiplied by a factor of $2^{-k}$, and adjust $y$ 
                accordingly.
            */
            z = ((x - 1) / two_to_the (k)) + 1; /* $z=\lceil x/2^k\rceil$ */
            while (x < fraction_four + z) {
                z = (z + 1)/2;
                k++;
            };
            y += mp_m_spec_log[k];
            x -= z;
        }
        ret->data.val = (y / 8);
    }
}

/*tex

    Conversely, the exponential routine calculates $\exp(x/2^8)$, when |x| is |scaled|. The result 
    is an integer approximation to $2^{16}\exp(x/2^{24})$, when |x| is regarded as an integer.

*/

static void mp_scaled_m_exp(MP mp, mp_number *ret, mp_number *x_orig)
{
    int y, z;  /* auxiliary registers */
    int x = x_orig->data.val;
    if (x > 174436200) {
        /* $2^{24}\ln((2^{31}-1)/2^{16})\approx 174436199.51$ */
        mp->arithmic_error = 1;
        ret->data.val = EL_GORDO;
    } else if (x < -197694359) {
        /* $2^{24}\ln(2^{-1}/2^{16})\approx-197694359.45$ */
        ret->data.val = 0;
    } else {
        if (x <= 0) {
            z = -8 * x;
            y = 0x100000; /* $y=2^{20}$ */
        } else {
            if (x <= 127919879) {
                z = 1023359037 - 8 * x;
                /* $2^{27}\ln((2^{31}-1)/2^{20})\approx 1023359037.125$ */
            } else {
                /* |z| is always nonnegative */
                z = 8 * (174436200 - x);
            }
            y = EL_GORDO;
        }
        /* Multiply |y| by $\exp(-z/2^{27})$: */
        {
            int k = 1;
            while (z > 0) {
                while (z >= mp_m_spec_log[k]) {
                    z -= mp_m_spec_log[k];
                    y = y - 1 - ((y - two_to_the(k - 1)) / two_to_the(k));
                }
                k++;
            }
        }
        if (x <= 127919879) {
            ret->data.val = ((y + 8) / 16);
        } else {
            ret->data.val = y;
        }
    }
}

/*tex

    The idea here is that subtracting |mp_m_spec_log[k]| from |z| corresponds to multiplying |y| 
    by $1-2^{-k}$.

    A subtle point (which had to be checked) was that if $x=127919879$, the value of~|y| will 
    decrease so that |y+8| doesn't overflow. In fact, $z$ will be 5 in this case, and |y| will 
    decrease by~64 when |k=25| and by~16 when |k=27|.The trigonometric subroutines use an auxiliary
    table such that |spec_atan[k]| contains an approximation to the |angle| whose tangent is~$1 /
    2^k$. $\arctan2^{-k}$ times $2^{20}\cdot180/\pi$Given integers |x| and |y|, not both zero, the 
    |n_arg| function returns the |angle| whose tangent points in the direction $(x,y)$. This 
    subroutine first determines the correct octant, then solves the problem for |0<=y<=x|, then
    converts the result appropriately to return an answer in the range |-one_eighty_deg <= 
    one_eighty_deg|. (The answer is |+one_eighty_deg| if |y=0| and |x<0|, but an answer of 
    |-one_eighty_deg| is possible if, for example, |y = -1| and $x = -2^{30}$.)

*/

static void mp_scaled_n_arg(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
{
    int octant; /* octant code */
    int x = x_orig->data.val;
    int y = y_orig->data.val;
    if (x >= 0) {
        octant = first_octant;
    } else {
        x = -x;
        octant = first_octant + negate_x;
    }
    if (y < 0) {
        y = -y;
        octant = octant + negate_y;
    }
    if (x < y) {
        int t = y;
        y = x;
        x = t;
        octant = octant + switch_x_and_y;
    }
    if (x == 0) {
        if (internal_value(mp_default_zero_angle_internal).data.val < 0) {
            mp_error(
                mp,
                "angle(0,0) is taken as zero",
                "The 'angle' between two identical points is undefined. I'm zeroing this one.\n"
                "Proceed, with fingers crossed."
            );
            ret->data.val = 0;
        } else { 
            ret->data.val = internal_value(mp_default_zero_angle_internal).data.val;
        }
    } else {
        int z = 0; /* auxiliary register */
        ret->type = mp_angle_type;
        /* Set variable |z| to the arg of $(x,y)$: */
        while (x >= fraction_two) {
            x = x/2;
            y = y/2;
        }
        if (y > 0) {
            int k= 0; /* loop counter */
            while (x < fraction_one) {
                x += x;
                y += y;
            };
            /* Increase |z| to the arg of $(x,y)$: */
            do {
                y += y;
                k++;
                if (y > x) {
                    int t = x; 
                    z = z + mp_m_spec_atan[k];
                    x = x + (y / two_to_the(k + k));
                    y = y - t;
                };
            } while (k != 15);
            do {
                y += y;
                k++;
                if (y > x) {
                    z = z + mp_m_spec_atan[k];
                    y = y - x;
                };
            } while (k != 26);
        }
        /* Return an appropriate answer based on |z| and |octant|: */
        switch (octant) {
            case first_octant:   ret->data.val =  z;                  break;
            case second_octant:  ret->data.val = -z + ninety_deg;     break;
            case third_octant:   ret->data.val =  z + ninety_deg;     break;
            case fourth_octant:  ret->data.val = -z + one_eighty_deg; break;
            case fifth_octant:   ret->data.val =  z - one_eighty_deg; break;
            case sixth_octant:   ret->data.val = -z - ninety_deg;     break;
            case seventh_octant: ret->data.val =  z - ninety_deg;     break;
            case eighth_octant:  ret->data.val = -z;                  break;
        }
// printf("I x=%f y=%f atan=%f\n",y_orig->data.val/65536.0,x_orig->data.val/65536.0,ret->data.val/65536.0);
    }
}

/*tex
    At this point we have |x>=y>=0|, and |x>0|. The numbers are scaled up or down until $2^{28}\L 
    x<2^{29}$, so that accurate fixed-point calculations will be made.
*/

/*tex

    During the calculations of this section, variables |x| and~|y| represent actual coordinates 
    $(x,2^{-k}y)$. We will maintain the condition |x>=y|, so that the tangent will be at most 
    $2^{-k}$. If $x<2y$, the tangent is greater than $2^{-k-1}$. The transformation $(a,b) \mapsto
    (a+b\tan\phi,b-a\tan\phi)$ replaces $(a,b)$ by coordinates whose angle has decreased by~$\phi$; 
    in the special case $a=x$, $b=2^{-k}y$, and $\tan\phi=2^{-k-1}$, this operation reduces to the 
    particularly simple iteration shown here. [Cf.~John E. Meggitt,  {\sl IBM Journal of Research 
    and Development \bf6} (1962), 210--226.]

    The initial value of |x| will be multiplied by at most $(1 + {1\over2})(1 + {1\over8}) (1 + 
    {1\over32})\cdots\approx 1.7584$; hence there is no chance of integer overflow. Conversely, the 
    |n_sin_cos| routine takes an |angle| and produces the sine and cosine of that angle. The 
    results of this routine are stored in global integer variables |n_sin| and |n_cos|.

    Given an integer |z| that is $2^{20}$ times an angle $\theta$ in degrees, the purpose of 
    |n_sin_cos(z)| is to set |x=| and |y=| (approximately), for some rather large number~|r|. The
    maximum of |x| and |y| will be between $2^{28}$ and $2^{30}$, so that there will be hardly any 
    loss of accuracy. Then |x| and~|y| are divided by~|r|.

*/

static void mp_scaled_n_sin_cos(MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin)
{
    int k;                    /* loop control variable */
    int q;                    /* specifies the quadrant */
    int x, y, t;              /* temporary registers */
    int z = z_orig->data.val; /* scaled */
    mp_number x_n, y_n, ret;
    mp_scaled_allocate_number(mp, &ret, mp_scaled_type);
    mp_scaled_allocate_number(mp, &x_n, mp_scaled_type);
    mp_scaled_allocate_number(mp, &y_n, mp_scaled_type);
    while (z < 0) {
        z = z + three_sixty_deg;
    }
    z = z % three_sixty_deg;
    /*tex Now |0 <= z < three_sixty_deg|. */
    q = z / forty_five_deg;
    z = z % forty_five_deg;
    x = fraction_one;
    y = x;
    if (! odd(q)) {
        z = forty_five_deg - z;
    }
    /*tex Subtract angle |z| from |(x,y)|: */
    k = 1;
    while (z > 0) {
        if (z >= mp_m_spec_atan[k]) {
            z = z - mp_m_spec_atan[k];
            t = x;
            x = t + y / two_to_the(k);
            y = y - t / two_to_the(k);
        }
        k++;
    }
    if (y < 0) {
        /*tex This precaution may never be needed */
        y = 0;
    }
    /*tex Convert |(x,y)| to the octant determined by~|q|: */
    switch (q) {
        case 0:                        break;
        case 1: t = x; x =  y; y =  t; break;
        case 2: t = x; x = -y; y =  t; break;
        case 3:        x = -x;         break;
        case 4:        x = -x; y = -y; break;
        case 5: t = x; x = -y; y = -t; break;
        case 6: t = x; x =  y; y = -t; break;
        case 7:                y = -y; break;
    }
    x_n.data.val = x;
    y_n.data.val = y;
    mp_scaled_pyth_add(mp, &ret, &x_n, &y_n);
    n_cos->data.val = mp_scaled_aux_make_fraction(mp, x, ret.data.val);
    n_sin->data.val = mp_scaled_aux_make_fraction(mp, y, ret.data.val);
    mp_scaled_free_number(mp, &ret);
    mp_scaled_free_number(mp, &x_n);
    mp_scaled_free_number(mp, &y_n);
}

/*tex

    In this case the octants are numbered sequentially.The main iteration of |n_sin_cos| is similar 
    to that of |n_arg| but applied in reverse. The values of |mp_m_spec_atan[k]| decrease slowly 
    enough that this loop is guaranteed to terminate before the (nonexistent) value 
    |mp_m_spec_atan[27]| would be required.To initialize the |randoms| table, we call the following 
    routine.

*/

static void mp_scaled_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;
        }
        mp->randoms[(i * 21) % 55].data.val = j;
    }
    /*tex We \quote {warm up} the array: */
    mp_new_randoms(mp);
    mp_new_randoms(mp);
    mp_new_randoms(mp);
}

static void mp_scaled_print(MP mp, mp_number *n)
{
    mp_print_e_str(mp, mp_string_scaled(mp, n->data.val));
}

static char *mp_scaled_tostring(MP mp, mp_number *n)
{
    return mp_string_scaled(mp, n->data.val);
}

static void mp_scaled_modulo(mp_number *a, mp_number *b)
{
    a->data.val = a->data.val % b->data.val;
}

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

/*tex

    To produce a uniform random number in the range |0 <= u < x| or |0 >= u > x| or |0 = u = x|, 
    given a |scaled| value~|x|, we proceed as shown here. Note that the call of |take_fraction| 
    will produce the values 0 and~|x| with about half the probability that it will produce any 
    other particular values between 0 and~|x|, because it rounds its answers.

*/

static void mp_scaled_m_unif_rand(MP mp, mp_number *ret, mp_number *x_orig)
{
    mp_number x, abs_x, u, y; /* |y| is trial value */
    mp_scaled_allocate_number(mp, &y, mp_fraction_type);
    mp_scaled_allocate_clone(mp, &x, mp_scaled_type, x_orig);
    mp_scaled_allocate_abs(mp, &abs_x, mp_scaled_type, &x);
    mp_scaled_allocate_number(mp, &u, mp_scaled_type);
    mp_next_random(mp, &u);
    /*|take_fraction (y, abs_x, u);|*/
    mp_scaled_take_fraction(mp, &y, &abs_x, &u);
    if (mp_scaled_equal(&y, &abs_x)) {
        /*|set_number_to_zero(*ret);|*/
        mp_scaled_clone(ret, &((math_data *)mp->math)->md_zero_t);
    } else if (mp_scaled_greater(&x, &((math_data *)mp->math)->md_zero_t)) {
        mp_scaled_clone(ret, &y);
    } else {
        mp_scaled_clone(ret, &y);
        mp_scaled_negate(ret);
    }
    mp_scaled_free_number(mp, &y);
    mp_scaled_free_number(mp, &abs_x);
    mp_scaled_free_number(mp, &x);
    mp_scaled_free_number(mp, &u);
}

/*tex

    Finally, a normal deviate with mean zero and unit standard deviation can readily be obtained 
    with the ratio method (Algorithm 3.4.1R in {\sl The Art of Computer Programming}).

*/

static void mp_scaled_m_norm_rand(MP mp, mp_number *ret)
{
    mp_number abs_x, u, r, la, xa;
    mp_scaled_allocate_number(mp, &la, mp_scaled_type);
    mp_scaled_allocate_number(mp, &xa, mp_scaled_type);
    mp_scaled_allocate_number(mp, &abs_x, mp_scaled_type);
    mp_scaled_allocate_number(mp, &u, mp_scaled_type);
    mp_scaled_allocate_number(mp, &r, mp_scaled_type);
    do {
        do {
            mp_number v;
            mp_scaled_allocate_number(mp, &v, mp_scaled_type);
            mp_next_random(mp, &v);
            mp_scaled_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t);
            mp_scaled_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v);
            mp_scaled_free_number(mp, &v);
            mp_next_random(mp, &u);
            mp_scaled_clone(&abs_x, &xa);
            mp_scaled_abs(&abs_x);
        } while (! mp_scaled_less(&abs_x, &u));
        mp_scaled_make_fraction(mp, &r, &xa, &u);
        mp_scaled_clone(&xa, &r);
        mp_scaled_m_log(mp, &la, &u);
        mp_scaled_set_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la);
    } while (mp_scaled_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0);
    mp_scaled_clone(ret, &xa);
    mp_scaled_free_number(mp, &r);
    mp_scaled_free_number(mp, &abs_x);
    mp_scaled_free_number(mp, &la);
    mp_scaled_free_number(mp, &xa);
    mp_scaled_free_number(mp, &u);
}

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

static void mp_scaled_free_math(MP mp)
{
    mp_scaled_free_number(mp, &(mp->math->md_epsilon_t));
    mp_scaled_free_number(mp, &(mp->math->md_inf_t));
    mp_scaled_free_number(mp, &(mp->math->md_negative_inf_t));
    mp_scaled_free_number(mp, &(mp->math->md_arc_tol_k));
    mp_scaled_free_number(mp, &(mp->math->md_three_sixty_deg_t));
    mp_scaled_free_number(mp, &(mp->math->md_one_eighty_deg_t));
    mp_scaled_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t));
    mp_scaled_free_number(mp, &(mp->math->md_fraction_one_t));
    mp_scaled_free_number(mp, &(mp->math->md_fraction_half_t));
    mp_scaled_free_number(mp, &(mp->math->md_fraction_three_t));
    mp_scaled_free_number(mp, &(mp->math->md_fraction_four_t));
    mp_scaled_free_number(mp, &(mp->math->md_zero_t));
    mp_scaled_free_number(mp, &(mp->math->md_half_unit_t));
    mp_scaled_free_number(mp, &(mp->math->md_three_quarter_unit_t));
    mp_scaled_free_number(mp, &(mp->math->md_unity_t));
    mp_scaled_free_number(mp, &(mp->math->md_two_t));
    mp_scaled_free_number(mp, &(mp->math->md_three_t));
    mp_scaled_free_number(mp, &(mp->math->md_one_third_inf_t));
    mp_scaled_free_number(mp, &(mp->math->md_warning_limit_t));
    mp_scaled_free_number(mp, &(mp->math->md_one_k));
    mp_scaled_free_number(mp, &(mp->math->md_sqrt_8_e_k));
    mp_scaled_free_number(mp, &(mp->math->md_twelve_ln_2_k));
    mp_scaled_free_number(mp, &(mp->math->md_coef_bound_k));
    mp_scaled_free_number(mp, &(mp->math->md_coef_bound_minus_1));
    mp_scaled_free_number(mp, &(mp->math->md_twelvebits_3));
    mp_scaled_free_number(mp, &(mp->math->md_twentysixbits_sqrt2_t));
    mp_scaled_free_number(mp, &(mp->math->md_twentyeightbits_d_t));
    mp_scaled_free_number(mp, &(mp->math->md_twentysevenbits_sqrt2_d_t));
    mp_scaled_free_number(mp, &(mp->math->md_fraction_threshold_t));
    mp_scaled_free_number(mp, &(mp->math->md_half_fraction_threshold_t));
    mp_scaled_free_number(mp, &(mp->math->md_scaled_threshold_t));
    mp_scaled_free_number(mp, &(mp->math->md_half_scaled_threshold_t));
    mp_scaled_free_number(mp, &(mp->math->md_near_zero_angle_t));
    mp_scaled_free_number(mp, &(mp->math->md_p_over_v_threshold_t));
    mp_scaled_free_number(mp, &(mp->math->md_equation_threshold_t));
    mp_memory_free(mp->math);
}

math_data *mp_initialize_scaled_math(MP mp)
{
    math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data));
    /* alloc */
    math->md_allocate        = mp_scaled_allocate_number;
    math->md_free            = mp_scaled_free_number;
    math->md_allocate_clone  = mp_scaled_allocate_clone;
    math->md_allocate_abs    = mp_scaled_allocate_abs;
    math->md_allocate_div    = mp_scaled_allocate_div;
    math->md_allocate_mul    = mp_scaled_allocate_mul;
    math->md_allocate_add    = mp_scaled_allocate_add;
    math->md_allocate_sub    = mp_scaled_allocate_sub;
    math->md_allocate_double = mp_scaled_allocate_double;
    /* precission */
    mp_scaled_allocate_number(mp, &math->md_precision_default, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_precision_max, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_precision_min, mp_scaled_type);
    /* here are the constants for |scaled| objects */
    mp_scaled_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_unity_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_two_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_three_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_zero_t, mp_scaled_type);
    /* |fractions| */
    mp_scaled_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
    /* |angles| */
    mp_scaled_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type);
    mp_scaled_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type);
    mp_scaled_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type);
    /* various approximations */
    mp_scaled_allocate_number(mp, &math->md_one_k, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type);
    /* thresholds */
    mp_scaled_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
    mp_scaled_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
    mp_scaled_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type);
    mp_scaled_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
    /* initializations */
    math->md_precision_default.data.val         = unity * 10;
    math->md_precision_max.data.val             = unity * 10;
    math->md_precision_min.data.val             = unity * 10;
    math->md_epsilon_t.data.val                 = 1;
    math->md_inf_t.data.val                     = EL_GORDO;
    math->md_negative_inf_t.data.val            = negative_EL_GORDO;
    math->md_one_third_inf_t.data.val           = one_third_EL_GORDO;
    math->md_warning_limit_t.data.val           = fraction_one;
    math->md_unity_t.data.val                   = unity;
    math->md_two_t.data.val                     = two;
    math->md_three_t.data.val                   = three;
    math->md_half_unit_t.data.val               = half_unit;
    math->md_three_quarter_unit_t.data.val      = three_quarter_unit;
    math->md_arc_tol_k.data.val                 = (unity/4096);
    math->md_fraction_one_t.data.val            = fraction_one;
    math->md_fraction_half_t.data.val           = fraction_half;
    math->md_fraction_three_t.data.val          = fraction_three;
    math->md_fraction_four_t.data.val           = fraction_four;
    math->md_three_sixty_deg_t.data.val         = three_sixty_deg;
    math->md_one_eighty_deg_t.data.val          = one_eighty_deg;
    math->md_negative_one_eighty_deg_t.data.val = negative_one_eighty_deg;
    math->md_one_k.data.val                     = 1024;
    math->md_sqrt_8_e_k.data.val                = 112429;         /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
    math->md_twelve_ln_2_k.data.val             = 139548960;      /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
    math->md_coef_bound_k.data.val              = coef_bound;
    math->md_coef_bound_minus_1.data.val        = coef_bound - 1;
    math->md_twelvebits_3.data.val              = 1365;           /* $1365\approx 2^{12}/3$ */
    math->md_twentysixbits_sqrt2_t.data.val     = 94906266;       /* $2^{26}\sqrt2\approx94906265.62$ */
    math->md_twentyeightbits_d_t.data.val       = 35596755;       /* $2^{28}d\approx35596754.69$ */
    math->md_twentysevenbits_sqrt2_d_t.data.val = 25170707;       /* $2^{27}\sqrt2\,d\approx25170706.63$ */
    math->md_fraction_threshold_t.data.val      = fraction_threshold;
    math->md_half_fraction_threshold_t.data.val = half_fraction_threshold;
    math->md_scaled_threshold_t.data.val        = scaled_threshold;
    math->md_half_scaled_threshold_t.data.val   = half_scaled_threshold;
    math->md_near_zero_angle_t.data.val         = near_zero_angle;
    math->md_p_over_v_threshold_t.data.val      = p_over_v_threshold;
    math->md_equation_threshold_t.data.val      = equation_threshold;
    /* functions */
    math->md_from_int                 = mp_scaled_set_from_int;
    math->md_from_boolean             = mp_scaled_set_from_boolean;
    math->md_from_scaled              = mp_scaled_set_from_scaled;
    math->md_from_double              = mp_scaled_set_from_double;
    math->md_from_addition            = mp_scaled_set_from_addition;
    math->md_half_from_addition       = mp_scaled_set_half_from_addition;
    math->md_from_subtraction         = mp_scaled_set_from_subtraction;
    math->md_half_from_subtraction    = mp_scaled_set_half_from_subtraction;
    math->md_from_of_the_way          = mp_scaled_set_from_of_the_way;
    math->md_from_div                 = mp_scaled_set_from_div;
    math->md_from_mul                 = mp_scaled_set_from_mul;
    math->md_from_int_div             = mp_scaled_set_from_int_div;
    math->md_from_int_mul             = mp_scaled_set_from_int_mul;
    math->md_negate                   = mp_scaled_negate;
    math->md_add                      = mp_scaled_add;
    math->md_subtract                 = mp_scaled_subtract;
    math->md_half                     = mp_scaled_half;
    math->md_do_double                = mp_scaled_double;
    math->md_abs                      = mp_scaled_abs;
    math->md_clone                    = mp_scaled_clone;
    math->md_negated_clone            = mp_scaled_negated_clone;
    math->md_abs_clone                = mp_scaled_abs_clone;
    math->md_swap                     = mp_scaled_swap;
    math->md_add_scaled               = mp_scaled_add_scaled;
    math->md_multiply_int             = mp_scaled_multiply_int;
    math->md_divide_int               = mp_scaled_divide_int;
    math->md_to_int                   = mp_scaled_to_int;
    math->md_to_boolean               = mp_scaled_to_boolean;
    math->md_to_scaled                = mp_scaled_to_scaled;
    math->md_to_double                = mp_scaled_to_double;
    math->md_odd                      = mp_scaled_odd;
    math->md_equal                    = mp_scaled_equal;
    math->md_less                     = mp_scaled_less;
    math->md_greater                  = mp_scaled_greater;
    math->md_non_equal_abs            = mp_scaled_non_equal_abs;
    math->md_round_unscaled           = mp_scaled_round_unscaled;
    math->md_floor_scaled             = mp_scaled_floor;
    math->md_fraction_to_round_scaled = mp_scaled_fraction_to_round_scaled;
    math->md_make_scaled              = mp_scaled_make_scaled;
    math->md_make_fraction            = mp_scaled_make_fraction;
    math->md_take_fraction            = mp_scaled_take_fraction;
    math->md_take_scaled              = mp_scaled_take_scaled;
    math->md_velocity                 = mp_scaled_velocity;
    math->md_n_arg                    = mp_scaled_n_arg;
    math->md_m_log                    = mp_scaled_m_log;
    math->md_m_exp                    = mp_scaled_m_exp;
    math->md_m_unif_rand              = mp_scaled_m_unif_rand;
    math->md_m_norm_rand              = mp_scaled_m_norm_rand;
    math->md_pyth_add                 = mp_scaled_pyth_add;
    math->md_pyth_sub                 = mp_scaled_pyth_sub;
    math->md_power_of                 = mp_scaled_power_of;
    math->md_fraction_to_scaled       = mp_scaled_fraction_to_scaled;
    math->md_scaled_to_fraction       = mp_scaled_scaled_to_fraction;
    math->md_scaled_to_angle          = mp_scaled_scaled_to_angle;
    math->md_angle_to_scaled          = mp_scaled_angle_to_scaled;
    math->md_init_randoms             = mp_scaled_init_randoms;
    math->md_sin_cos                  = mp_scaled_n_sin_cos;
    math->md_slow_add                 = mp_scaled_slow_add;
    math->md_sqrt                     = mp_scaled_sqrt;
    math->md_print                    = mp_scaled_print;
    math->md_tostring                 = mp_scaled_tostring;
    math->md_modulo                   = mp_scaled_modulo;
    math->md_ab_vs_cd                 = mp_scaled_ab_vs_cd;
    math->md_crossing_point           = mp_scaled_crossing_point;
    math->md_scan_numeric             = mp_scaled_scan_numeric_token;
    math->md_scan_fractional          = mp_scaled_scan_fractional_token;
    /* housekeeping */
    math->md_free_math                = mp_scaled_free_math;
    math->md_set_precision            = mp_scaled_set_precision;
    return math;
}