path: root/source/luametatex/source/mp/mpc/mpmath.c

/* This file is generated by "mtxrun --script "mtx-wtoc.lua" from the metapost cweb files. */


# include "mpconfig.h"
# include "mpmath.h"
# include "mpstrings.h"

# define  coef_bound               04525252525
# 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
# 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       05252525252
# 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
# define  fraction_half            0x08000000
# define  fraction_one             0x10000000
# define  fraction_two             0x20000000
# define  fraction_three           0x30000000
# define  fraction_four            0x40000000
# 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 int     mp_ab_vs_cd                        (mp_number *a, mp_number *b, mp_number *c, mp_number *d);
static void    mp_allocate_abs                    (MP mp, mp_number *n, mp_number_type t, mp_number *B);
static void    mp_allocate_clone                  (MP mp, mp_number *n, mp_number_type t, mp_number *B);
static void    mp_allocate_double                 (MP mp, mp_number *n, double v);
static void    mp_allocate_number                 (MP mp, mp_number *n, mp_number_type t);
static void    mp_crossing_point                  (MP mp, mp_number *ret, mp_number *a, mp_number *b, mp_number *c);
static void    mp_fraction_to_round_scaled        (mp_number *x);
static void    mp_free_number                     (MP mp, mp_number *n);
static void    mp_free_scaled_math                (MP mp);
static void    mp_init_randoms                    (MP mp, int seed);
static void    mp_m_exp                           (MP mp, mp_number *ret, mp_number *x_orig);
static void    mp_m_log                           (MP mp, mp_number *ret, mp_number *x_orig);
static void    mp_m_norm_rand                     (MP mp, mp_number *ret);
static void    mp_m_unif_rand                     (MP mp, mp_number *ret, mp_number *x_orig);
static int     mp_make_scaled                     (MP mp, int p, int q);
static void    mp_n_arg                           (MP mp, mp_number *ret, mp_number *x, mp_number *y);
static void    mp_n_sin_cos                       (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin);
static void    mp_number_abs                      (mp_number *A);
static void    mp_number_abs_clone                (mp_number *A, mp_number *B);
static void    mp_number_add                      (mp_number *A, mp_number *B);
static void    mp_number_add_scaled               (mp_number *A, int B);
static void    mp_number_angle_to_scaled          (mp_number *A);
static void    mp_number_clone                    (mp_number *A, mp_number *B);
static void    mp_number_divide_int               (mp_number *A, int B);
static void    mp_number_double                   (mp_number *A);
static int     mp_number_equal                    (mp_number *A, mp_number *B);
static void    mp_number_floor                    (mp_number *i);
static void    mp_number_fraction_to_scaled       (mp_number *A);
static int     mp_number_greater                  (mp_number *A, mp_number *B);
static void    mp_number_half                     (mp_number *A);
static int     mp_number_less                     (mp_number *A, mp_number *B);
static void    mp_number_make_fraction            (MP mp, mp_number *r, mp_number *p, mp_number *q);
static void    mp_number_make_scaled              (MP mp, mp_number *r, mp_number *p, mp_number *q);
static void    mp_number_modulo                   (mp_number *a, mp_number *b);
static void    mp_number_multiply_int             (mp_number *A, int B);
static void    mp_number_negate                   (mp_number *A);
static void    mp_number_negated_clone            (mp_number *A, mp_number *B);
static int     mp_number_nonequalabs              (mp_number *A, mp_number *B);
static int     mp_number_odd                      (mp_number *A);
static void    mp_number_scaled_to_angle          (mp_number *A);
static void    mp_number_scaled_to_fraction       (mp_number *A);
static void    mp_number_subtract                 (mp_number *A, mp_number *B);
static void    mp_number_swap                     (mp_number *A, mp_number *B);
static void    mp_number_take_fraction            (MP mp, mp_number *r, mp_number *p, mp_number *q);
static void    mp_number_take_scaled              (MP mp, mp_number *r, mp_number *p, mp_number *q);
static int     mp_number_to_boolean               (mp_number *A);
static double  mp_number_to_double                (mp_number *A);
static int     mp_number_to_int                   (mp_number *A);
static int     mp_number_to_scaled                (mp_number *A);
static void    mp_power_of                        (MP mp, mp_number *r, mp_number *a, mp_number *b);
static void    mp_print_number                    (MP mp, mp_number *n);
static void    mp_pyth_add                        (MP mp, mp_number *r, mp_number *a, mp_number *b);
static void    mp_pyth_sub                        (MP mp, mp_number *r, mp_number *a, mp_number *b);
static int     mp_round_decimals                  (MP mp, unsigned char *b, int k);
static int     mp_round_unscaled                  (mp_number *x_orig);
static void    mp_scaled_set_precision            (MP mp);
static void    mp_scan_fractional_token           (MP mp, int n);
static void    mp_scan_numeric_token              (MP mp, int n);
static void    mp_set_number_from_addition        (mp_number *A, mp_number *B, mp_number *C);
static void    mp_set_number_from_boolean         (mp_number *A, int B);
static void    mp_set_number_from_div             (mp_number *A, mp_number *B, mp_number *C);
static void    mp_set_number_from_double          (mp_number *A, double B);
static void    mp_set_number_from_int             (mp_number *A, int B);
static void    mp_set_number_from_int_div         (mp_number *A, mp_number *B, int C);
static void    mp_set_number_from_int_mul         (mp_number *A, mp_number *B, int C);
static void    mp_set_number_from_mul             (mp_number *A, mp_number *B, mp_number *C);
static void    mp_set_number_from_of_the_way      (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C);
static void    mp_set_number_from_scaled          (mp_number *A, int B);
static void    mp_set_number_from_subtraction     (mp_number *A, mp_number *B, mp_number *C);
static void    mp_set_number_half_from_addition   (mp_number *A, mp_number *B, mp_number *C);
static void    mp_set_number_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C);
static void    mp_slow_add                        (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig);
static void    mp_square_rt                       (MP mp, mp_number *ret, mp_number *x_orig);
static int     mp_take_fraction                   (MP mp, int q, int f);
static int     mp_take_scaled                     (MP mp, int q, int f);
static void    mp_velocity                        (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t);
static void    mp_wrapup_numeric_token            (MP mp, int n, int f);
static char   *mp_number_tostring                 (MP mp, mp_number *n);
static char   *mp_string_scaled                   (MP mp, int s);
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
};

