SourceBrowser

mpmathdouble.c /size: 41 Kb last modification: 2025-02-21 11:03
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/* 
7    todo: check the random code, we might as well use the lua randomizer for all number models. 
8*/
9
10# include "mpmathdouble.h"
11
12# define  PI                      3.1415926535897932384626433832795028841971
13# define  fraction_multiplier     4096.0
14# define  angle_multiplier        16.0
15# define  coef_bound              ((7.0/3.0)*fraction_multiplier)
16# define  fraction_threshold      0.04096
17# define  half_fraction_threshold (fraction_threshold/2)
18# define  scaled_threshold        0.000122
19# define  half_scaled_threshold   (scaled_threshold/2)
20# define  near_zero_angle         (0.0256*angle_multiplier)
21# define  p_over_v_threshold      0x80000
22# define  equation_threshold      0.001
23# define  warning_limit           pow(2.0,52.0)
24# define  epsilon                 pow(2.0,-52.0)
25# define  unity                   1.0
26# define  two                     2.0
27# define  three                   3.0
28# define  half_unit               0.5
29# define  three_quarter_unit      0.75
30# define  EL_GORDO                (DBL_MAX / 2.0 - 1.0)
31# define  negative_EL_GORDO       (-EL_GORDO)
32# define  one_third_EL_GORDO      (EL_GORDO / 3.0)
33# define  fraction_half           (0.5*fraction_multiplier)
34# define  fraction_one            (1.0*fraction_multiplier)
35# define  fraction_two            (2.0*fraction_multiplier)
36# define  fraction_three          (3.0*fraction_multiplier)
37# define  fraction_four           (4.0*fraction_multiplier)
38# define  no_crossing             (fraction_one + 1)
39# define  one_crossing            fraction_one
40# define  zero_crossing           0
41# define  one_eighty_deg          (180.0*angle_multiplier)
42# define  negative_one_eighty_deg (-180.0*angle_multiplier)
43# define  three_sixty_deg         (360.0*angle_multiplier)
44# define  odd(A)                  (abs(A) % 2 == 1)
45# define  two_to_the(A)           (1<<(unsigned)(A))
46# define  set_cur_cmd(A)          mp->cur_mod_->command = (A)
47# define  set_cur_mod(A)          mp->cur_mod_->data.n.data.dval = (A)
48
49/*tex 
50
51Apart from it looking weird in results a |-0.0| can also give wrong results, for instance in a 
52|tan2|, see there for a comment. This is why we now do some checking on zero which in some cases
53also might be a bit more efficient as it avoids a multiplication or division, so likely we break 
54even. 
55
56*/
57
58static inline double mp_double_make_fraction (double p, double q) { 
59 // return (p / q) * fraction_multiplier; 
60    return p == 0.0 ? 0.0 : (p / q) * fraction_multiplier; 
61}
62
63// printf("%f %f %f %f\n",p,q,p*q,(p * q) / fraction_multiplier);
64// -4096.000000 0.000000 -0.000000 -0.000000
65
66static inline double mp_double_take_fraction(double p, double q) { 
67 // return (p * q) / fraction_multiplier; 
68    return p == 0.0 ? 0.0 : q == 0.0 ? 0.0 : (p * q) / fraction_multiplier; 
69}
70
71static inline double mp_double_make_scaled(double p, double q) { 
72 // return p / q; 
73    return p == 0.0 ? 0.0 : p / q; 
74}
75
76static void mp_double_allocate_number(MP mp, mp_number *n, mp_number_type t)
77{
78    (void) mp;
79    n->data.dval = 0.0;
80    n->type = t;
81}
82
83static void mp_double_allocate_clone(MP mp, mp_number *n, mp_number_type t, mp_number *v)
84{
85    (void) mp;
86    n->type = t;
87    n->data.dval = v->data.dval;
88}
89
90static void mp_double_allocate_abs(MP mp, mp_number *n, mp_number_type t, mp_number *v)
91{
92    (void) mp;
93    n->type = t;
94    n->data.dval = fabs(v->data.dval);
95}
96
97static void mp_double_allocate_double(MP mp, mp_number *n, double v)
98{
99    (void) mp;
100    n->type = mp_scaled_type;
101    n->data.dval = v;
102}
103
104static void mp_double_free_number(MP mp, mp_number *n)
105{
106    (void) mp;
107    n->type = mp_nan_type;
108}
109
110static void mp_double_set_from_int(mp_number *A, int B)
111{
112    A->data.dval = B;
113}
114
115static void mp_double_set_from_boolean(mp_number *A, int B)
116{
117    A->data.dval = B;
118}
119
120static void mp_double_set_from_scaled(mp_number *A, int B)
121{
122    A->data.dval = B == 0 ? 0.0 : B / 65536.0;
123}
124
125static void mp_double_set_from_double(mp_number *A, double B)
126{
127    A->data.dval = B;
128}
129
130static void mp_double_set_from_addition(mp_number *A, mp_number *B, mp_number *C)
131{
132    A->data.dval = B->data.dval + C->data.dval;
133}
134
135static void mp_double_set_half_from_addition(mp_number *A, mp_number *B, mp_number *C)
136{
137    A->data.dval = (B->data.dval + C->data.dval) / 2.0;
138}
139
140static void mp_double_set_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
141{
142    A->data.dval = B->data.dval - C->data.dval;
143}
144
145static void mp_double_set_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
146{
147    A->data.dval = (B->data.dval - C->data.dval) / 2.0;
148}
149
150static void mp_double_set_from_div(mp_number *A, mp_number *B, mp_number *C)
151{
152    A->data.dval = B->data.dval == 0.0 ? 0.0 : B->data.dval / C->data.dval;
153}
154
155static void mp_double_set_from_mul(mp_number *A, mp_number *B, mp_number *C)
156{
157    A->data.dval = B->data.dval == 0.0 || C->data.dval == 0.0 ? 0.0 : B->data.dval * C->data.dval;
158}
159
160static void mp_double_set_from_int_div(mp_number *A, mp_number *B, int C)
161{
162    A->data.dval = B->data.dval == 0.0 ? 0.0 : B->data.dval / C;
163}
164
165static void mp_double_set_from_int_mul(mp_number *A, mp_number *B, int C)
166{
167    A->data.dval = B->data.dval == 0.0 || C == 0 ? 0.0 : B->data.dval * C;
168}
169
170static void mp_double_set_from_of_the_way(MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C)
171{
172    (void) mp;
173    A->data.dval = B->data.dval - mp_double_take_fraction(B->data.dval - C->data.dval, t->data.dval);
174}
175
176static void mp_double_negate(mp_number *A)
177{
178    if (A->data.dval == -0.0) { /* already checked */
179        A->data.dval = 0.0;
180    } else if (A->data.dval != 0.0) {
181        A->data.dval = -A->data.dval;
182    }
183}
184
185static void mp_double_add(mp_number *A, mp_number *B)
186{
187    A->data.dval += B->data.dval;
188}
189
190static void mp_double_subtract(mp_number *A, mp_number *B)
191{
192    A->data.dval -= B->data.dval;
193}
194
195static void mp_double_half(mp_number *A)
196{
197    if (A->data.dval != 0.0) {
198        A->data.dval /= 2.0;
199    }
200}
201
202static void mp_double_double(mp_number *A)
203{
204    if (A->data.dval != 0.0) {
205        A->data.dval *= 2.0;
206    }
207}
208
209static void mp_double_add_scaled(mp_number *A, int B)
210{
211    /* also for negative B */
212    A->data.dval += (B / 65536.0);
213}
214
215static void mp_double_multiply_int(mp_number *A, int B)
216{
217    if (A->data.dval != 0.0) {
218        A->data.dval *= (double) B;
219    }
220}
221
222static void mp_double_divide_int(mp_number *A, int B)
223{
224    if (A->data.dval != 0.0) {
225        A->data.dval /= (double) B;
226    }
227}
228
229static void mp_double_abs(mp_number *A)
230{
231    A->data.dval = fabs(A->data.dval);
232}
233
234static void mp_double_clone(mp_number *A, mp_number *B)
235{
236    A->data.dval = B->data.dval;
237}
238
239static void mp_double_negated_clone(mp_number *A, mp_number *B)
240{
241    A->data.dval = -B->data.dval;
242    if (A->data.dval == -0.0) {
243        A->data.dval = 0.0;
244    }
245}
246
247static void mp_double_abs_clone(mp_number *A, mp_number *B)
248{
249    A->data.dval = fabs(B->data.dval);
250}
251
252static void mp_double_swap(mp_number *A, mp_number *B)
253{
254    double swap_tmp = A->data.dval;
255    A->data.dval = B->data.dval;
256    B->data.dval = swap_tmp;
257}
258
259static void mp_double_fraction_to_scaled(mp_number *A)
260{
261    A->type = mp_scaled_type;
262    if (A->data.dval != 0.0) {
263        A->data.dval /= fraction_multiplier;
264    }
265}
266
267static void mp_double_angle_to_scaled(mp_number *A)
268{
269    A->type = mp_scaled_type;
270    if (A->data.dval != 0.0) {
271        A->data.dval /= angle_multiplier;
272    }
273}
274
275static void mp_double_scaled_to_fraction(mp_number *A)
276{
277    A->type = mp_fraction_type;
278    if (A->data.dval != 0.0) {
279        A->data.dval *= fraction_multiplier;
280    }
281}
282
283static void mp_double_scaled_to_angle(mp_number *A)
284{
285    A->type = mp_angle_type;
286    if (A->data.dval != 0.0) {
287        A->data.dval *= angle_multiplier;
288    }
289}
290
291static int mp_double_to_scaled(mp_number *A)
292{
293    return (int) lround(A->data.dval * 65536.0);
294}
295
296static int mp_double_to_int(mp_number *A)
297{
298    return (int) (A->data.dval);
299}
300
301static int mp_double_to_boolean(mp_number *A)
302{
303    return (int) (A->data.dval);
304}
305
306static double mp_double_to_double(mp_number *A)
307{
308    return A->data.dval;
309}
310
311static int mp_double_odd(mp_number *A)
312{
313    return odd((int) lround(A->data.dval));
314}
315
316static int mp_double_equal(mp_number *A, mp_number *B)
317{
318    return A->data.dval == B->data.dval;
319}
320
321static int mp_double_greater(mp_number *A, mp_number *B)
322{
323    return A->data.dval > B->data.dval;
324}
325
326static int mp_double_less(mp_number *A, mp_number *B)
327{
328    return A->data.dval < B->data.dval;
329}
330
331static int mp_double_non_equal_abs(mp_number *A, mp_number *B)
332{
333    return fabs(A->data.dval) != fabs(B->data.dval);
334}
335
336static char *mp_double_number_tostring(MP mp, mp_number *n)
337{
338    static char set[64];
339    int l = 0;
340    char *ret = mp_memory_allocate(64);
341    (void) mp;
342    snprintf(set, 64, mp->less_digits ? "%.3g" : "%.17g", n->data.dval);
343    while (set[l] == ' ') {
344        l++;
345    }
346    strcpy(ret, set+l);
347    return ret;
348}
349
350static void mp_double_print_number(MP mp, mp_number *n)
351{
352    char *str = mp_double_number_tostring(mp, n);
353    mp_print_e_str(mp, str);
354    mp_memory_free(str);
355}
356
357static void mp_double_slow_add(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
358{
359    double x = x_orig->data.dval;
360    double y = y_orig->data.dval;
361    if (x >= 0.0) {
362        if (y <= EL_GORDO - x) {
363            ret->data.dval = x + y;
364        } else {
365            mp->arith_error = 1;
366            ret->data.dval = EL_GORDO;
367        }
368    } else if (-y <= EL_GORDO + x) {
369        ret->data.dval = x + y;
370    } else {
371        mp->arith_error = 1;
372        ret->data.dval = negative_EL_GORDO;
373    }
374}
375
376static void mp_double_number_make_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q) {
377    (void) mp;
378    ret->data.dval = mp_double_make_fraction(p->data.dval, q->data.dval);
379}
380
381static void mp_double_number_take_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q) {
382   (void) mp;
383   ret->data.dval = mp_double_take_fraction(p->data.dval, q->data.dval);
384}
385
386static void mp_double_number_take_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
387{
388    (void) mp;
389    ret->data.dval = p_orig->data.dval * q_orig->data.dval;
390}
391
392static void mp_double_number_make_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
393{
394    (void) mp;
395    ret->data.dval = p_orig->data.dval / q_orig->data.dval;
396}
397
398static void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop)
399{
400    double result;
401    char *end = (char *) stop;
402    errno = 0;
403    result = strtod((char *) start, &end);
404    if (errno == 0) {
405        set_cur_mod(result);
406        if (result >= warning_limit) {
407            if (internal_value(mp_warning_check_internal).data.dval > 0 && (mp->scanner_status != mp_tex_flushing_state)) {
408                char msg[256];
409                snprintf(msg, 256, "Number is too large (%g)", result);
410                mp_error(
411                    mp,
412                    msg,
413                    "Continue and I'll try to cope with that big value; but it might be dangerous."
414                    "(Set warningcheck := 0 to suppress this message.)"
415                );
416            }
417        }
418    } else if (mp->scanner_status != mp_tex_flushing_state) {
419        mp_error(
420            mp,
421            "Enormous number has been reduced.",
422            "I could not handle this number specification probably because it is out of"
423            "range."
424        );
425        set_cur_mod(EL_GORDO);
426    }
427    set_cur_cmd(mp_numeric_command);
428}
429
430static void mp_double_aux_find_exponent(MP mp)
431{
432    if (mp->buffer[mp->cur_input.loc_field] == 'e' || mp->buffer[mp->cur_input.loc_field] == 'E') {
433        mp->cur_input.loc_field++;
434        if (!(mp->buffer[mp->cur_input.loc_field] == '+'
435           || mp->buffer[mp->cur_input.loc_field] == '-'
436           || mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class)) {
437            mp->cur_input.loc_field--;
438            return;
439        }
440        if (mp->buffer[mp->cur_input.loc_field] == '+'
441         || mp->buffer[mp->cur_input.loc_field] == '-') {
442            mp->cur_input.loc_field++;
443        }
444        while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
445            mp->cur_input.loc_field++;
446        }
447    }
448}
449
450static void mp_double_scan_fractional_token(MP mp, int n) /* n is scaled */
451{
452    unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
453    unsigned char *stop;
454    (void) n;
455    while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
456        mp->cur_input.loc_field++;
457    }
458    mp_double_aux_find_exponent(mp);
459    stop = &mp->buffer[mp->cur_input.loc_field-1];
460    mp_wrapup_numeric_token(mp, start, stop);
461}
462
463/*tex
464    The input format is the same as for the C language, so we just collect valid bytes in the 
465    buffer, then call |strtod()|. It looks like we have no buffer overflow check here. 
466*/
467
468static void mp_double_scan_numeric_token(MP mp, int n) /* n is scaled */
469{
470    unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
471    unsigned char *stop;
472    (void) n;
473    while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
474        mp->cur_input.loc_field++;
475    }
476    if (mp->buffer[mp->cur_input.loc_field] == '.' && mp->buffer[mp->cur_input.loc_field+1] != '.') {
477        mp->cur_input.loc_field++;
478        while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
479            mp->cur_input.loc_field++;
480        }
481    }
482    mp_double_aux_find_exponent(mp);
483    stop = &mp->buffer[mp->cur_input.loc_field-1];
484    mp_wrapup_numeric_token(mp, start, stop);
485}
486
487static void mp_double_velocity(MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t)
488{
489    double acc, num, denom; /* registers for intermediate calculations */
490    (void) mp;
491    acc = mp_double_take_fraction(st->data.dval - (sf->data.dval / 16.0), sf->data.dval - (st->data.dval / 16.0));
492    acc = mp_double_take_fraction(acc, ct->data.dval - cf->data.dval);
493    num = fraction_two + mp_double_take_fraction(acc, sqrt(2)*fraction_one);
494    denom = fraction_three
495          + mp_double_take_fraction(ct->data.dval, 3*fraction_half*(sqrt(5.0)-1.0))
496          + mp_double_take_fraction(cf->data.dval, 3*fraction_half*(3.0-sqrt(5.0)));
497    if (t->data.dval != unity) {
498        num = mp_double_make_scaled(num, t->data.dval);
499    }
500    if (num / 4 >= denom) {
501        ret->data.dval = fraction_four;
502    } else {
503        ret->data.dval = mp_double_make_fraction(num, denom);
504    }
505}
506
507static int mp_double_ab_vs_cd(mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig)
508{
509    double ab = a_orig->data.dval * b_orig->data.dval;
510    double cd = c_orig->data.dval * d_orig->data.dval;
511    if (ab > cd) {
512        return 1;
513    } else if (ab < cd) {
514        return -1;
515    } else {
516        return 0;
517    }
518}
519
520static void mp_double_crossing_point(MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc)
521{
522    double a = aa->data.dval;
523    double b = bb->data.dval;
524    double c = cc->data.dval;
525    (void) mp;
526    if (a < 0.0) {
527        ret->data.dval = zero_crossing;
528        return;
529    } else if (c >= 0.0) {
530        if (b >= 0.0) {
531            if (c > 0.0) {
532                ret->data.dval = no_crossing;
533            } else if ((a == 0.0) && (b == 0.0)) {
534                ret->data.dval = no_crossing;
535            } else {
536                ret->data.dval = one_crossing;
537            }
538            return;
539        } else if (a == 0.0) {
540            ret->data.dval = zero_crossing;
541            return;
542        }
543    } else if ((a == 0.0) && (b <= 0.0)) {
544        ret->data.dval = zero_crossing;
545        return;
546    }
547    /* Use bisection to find the crossing point... */
548    {
549        double d = epsilon;
550        double x0 = a;
551        double x1 = a - b;
552        double x2 = b - c;
553        do {
554            /* not sure why the error correction has to be >= 1E-12 */
555            double x = (x1 + x2) / 2 + 1E-12;
556            if (x1 - x0 > x0) {
557                x2 = x;
558                x0 += x0;
559                d += d;
560            } else {
561                double xx = x1 + x - x0;
562                if (xx > x0) {
563                    x2 = x;
564                    x0 += x0;
565                    d += d;
566                } else {
567                    x0 = x0 - xx;
568                    if ((x <= x0) && (x + x2 <= x0)) {
569                        ret->data.dval = no_crossing;
570                        return;
571                    }
572                    x1 = x;
573                    d = d + d + epsilon;
574                }
575            }
576        } while (d < fraction_one);
577        ret->data.dval = (d - fraction_one);
578    }
579}
580
581/*tex 
582
583    This is the same one as in |mpmath.c|:
584
585*/
586
587static int mp_double_unscaled(mp_number *x_orig)
588{
589    return (int) lround(x_orig->data.dval);
590}
591
592static void mp_double_floor(mp_number *i)
593{
594    i->data.dval = floor(i->data.dval);
595}
596
597static void mp_double_fraction_to_round_scaled(mp_number *x_orig)
598{
599    double x = x_orig->data.dval;
600    x_orig->type = mp_scaled_type;
601    x_orig->data.dval = x / fraction_multiplier;
602}
603
604static void mp_double_square_rt(MP mp, mp_number *ret, mp_number *x_orig) /* return, x: scaled */
605{
606    double x = x_orig->data.dval;
607    if (x > 0) {
608        ret->data.dval = sqrt(x);
609    } else {
610        if (x < 0) {
611            char msg[256];
612            char *xstr = mp_double_number_tostring(mp, x_orig);
613            snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr);
614            mp_memory_free(xstr);
615            mp_error(
616                mp,
617                msg,
618                "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
619                "Proceed, with fingers crossed."
620            );
621        }
622        ret->data.dval = 0;
623    }
624}
625
626static void mp_double_pyth_add(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
627{
628    double a = fabs(a_orig->data.dval);
629    double b = fabs(b_orig->data.dval);
630    errno = 0;
631    ret->data.dval = sqrt(a*a + b*b);
632    if (errno) {
633        mp->arith_error = 1;
634        ret->data.dval = EL_GORDO;
635    }
636}
637
638static void mp_double_pyth_sub(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
639{
640    double a = fabs(a_orig->data.dval);
641    double b = fabs(b_orig->data.dval);
642    if (a > b) {
643        a = sqrt(a*a - b*b);
644    } else {
645        if (a < b) {
646            char msg[256];
647            char *astr = mp_double_number_tostring(mp, a_orig);
648            char *bstr = mp_double_number_tostring(mp, b_orig);
649            snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, bstr);
650            mp_memory_free(astr);
651            mp_memory_free(bstr);
652            mp_error(
653                mp,
654                msg,
655                "Since I don't take square roots of negative numbers, Im zeroing this one.\n"
656                "Proceed, with fingers crossed."
657            );
658        }
659        a = 0;
660    }
661    ret->data.dval = a;
662}
663
664static void mp_double_power_of(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
665{
666    errno = 0;
667    ret->data.dval = pow(a_orig->data.dval, b_orig->data.dval);
668    if (errno) {
669        mp->arith_error = 1;
670        ret->data.dval = EL_GORDO;
671    }
672}
673
674static void mp_double_m_log(MP mp, mp_number *ret, mp_number *x_orig)
675{
676    if (x_orig->data.dval > 0) {
677        ret->data.dval = log(x_orig->data.dval)*256.0;
678    } else {
679        char msg[256];
680        char *xstr = mp_double_number_tostring(mp, x_orig);
681        snprintf(msg, 256, "Logarithm of %s has been replaced by 0", xstr);
682        mp_memory_free(xstr);
683        mp_error(
684            mp,
685            msg,
686            "Since I don't take logs of non-positive numbers, I'm zeroing this one.\n"
687            "Proceed, with fingers crossed."
688        );
689        ret->data.dval = 0;
690    }
691}
692
693static void mp_double_m_exp(MP mp, mp_number *ret, mp_number *x_orig)
694{
695    errno = 0;
696    ret->data.dval = exp(x_orig->data.dval / 256.0);
697    if (errno) {
698        if (x_orig->data.dval > 0) {
699            mp->arith_error = 1;
700            ret->data.dval = EL_GORDO;
701        } else {
702            ret->data.dval = 0;
703        }
704    }
705}
706
707/*tex
708
709    We (MS/HH) ran into an issue @ 2024-03-09 that showed up because a file from early 2023 gave 
710    different results: |draw (50,0) {dir 120} .. (-50,0)  ..{dir 60} (50,0);| and in the end the 
711    issue was due to a |-0| not making |atan2| happy. It was ok in scaled (what was the default
712    in \CONTEXT\ before 2024), posit and decimal, but failed in double (which became default in 
713    mid 2023) and in binary (both also tested in \LUATEX). 
714
715    We now try to catch the |-0.0| upstream but we keep the patch commented. However, it was 
716    decided that in mplib (in \LUAMETATEX) we patch the |n_arg| function, so there the commented 
717    code below is used. 
718
719*/
720
721static void mp_double_n_arg(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
722{
723    /*        
724        if (x_orig->data.dval == -0.0) {
725            x_orig->data.dval = 0.0;
726        }
727        if (y_orig->data.dval == -0.0) {
728            y_orig->data.dval = 0.0;
729        }
730    */
731    if (x_orig->data.dval == 0.0 && y_orig->data.dval == 0.0) {
732        if (internal_value(mp_default_zero_angle_internal).data.dval < 0) {
733            mp_error(
734                mp,
735                "angle(0,0) is taken as zero",
736                "The 'angle' between two identical points is undefined. I'm zeroing this one.\n"
737                "Proceed, with fingers crossed."
738            );
739            ret->data.dval = internal_value(mp_default_zero_angle_internal).data.dval;
740        } else { 
741            ret->data.dval = 0;
742        }
743    } else {
744        ret->type = mp_angle_type;
745        ret->data.dval = atan2(y_orig->data.dval, x_orig->data.dval) * (180.0 / PI) * angle_multiplier;
746        if (ret->data.dval == -0.0) {
747            ret->data.dval = 0.0;
748        }
749// printf("D x=%f y=%f atan=%f\n",y_orig->data.dval,x_orig->data.dval,ret->data.dval);
750    }
751}
752
753static void mp_double_sin_cos(MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin)
754{
755    double rad = (z_orig->data.dval / angle_multiplier); /* still degrees */
756    (void) mp;
757    if ((rad == 90.0) || (rad == -270)){
758        n_cos->data.dval = 0.0;
759        n_sin->data.dval = fraction_multiplier;
760    } else if ((rad == -90.0) || (rad == 270.0)) {
761        n_cos->data.dval = 0.0;
762        n_sin->data.dval = -fraction_multiplier;
763    } else if ((rad == 180.0) || (rad == -180.0)) {
764        n_cos->data.dval = -fraction_multiplier;
765        n_sin->data.dval = 0.0;
766    } else {
767        rad = rad * PI / 180.0;
768        n_cos->data.dval = cos(rad) * fraction_multiplier;
769        n_sin->data.dval = sin(rad) * fraction_multiplier;
770    }
771}
772
773/*tex
774    This is the |http://www-cs-faculty.stanford.edu/~uno/programs/rng.c| with small cosmetic 
775    modifications. The code only documented here as the other non scaled number models use the 
776    same method.
777*/
778
779# define KK            100                /* the long lag  */
780# define LL            37                 /* the short lag */
781# define MM            (1L<<30)           /* the modulus   */
782# define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
783# define TT            70                 /* guaranteed separation between streams */
784# define is_odd(x)     ((x)&1)            /* units bit of x */
785# define QUALITY       1009               /* recommended quality level for high-res use */
786
787/*tex 
788    The destination array length (must be at least KK).
789*/
790
791typedef struct mp_double_random_info {
792    long  x[KK];
793    long  buf[QUALITY];
794    long  dummy;
795    long  started;
796    long *ptr;
797} mp_double_random_info;
798
799static mp_double_random_info mp_double_random_data = {
800    .dummy   = -1,
801    .started = -1,
802    .ptr     = &mp_double_random_data.dummy
803};
804
805static void mp_double_aux_ran_array(long aa[], int n)
806{
807    int i, j;
808    for (j = 0; j < KK; j++) {
809        aa[j] = mp_double_random_data.x[j];
810    }
811    for (; j < n; j++) {
812        aa[j] = mod_diff(aa[j - KK], aa[j - LL]);
813    }
814    for (i = 0; i < LL; i++, j++) {
815        mp_double_random_data.x[i] = mod_diff(aa[j - KK], aa[j - LL]);
816    }
817    for (; i < KK; i++, j++) {
818        mp_double_random_data.x[i] = mod_diff(aa[j - KK], mp_double_random_data.x[i - LL]);
819    }
820}
821
822static void mp_double_aux_ran_start(long seed)
823{
824    int t, j;
825    long x[KK + KK - 1]; /* the preparation buffer */
826    long ss = (seed+2) & (MM - 2);
827    for (j = 0; j < KK; j++) {
828        /* bootstrap the buffer */
829        x[j] = ss;
830        /* cyclic shift 29 bits */
831        ss <<= 1;
832        if (ss >= MM) {
833            ss -= MM - 2;
834        }
835    }
836    /* make x[1] (and only x[1]) odd */
837    x[1]++;
838    for (ss = seed & (MM - 1), t = TT - 1; t;) {
839        for (j = KK - 1; j > 0; j--) {
840            /* "square" */
841            x[j + j] = x[j];
842            x[j + j - 1] = 0;
843        }
844        for (j = KK + KK - 2; j >= KK; j--) {
845            x[j - (KK -LL)] = mod_diff(x[j - (KK - LL)], x[j]);
846            x[j - KK] = mod_diff(x[j - KK], x[j]);
847        }
848        if (is_odd(ss)) {
849            /* "multiply by z" */
850            for (j = KK; j>0; j--) {
851                x[j] = x[j-1];
852            }
853            x[0] = x[KK];
854            /* shift the buffer cyclically */
855            x[LL] = mod_diff(x[LL], x[KK]);
856        }
857        if (ss) {
858            ss >>= 1;
859        } else {
860            t--;
861        }
862    }
863    for (j = 0; j < LL; j++) {
864        mp_double_random_data.x[j + KK - LL] = x[j];
865    }
866    for (;j < KK; j++) {
867        mp_double_random_data.x[j - LL] = x[j];
868    }
869    for (j = 0; j < 10; j++) {
870        /* warm things up */
871        mp_double_aux_ran_array(x, KK + KK - 1);
872    }
873    mp_double_random_data.ptr = &mp_double_random_data.started;
874}
875
876static long mp_double_aux_ran_arr_cycle(void)
877{
878    if (mp_double_random_data.ptr == &mp_double_random_data.dummy) {
879        /* the user forgot to initialize */
880        mp_double_aux_ran_start(314159L);
881    }
882    mp_double_aux_ran_array(mp_double_random_data.buf, QUALITY);
883    mp_double_random_data.buf[KK] = -1;
884    mp_double_random_data.ptr = mp_double_random_data.buf + 1;
885    return mp_double_random_data.buf[0];
886}
887
888static void mp_double_init_randoms(MP mp, int seed)
889{
890    int k = 1;
891    int j = abs(seed);
892    int f = (int) fraction_one; /* avoid warnings */
893    while (j >= f) {
894        j = j / 2;
895    }
896    for (int i = 0; i <= 54; i++) {
897        int jj = k;
898        k = j - k;
899        j = jj;
900        if (k < 0) {
901            k += f;
902        }
903        mp->randoms[(i * 21) % 55].data.dval = j;
904    }
905    mp_new_randoms(mp);
906    mp_new_randoms(mp);
907    mp_new_randoms(mp);
908    /* warm up the array */
909    mp_double_aux_ran_start((unsigned long) seed);
910}
911
912static void mp_double_modulo(mp_number *a, mp_number *b)
913{
914    double tmp;
915    a->data.dval = modf((double) a->data.dval / (double) b->data.dval, &tmp) * (double) b->data.dval;
916}
917
918static void mp_double_aux_next_unif_random(MP mp, mp_number *ret)
919{
920    unsigned long int op = (unsigned) (*mp_double_random_data.ptr>=0? *mp_double_random_data.ptr++: mp_double_aux_ran_arr_cycle());
921    double a = op / (MM * 1.0);
922    (void) mp;
923    ret->data.dval = a;
924}
925
926static void mp_double_aux_next_random(MP mp, mp_number *ret)
927{
928    if ( mp->j_random==0) {
929        mp_new_randoms(mp);
930    } else {
931        mp->j_random = mp->j_random-1;
932    }
933    mp_double_clone(ret, &(mp->randoms[mp->j_random]));
934}
935
936static void mp_double_m_unif_rand(MP mp, mp_number *ret, mp_number *x_orig)
937{
938    mp_number x, abs_x, u, y; /* |y| is trial value */
939    mp_double_allocate_number(mp, &y, mp_fraction_type);
940    mp_double_allocate_clone(mp, &x, mp_scaled_type, x_orig);
941    mp_double_allocate_abs(mp, &abs_x, mp_scaled_type, &x);
942    mp_double_allocate_number(mp, &u, mp_scaled_type);
943    mp_double_aux_next_unif_random(mp, &u);
944    y.data.dval = abs_x.data.dval * u.data.dval;
945    mp_double_free_number(mp, &u);
946    if (mp_double_equal(&y, &abs_x)) {
947        mp_double_clone(ret, &((math_data *)mp->math)->md_zero_t);
948    } else if (mp_double_greater(&x, &((math_data *)mp->math)->md_zero_t)) {
949        mp_double_clone(ret, &y);
950    } else {
951        mp_double_negated_clone(ret, &y);
952    }
953    mp_double_free_number(mp, &abs_x);
954    mp_double_free_number(mp, &x);
955    mp_double_free_number(mp, &y);
956}
957
958static void mp_double_m_norm_rand(MP mp, mp_number *ret)
959{
960    mp_number abs_x, u, r, la, xa;
961    mp_double_allocate_number(mp, &la, mp_scaled_type);
962    mp_double_allocate_number(mp, &xa, mp_scaled_type);
963    mp_double_allocate_number(mp, &abs_x, mp_scaled_type);
964    mp_double_allocate_number(mp, &u, mp_scaled_type);
965    mp_double_allocate_number(mp, &r, mp_scaled_type);
966    do {
967        do {
968            mp_number v;
969            mp_double_allocate_number(mp, &v, mp_scaled_type);
970            mp_double_aux_next_random(mp, &v);
971            mp_double_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t);
972            mp_double_number_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v);
973            mp_double_free_number(mp, &v);
974            mp_double_aux_next_random(mp, &u);
975            mp_double_clone(&abs_x, &xa);
976            mp_double_abs(&abs_x);
977        } while (! mp_double_less(&abs_x, &u));
978        mp_double_number_make_fraction(mp, &r, &xa, &u);
979        mp_double_clone(&xa, &r);
980        mp_double_m_log(mp, &la, &u);
981        mp_double_set_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la);
982    } while (mp_double_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0);
983    mp_double_clone(ret, &xa);
984    mp_double_free_number(mp, &r);
985    mp_double_free_number(mp, &abs_x);
986    mp_double_free_number(mp, &la);
987    mp_double_free_number(mp, &xa);
988    mp_double_free_number(mp, &u);
989}
990
991static void mp_double_set_precision(MP mp)
992{
993    (void) mp;
994}
995
996static void mp_double_free_math(MP mp)
997{
998    mp_double_free_number(mp, &(mp->math->md_three_sixty_deg_t));
999    mp_double_free_number(mp, &(mp->math->md_one_eighty_deg_t));
1000    mp_double_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t));
1001    mp_double_free_number(mp, &(mp->math->md_fraction_one_t));
1002    mp_double_free_number(mp, &(mp->math->md_zero_t));
1003    mp_double_free_number(mp, &(mp->math->md_half_unit_t));
1004    mp_double_free_number(mp, &(mp->math->md_three_quarter_unit_t));
1005    mp_double_free_number(mp, &(mp->math->md_unity_t));
1006    mp_double_free_number(mp, &(mp->math->md_two_t));
1007    mp_double_free_number(mp, &(mp->math->md_three_t));
1008    mp_double_free_number(mp, &(mp->math->md_one_third_inf_t));
1009    mp_double_free_number(mp, &(mp->math->md_inf_t));
1010    mp_double_free_number(mp, &(mp->math->md_negative_inf_t));
1011    mp_double_free_number(mp, &(mp->math->md_warning_limit_t));
1012    mp_double_free_number(mp, &(mp->math->md_one_k));
1013    mp_double_free_number(mp, &(mp->math->md_sqrt_8_e_k));
1014    mp_double_free_number(mp, &(mp->math->md_twelve_ln_2_k));
1015    mp_double_free_number(mp, &(mp->math->md_coef_bound_k));
1016    mp_double_free_number(mp, &(mp->math->md_coef_bound_minus_1));
1017    mp_double_free_number(mp, &(mp->math->md_fraction_threshold_t));
1018    mp_double_free_number(mp, &(mp->math->md_half_fraction_threshold_t));
1019    mp_double_free_number(mp, &(mp->math->md_scaled_threshold_t));
1020    mp_double_free_number(mp, &(mp->math->md_half_scaled_threshold_t));
1021    mp_double_free_number(mp, &(mp->math->md_near_zero_angle_t));
1022    mp_double_free_number(mp, &(mp->math->md_p_over_v_threshold_t));
1023    mp_double_free_number(mp, &(mp->math->md_equation_threshold_t));
1024    mp_memory_free(mp->math);
1025}
1026
1027math_data *mp_initialize_double_math(MP mp)
1028{
1029    math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data));
1030    /* alloc */
1031    math->md_allocate        = mp_double_allocate_number;
1032    math->md_free            = mp_double_free_number;
1033    math->md_allocate_clone  = mp_double_allocate_clone;
1034    math->md_allocate_abs    = mp_double_allocate_abs;
1035    math->md_allocate_double = mp_double_allocate_double;
1036    /* precission */
1037    mp_double_allocate_number(mp, &math->md_precision_default, mp_scaled_type);
1038    mp_double_allocate_number(mp, &math->md_precision_max, mp_scaled_type);
1039    mp_double_allocate_number(mp, &math->md_precision_min, mp_scaled_type);
1040    /* here are the constants for |scaled| objects */
1041    mp_double_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type);
1042    mp_double_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
1043    mp_double_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type);
1044    mp_double_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
1045    mp_double_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type);
1046    mp_double_allocate_number(mp, &math->md_unity_t, mp_scaled_type);
1047    mp_double_allocate_number(mp, &math->md_two_t, mp_scaled_type);
1048    mp_double_allocate_number(mp, &math->md_three_t, mp_scaled_type);
1049    mp_double_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type);
1050    mp_double_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type);
1051    mp_double_allocate_number(mp, &math->md_zero_t, mp_scaled_type);
1052    /* |fractions| */
1053    mp_double_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type);
1054    mp_double_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type);
1055    mp_double_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
1056    mp_double_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
1057    mp_double_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
1058    /* |angles| */
1059    mp_double_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type);
1060    mp_double_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type);
1061    mp_double_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type);
1062    /* various approximations */
1063    mp_double_allocate_number(mp, &math->md_one_k, mp_scaled_type);
1064    mp_double_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type);
1065    mp_double_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type);
1066    mp_double_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type);
1067    mp_double_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type);
1068    mp_double_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type);
1069    mp_double_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type);
1070    mp_double_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type);
1071    mp_double_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type);
1072    /* thresholds */
1073    mp_double_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type);
1074    mp_double_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
1075    mp_double_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
1076    mp_double_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
1077    mp_double_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
1078    mp_double_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type);
1079    mp_double_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
1080    /* initializations */
1081    math->md_precision_default.data.dval         = 16 * unity; 
1082    math->md_precision_max.data.dval             = 16 * unity;
1083    math->md_precision_min.data.dval             = 16 * unity;
1084    math->md_epsilon_t.data.dval                 = epsilon;
1085    math->md_inf_t.data.dval                     = EL_GORDO;
1086    math->md_negative_inf_t.data.dval            = negative_EL_GORDO;
1087    math->md_one_third_inf_t.data.dval           = one_third_EL_GORDO;
1088    math->md_warning_limit_t.data.dval           = warning_limit;
1089    math->md_unity_t.data.dval                   = unity;
1090    math->md_two_t.data.dval                     = two;
1091    math->md_three_t.data.dval                   = three;
1092    math->md_half_unit_t.data.dval               = half_unit;
1093    math->md_three_quarter_unit_t.data.dval      = three_quarter_unit;
1094    math->md_arc_tol_k.data.dval                 = (unity/4096);                /* quit when change in arc length estimate reaches this */
1095    math->md_fraction_one_t.data.dval            = fraction_one;
1096    math->md_fraction_half_t.data.dval           = fraction_half;
1097    math->md_fraction_three_t.data.dval          = fraction_three;
1098    math->md_fraction_four_t.data.dval           = fraction_four;
1099    math->md_three_sixty_deg_t.data.dval         = three_sixty_deg;
1100    math->md_one_eighty_deg_t.data.dval          = one_eighty_deg;
1101    math->md_negative_one_eighty_deg_t.data.dval = negative_one_eighty_deg;
1102    math->md_one_k.data.dval                     = 1.0/64 ;
1103    math->md_sqrt_8_e_k.data.dval                = 1.71552776992141359295;       /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
1104    math->md_twelve_ln_2_k.data.dval             = 8.31776616671934371292 * 256; /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
1105    math->md_coef_bound_k.data.dval              = coef_bound;
1106    math->md_coef_bound_minus_1.data.dval        = coef_bound - 1/65536.0;
1107    math->md_twelvebits_3.data.dval              =     1365 / 65536.0;           /* $1365\approx 2^{12}/3$ */
1108    math->md_twentysixbits_sqrt2_t.data.dval     = 94906266 / 65536.0;           /* $2^{26}\sqrt2\approx94906265.62$ */
1109    math->md_twentyeightbits_d_t.data.dval       = 35596755 / 65536.0;           /* $2^{28}d\approx35596754.69$ */
1110    math->md_twentysevenbits_sqrt2_d_t.data.dval = 25170707 / 65536.0;           /* $2^{27}\sqrt2\,d\approx25170706.63$ */
1111    math->md_fraction_threshold_t.data.dval      = fraction_threshold;
1112    math->md_half_fraction_threshold_t.data.dval = half_fraction_threshold;
1113    math->md_scaled_threshold_t.data.dval        = scaled_threshold;
1114    math->md_half_scaled_threshold_t.data.dval   = half_scaled_threshold;
1115    math->md_near_zero_angle_t.data.dval         = near_zero_angle;
1116    math->md_p_over_v_threshold_t.data.dval      = p_over_v_threshold;
1117    math->md_equation_threshold_t.data.dval      = equation_threshold;
1118    /* functions */
1119    math->md_from_int                 = mp_double_set_from_int;
1120    math->md_from_boolean             = mp_double_set_from_boolean;
1121    math->md_from_scaled              = mp_double_set_from_scaled;
1122    math->md_from_double              = mp_double_set_from_double;
1123    math->md_from_addition            = mp_double_set_from_addition;
1124    math->md_half_from_addition       = mp_double_set_half_from_addition;
1125    math->md_from_subtraction         = mp_double_set_from_subtraction;
1126    math->md_half_from_subtraction    = mp_double_set_half_from_subtraction;
1127    math->md_from_of_the_way          = mp_double_set_from_of_the_way;
1128    math->md_from_div                 = mp_double_set_from_div;
1129    math->md_from_mul                 = mp_double_set_from_mul;
1130    math->md_from_int_div             = mp_double_set_from_int_div;
1131    math->md_from_int_mul             = mp_double_set_from_int_mul;
1132    math->md_negate                   = mp_double_negate;
1133    math->md_add                      = mp_double_add;
1134    math->md_subtract                 = mp_double_subtract;
1135    math->md_half                     = mp_double_half;
1136    math->md_do_double                = mp_double_double;
1137    math->md_abs                      = mp_double_abs;
1138    math->md_clone                    = mp_double_clone;
1139    math->md_negated_clone            = mp_double_negated_clone;
1140    math->md_abs_clone                = mp_double_abs_clone;
1141    math->md_swap                     = mp_double_swap;
1142    math->md_add_scaled               = mp_double_add_scaled;
1143    math->md_multiply_int             = mp_double_multiply_int;
1144    math->md_divide_int               = mp_double_divide_int;
1145    math->md_to_int                   = mp_double_to_int;
1146    math->md_to_boolean               = mp_double_to_boolean;
1147    math->md_to_scaled                = mp_double_to_scaled;
1148    math->md_to_double                = mp_double_to_double;
1149    math->md_odd                      = mp_double_odd;
1150    math->md_equal                    = mp_double_equal;
1151    math->md_less                     = mp_double_less;
1152    math->md_greater                  = mp_double_greater;
1153    math->md_non_equal_abs            = mp_double_non_equal_abs;
1154    math->md_round_unscaled           = mp_double_unscaled;
1155    math->md_floor_scaled             = mp_double_floor;
1156    math->md_fraction_to_round_scaled = mp_double_fraction_to_round_scaled;
1157    math->md_make_scaled              = mp_double_number_make_scaled;
1158    math->md_make_fraction            = mp_double_number_make_fraction;
1159    math->md_take_fraction            = mp_double_number_take_fraction;
1160    math->md_take_scaled              = mp_double_number_take_scaled;
1161    math->md_velocity                 = mp_double_velocity;
1162    math->md_n_arg                    = mp_double_n_arg;
1163    math->md_m_log                    = mp_double_m_log;
1164    math->md_m_exp                    = mp_double_m_exp;
1165    math->md_m_unif_rand              = mp_double_m_unif_rand;
1166    math->md_m_norm_rand              = mp_double_m_norm_rand;
1167    math->md_pyth_add                 = mp_double_pyth_add;
1168    math->md_pyth_sub                 = mp_double_pyth_sub;
1169    math->md_power_of                 = mp_double_power_of;
1170    math->md_fraction_to_scaled       = mp_double_fraction_to_scaled;
1171    math->md_scaled_to_fraction       = mp_double_scaled_to_fraction;
1172    math->md_scaled_to_angle          = mp_double_scaled_to_angle;
1173    math->md_angle_to_scaled          = mp_double_angle_to_scaled;
1174    math->md_init_randoms             = mp_double_init_randoms;
1175    math->md_sin_cos                  = mp_double_sin_cos;
1176    math->md_slow_add                 = mp_double_slow_add;
1177    math->md_sqrt                     = mp_double_square_rt;
1178    math->md_print                    = mp_double_print_number;
1179    math->md_tostring                 = mp_double_number_tostring;
1180    math->md_modulo                   = mp_double_modulo;
1181    math->md_ab_vs_cd                 = mp_double_ab_vs_cd;
1182    math->md_crossing_point           = mp_double_crossing_point;
1183    math->md_scan_numeric             = mp_double_scan_numeric_token;
1184    math->md_scan_fractional          = mp_double_scan_fractional_token;
1185    /* housekeeping */
1186    math->md_free_math                = mp_double_free_math;
1187    math->md_set_precision            = mp_double_set_precision;
1188    return math;
1189}
1190
Source Browser ?
