Source Browser ?

1/* 
2    This file was generated by "mtxrun --script "mtx-wtoc.lua" from the metapost cweb 
3    files but now maintained as C file.
4*/
5
6# include "mpmathscaled.h"
7# include "mpstrings.h"
8
9# define  coef_bound               0x25555555
10# define  fraction_threshold       2685
11# define  half_fraction_threshold  1342
12# define  scaled_threshold         8
13# define  half_scaled_threshold    4
14# define  near_zero_angle          26844
15# define  p_over_v_threshold       0x80000
16# define  equation_threshold       64
17
18/*tex
19    Fixed-point arithmetic is done on {\sl scaled integers} that are multiples of $2^{-16}$. In 
20    other words, a binary point is assumed to be sixteen bit positions from the right end of a 
21    binary computer word.
22*/
23
24# define  unity                    0x10000
25# define  two                      (2*unity)
26# define  three                    (3*unity)
27# define  half_unit                (unity/2)
28# define  three_quarter_unit       (3*(unity/4))
29# define  EL_GORDO                 0x7FFFFFFF
30# define  negative_EL_GORDO        (-EL_GORDO)
31# define  one_third_EL_GORDO       0x2AAAAAAA
32
33/*tex We need these preprocessor values */
34
35# define  TWEXP31                  2147483648.0
36# define  TWEXP28                  268435456.0
37# define  TWEXP16                  65536.0
38# define  TWEXP_16                 (1.0/65536.0)
39# define  TWEXP_28                 (1.0/268435456.0)
40# define  no_crossing              (fraction_one + 1)
41# define  one_crossing             fraction_one
42# define  zero_crossing            0
43
44/*tex
45    The |scaled| quantities in \MP\ programs are generally supposed to be less than $2^{12}$ in 
46    absolute value, so \MP\ does much of its internal arithmetic with 28~significant bits of 
47    precision. A |fraction| denotes a scaled integer whose binary point is assumed to be 28 bit 
48    positions from the right.
49*/
50
51# define  fraction_half            0x08000000
52# define  fraction_one             0x10000000
53# define  fraction_two             0x20000000
54# define  fraction_three           0x30000000
55# define  fraction_four            0x40000000
56
57/*tex
58    Octants are represented in a \quote {Gray code,} since that turns out to be computationally 
59    simplest.
60*/
61
62# define  negate_x                 1
63# define  negate_y                 2
64# define  switch_x_and_y           4
65# define  first_octant             1
66# define  second_octant            (first_octant + switch_x_and_y)
67# define  third_octant             (first_octant + switch_x_and_y + negate_x)
68# define  fourth_octant            (first_octant                  + negate_x)
69# define  fifth_octant             (first_octant                  + negate_x + negate_y)
70# define  sixth_octant             (first_octant + switch_x_and_y + negate_x + negate_y)
71# define  seventh_octant           (first_octant + switch_x_and_y + negate_y)
72# define  eighth_octant            (first_octant                  + negate_y)
73# define  forty_five_deg           0x02D00000
74# define  ninety_deg               0x05A00000
75# define  one_eighty_deg           0x0B400000
76# define  negative_one_eighty_deg -0x0B400000
77# define  three_sixty_deg          0x16800000
78# define  odd(A)                   (abs(A)%2==1)
79# define  two_to_the(A)            (1<<(unsigned)(A))
80# define  set_cur_cmd(A)           mp->cur_mod_->command = (A)
81# define  set_cur_mod(A)           mp->cur_mod_->data.n.data.val = (A)
82
83static const int mp_m_spec_log[29] = {
84    0, 93032640, 38612034, 17922280, 8662214, 4261238, 2113709, 1052693, 525315,
85    262400, 131136, 65552, 32772, 16385, 8192, 4096, 2048, 1024, 512, 256, 128,
86    64, 32, 16, 8, 4, 2, 1, 1
87};
88
89static const int mp_m_spec_atan[27] = {
90    0, 27855475, 14718068, 7471121, 3750058, 1876857, 938658, 469357, 234682,
91    117342, 58671, 29335, 14668, 7334, 3667, 1833, 917, 458, 229, 115, 57, 29,
92    14, 7, 4, 2, 1
93};
94
95static int mp_scaled_aux_take_fraction (MP mp, int p, int q);
96static int mp_scaled_aux_make_fraction (MP mp, int p, int q);
97static int mp_scaled_round_unscaled    (mp_number *x_orig);
98
99static void mp_scaled_allocate_number(MP mp, mp_number *n, mp_number_type t)
100{
101    (void) mp;
102    n->data.val = 0;
103    n->type = t;
104}
105
106static void mp_scaled_allocate_clone(MP mp, mp_number *n, mp_number_type t, mp_number *v)
107{
108    (void) mp;
109    n->type = t;
110    n->data.val = v->data.val;
111}
112
113static void mp_scaled_allocate_abs(MP mp, mp_number *n, mp_number_type t, mp_number *v)
114{
115    (void) mp;
116    n->type = t;
117    n->data.val = abs(v->data.val);
118}
119
120static void mp_scaled_allocate_double(MP mp, mp_number *n, double v)
121{
122    (void) mp;
123    n->type = mp_scaled_type;
124    n->data.val = (int) (v * 65536.0);
125}
126
127static void mp_scaled_free_number(MP mp, mp_number *n)
128{
129    (void) mp;
130    n->type = mp_nan_type;
131}
132
133static void mp_scaled_set_from_int(mp_number *A, int B)
134{
135    A->data.val = B * 65536;
136}
137
138static void mp_scaled_set_from_boolean(mp_number *A, int B)
139{
140    A->data.val = B;
141}
142
143static void mp_scaled_set_from_scaled(mp_number *A, int B)
144{
145    A->data.val = B;
146}
147
148static void mp_scaled_set_from_double(mp_number *A, double B)
149{
150    A->data.val = (int) (B * 65536.0);
151}
152
153static void mp_scaled_set_from_addition(mp_number *A, mp_number *B, mp_number *C)
154{
155    A->data.val = B->data.val + C->data.val;
156}
157
158static void mp_scaled_set_half_from_addition(mp_number *A, mp_number *B, mp_number *C)
159{
160    A->data.val = (B->data.val + C->data.val) / 2;
161}
162
163static void mp_scaled_set_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
164{
165    A->data.val = B->data.val - C->data.val;
166}
167
168static void mp_scaled_set_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
169{
170    A->data.val = (B->data.val - C->data.val) / 2;
171}
172
173static void mp_scaled_set_from_div(mp_number *A, mp_number *B, mp_number *C)
174{
175    A->data.val = B->data.val / C->data.val;
176}
177
178static void mp_scaled_set_from_mul(mp_number *A, mp_number *B, mp_number *C)
179{
180    A->data.val = B->data.val * C->data.val;
181}
182
183static void mp_scaled_set_from_int_div(mp_number *A, mp_number *B, int C)
184{
185    A->data.val = B->data.val / C;
186}
187
188static void mp_scaled_set_from_int_mul(mp_number *A, mp_number *B, int C)
189{
190    A->data.val = B->data.val * C;
191}
192
193static void mp_scaled_set_from_of_the_way(MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C)
194{
195    (void) mp;
196    A->data.val = B->data.val - mp_scaled_aux_take_fraction(mp, (B->data.val - C->data.val), t->data.val);
197}
198
199static void mp_scaled_negate(mp_number *A)
200{
201    A->data.val = -A->data.val;
202}
203
204static void mp_scaled_add(mp_number *A, mp_number *B)
205{
206    A->data.val = A->data.val + B->data.val;
207}
208
209static void mp_scaled_subtract(mp_number *A, mp_number *B)
210{
211    A->data.val = A->data.val - B->data.val;
212}
213
214static void mp_scaled_half(mp_number *A)
215{
216    A->data.val = A->data.val / 2;
217}
218
219static void mp_scaled_double(mp_number *A)
220{
221    A->data.val = A->data.val + A->data.val;
222}
223
224static void mp_scaled_add_scaled(mp_number *A, int B)
225{
226    /* also for negative B */
227    A->data.val = A->data.val + B;
228}
229
230static void mp_scaled_multiply_int(mp_number *A, int B)
231{
232    A->data.val = B * A->data.val;
233}
234
235static void mp_scaled_divide_int(mp_number *A, int B)
236{
237    A->data.val = A->data.val / B;
238}
239
240static void mp_scaled_abs(mp_number *A)
241{
242    A->data.val = abs(A->data.val);
243}
244
245static void mp_scaled_clone(mp_number *A, mp_number *B)
246{
247    A->data.val = B->data.val;
248}
249
250static void mp_scaled_negated_clone(mp_number *A, mp_number *B)
251{
252    A->data.val = -B->data.val;
253}
254
255static void mp_scaled_abs_clone(mp_number *A, mp_number *B)
256{
257    A->data.val = abs(B->data.val);
258}
259
260static void mp_scaled_swap(mp_number *A, mp_number *B)
261{
262    int swap_tmp = A->data.val;
263    A->data.val = B->data.val;
264    B->data.val = swap_tmp;
265}
266
267static void mp_scaled_fraction_to_scaled(mp_number *A)
268{
269    A->type = mp_scaled_type;
270    A->data.val = A->data.val / 4096;
271}
272
273static void mp_scaled_angle_to_scaled(mp_number *A)
274{
275    A->type = mp_scaled_type;
276    if (A->data.val >= 0) {
277        A->data.val =   ( A->data.val + 8) / 16;
278    } else {
279        A->data.val = -((-A->data.val + 8) / 16);
280    }
281}
282
283static void mp_scaled_scaled_to_fraction(mp_number *A)
284{
285    A->type = mp_fraction_type;
286    A->data.val = A->data.val * 4096;
287}
288
289static void mp_scaled_scaled_to_angle(mp_number *A)
290{
291    A->type = mp_angle_type;
292    A->data.val = A->data.val * 16;
293}
294
295static int mp_scaled_to_int(mp_number *A)
296{
297    return A->data.val;
298}
299
300static int mp_scaled_to_scaled(mp_number *A)
301{
302    return A->data.val;
303}
304
305static int mp_scaled_to_boolean(mp_number *A)
306{
307    return A->data.val;
308}
309
310static double mp_scaled_to_double(mp_number *A)
311{
312    return A->data.val / 65536.0;
313}
314
315static int mp_scaled_odd(mp_number *A)
316{
317    return odd(mp_scaled_round_unscaled(A));
318}
319
320static int mp_scaled_equal(mp_number *A, mp_number *B) {
321    return A->data.val == B->data.val;
322}
323
324static int mp_scaled_greater(mp_number *A, mp_number *B)
325{
326    return A->data.val > B->data.val;
327}
328
329static int mp_scaled_less(mp_number *A, mp_number *B)
330{
331    return A->data.val < B->data.val;
332}
333
334static int mp_scaled_non_equal_abs(mp_number *A, mp_number *B)
335{
336    return abs(A->data.val) != abs(B->data.val);
337}
338
339/*tex 
340
341    One of \MP's most common operations is the calculation of $\lfloor {a + b \over 2}\rfloor$, 
342    the midpoint of two given integers |a| and~|b|. The most decent way to do this is to write 
343    |(a + b) / 2|; but on many machines it is more efficient to calculate |(a + b) >> 1|.
344
345    Therefore the midpoint operation will always be denoted by |half(a + b)| in this program. If 
346    \MP\ is being implemented with languages that permit binary shifting, the |half| macro should 
347    be changed to make this operation as efficient as possible. Since some systems have shift 
348    operators that can only be trusted to work on positive numbers, there is also a macro |halfp| 
349    that is used only when the quantity being halved is known to be positive or zero.
350
351    Here is a procedure analogous to |print_int|. If the output of this procedure is subsequently 
352    read by \MP\ and converted by the |round_decimals| routine above, it turns out that the 
353    original value will be reproduced exactly. A decimal point is printed only if the value is not 
354    an integer. If there is more than one way to print the result with the optimum number of digits
355    following the decimal point, the closest possible value is given.
356
357    The invariant relation in the |repeat| loop is that a sequence of decimal digits yet to be 
358    printed will yield the original number if and only if they form a fraction~$f$ in the range 
359    $s - \delta \L 10 \cdot 2^{16} f < s$. We can stop if and only if $f = 0$ satisfies this 
360    condition; the loop will terminate before $s$ can possibly become zero. We round to five 
361    digits.
362
363*/
364
365/*tex Todo: check if we can replace this: */
366
367static char *mp_string_scaled(MP mp, int s)
368{
369    static char scaled_string[32];
370    int i = 0;
371    (void) mp;
372    if (s < 0) {
373        scaled_string[i++] = '-';
374        s = -s;
375    }
376    /*tex The integer part: */
377    snprintf ((scaled_string+i), 12, "%d", (int) (s / unity));
378    while (*(scaled_string+i)) {
379        i++;
380    }
381    s = 10 * (s % unity) + 5;
382    if (s != 5) {
383        /*tex The amount of allowable inaccuracy, scaled: */
384        int delta = 10;
385        scaled_string[i++] = '.';
386        do {
387            /*tex Round the final digit: */
388            if (delta > unity) {
389                s = s + 0x8000 - (delta / 2);
390            }
391            scaled_string[i++] = (char) ((s / unity) + '0');
392            s = 10 * (s % unity);
393            delta = delta * 10;
394        } while (s > delta);
395    }
396    scaled_string[i] = '\0';
397    return scaled_string;
398}
399
400/*tex
401    Addition is not always checked to make sure that it doesn't overflow, but in places where 
402    overflow isn't too unlikely the |slow_add| routine is used.
403*/
404
405static void mp_scaled_slow_add(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
406{
407    int x = x_orig->data.val;
408    int y = y_orig->data.val;
409    if (x >= 0) {
410        if (y <= EL_GORDO - x) {
411            ret->data.val = x + y;
412        } else {
413            mp->arith_error = 1;
414            ret->data.val =  EL_GORDO;
415        }
416    } else if (-y <= EL_GORDO + x) {
417        ret->data.val = x + y;
418    } else {
419        mp->arith_error = 1;
420        ret->data.val = negative_EL_GORDO;
421    }
422}
423
424/*tex
425
426    The |make_fraction| routine produces the |fraction| equivalent of |p/q|, given integers |p| 
427    and~|q|; it computes the integer $f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are
428    positive. If |p| and |q| are both of the same scaled type |t|, the \quote {type relation}
429    |make_fraction(t,t) = fraction| is valid; and it's also possible to use the subroutine \quote 
430    {backwards,} using the relation |make_fraction(t,fraction) = t| between scaled types.
431
432    If the result would have magnitude $2^{31}$ or more, |make_fraction| sets |arith_error := 
433    true|. Most of \MP's internal computations have been designed to avoid this sort of error.
434
435    If this subroutine were programmed in assembly language on a typical machine, we could simply
436    compute |(*p)div q|, since a double-precision product can often be input to a fixed-point 
437    division instruction. But when we are restricted to integer arithmetic it is necessary either
438    to resort to multiple-precision maneuvering or to use a simple but slow iteration. The 
439    multiple-precision technique would be about three times faster than the code adopted here, but 
440    it would be comparatively long and tricky, involving about sixteen additional multiplications 
441    and divisions.
442
443    This operation is part of \MP's \quote {inner loop}; indeed, it will consume nearly 10\pct! of 
444    the running time (exclusive of input and output) if the code below is left unchanged. A 
445    machine-dependent recoding will therefore make \MP\ run faster. The present implementation is 
446    highly portable, but slow; it avoids multiplication and division except in the initial stage. 
447    System wizards should be careful to replace it with a routine that is guaranteed to produce 
448    identical results in all cases. 
449
450    As noted below, a few more routines should also be replaced by machine-dependent code, for 
451    efficiency. But when a procedure is not part of the \quote {inner loop,} such changes aren't 
452    advisable; simplicity and robustness are preferable to trickery, unless the cost is too high. 
453
454*/
455
456/* 
457    Todo: check if we can replace this (see mpmathdouble): 
458*/
459
460static int mp_scaled_aux_make_fraction(MP mp, int p, int q)
461{
462    if (q == 0) {
463        mp_confusion (mp, "division by zero");
464        return 0;
465    } else {
466        double d = TWEXP28 * (double) p / (double) q;
467        if ((p ^ q) >= 0) {
468            d += 0.5;
469            if (d >= TWEXP31) {
470                mp->arith_error = 1;
471                return EL_GORDO;
472            } else {
473                int i = (int) d;
474                if (d == (double) i && (((q > 0 ? -q : q) & 0x7FFF) * (((i & 0x3FFF) << 1) - 1) & 0x800) != 0) {
475                    --i;
476                }
477                return i;
478            }
479        } else {
480            d -= 0.5;
481            if (d <= -TWEXP31) {
482                mp->arith_error = 1;
483                return -negative_EL_GORDO;
484            } else {
485                int i = (int) d;
486                if (d == (double) i && (((q > 0 ? q : -q) & 0x7FFF) * (((i & 0x3FFF) << 1) + 1) & 0x800) != 0) {
487                    ++i;
488                }
489                return i;
490            }
491        }
492    }
493}
494
495static void mp_scaled_make_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q)
496{
497    ret->data.val = mp_scaled_aux_make_fraction (mp, p->data.val, q->data.val);
498}
499
500/*tex
501    The dual of |make_fraction| is |take_fraction|, which multiplies a given integer~|q| by a 
502    fraction~|f|. When the operands are positive, it computes $p = \lfloor qf/2^{28} + { 1\over 2}
503    \rfloor$, a symmetric function of |q| and~|f|.
504
505    This routine is even more \quote {inner loopy} than |make_fraction|; the present implementation 
506    consumes almost 20\pct! of \MP's computation time during typical jobs, so a machine-language 
507    substitute is advisable.
508*/
509
510static int mp_scaled_aux_take_fraction(MP mp, int p, int q)
511{
512    double d = (double) p *(double) q *TWEXP_28;
513    if ((p ^ q) >= 0) {
514        d += 0.5;
515        if (d >= TWEXP31) {
516            if (d != TWEXP31 || (((p & 0x7FFF) * (q & 0x7FFF)) & 040000) == 0) {
517                mp->arith_error = 1;
518            }
519            return EL_GORDO;
520        } else {
521            int i = (int) d;
522            if (d == (double) i && (((p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
523                --i;
524            }
525            return i;
526        }
527    } else {
528        d -= 0.5;
529        if (d <= -TWEXP31) {
530            if (d != -TWEXP31 || ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) == 0) {
531                mp->arith_error = 1;
532            }
533            return -negative_EL_GORDO;
534        } else {
535            int i = (int) d;
536            if (d == (double) i && ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
537                ++i;
538            }
539            return i;
540        }
541    }
542}
543
544static void mp_scaled_take_fraction(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
545{
546    ret->data.val = mp_scaled_aux_take_fraction(mp, p_orig->data.val, q_orig->data.val);
547}
548
549/*tex
550
551    When we want to multiply something by a |scaled| quantity, we use a scheme analogous to 
552    |take_fraction| but with a different scaling. Given positive operands, |take_scaled| computes 
553    the quantity $p=\lfloor qf/2^{16}+{1\over2}\rfloor$.
554
555    Once again it is a good idea to use a machine-language replacement if possible; otherwise 
556    |take_scaled| will use more than 2\pct! of the running time when the Computer Modern fonts are
557    being generated. 
558
559*/
560
561static int mp_take_scaled(MP mp, int p, int q)
562{ 
563    /* q = scaled */
564    double d = (double) p * (double) q * TWEXP_16;
565    if ((p ^ q) >= 0) {
566        d += 0.5;
567        if (d >= TWEXP31) {
568            if (d != TWEXP31 || (((p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) == 0) {
569                mp->arith_error = 1;
570            }
571            return EL_GORDO;
572        } else {
573            int i = (int) d;
574            if (d == (double) i && (((p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
575                --i;
576            }
577            return i;
578        }
579    } else {
580        d -= 0.5;
581        if (d <= -TWEXP31) {
582            if (d != -TWEXP31 || ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) == 0) {
583                mp->arith_error = 1;
584            }
585            return -negative_EL_GORDO;
586        } else {
587            int i = (int) d;
588            if (d == (double) i && ((-(p & 0x7FFF) * (q & 0x7FFF)) & 0x4000) != 0) {
589                ++i;
590            }
591            return i;
592        }
593    }
594}
595
596static void mp_scaled_take_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
597{
598    ret->data.val = mp_take_scaled(mp, p_orig->data.val, q_orig->data.val);
599}
600
601/*tex 
602
603    For completeness, there's also |make_scaled|, which computes a quotient as a |scaled| number 
604    instead of as a |fraction|. In other words, the result is $\lfloor2^{16}p/q+{1\over2}\rfloor$, 
605    if the operands are positive. \ (This procedure is not used especially often, so it is not 
606    part of \MP's inner loop.)
607
608*/
609
610static int mp_make_scaled(MP mp, int p, int q)
611{
612    if (q == 0) {
613        mp_confusion(mp, "division by zero");
614        return 0;
615    } else {
616        double d = TWEXP16 * (double) p / (double) q;
617        if ((p ^ q) >= 0) {
618            d += 0.5;
619            if (d >= TWEXP31) {
620                mp->arith_error = 1;
621                return EL_GORDO;
622            } else {
623                int i = (int) d;
624                if (d == (double) i && (((q > 0 ? -q : q) & 0x7FFF) * (((i & 0x3FFF) << 1) - 1) & 0x800) != 0) {
625                    --i;
626                }
627                return i;
628            }
629        } else {
630            d -= 0.5;
631            if (d <= -TWEXP31) {
632                mp->arith_error = 1;
633                return -negative_EL_GORDO;
634            } else {
635                int i = (int) d;
636                if (d == (double) i && (((q > 0 ? q : -q) & 0x7FFF) * (((i & 0x3FFF) << 1) + 1) & 0x800) != 0) {
637                    ++i;
638                }
639                return i;
640            }
641        }
642    }
643}
644
645static void mp_scaled_make_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
646{
647    ret->data.val = mp_make_scaled(mp, p_orig->data.val, q_orig->data.val);
648}
649
650/*tex
651    The following function is used to create a scaled integer from a given decimal fraction 
652    $(.d_0d_1\ldots d_{k-1})$, where |0<=k<=17|. This converts a decimal fraction.
653*/
654
655static int mp_round_decimals(MP mp, unsigned char *b, int k)
656{
657    unsigned a = 0;
658    (void) mp; 
659    for (int l = k-1; l >= 0; l-- ) {
660        if (l < 16) {
661            /*tex digits for |k>=17| cannot affect the result */
662            a = (a + (unsigned) (*(b+l) - '0') * two) / 10;
663        }
664    }
665    return (int) (a + 1)/2;
666}
667
668/*tex 
669
670    We no longer have these character mapping so we can just as well do the next without class 
671    checking. Also because signs are hard checked.
672*/
673
674static void mp_scaled_wrapup_numeric_token(MP mp, int n, int f)
675{
676    if (n < 32768) {
677        int mod = (n * unity + f); /* scaled */
678        set_cur_mod(mod);
679        if (mod >= fraction_one) {
680            if (internal_value(mp_warning_check_internal).data.val > 0 && (mp->scanner_status != mp_tex_flushing_state)) {
681                char msg[256];
682                snprintf(msg, 256, "Number is too large (%s)", mp_string_scaled(mp,mod));
683                mp_error(
684                    mp,
685                    msg,
686                    "It is at least 4096. Continue and I'll try to cope with that big value;\n"
687                    "but it might be dangerous. (Set warningcheck:=0 to suppress this message.)"
688                );
689            }
690        }
691    } else if (mp->scanner_status != mp_tex_flushing_state) {
692        mp_error(
693            mp,
694            "Enormous number has been reduced",
695            "I can\'t handle numbers bigger than 32767.99998; so I've changed your constant\n"
696            "to that maximum amount."
697        );
698        set_cur_mod(EL_GORDO);
699    }
700    set_cur_cmd(mp_numeric_command);
701}
702
703static void mp_scaled_scan_fractional_token(MP mp, int n)
704{ 
705    /* n: scaled */
706    int f; /* scaled */
707    int k = 0;
708    do {
709        k++;
710        mp->cur_input.loc_field++;
711    } while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class);
712    f = mp_round_decimals(mp, (unsigned char *)(mp->buffer+mp->cur_input.loc_field-k), (int) k);
713    if (f == unity) {
714        n++;
715        f = 0;
716    }
717    mp_scaled_wrapup_numeric_token(mp, n, f);
718}
719
720static void mp_scaled_scan_numeric_token(MP mp, int n)
721{
722    while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
723        if (n < 32768) {
724            n = 10 * n + mp->buffer[mp->cur_input.loc_field] - '0';
725        }
726        mp->cur_input.loc_field++;
727    }
728    if (! (mp->buffer[mp->cur_input.loc_field] == '.' && mp->char_class[mp->buffer[mp->cur_input.loc_field + 1]] == mp_digit_class)) {
729        mp_scaled_wrapup_numeric_token(mp, n, 0);
730    } else {
731        mp->cur_input.loc_field++;
732        mp_scaled_scan_fractional_token(mp, n);
733    }
734}
735
736/*tex
737
738    Here is a typical example of how the routines above can be used. It computes the function 
739
740    $$
741        {1 \over 3 \tau}f(\theta,\phi)= {\tau^{-1} \bigl( 2 + \sqrt2\, (\sin \theta -{1 \over 16}
742        \sin \phi) (\sin \phi - {1 \over 16} \sin \theta) (\cos \theta - \cos \phi) \bigr) \over
743        3\, \bigl( 1 + {1 \over 2} (\sqrt5 - 1) \cos \theta + {1 \over 2} (3 - \sqrt 5 \,) \cos 
744        \phi \bigr)},
745    $$
746
747    where $\tau$ is a |scaled| \quote {tension} parameter. This is \MP's magic fudge factor for 
748    placing the first control point of a curve that starts at an angle $\theta$ and ends at an 
749    angle $\phi$ from the straight path. (Actually, if the stated quantity exceeds 4, \MP\ reduces 
750    it to~4.)
751
752    The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$. (It's a sum of 
753    eight terms whose absolute values can be bounded using relations such as $\sin \theta \cos 
754    \theta |1 \over 2|$.) Thus the numerator is positive; and since the tension $\tau$ is 
755    constrained to be at least $3\over4$, the numerator is less than $16 \over 3$. The denominator 
756    is nonnegative and at most~6. Hence the fixed-point calculations below are guaranteed to stay 
757    within the bounds of a 32-bit computer word.
758
759    The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction| arguments |st|, 
760    |ct|, |sf|, and |cf|, representing $\sin \theta$, $\cos \theta$, $\sin \phi$, and $\cos \phi$, 
761    respectively.
762
763*/
764
765static void mp_scaled_velocity(MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t)
766{
767    int acc, num, denom; 
768    acc = mp_scaled_aux_take_fraction(mp, st->data.val - (sf->data.val / 16), sf->data.val - (st->data.val / 16));
769    acc = mp_scaled_aux_take_fraction(mp, acc, ct->data.val - cf->data.val);
770    num = fraction_two + mp_scaled_aux_take_fraction(mp, acc, 379625062);
771    /*tex
772        $2^{28}\sqrt2\approx379625062.497$
773    */
774    denom = fraction_three + mp_scaled_aux_take_fraction(mp, ct->data.val, 497706707) + mp_scaled_aux_take_fraction(mp, cf->data.val, 307599661);
775    /*tex
776        $3\cdot2^{27}\cdot(\sqrt5-1)\approx497706706.78$ and
777        $3\cdot2^{27}\cdot(3-\sqrt5\,)\approx307599661.22$
778    */
779    if (t->data.val != unity) {
780        /*tex 
781            |make_scaled(fraction,scaled) = fraction| 
782        */
783        num = mp_make_scaled(mp, num, t->data.val);
784    }
785    if (num / 4 >= denom) {
786        ret->data.val = fraction_four;
787    } else {
788        ret->data.val = mp_scaled_aux_make_fraction(mp, num, denom);
789    }
790}
791
792/*tex
793    The following somewhat different subroutine tests rigorously if $ab$ is greater than, equal to,
794    or less than~$cd$, given integers $(a,b,c,d)$. In most cases a quick decision is reached. The 
795    result is $+1$, 0, or~$-1$ in the three respective cases.
796*/
797
798static int mp_scaled_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig)
799{
800    int a = a_orig->data.val;
801    int b = b_orig->data.val;
802    int c = c_orig->data.val;
803    int d = d_orig->data.val;
804    if (a < 0) {
805        a = -a;
806        b = -b;
807    }
808    if (c < 0) {
809        c = -c;
810        d = -d;
811    }
812    if (d <= 0) {
813        if (b >= 0) {
814            if ((a == 0 || b == 0) && (c == 0 || d == 0)) {
815                return 0;
816            } else {
817                return 1;
818            }
819        } else if (d == 0) {
820            return a == 0 ? 0 : -1;
821        } else {
822            int q = a;
823            a = c;
824            c = q;
825            q = -b;
826            b = -d;
827            d = q;
828        }
829    } else if (b <= 0) {
830        if (b < 0 && a > 0) {
831            return -1;
832        } else {
833            return c == 0 ? 0 : -1;
834        }
835    }
836    while (1) {
837        int q = a / d;
838        int r = c / b;
839        if (q != r) {
840            return q > r ? 1 : -1;
841        } else {
842            q = a % d;
843            r = c % b;
844            if (r == 0) {
845                return q ? 1 : 0;
846            } else if (q == 0) {
847                return -1;
848            } else {
849                a = b;
850                b = q;
851                c = d;
852                d = r;
853            }
854        }
855    }
856    /*tex Now |a > d > 0| and |c > b > 0|. */
857}
858
859/* 
860    Now here's a subroutine that's handy for all sorts of path computations: Given a quadratic
861    polynomial $B(a,b,c;t)$, the |crossing_point| function returns the unique |fraction| value |t|
862    between 0 and~1 at which $B(a,b,c;t)$ changes from positive to negative, or returns |t = 
863    fraction_one + 1| if no such value exists. If |a < 0| (so that $B(a,b,c;t)$ is already negative
864    at |t = 0|), |crossing_point| returns the value zero.
865
866    The general bisection method is quite simple when $n=2$, hence |crossing_point| does not take 
867    much time. At each stage in the recursion we have a subinterval defined by |l| and~|j| such 
868    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 
869    subinterval where $x_0\G0$ and $\min(x_1,x_2)<0$.
870
871    It is convenient for purposes of calculation to combine the values of |l| and~|j| in a single 
872    variable $d=2^l+j$, because the operation of bisection then corresponds simply to doubling 
873    $d$ and possibly adding~1. Furthermore it proves to be convenient to modify our previous 
874    conventions for bisection slightly, maintaining the variables $X_0 = 2^lx_0$, $X_1 = 2^l(x_0 -
875    x_1)$, and $X_2 = 2^l(x_1 - x_2)$. With these variables the conditions $x_0 \ge 0$ and
876    $\min(x_1,x_2) < 0$ are equivalent to $\max(X_1,X_1+X_2) > X_0 \ge 0$.
877
878    The following code maintains the invariant relations $0\L|x0| < \max(|x1|,|x1| + |x2|)$, 
879    $\vert|x1|\vert < 2^{30}$, $\vert|x2|\vert < 2^{30}$; it has been constructed in such a way 
880    that no arithmetic overflow will occur if the inputs satisfy $a < 2^{30}$, $\vert a - b\vert <
881    2^{30}$, and $\vert b - c\vert < 2^{30}$.
882*/
883
884static void mp_scaled_crossing_point(MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc)
885{
886    int x, xx, x0, x1, x2; /* temporary registers for bisection */
887    int a = aa->data.val;
888    int b = bb->data.val;
889    int c = cc->data.val;
890    int d;                 /* recursive counter */
891    (void) mp;
892    if (a < 0) {
893        ret->data.val = zero_crossing;
894        return;
895    } else if (c >= 0) {
896        if (b >= 0) {
897            if (c > 0) {
898                ret->data.val = no_crossing;
899            } else if ((a == 0) && (b == 0)) {
900                ret->data.val = no_crossing;
901            } else {
902                ret->data.val = one_crossing;
903            }
904            return;
905        } else if (a == 0) {
906            ret->data.val = zero_crossing;
907            return;
908        }
909    } else if (a == 0) {
910        if (b <= 0) {
911            ret->data.val = zero_crossing;
912            return;
913        }
914    }
915    /* Use bisection to find the crossing point... */
916    d = 1;
917    x0 = a;
918    x1 = a - b;
919    x2 = b - c;
920    do {
921        x = (x1 + x2) / 2;
922        if (x1 - x0 > x0) {
923            x2 = x;
924            x0 += x0;
925            d += d;
926        } else {
927            xx = x1 + x - x0;
928            if (xx > x0) {
929                x2 = x;
930                x0 += x0;
931                d += d;
932            } else {
933                x0 = x0 - xx;
934                if ((x <= x0) && (x + x2 <= x0)) {
935                    ret->data.val = no_crossing;
936                    return;
937                } else {
938                    x1 = x;
939                    d = d + d + 1;
940                }
941            }
942        }
943    } while (d < fraction_one);
944    ret->data.val = d - fraction_one;
945}
946
947/*tex
948    We conclude this set of elementary routines with some simple rounding and truncation 
949    operations; |round_unscaled| rounds a |scaled| and converts it to |int|
950*/
951
952static int mp_scaled_round_unscaled(mp_number *x_orig)
953{
954    int x = x_orig->data.val;
955    if (x >= 32768) {
956        return 1 + ((x-32768) / 65536);
957    } else if (x >= -32768) {
958        return 0;
959    } else {
960        return -(1+((-(x+1)-32768) / 65536));
961    }
962}
963static void mp_scaled_floor(mp_number *i)
964{
965    i->data.val = i->data.val & -65536;
966}
967
968static void mp_scaled_fraction_to_round_scaled(mp_number *x_orig)
969{
970    int x = x_orig->data.val;
971    x_orig->type = mp_scaled_type;
972    x_orig->data.val = (x>=2048 ? 1+((x-2048) / 4096)  : ( x>=-2048 ? 0 : -(1+((-(x+1)-2048) / 4096))));
973}
974
975/*tex
976    \MP\ computes all of the necessary special functions from scratch, without relying on |real| 
977    arithmetic or system subroutines for sines, cosines, etc.
978    
979    To get the square root of a |scaled| number |x|, we want to calculate $s = \lfloor 2^8\! \sqrt 
980    x + {1\over2} \rfloor$. If $x > 0$, this is the unique integer such that $2^{16}x - s\L s^2 < 
981    2^{16}x + s$. The following subroutine determines $s$ by an iterative method that maintains the
982    invariant relations $x=2^{46 - 2k}x_0 \bmod 2^{30}$, $0 < y = \lfloor 2^{16 - 2k}x_0\rfloor -
983    s^2 + s\L q = 2s$, where $x_0$ is the initial value of $x$. The value of~$y$ might, however, be 
984    zero at the start of the first iteration.
985*/
986
987static void mp_scaled_sqrt(MP mp, mp_number *ret, mp_number *x_orig)
988{
989    int x = x_orig->data.val;
990    if (x <= 0) {
991        if (x < 0) {
992            char msg[256];
993            snprintf(msg, 256, "Square root of %s has been replaced by 0", mp_string_scaled(mp, x));
994            mp_error(
995                mp,
996                msg,
997                "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
998                "Proceed, with fingers crossed."
999            );
1000        }
1001        ret->data.val = 0;
1002    } else {
1003        int k = 23; /* iteration control counter */
1004        int y;
1005        int q = 2;
1006        while (x < fraction_two) {  /* i.e., |while x<|\unskip */
1007            k--;
1008            x = x + x + x + x;
1009        }
1010        if (x < fraction_four)
1011            y = 0;
1012        else {
1013            x = x - fraction_four;
1014            y = 1;
1015        }
1016        do {
1017            /*tex Decrease |k| by 1, maintaining the invariant relations between |x|, |y|, and~|q|: */
1018            x += x;
1019            y += y;
1020            if (x >= fraction_four) {
1021                /* note that |fraction_four=| */
1022                x = x - fraction_four;
1023                y++;
1024            };
1025            x += x;
1026            y = y + y - q;
1027            q += q;
1028            if (x >= fraction_four) {
1029                x = x - fraction_four;
1030                y++;
1031            };
1032            if (y > (int) q) {
1033                y -= q;
1034                q += 2;
1035            } else if (y <= 0) {
1036                q -= 2;
1037                y += q;
1038            };
1039            k--;
1040        } while (k != 0);
1041        ret->data.val = (int) (q/2);
1042    }
1043}
1044
1045/*tex
1046    Pythagorean addition $\psqrt{a^2 + b^2}$ is implemented by an elegant iterative scheme due to 
1047    Cleve Moler and Donald Morrison [{\sl IBM Journal   of Research and Development \bf27} (1983),
1048    577--581]. It modifies |a| and~|b| in such a way that their Pythagorean sum remains invariant, 
1049    while the smaller argument decreases.
1050*/
1051
1052static void mp_scaled_pyth_add(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
1053{
1054    int a = abs(a_orig->data.val);
1055    int b = abs(b_orig->data.val);
1056    if (a < b) {
1057        int r = b;
1058        b = a;
1059        a = r;
1060    }
1061    /* now |0<=b<=a| */
1062    if (b > 0) {
1063        int big;  
1064        /*tex Is the result dangerously near $2^{31}$? */
1065        if (a < fraction_two) {
1066            big = 0;
1067        } else {
1068            a = a / 4;
1069            b = b / 4;
1070            big = 1;
1071        }
1072        /*tex 
1073            We reduced the precision to avoid arithmetic overflow. Replace |a| by an approximation 
1074            to $\psqrt{a^2+b^2}:$ 
1075        */
1076        while (1) {
1077            int r = mp_scaled_aux_make_fraction(mp, b, a);
1078            r = mp_scaled_aux_take_fraction(mp, r, r);
1079            /*tex Now $r\approx b^2/a^2$: */
1080            if (r == 0) {
1081                break;
1082            } else {
1083                r = mp_scaled_aux_make_fraction(mp, r, fraction_four + r);
1084                a = a + mp_scaled_aux_take_fraction(mp, a + a, r);
1085                b = mp_scaled_aux_take_fraction(mp, b, r);
1086            }
1087        }
1088        if (big) {
1089            if (a < fraction_two) {
1090                a = a + a + a + a;
1091            } else {
1092                mp->arith_error = 1;
1093                a = EL_GORDO;
1094            }
1095        }
1096    }
1097    ret->data.val = a;
1098}
1099
1100/*tex
1101    The key idea here is to reflect the vector $(a,b)$ about the line through $(a,b/2)$. Here is a
1102    similar algorithm for $\psqrt{a^2-b^2}$. It converges slowly when $b$ is near $a$, but 
1103    otherwise it works fine.
1104*/
1105
1106static void mp_scaled_pyth_sub(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
1107{
1108    int a = abs(a_orig->data.val);
1109    int b = abs(b_orig->data.val);
1110    if (a <= b) {
1111        /*tex Handle erroneous |pyth_sub| and set |a:=0|: */
1112        if (a < b) {
1113            char msg[256];
1114            char *astr = mp_strdup(mp_string_scaled(mp, a));
1115            snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, mp_string_scaled(mp, b));
1116            mp_memory_free(astr);
1117            mp_error(
1118                mp,
1119                msg,
1120                "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
1121                "Proceed, with fingers crossed."
1122            );
1123        }
1124        a = 0;
1125    } else {
1126        int big;  
1127        /*tex Is the result dangerously near $2^{31}$? */
1128        if (a < fraction_four) {
1129            big = 0;
1130        } else {
1131            a = (int) a/2;
1132            b = (int) b/2;
1133            big = 1;
1134        }
1135        /*tex Replace |a| by an approximation to $\psqrt{a^2-b^2}$ */
1136        while (1) {
1137            int r = mp_scaled_aux_make_fraction(mp, b, a);
1138            r = mp_scaled_aux_take_fraction(mp, r, r);
1139            /*tex Now $r\approx b^2/a^2$. */
1140            if (r == 0) {
1141                break;
1142            } else {
1143                r = mp_scaled_aux_make_fraction(mp, r, fraction_four - r);
1144                a = a - mp_scaled_aux_take_fraction(mp, a + a, r);
1145                b = mp_scaled_aux_take_fraction(mp, b, r);
1146            }
1147        }
1148        if (big) {
1149            a *= 2;
1150        }
1151    }
1152    ret->data.val = a;
1153}
1154
1155/*tex 
1156    We just abuse doubles here.
1157*/
1158
1159static void mp_scaled_power_of(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
1160{
1161    double p = pow(mp_scaled_to_double(a_orig), mp_scaled_to_double(b_orig));
1162    long r = lround(p * 65536.0);
1163    if (r > 0) {
1164        if (r >= EL_GORDO) {
1165            mp->arith_error = 1;
1166            r = EL_GORDO;
1167        }
1168    } else if (r < 0) {
1169        if (r <= - EL_GORDO) {
1170            mp->arith_error = 1;
1171            r = - EL_GORDO;
1172        }
1173    }
1174    ret->data.val = r;
1175}
1176
1177/*tex
1178
1179    The subroutines for logarithm and exponential involve two tables. The first is simple: 
1180    |two_to_the[k]| equals $2^k$. The second involves a bit more calculation, which the author 
1181    claims to have done correctly: |mp_m_spec_log[k]| is $2^{27}$ times $\ln\bigl(1/(1-2^{-k})
1182    \bigr)= 2^{-k}+{1\over2}2^{-2k}+{1\over3}2^{-3k}+\cdots\,$, rounded to the nearest integer.
1183    Here is the routine that calculates $2^8$ times the natural logarithm of a |scaled| quantity; 
1184    it is an integer approximation to $2^{24}\ln(x/2^{16})$, when |x| is a given positive integer.
1185
1186    The method is based on exercise 1.2.2--25 in {\sl The Art of Computer Programming}: During the
1187    main iteration we have $1\L 2^{-30}x<1/(1-2^{1-k})$, and the logarithm of $2^{30}x$ remains to 
1188    be added to an accumulator register called~$y$. Three auxiliary bits of accuracy are retained 
1189    in~$y$ during the calculation, and sixteen auxiliary bits to extend |y| are kept in~|z| during 
1190    the initial argument reduction. (We add $100\cdot2^{16}=6553600$ to~|z| and subtract 100 
1191    from~|y| so that |z| will not become negative; also, the actual amount subtracted from~|y| 
1192    is~96, not~100, because we want to add~4 for rounding before the final division by~8.)
1193
1194*/
1195
1196static void mp_scaled_m_log(MP mp, mp_number *ret, mp_number *x_orig)
1197{
1198    int x = x_orig->data.val;
1199    if (x <= 0) {
1200        /* Handle non-positive logarithm: */
1201        char msg[256];
1202        snprintf(msg, 256, "Logarithm of %s has been replaced by 0", mp_string_scaled(mp, x));
1203        mp_error(
1204            mp,
1205            msg,
1206            "Since I don't take logs of non-positive numbers, I'm zeroing this one.\n"
1207            "Proceed, with fingers crossed."
1208        );
1209        ret->data.val = 0;
1210    } else {
1211        int k = 2;                    /* iteration counter, starts at 2 */
1212        int y = 1302456956 + 4 - 100; /* $14\times2^{27}\ln2\approx1302456956.421063$ */
1213        int z = 27595 + 6553600;      /* and $2^{16}\times .421063\approx 27595$ */
1214        /* $2^{27}\ln2\approx 93032639.74436163$ and $2^{16}\times.74436163\approx 48782$ */
1215        while (x < fraction_four) {
1216            x = 2*x;
1217            y -= 93032639;
1218            z -= 48782;
1219        }
1220        y = y + (z / unity);
1221        while (x > fraction_four + 4) {
1222            /*tex
1223                Increase |k| until |x| can be multiplied by a factor of $2^{-k}$, and adjust $y$ 
1224                accordingly.
1225            */
1226            z = ((x - 1) / two_to_the (k)) + 1; /* $z=\lceil x/2^k\rceil$ */
1227            while (x < fraction_four + z) {
1228                z = (z + 1)/2;
1229                k++;
1230            };
1231            y += mp_m_spec_log[k];
1232            x -= z;
1233        }
1234        ret->data.val = (y / 8);
1235    }
1236}
1237
1238/*tex
1239
1240    Conversely, the exponential routine calculates $\exp(x/2^8)$, when |x| is |scaled|. The result 
1241    is an integer approximation to $2^{16}\exp(x/2^{24})$, when |x| is regarded as an integer.
1242
1243*/
1244
1245static void mp_scaled_m_exp(MP mp, mp_number *ret, mp_number *x_orig)
1246{
1247    int y, z;  /* auxiliary registers */
1248    int x = x_orig->data.val;
1249    if (x > 174436200) {
1250        /* $2^{24}\ln((2^{31}-1)/2^{16})\approx 174436199.51$ */
1251        mp->arith_error = 1;
1252        ret->data.val = EL_GORDO;
1253    } else if (x < -197694359) {
1254        /* $2^{24}\ln(2^{-1}/2^{16})\approx-197694359.45$ */
1255        ret->data.val = 0;
1256    } else {
1257        if (x <= 0) {
1258            z = -8 * x;
1259            y = 0x100000; /* $y=2^{20}$ */
1260        } else {
1261            if (x <= 127919879) {
1262                z = 1023359037 - 8 * x;
1263                /* $2^{27}\ln((2^{31}-1)/2^{20})\approx 1023359037.125$ */
1264            } else {
1265                /* |z| is always nonnegative */
1266                z = 8 * (174436200 - x);
1267            }
1268            y = EL_GORDO;
1269        }
1270        /* Multiply |y| by $\exp(-z/2^{27})$: */
1271        {
1272            int k = 1;
1273            while (z > 0) {
1274                while (z >= mp_m_spec_log[k]) {
1275                    z -= mp_m_spec_log[k];
1276                    y = y - 1 - ((y - two_to_the(k - 1)) / two_to_the(k));
1277                }
1278                k++;
1279            }
1280        }
1281        if (x <= 127919879) {
1282            ret->data.val = ((y + 8) / 16);
1283        } else {
1284            ret->data.val = y;
1285        }
1286    }
1287}
1288
1289/*tex
1290
1291    The idea here is that subtracting |mp_m_spec_log[k]| from |z| corresponds to multiplying |y| 
1292    by $1-2^{-k}$.
1293
1294    A subtle point (which had to be checked) was that if $x=127919879$, the value of~|y| will 
1295    decrease so that |y+8| doesn't overflow. In fact, $z$ will be 5 in this case, and |y| will 
1296    decrease by~64 when |k=25| and by~16 when |k=27|.The trigonometric subroutines use an auxiliary
1297    table such that |spec_atan[k]| contains an approximation to the |angle| whose tangent is~$1 /
1298    2^k$. $\arctan2^{-k}$ times $2^{20}\cdot180/\pi$Given integers |x| and |y|, not both zero, the 
1299    |n_arg| function returns the |angle| whose tangent points in the direction $(x,y)$. This 
1300    subroutine first determines the correct octant, then solves the problem for |0<=y<=x|, then
1301    converts the result appropriately to return an answer in the range |-one_eighty_deg <= 
1302    one_eighty_deg|. (The answer is |+one_eighty_deg| if |y=0| and |x<0|, but an answer of 
1303    |-one_eighty_deg| is possible if, for example, |y = -1| and $x = -2^{30}$.)
1304
1305*/
1306
1307static void mp_scaled_n_arg(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
1308{
1309    int octant; /* octant code */
1310    int x = x_orig->data.val;
1311    int y = y_orig->data.val;
1312    if (x >= 0) {
1313        octant = first_octant;
1314    } else {
1315        x = -x;
1316        octant = first_octant + negate_x;
1317    }
1318    if (y < 0) {
1319        y = -y;
1320        octant = octant + negate_y;
1321    }
1322    if (x < y) {
1323        int t = y;
1324        y = x;
1325        x = t;
1326        octant = octant + switch_x_and_y;
1327    }
1328    if (x == 0) {
1329        if (internal_value(mp_default_zero_angle_internal).data.val < 0) {
1330            mp_error(
1331                mp,
1332                "angle(0,0) is taken as zero",
1333                "The 'angle' between two identical points is undefined. I'm zeroing this one.\n"
1334                "Proceed, with fingers crossed."
1335            );
1336            ret->data.val = 0;
1337        } else { 
1338            ret->data.val = internal_value(mp_default_zero_angle_internal).data.val;
1339        }
1340    } else {
1341        int z = 0; /* auxiliary register */
1342        ret->type = mp_angle_type;
1343        /* Set variable |z| to the arg of $(x,y)$: */
1344        while (x >= fraction_two) {
1345            x = x/2;
1346            y = y/2;
1347        }
1348        if (y > 0) {
1349            int k= 0; /* loop counter */
1350            while (x < fraction_one) {
1351                x += x;
1352                y += y;
1353            };
1354            /* Increase |z| to the arg of $(x,y)$: */
1355            do {
1356                y += y;
1357                k++;
1358                if (y > x) {
1359                    int t = x; 
1360                    z = z + mp_m_spec_atan[k];
1361                    x = x + (y / two_to_the(k + k));
1362                    y = y - t;
1363                };
1364            } while (k != 15);
1365            do {
1366                y += y;
1367                k++;
1368                if (y > x) {
1369                    z = z + mp_m_spec_atan[k];
1370                    y = y - x;
1371                };
1372            } while (k != 26);
1373        }
1374        /* Return an appropriate answer based on |z| and |octant|: */
1375        switch (octant) {
1376            case first_octant:   ret->data.val =  z;                  break;
1377            case second_octant:  ret->data.val = -z + ninety_deg;     break;
1378            case third_octant:   ret->data.val =  z + ninety_deg;     break;
1379            case fourth_octant:  ret->data.val = -z + one_eighty_deg; break;
1380            case fifth_octant:   ret->data.val =  z - one_eighty_deg; break;
1381            case sixth_octant:   ret->data.val = -z - ninety_deg;     break;
1382            case seventh_octant: ret->data.val =  z - ninety_deg;     break;
1383            case eighth_octant:  ret->data.val = -z;                  break;
1384        }
1385// 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);
1386    }
1387}
1388
1389/*tex
1390    At this point we have |x>=y>=0|, and |x>0|. The numbers are scaled up or down until $2^{28}\L 
1391    x<2^{29}$, so that accurate fixed-point calculations will be made.
1392*/
1393
1394/*tex
1395
1396    During the calculations of this section, variables |x| and~|y| represent actual coordinates 
1397    $(x,2^{-k}y)$. We will maintain the condition |x>=y|, so that the tangent will be at most 
1398    $2^{-k}$. If $x<2y$, the tangent is greater than $2^{-k-1}$. The transformation $(a,b) \mapsto
1399    (a+b\tan\phi,b-a\tan\phi)$ replaces $(a,b)$ by coordinates whose angle has decreased by~$\phi$; 
1400    in the special case $a=x$, $b=2^{-k}y$, and $\tan\phi=2^{-k-1}$, this operation reduces to the 
1401    particularly simple iteration shown here. [Cf.~John E. Meggitt,  {\sl IBM Journal of Research 
1402    and Development \bf6} (1962), 210--226.]
1403
1404    The initial value of |x| will be multiplied by at most $(1 + {1\over2})(1 + {1\over8}) (1 + 
1405    {1\over32})\cdots\approx 1.7584$; hence there is no chance of integer overflow. Conversely, the 
1406    |n_sin_cos| routine takes an |angle| and produces the sine and cosine of that angle. The 
1407    results of this routine are stored in global integer variables |n_sin| and |n_cos|.
1408
1409    Given an integer |z| that is $2^{20}$ times an angle $\theta$ in degrees, the purpose of 
1410    |n_sin_cos(z)| is to set |x=| and |y=| (approximately), for some rather large number~|r|. The
1411    maximum of |x| and |y| will be between $2^{28}$ and $2^{30}$, so that there will be hardly any 
1412    loss of accuracy. Then |x| and~|y| are divided by~|r|.
1413
1414*/
1415
1416static void mp_scaled_n_sin_cos(MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin)
1417{
1418    int k;                    /* loop control variable */
1419    int q;                    /* specifies the quadrant */
1420    int x, y, t;              /* temporary registers */
1421    int z = z_orig->data.val; /* scaled */
1422    mp_number x_n, y_n, ret;
1423    mp_scaled_allocate_number(mp, &ret, mp_scaled_type);
1424    mp_scaled_allocate_number(mp, &x_n, mp_scaled_type);
1425    mp_scaled_allocate_number(mp, &y_n, mp_scaled_type);
1426    while (z < 0) {
1427        z = z + three_sixty_deg;
1428    }
1429    z = z % three_sixty_deg;
1430    /*tex Now |0 <= z < three_sixty_deg|. */
1431    q = z / forty_five_deg;
1432    z = z % forty_five_deg;
1433    x = fraction_one;
1434    y = x;
1435    if (! odd(q)) {
1436        z = forty_five_deg - z;
1437    }
1438    /*tex Subtract angle |z| from |(x,y)|: */
1439    k = 1;
1440    while (z > 0) {
1441        if (z >= mp_m_spec_atan[k]) {
1442            z = z - mp_m_spec_atan[k];
1443            t = x;
1444            x = t + y / two_to_the(k);
1445            y = y - t / two_to_the(k);
1446        }
1447        k++;
1448    }
1449    if (y < 0) {
1450        /*tex This precaution may never be needed */
1451        y = 0;
1452    }
1453    /*tex Convert |(x,y)| to the octant determined by~|q|: */
1454    switch (q) {
1455        case 0:                        break;
1456        case 1: t = x; x =  y; y =  t; break;
1457        case 2: t = x; x = -y; y =  t; break;
1458        case 3:        x = -x;         break;
1459        case 4:        x = -x; y = -y; break;
1460        case 5: t = x; x = -y; y = -t; break;
1461        case 6: t = x; x =  y; y = -t; break;
1462        case 7:                y = -y; break;
1463    }
1464    x_n.data.val = x;
1465    y_n.data.val = y;
1466    mp_scaled_pyth_add(mp, &ret, &x_n, &y_n);
1467    n_cos->data.val = mp_scaled_aux_make_fraction(mp, x, ret.data.val);
1468    n_sin->data.val = mp_scaled_aux_make_fraction(mp, y, ret.data.val);
1469    mp_scaled_free_number(mp, &ret);
1470    mp_scaled_free_number(mp, &x_n);
1471    mp_scaled_free_number(mp, &y_n);
1472}
1473
1474/*tex
1475
1476    In this case the octants are numbered sequentially.The main iteration of |n_sin_cos| is similar 
1477    to that of |n_arg| but applied in reverse. The values of |mp_m_spec_atan[k]| decrease slowly 
1478    enough that this loop is guaranteed to terminate before the (nonexistent) value 
1479    |mp_m_spec_atan[27]| would be required.To initialize the |randoms| table, we call the following 
1480    routine.
1481
1482*/
1483
1484static void mp_scaled_init_randoms(MP mp, int seed)
1485{
1486    int k = 1; /* more or less random integers */
1487    int j = abs(seed);
1488    while (j >= fraction_one) {
1489        j = j/2;
1490    }
1491    for (int i = 0; i <= 54; i++) {
1492        int jj = k;
1493        k = j - k;
1494        j = jj;
1495        if (k < 0) {
1496            k += fraction_one;
1497        }
1498        mp->randoms[(i * 21) % 55].data.val = j;
1499    }
1500    /*tex We \quote {warm up} the array: */
1501    mp_new_randoms(mp);
1502    mp_new_randoms(mp);
1503    mp_new_randoms(mp);
1504}
1505
1506static void mp_scaled_print(MP mp, mp_number *n)
1507{
1508    mp_print_e_str(mp, mp_string_scaled(mp, n->data.val));
1509}
1510
1511static char *mp_scaled_tostring(MP mp, mp_number *n)
1512{
1513    return mp_string_scaled(mp, n->data.val);
1514}
1515
1516static void mp_scaled_modulo(mp_number *a, mp_number *b)
1517{
1518    a->data.val = a->data.val % b->data.val;
1519}
1520
1521static void mp_next_random(MP mp, mp_number *ret)
1522{
1523    if ( mp->j_random == 0) {
1524        mp_new_randoms(mp);
1525    } else {
1526        mp->j_random = mp->j_random-1;
1527    }
1528    mp_scaled_clone(ret, &(mp->randoms[mp->j_random]));
1529}
1530
1531/*tex
1532
1533    To produce a uniform random number in the range |0 <= u < x| or |0 >= u > x| or |0 = u = x|, 
1534    given a |scaled| value~|x|, we proceed as shown here. Note that the call of |take_fraction| 
1535    will produce the values 0 and~|x| with about half the probability that it will produce any 
1536    other particular values between 0 and~|x|, because it rounds its answers.
1537
1538*/
1539
1540static void mp_scaled_m_unif_rand(MP mp, mp_number *ret, mp_number *x_orig)
1541{
1542    mp_number x, abs_x, u, y; /* |y| is trial value */
1543    mp_scaled_allocate_number(mp, &y, mp_fraction_type);
1544    mp_scaled_allocate_clone(mp, &x, mp_scaled_type, x_orig);
1545    mp_scaled_allocate_abs(mp, &abs_x, mp_scaled_type, &x);
1546    mp_scaled_allocate_number(mp, &u, mp_scaled_type);
1547    mp_next_random(mp, &u);
1548    /*|take_fraction (y, abs_x, u);|*/
1549    mp_scaled_take_fraction(mp, &y, &abs_x, &u);
1550    if (mp_scaled_equal(&y, &abs_x)) {
1551        /*|set_number_to_zero(*ret);|*/
1552        mp_scaled_clone(ret, &((math_data *)mp->math)->md_zero_t);
1553    } else if (mp_scaled_greater(&x, &((math_data *)mp->math)->md_zero_t)) {
1554        mp_scaled_clone(ret, &y);
1555    } else {
1556        mp_scaled_clone(ret, &y);
1557        mp_scaled_negate(ret);
1558    }
1559    mp_scaled_free_number(mp, &y);
1560    mp_scaled_free_number(mp, &abs_x);
1561    mp_scaled_free_number(mp, &x);
1562    mp_scaled_free_number(mp, &u);
1563}
1564
1565/*tex
1566
1567    Finally, a normal deviate with mean zero and unit standard deviation can readily be obtained 
1568    with the ratio method (Algorithm 3.4.1R in {\sl The Art of Computer Programming}).
1569
1570*/
1571
1572static void mp_scaled_m_norm_rand(MP mp, mp_number *ret)
1573{
1574    mp_number abs_x, u, r, la, xa;
1575    mp_scaled_allocate_number(mp, &la, mp_scaled_type);
1576    mp_scaled_allocate_number(mp, &xa, mp_scaled_type);
1577    mp_scaled_allocate_number(mp, &abs_x, mp_scaled_type);
1578    mp_scaled_allocate_number(mp, &u, mp_scaled_type);
1579    mp_scaled_allocate_number(mp, &r, mp_scaled_type);
1580    do {
1581        do {
1582            mp_number v;
1583            mp_scaled_allocate_number(mp, &v, mp_scaled_type);
1584            mp_next_random(mp, &v);
1585            mp_scaled_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t);
1586            mp_scaled_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v);
1587            mp_scaled_free_number(mp, &v);
1588            mp_next_random(mp, &u);
1589            mp_scaled_clone(&abs_x, &xa);
1590            mp_scaled_abs(&abs_x);
1591        } while (! mp_scaled_less(&abs_x, &u));
1592        mp_scaled_make_fraction(mp, &r, &xa, &u);
1593        mp_scaled_clone(&xa, &r);
1594        mp_scaled_m_log(mp, &la, &u);
1595        mp_scaled_set_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la);
1596    } while (mp_scaled_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0);
1597    mp_scaled_clone(ret, &xa);
1598    mp_scaled_free_number(mp, &r);
1599    mp_scaled_free_number(mp, &abs_x);
1600    mp_scaled_free_number(mp, &la);
1601    mp_scaled_free_number(mp, &xa);
1602    mp_scaled_free_number(mp, &u);
1603}
1604
1605static void mp_scaled_set_precision(MP mp)
1606{
1607    (void) mp;
1608}
1609
1610static void mp_scaled_free_math(MP mp)
1611{
1612    mp_scaled_free_number(mp, &(mp->math->md_epsilon_t));
1613    mp_scaled_free_number(mp, &(mp->math->md_inf_t));
1614    mp_scaled_free_number(mp, &(mp->math->md_negative_inf_t));
1615    mp_scaled_free_number(mp, &(mp->math->md_arc_tol_k));
1616    mp_scaled_free_number(mp, &(mp->math->md_three_sixty_deg_t));
1617    mp_scaled_free_number(mp, &(mp->math->md_one_eighty_deg_t));
1618    mp_scaled_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t));
1619    mp_scaled_free_number(mp, &(mp->math->md_fraction_one_t));
1620    mp_scaled_free_number(mp, &(mp->math->md_fraction_half_t));
1621    mp_scaled_free_number(mp, &(mp->math->md_fraction_three_t));
1622    mp_scaled_free_number(mp, &(mp->math->md_fraction_four_t));
1623    mp_scaled_free_number(mp, &(mp->math->md_zero_t));
1624    mp_scaled_free_number(mp, &(mp->math->md_half_unit_t));
1625    mp_scaled_free_number(mp, &(mp->math->md_three_quarter_unit_t));
1626    mp_scaled_free_number(mp, &(mp->math->md_unity_t));
1627    mp_scaled_free_number(mp, &(mp->math->md_two_t));
1628    mp_scaled_free_number(mp, &(mp->math->md_three_t));
1629    mp_scaled_free_number(mp, &(mp->math->md_one_third_inf_t));
1630    mp_scaled_free_number(mp, &(mp->math->md_warning_limit_t));
1631    mp_scaled_free_number(mp, &(mp->math->md_one_k));
1632    mp_scaled_free_number(mp, &(mp->math->md_sqrt_8_e_k));
1633    mp_scaled_free_number(mp, &(mp->math->md_twelve_ln_2_k));
1634    mp_scaled_free_number(mp, &(mp->math->md_coef_bound_k));
1635    mp_scaled_free_number(mp, &(mp->math->md_coef_bound_minus_1));
1636    mp_scaled_free_number(mp, &(mp->math->md_twelvebits_3));
1637    mp_scaled_free_number(mp, &(mp->math->md_twentysixbits_sqrt2_t));
1638    mp_scaled_free_number(mp, &(mp->math->md_twentyeightbits_d_t));
1639    mp_scaled_free_number(mp, &(mp->math->md_twentysevenbits_sqrt2_d_t));
1640    mp_scaled_free_number(mp, &(mp->math->md_fraction_threshold_t));
1641    mp_scaled_free_number(mp, &(mp->math->md_half_fraction_threshold_t));
1642    mp_scaled_free_number(mp, &(mp->math->md_scaled_threshold_t));
1643    mp_scaled_free_number(mp, &(mp->math->md_half_scaled_threshold_t));
1644    mp_scaled_free_number(mp, &(mp->math->md_near_zero_angle_t));
1645    mp_scaled_free_number(mp, &(mp->math->md_p_over_v_threshold_t));
1646    mp_scaled_free_number(mp, &(mp->math->md_equation_threshold_t));
1647    mp_memory_free(mp->math);
1648}
1649
1650math_data *mp_initialize_scaled_math(MP mp)
1651{
1652    math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data));
1653    /* alloc */
1654    math->md_allocate        = mp_scaled_allocate_number;
1655    math->md_free            = mp_scaled_free_number;
1656    math->md_allocate_clone  = mp_scaled_allocate_clone;
1657    math->md_allocate_abs    = mp_scaled_allocate_abs;
1658    math->md_allocate_double = mp_scaled_allocate_double;
1659    /* precission */
1660    mp_scaled_allocate_number(mp, &math->md_precision_default, mp_scaled_type);
1661    mp_scaled_allocate_number(mp, &math->md_precision_max, mp_scaled_type);
1662    mp_scaled_allocate_number(mp, &math->md_precision_min, mp_scaled_type);
1663    /* here are the constants for |scaled| objects */
1664    mp_scaled_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type);
1665    mp_scaled_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
1666    mp_scaled_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type);
1667    mp_scaled_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
1668    mp_scaled_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type);
1669    mp_scaled_allocate_number(mp, &math->md_unity_t, mp_scaled_type);
1670    mp_scaled_allocate_number(mp, &math->md_two_t, mp_scaled_type);
1671    mp_scaled_allocate_number(mp, &math->md_three_t, mp_scaled_type);
1672    mp_scaled_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type);
1673    mp_scaled_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type);
1674    mp_scaled_allocate_number(mp, &math->md_zero_t, mp_scaled_type);
1675    /* |fractions| */
1676    mp_scaled_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type);
1677    mp_scaled_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type);
1678    mp_scaled_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
1679    mp_scaled_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
1680    mp_scaled_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
1681    /* |angles| */
1682    mp_scaled_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type);
1683    mp_scaled_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type);
1684    mp_scaled_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type);
1685    /* various approximations */
1686    mp_scaled_allocate_number(mp, &math->md_one_k, mp_scaled_type);
1687    mp_scaled_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type);
1688    mp_scaled_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type);
1689    mp_scaled_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type);
1690    mp_scaled_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type);
1691    mp_scaled_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type);
1692    mp_scaled_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type);
1693    mp_scaled_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type);
1694    mp_scaled_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type);
1695    /* thresholds */
1696    mp_scaled_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type);
1697    mp_scaled_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
1698    mp_scaled_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
1699    mp_scaled_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
1700    mp_scaled_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
1701    mp_scaled_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type);
1702    mp_scaled_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
1703    /* initializations */
1704    math->md_precision_default.data.val         = unity * 10;
1705    math->md_precision_max.data.val             = unity * 10;
1706    math->md_precision_min.data.val             = unity * 10;
1707    math->md_epsilon_t.data.val                 = 1;
1708    math->md_inf_t.data.val                     = EL_GORDO;
1709    math->md_negative_inf_t.data.val            = negative_EL_GORDO;
1710    math->md_one_third_inf_t.data.val           = one_third_EL_GORDO;
1711    math->md_warning_limit_t.data.val           = fraction_one;
1712    math->md_unity_t.data.val                   = unity;
1713    math->md_two_t.data.val                     = two;
1714    math->md_three_t.data.val                   = three;
1715    math->md_half_unit_t.data.val               = half_unit;
1716    math->md_three_quarter_unit_t.data.val      = three_quarter_unit;
1717    math->md_arc_tol_k.data.val                 = (unity/4096);
1718    math->md_fraction_one_t.data.val            = fraction_one;
1719    math->md_fraction_half_t.data.val           = fraction_half;
1720    math->md_fraction_three_t.data.val          = fraction_three;
1721    math->md_fraction_four_t.data.val           = fraction_four;
1722    math->md_three_sixty_deg_t.data.val         = three_sixty_deg;
1723    math->md_one_eighty_deg_t.data.val          = one_eighty_deg;
1724    math->md_negative_one_eighty_deg_t.data.val = negative_one_eighty_deg;
1725    math->md_one_k.data.val                     = 1024;
1726    math->md_sqrt_8_e_k.data.val                = 112429;         /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
1727    math->md_twelve_ln_2_k.data.val             = 139548960;      /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
1728    math->md_coef_bound_k.data.val              = coef_bound;
1729    math->md_coef_bound_minus_1.data.val        = coef_bound - 1;
1730    math->md_twelvebits_3.data.val              = 1365;           /* $1365\approx 2^{12}/3$ */
1731    math->md_twentysixbits_sqrt2_t.data.val     = 94906266;       /* $2^{26}\sqrt2\approx94906265.62$ */
1732    math->md_twentyeightbits_d_t.data.val       = 35596755;       /* $2^{28}d\approx35596754.69$ */
1733    math->md_twentysevenbits_sqrt2_d_t.data.val = 25170707;       /* $2^{27}\sqrt2\,d\approx25170706.63$ */
1734    math->md_fraction_threshold_t.data.val      = fraction_threshold;
1735    math->md_half_fraction_threshold_t.data.val = half_fraction_threshold;
1736    math->md_scaled_threshold_t.data.val        = scaled_threshold;
1737    math->md_half_scaled_threshold_t.data.val   = half_scaled_threshold;
1738    math->md_near_zero_angle_t.data.val         = near_zero_angle;
1739    math->md_p_over_v_threshold_t.data.val      = p_over_v_threshold;
1740    math->md_equation_threshold_t.data.val      = equation_threshold;
1741    /* functions */
1742    math->md_from_int                 = mp_scaled_set_from_int;
1743    math->md_from_boolean             = mp_scaled_set_from_boolean;
1744    math->md_from_scaled              = mp_scaled_set_from_scaled;
1745    math->md_from_double              = mp_scaled_set_from_double;
1746    math->md_from_addition            = mp_scaled_set_from_addition;
1747    math->md_half_from_addition       = mp_scaled_set_half_from_addition;
1748    math->md_from_subtraction         = mp_scaled_set_from_subtraction;
1749    math->md_half_from_subtraction    = mp_scaled_set_half_from_subtraction;
1750    math->md_from_of_the_way          = mp_scaled_set_from_of_the_way;
1751    math->md_from_div                 = mp_scaled_set_from_div;
1752    math->md_from_mul                 = mp_scaled_set_from_mul;
1753    math->md_from_int_div             = mp_scaled_set_from_int_div;
1754    math->md_from_int_mul             = mp_scaled_set_from_int_mul;
1755    math->md_negate                   = mp_scaled_negate;
1756    math->md_add                      = mp_scaled_add;
1757    math->md_subtract                 = mp_scaled_subtract;
1758    math->md_half                     = mp_scaled_half;
1759    math->md_do_double                = mp_scaled_double;
1760    math->md_abs                      = mp_scaled_abs;
1761    math->md_clone                    = mp_scaled_clone;
1762    math->md_negated_clone            = mp_scaled_negated_clone;
1763    math->md_abs_clone                = mp_scaled_abs_clone;
1764    math->md_swap                     = mp_scaled_swap;
1765    math->md_add_scaled               = mp_scaled_add_scaled;
1766    math->md_multiply_int             = mp_scaled_multiply_int;
1767    math->md_divide_int               = mp_scaled_divide_int;
1768    math->md_to_int                   = mp_scaled_to_int;
1769    math->md_to_boolean               = mp_scaled_to_boolean;
1770    math->md_to_scaled                = mp_scaled_to_scaled;
1771    math->md_to_double                = mp_scaled_to_double;
1772    math->md_odd                      = mp_scaled_odd;
1773    math->md_equal                    = mp_scaled_equal;
1774    math->md_less                     = mp_scaled_less;
1775    math->md_greater                  = mp_scaled_greater;
1776    math->md_non_equal_abs            = mp_scaled_non_equal_abs;
1777    math->md_round_unscaled           = mp_scaled_round_unscaled;
1778    math->md_floor_scaled             = mp_scaled_floor;
1779    math->md_fraction_to_round_scaled = mp_scaled_fraction_to_round_scaled;
1780    math->md_make_scaled              = mp_scaled_make_scaled;
1781    math->md_make_fraction            = mp_scaled_make_fraction;
1782    math->md_take_fraction            = mp_scaled_take_fraction;
1783    math->md_take_scaled              = mp_scaled_take_scaled;
1784    math->md_velocity                 = mp_scaled_velocity;
1785    math->md_n_arg                    = mp_scaled_n_arg;
1786    math->md_m_log                    = mp_scaled_m_log;
1787    math->md_m_exp                    = mp_scaled_m_exp;
1788    math->md_m_unif_rand              = mp_scaled_m_unif_rand;
1789    math->md_m_norm_rand              = mp_scaled_m_norm_rand;
1790    math->md_pyth_add                 = mp_scaled_pyth_add;
1791    math->md_pyth_sub                 = mp_scaled_pyth_sub;
1792    math->md_power_of                 = mp_scaled_power_of;
1793    math->md_fraction_to_scaled       = mp_scaled_fraction_to_scaled;
1794    math->md_scaled_to_fraction       = mp_scaled_scaled_to_fraction;
1795    math->md_scaled_to_angle          = mp_scaled_scaled_to_angle;
1796    math->md_angle_to_scaled          = mp_scaled_angle_to_scaled;
1797    math->md_init_randoms             = mp_scaled_init_randoms;
1798    math->md_sin_cos                  = mp_scaled_n_sin_cos;
1799    math->md_slow_add                 = mp_scaled_slow_add;
1800    math->md_sqrt                     = mp_scaled_sqrt;
1801    math->md_print                    = mp_scaled_print;
1802    math->md_tostring                 = mp_scaled_tostring;
1803    math->md_modulo                   = mp_scaled_modulo;
1804    math->md_ab_vs_cd                 = mp_scaled_ab_vs_cd;
1805    math->md_crossing_point           = mp_scaled_crossing_point;
1806    math->md_scan_numeric             = mp_scaled_scan_numeric_token;
1807    math->md_scan_fractional          = mp_scaled_scan_fractional_token;
1808    /* housekeeping */
1809    math->md_free_math                = mp_scaled_free_math;
1810    math->md_set_precision            = mp_scaled_set_precision;
1811    return math;
1812}
1813