math_data *mp_initialize_scaled_math(MP mp)
{
    math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data));
    math->md_allocate        = mp_allocate_number;
    math->md_free            = mp_free_number;
    math->md_allocate_clone  = mp_allocate_clone;
    math->md_allocate_abs    = mp_allocate_abs;
    math->md_allocate_double = mp_allocate_double;
    mp_allocate_number(mp, &math->md_precision_default, mp_scaled_type);
    mp_allocate_number(mp, &math->md_precision_max, mp_scaled_type);
    mp_allocate_number(mp, &math->md_precision_min, mp_scaled_type);
    mp_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_unity_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_two_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_three_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_zero_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type);
    mp_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type);
    mp_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type);
    mp_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type);
    mp_allocate_number(mp, &math->md_one_k, mp_scaled_type);
    mp_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type);
    mp_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type);
    mp_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type);
    mp_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type);
    mp_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type);
    mp_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
    mp_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
    mp_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type);
    mp_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
    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;
    math->md_twelve_ln_2_k.data.val             = 139548960;
    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;
    math->md_twentysixbits_sqrt2_t.data.val     = 94906266;
    math->md_twentyeightbits_d_t.data.val       = 35596755;
    math->md_twentysevenbits_sqrt2_d_t.data.val = 25170707;
    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;
    math->md_from_int                 = mp_set_number_from_int;
    math->md_from_boolean             = mp_set_number_from_boolean;
    math->md_from_scaled              = mp_set_number_from_scaled;
    math->md_from_double              = mp_set_number_from_double;
    math->md_from_addition            = mp_set_number_from_addition;
    math->md_half_from_addition       = mp_set_number_half_from_addition;
    math->md_from_subtraction         = mp_set_number_from_subtraction;
    math->md_half_from_subtraction    = mp_set_number_half_from_subtraction;
    math->md_from_oftheway            = mp_set_number_from_of_the_way;
    math->md_from_div                 = mp_set_number_from_div;
    math->md_from_mul                 = mp_set_number_from_mul;
    math->md_from_int_div             = mp_set_number_from_int_div;
    math->md_from_int_mul             = mp_set_number_from_int_mul;
    math->md_negate                   = mp_number_negate;
    math->md_add                      = mp_number_add;
    math->md_subtract                 = mp_number_subtract;
    math->md_half                     = mp_number_half;
    math->md_do_double                = mp_number_double;
    math->md_abs                      = mp_number_abs;
    math->md_clone                    = mp_number_clone;
    math->md_negated_clone            = mp_number_negated_clone;
    math->md_abs_clone                = mp_number_abs_clone;
    math->md_swap                     = mp_number_swap;
    math->md_add_scaled               = mp_number_add_scaled;
    math->md_multiply_int             = mp_number_multiply_int;
    math->md_divide_int               = mp_number_divide_int;
    math->md_to_int                   = mp_number_to_int;
    math->md_to_boolean               = mp_number_to_boolean;
    math->md_to_scaled                = mp_number_to_scaled;
    math->md_to_double                = mp_number_to_double;
    math->md_odd                      = mp_number_odd;
    math->md_equal                    = mp_number_equal;
    math->md_less                     = mp_number_less;
    math->md_greater                  = mp_number_greater;
    math->md_nonequalabs              = mp_number_nonequalabs;
    math->md_round_unscaled           = mp_round_unscaled;
    math->md_floor_scaled             = mp_number_floor;
    math->md_fraction_to_round_scaled = mp_fraction_to_round_scaled;
    math->md_make_scaled              = mp_number_make_scaled;
    math->md_make_fraction            = mp_number_make_fraction;
    math->md_take_fraction            = mp_number_take_fraction;
    math->md_take_scaled              = mp_number_take_scaled;
    math->md_velocity                 = mp_velocity;
    math->md_n_arg                    = mp_n_arg;
    math->md_m_log                    = mp_m_log;
    math->md_m_exp                    = mp_m_exp;
    math->md_m_unif_rand              = mp_m_unif_rand;
    math->md_m_norm_rand              = mp_m_norm_rand;
    math->md_pyth_add                 = mp_pyth_add;
    math->md_pyth_sub                 = mp_pyth_sub;
    math->md_power_of                 = mp_power_of;
    math->md_fraction_to_scaled       = mp_number_fraction_to_scaled;
    math->md_scaled_to_fraction       = mp_number_scaled_to_fraction;
    math->md_scaled_to_angle          = mp_number_scaled_to_angle;
    math->md_angle_to_scaled          = mp_number_angle_to_scaled;
    math->md_init_randoms             = mp_init_randoms;
    math->md_sin_cos                  = mp_n_sin_cos;
    math->md_slow_add                 = mp_slow_add;
    math->md_sqrt                     = mp_square_rt;
    math->md_print                    = mp_print_number;
    math->md_tostring                 = mp_number_tostring;
    math->md_modulo                   = mp_number_modulo;
    math->md_ab_vs_cd                 = mp_ab_vs_cd;
    math->md_crossing_point           = mp_crossing_point;
    math->md_scan_numeric             = mp_scan_numeric_token;
    math->md_scan_fractional          = mp_scan_fractional_token;
    math->md_free_math                = mp_free_scaled_math;
    math->md_set_precision            = mp_scaled_set_precision;
    return math;
}

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

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

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

void mp_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;
}

void mp_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);
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

void mp_set_number_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_take_fraction(mp, (B->data.val - C->data.val), t->data.val);
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

void mp_number_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);
    }
}

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

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

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

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

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

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

int mp_number_odd(mp_number *A)
{
    return odd(mp_round_unscaled(A));
}

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

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

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

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

static char *mp_string_scaled (MP mp, int s)
{
    (void) mp;
    static char scaled_string[32];
    int i = 0;
    if (s < 0) {
        scaled_string[i++] = '-';
        s = -s;
    }
    mp_snprintf ((scaled_string+i), 12, "%d", (int) (s / unity));
    while (*(scaled_string+i)) {
        i++;
    }
    s = 10 * (s % unity) + 5;
    if (s != 5) {
        int delta = 10;
        scaled_string[i++] = '.';
        do {
            if (delta > unity) {
                s = s + 0100000 - (delta / 2);
            }
            scaled_string[i++] = '0' + (s / unity);
            s = 10 * (s % unity);
            delta = delta * 10;
        } while (s > delta);
    }
    scaled_string[i] = '\0';
    return scaled_string;
}

void mp_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->arith_error = 1;
            ret->data.val =  EL_GORDO;
        }
    } else if (-y <= EL_GORDO + x) {
        ret->data.val = x + y;
    } else {
        mp->arith_error = 1;
        ret->data.val = negative_EL_GORDO;
    }
}

static int mp_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->arith_error = 1;
                return EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? -q : q) & 077777) * (((i & 037777) << 1) - 1) & 04000) != 0) {
                    --i;
                }
                return i;
            }
        } else {
            d -= 0.5;
            if (d <= -TWEXP31) {
                mp->arith_error = 1;
                return -negative_EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? q : -q) & 077777) * (((i & 037777) << 1) + 1) & 04000) != 0) {
                    ++i;
                }
                return i;
            }
        }
    }
}

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

int mp_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 & 077777) * (q & 077777)) & 040000) == 0) {
                mp->arith_error = 1;
            }
            return EL_GORDO;
        } else {
            int i = (int) d;
            if (d == (double) i && (((p & 077777) * (q & 077777)) & 040000) != 0) {
                --i;
            }
            return i;
        }
    } else {
        d -= 0.5;
        if (d <= -TWEXP31) {
            if (d != -TWEXP31 || ((-(p & 077777) * (q & 077777)) & 040000) == 0) {
                mp->arith_error = 1;
            }
            return -negative_EL_GORDO;
        } else {
            int i = (int) d;
            if (d == (double) i && ((-(p & 077777) * (q & 077777)) & 040000) != 0) {
                ++i;
            }
            return i;
        }
    }
}

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

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

void mp_number_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);
}

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->arith_error = 1;
                return EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? -q : q) & 077777) * (((i & 037777) << 1) - 1) & 04000) != 0) {
                    --i;
                }
                return i;
            }
        } else {
            d -= 0.5;
            if (d <= -TWEXP31) {
                mp->arith_error = 1;
                return -negative_EL_GORDO;
            } else {
                int i = (int) d;
                if (d == (double) i && (((q > 0 ? q : -q) & 077777) * (((i & 037777) << 1) + 1) & 04000) != 0) {
                    ++i;
                }
                return i;
            }
        }
    }
}

void mp_number_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);
}

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

static void mp_wrapup_numeric_token (MP mp, int n, int f)
{
    if (n < 32768) {
        int mod = (n * unity + f);
        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];
                mp_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);
}

void mp_scan_fractional_token (MP mp, int n)
{
    int f;
    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_wrapup_numeric_token(mp, n, f);
}

void mp_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_wrapup_numeric_token(mp, n, 0);
    } else {
        mp->cur_input.loc_field++;
        mp_scan_fractional_token(mp, n);
    }
}

void mp_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_take_fraction(mp, st->data.val - (sf->data.val / 16), sf->data.val - (st->data.val / 16));
    acc = mp_take_fraction(mp, acc, ct->data.val - cf->data.val);
    num = fraction_two + mp_take_fraction(mp, acc, 379625062);
    denom = fraction_three + mp_take_fraction(mp, ct->data.val, 497706707) + mp_take_fraction (mp, cf->data.val, 307599661);
    if (t->data.val != unity) {
        num = mp_make_scaled (mp, num, t->data.val);
    }
    if (num / 4 >= denom) {
        ret->data.val = fraction_four;
    } else {
        ret->data.val = mp_make_fraction(mp, num, denom);
    }
}

static int mp_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;
            }
        }
    }
}

static void mp_crossing_point (MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc)
{
    int x, xx, x0, x1, x2;
    int a = aa->data.val;
    int b = bb->data.val;
    int c = cc->data.val;
    int d;
    (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;
        }
    }
    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;
}

int mp_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));
    }
}

void mp_number_floor(mp_number *i)
{
    i->data.val = i->data.val&-65536;
}

void mp_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))));
}

void mp_square_rt (MP mp, mp_number *ret, mp_number *x_orig)
{
    int x = x_orig->data.val;
    if (x <= 0) {
        if (x < 0) {
            char msg[256];
            mp_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;
        int y;
        int q = 2;
        while (x < fraction_two) {
            k--;
            x = x + x + x + x;
        }
        if (x < fraction_four)
            y = 0;
        else {
            x = x - fraction_four;
            y = 1;
        }
        do {
            x += x;
            y += y;
            if (x >= 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);
    }
}

void mp_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;
    }
    if (b > 0) {
        int big;
        if (a < fraction_two) {
            big = 0;
        } else {
            a = a / 4;
            b = b / 4;
            big = 1;
        }
        while (1) {
            int r = mp_make_fraction(mp, b, a);
            r = mp_take_fraction(mp, r, r);
            if (r == 0) {
                break;
            } else {
                r = mp_make_fraction(mp, r, fraction_four + r);
                a = a + mp_take_fraction(mp, a + a, r);
                b = mp_take_fraction(mp, b, r);
            }
        }
        if (big) {
            if (a < fraction_two) {
                a = a + a + a + a;
            } else {
                mp->arith_error = 1;
                a = EL_GORDO;
            }
        }
    }
    ret->data.val = a;
}

void mp_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) {
        if (a < b) {
            char msg[256];
            char *astr = mp_strdup(mp_string_scaled(mp, a));
            mp_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;
        if (a < fraction_four) {
            big = 0;
        } else {
            a = (int) a/2;
            b = (int) b/2;
            big = 1;
        }
        while (1) {
            int r = mp_make_fraction(mp, b, a);
            r = mp_take_fraction(mp, r, r);
            if (r == 0) {
                break;
            } else {
                r = mp_make_fraction(mp, r, fraction_four - r);
                a = a - mp_take_fraction(mp, a + a, r);
                b = mp_take_fraction(mp, b, r);
            }
        }
        if (big) {
            a *= 2;
        }
    }
    ret->data.val = a;
}

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

void mp_m_log (MP mp, mp_number *ret, mp_number *x_orig)
{
    int x = x_orig->data.val;
    if (x <= 0) {
        {
            char msg[256];
            mp_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;
        int y = 1302456956 + 4 - 100;
        int z = 27595 + 6553600;
        while (x < fraction_four) {
            x = 2*x;
            y -= 93032639;
            z -= 48782;
        }
        y = y + (z / unity);
        while (x > fraction_four + 4) {
            {
                z = ((x - 1) / two_to_the (k)) + 1;
                while (x < fraction_four + z) {
                    z = (z + 1)/2;
                    k++;
                };
                y += mp_m_spec_log[k];
                x -= z;
            }
        }
        ret->data.val = (y / 8);
    }
}

void mp_m_exp (MP mp, mp_number *ret, mp_number *x_orig)
{
    int y, z;
    int x = x_orig->data.val;
    if (x > 174436200) {
        mp->arith_error = 1;
        ret->data.val = EL_GORDO;
    } else if (x < -197694359) {
        ret->data.val = 0;
    } else {
        if (x <= 0) {
            z = -8 * x;
            y = 04000000;
        } else {
            if (x <= 127919879) {
                z = 1023359037 - 8 * x;
            } else {
                z = 8 * (174436200 - x);
            }
            y = EL_GORDO;
        }
        {
            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;
        }
    }
}

void mp_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
{
    int z;
    int t;
    int k;
    int octant;
    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) {
        t = y;
        y = x;
        x = t;
        octant = octant + switch_x_and_y;
    }
    if (x == 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->type = mp_angle_type;
        while (x >= fraction_two) {
            x = x/2;
            y = y/2;
        }
        z = 0;
        if (y > 0) {
            while (x < fraction_one) {
                x += x;
                y += y;
            };
            k = 0;
            do {
                y += y;
                k++;
                if (y > x) {
                    z = z + mp_m_spec_atan[k];
                    t = x;
                    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);
        }

        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;
        }
    }
}

void mp_n_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin)
{
    int k;
    int q;
    int x, y, t;
    int z = z_orig->data.val;
    mp_number x_n, y_n, ret;
    mp_allocate_number(mp, &ret, mp_scaled_type);
    mp_allocate_number(mp, &x_n, mp_scaled_type);
    mp_allocate_number(mp, &y_n, mp_scaled_type);
    while (z < 0) {
        z = z + three_sixty_deg;
    }
    z = 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;
    }
    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) {
        y = 0;
    }
    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_pyth_add(mp, &ret, &x_n, &y_n);
    n_cos->data.val = mp_make_fraction(mp, x, ret.data.val);
    n_sin->data.val = mp_make_fraction(mp, y, ret.data.val);
    mp_free_number(mp, &ret);
    mp_free_number(mp, &x_n);
    mp_free_number(mp, &y_n);
}

void mp_init_randoms (MP mp, int seed)
{
    int k = 1;
    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;
    }
    mp_new_randoms(mp);
    mp_new_randoms(mp);
    mp_new_randoms(mp);
}

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

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

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

static void mp_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig)
{
    mp_number x, abs_x, u, y;
    mp_allocate_number(mp, &y, mp_fraction_type);
    mp_allocate_clone(mp, &x, mp_scaled_type, x_orig);
    mp_allocate_abs(mp, &abs_x, mp_scaled_type, &x);
    mp_allocate_number(mp, &u, mp_scaled_type);
    mp_next_random(mp, &u);
    mp_number_take_fraction(mp, &y, &abs_x, &u);
    if (mp_number_equal(&y, &abs_x)) {
        mp_number_clone(ret, &((math_data *)mp->math)->md_zero_t);
    } else if (mp_number_greater(&x, &((math_data *)mp->math)->md_zero_t)) {
        mp_number_clone(ret, &y);
    } else {
        mp_number_clone(ret, &y);
        mp_number_negate(ret);
    }
    mp_free_number(mp, &y);
    mp_free_number(mp, &abs_x);
    mp_free_number(mp, &x);
    mp_free_number(mp, &u);
}

static void mp_m_norm_rand (MP mp, mp_number *ret)
{
    mp_number abs_x, u, r, la, xa;
    mp_allocate_number(mp, &la, mp_scaled_type);
    mp_allocate_number(mp, &xa, mp_scaled_type);
    mp_allocate_number(mp, &abs_x, mp_scaled_type);
    mp_allocate_number(mp, &u, mp_scaled_type);
    mp_allocate_number(mp, &r, mp_scaled_type);
    do {
        do {
            mp_number v;
            mp_allocate_number(mp, &v, mp_scaled_type);
            mp_next_random(mp, &v);
            mp_number_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t);
            mp_number_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v);
            mp_free_number(mp, &v);
            mp_next_random(mp, &u);
            mp_number_clone(&abs_x, &xa);
            mp_number_abs(&abs_x);
        } while (! mp_number_less(&abs_x, &u));
        mp_number_make_fraction(mp, &r, &xa, &u);
        mp_number_clone(&xa, &r);
        mp_m_log(mp, &la, &u);
        mp_set_number_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la);
    } while (mp_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0);
    mp_number_clone(ret, &xa);
    mp_free_number(mp, &r);
    mp_free_number(mp, &abs_x);
    mp_free_number(mp, &la);
    mp_free_number(mp, &xa);
    mp_free_number(mp, &u);
}