1 /*
2 * mpfr.c - routines for arbitrary-precision number support in gawk.
3 */
4
5 /*
6 * Copyright (C) 2012, 2013, 2015, 2017, 2018, 2019, 2021, 2022,
7 * the Free Software Foundation, Inc.
8 *
9 * This file is part of GAWK, the GNU implementation of the
10 * AWK Programming Language.
11 *
12 * GAWK is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 3 of the License, or
15 * (at your option) any later version.
16 *
17 * GAWK is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with this program; if not, write to the Free Software
24 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
25 */
26
27 #include "awk.h"
28
29 #ifdef HAVE_MPFR
30
31 int MPFR_round_mode = 'N'; // default value
32
33 #if !defined(MPFR_VERSION_MAJOR) || MPFR_VERSION_MAJOR < 3
34 typedef mp_exp_t mpfr_exp_t;
35 #endif
36
37 extern NODE **fmt_list; /* declared in eval.c */
38
39 mpz_t mpzval; /* GMP integer type, used as temporary in few places */
40 mpz_t MNR;
41 mpz_t MFNR;
42 bool do_ieee_fmt; /* IEEE-754 floating-point emulation */
43 mpfr_rnd_t ROUND_MODE;
44
45 static mpfr_prec_t default_prec;
46
47 static mpfr_rnd_t get_rnd_mode(const char rmode);
48 static NODE *mpg_force_number(NODE *n);
49 static NODE *mpg_make_number(double);
50 static NODE *mpg_format_val(const char *format, int index, NODE *s);
51 static int mpg_interpret(INSTRUCTION **cp);
52
53 static mpfr_exp_t min_exp = MPFR_EMIN_DEFAULT;
54 static mpfr_exp_t max_exp = MPFR_EMAX_DEFAULT;
55
56 /* temporary MPFR floats used to hold converted GMP integer operands */
57 static mpfr_t _mpf_t1;
58 static mpfr_t _mpf_t2;
59
60 /*
61 * PRECISION_MIN is the precision used to initialize _mpf_t1 and _mpf_t2.
62 * 64 bits should be enough for exact conversion of most integers to floats.
63 */
64
65 #define PRECISION_MIN 64
66
67 /* mf = { _mpf_t1, _mpf_t2 } */
68 static inline mpfr_ptr mpg_tofloat(mpfr_ptr mf, mpz_ptr mz);
69 /* T = {t1, t2} */
70 #define MP_FLOAT(T) is_mpg_integer(T) ? mpg_tofloat(_mpf_##T, (T)->mpg_i) : (T)->mpg_numbr
71
72
73 /* init_mpfr --- set up MPFR related variables */
74
75 void
76 init_mpfr(mpfr_prec_t prec, const char *rmode)
77 {
78 mpfr_set_default_prec(default_prec = prec);
79 ROUND_MODE = get_rnd_mode(rmode[0]);
80 mpfr_set_default_rounding_mode(ROUND_MODE);
81 make_number = mpg_make_number;
82 str2number = mpg_force_number;
83 format_val = mpg_format_val;
84 cmp_numbers = mpg_cmp;
85
86 mpz_init(MNR);
87 mpz_init(MFNR);
88 do_ieee_fmt = false;
89
90 mpfr_init2(_mpf_t1, PRECISION_MIN);
91 mpfr_init2(_mpf_t2, PRECISION_MIN);
92 mpz_init(mpzval);
93
94 register_exec_hook(mpg_interpret, 0);
95 }
96
97 /* cleanup_mpfr --- clean stuff up, mainly for valgrind */
98
99 void
100 cleanup_mpfr(void)
101 {
102 mpfr_clear(_mpf_t1);
103 mpfr_clear(_mpf_t2);
104 }
105
106 /* mpg_node --- allocate a node to store MPFR float or GMP integer */
107
108 NODE *
109 mpg_node(unsigned int flags)
110 {
111 NODE *r = make_number_node(flags);
112
113 if (flags == MPFN)
114 /* Initialize, set precision to the default precision, and value to NaN */
115 mpfr_init(r->mpg_numbr);
116 else
117 /* Initialize and set value to 0 */
118 mpz_init(r->mpg_i);
119 return r;
120 }
121
122 /*
123 * mpg_make_number --- make a arbitrary-precision number node
124 * and initialize with a C double
125 */
126
127 static NODE *
128 mpg_make_number(double x)
129 {
130 NODE *r;
131 double ival;
132
133 if ((ival = double_to_int(x)) != x) {
134 int tval;
135 r = mpg_float();
136 tval = mpfr_set_d(r->mpg_numbr, x, ROUND_MODE);
137 IEEE_FMT(r->mpg_numbr, tval);
138 } else {
139 r = mpg_integer();
140 mpz_set_d(r->mpg_i, ival);
141 }
142 return r;
143 }
144
145 /* mpg_strtoui --- assign arbitrary-precision integral value from a string */
146
147 int
148 mpg_strtoui(mpz_ptr zi, char *str, size_t len, char **end, int base)
149 {
150 char *s = str;
151 char *start;
152 int ret = -1;
153
154 /*
155 * mpz_set_str does not like leading 0x or 0X for hex (or 0 for octal)
156 * with a non-zero base argument.
157 */
158 if (base == 16 && len >= 2 && *s == '0' && (s[1] == 'x' || s[1] == 'X')) {
159 s += 2; len -= 2;
160 } else if (base == 8 && len >= 1 && *s == '0') {
161 s++; len--;
162 }
163 start = s;
164
165 while (len > 0) {
166 switch (*s) {
167 case '0':
168 case '1':
169 case '2':
170 case '3':
171 case '4':
172 case '5':
173 case '6':
174 case '7':
175 break;
176 case '8':
177 case '9':
178 if (base == 8)
179 base = 10;
180 break;
181 case 'a':
182 case 'b':
183 case 'c':
184 case 'd':
185 case 'e':
186 case 'f':
187 case 'A':
188 case 'B':
189 case 'C':
190 case 'D':
191 case 'E':
192 case 'F':
193 if (base == 16)
194 break;
195 default:
196 goto done;
197 }
198 s++; len--;
199 }
200 done:
201 if (s > start) {
202 char save = *s;
203 *s = '\0';
204 ret = mpz_set_str(zi, start, base);
205 *s = save;
206 }
207 if (end != NULL)
208 *end = s;
209 return ret;
210 }
211
212
213 /* mpg_maybe_float --- test if a string may contain arbitrary-precision float */
214
215 static int
216 mpg_maybe_float(const char *str, int use_locale)
217 {
218 int dec_point = '.';
219 const char *s = str;
220
221 #if defined(HAVE_LOCALE_H)
222 /*
223 * loc.decimal_point may not have been initialized yet,
224 * so double check it before using it.
225 */
226 if (use_locale && loc.decimal_point != NULL && loc.decimal_point[0] != '\0')
227 dec_point = loc.decimal_point[0]; /* XXX --- assumes one char */
228 #endif
229
230 if (strlen(s) >= 3
231 && ( ( (s[0] == 'i' || s[0] == 'I')
232 && (s[1] == 'n' || s[1] == 'N')
233 && (s[2] == 'f' || s[2] == 'F'))
234 || ( (s[0] == 'n' || s[0] == 'N')
235 && (s[1] == 'a' || s[1] == 'A')
236 && (s[2] == 'n' || s[2] == 'N'))))
237 return true;
238
239 for (; *s != '\0'; s++) {
240 if (*s == dec_point || *s == 'e' || *s == 'E')
241 return true;
242 }
243
244 return false;
245 }
246
247
248 /* mpg_zero --- initialize with arbitrary-precision integer(GMP) and set value to zero */
249
250 void
251 mpg_zero(NODE *n)
252 {
253 if (is_mpg_float(n)) {
254 mpfr_clear(n->mpg_numbr);
255 n->flags &= ~MPFN;
256 }
257 if (! is_mpg_integer(n)) {
258 mpz_init(n->mpg_i); /* this also sets its value to 0 */
259 n->flags |= MPZN;
260 } else
261 mpz_set_si(n->mpg_i, 0);
262 }
263
264
265 /* force_mpnum --- force a value to be a GMP integer or MPFR float */
266
267 static int
268 force_mpnum(NODE *n, int do_nondec, int use_locale)
269 {
270 char *cp, *cpend, *ptr, *cp1;
271 char save;
272 int tval, base = 10;
273
274 if (n->stlen == 0 || (n->flags & REGEX) != 0) {
275 mpg_zero(n);
276 return false;
277 }
278
279 cp = n->stptr;
280 cpend = n->stptr + n->stlen;
281 while (cp < cpend && isspace((unsigned char) *cp))
282 cp++;
283 if (cp == cpend) { /* only spaces */
284 mpg_zero(n);
285 return false;
286 }
287
288 save = *cpend;
289 *cpend = '\0';
290
291 if (*cp == '+' || *cp == '-')
292 cp1 = cp + 1;
293 else
294 cp1 = cp;
295
296 /*
297 * Maybe "+" or "-" was the field. mpg_strtoui
298 * won't check for that and set errno, so we have
299 * to check manually.
300 */
301 if (*cp1 == '\0') {
302 *cpend = save;
303 mpg_zero(n);
304 return false;
305 }
306
307 if (do_nondec)
308 base = get_numbase(cp1, cpend - cp1, use_locale);
309
310 if (base != 10 || ! mpg_maybe_float(cp1, use_locale)) {
311 mpg_zero(n);
312 errno = 0;
313 mpg_strtoui(n->mpg_i, cp1, cpend - cp1, & ptr, base);
314 if (*cp == '-')
315 mpz_neg(n->mpg_i, n->mpg_i);
316 goto done;
317 }
318
319 if (is_mpg_integer(n)) {
320 mpz_clear(n->mpg_i);
321 n->flags &= ~MPZN;
322 }
323
324 if (! is_mpg_float(n)) {
325 mpfr_init(n->mpg_numbr);
326 n->flags |= MPFN;
327 }
328
329 errno = 0;
330 tval = mpfr_strtofr(n->mpg_numbr, cp, & ptr, base, ROUND_MODE);
331 if (mpfr_nan_p(n->mpg_numbr) && *cp == '-')
332 tval = mpfr_setsign(n->mpg_numbr, n->mpg_numbr, 1, ROUND_MODE);
333 IEEE_FMT(n->mpg_numbr, tval);
334 done:
335 /* trailing space is OK for NUMBER */
336 while (ptr < cpend && isspace((unsigned char) *ptr))
337 ptr++;
338 *cpend = save;
339 if (errno == 0 && ptr == cpend)
340 return true;
341 errno = 0;
342 return false;
343 }
344
345 /* mpg_force_number --- force a value to be a multiple-precision number */
346
347 static NODE *
348 mpg_force_number(NODE *n)
349 {
350 char *cp, *cpend;
351
352 if (n->type == Node_elem_new) {
353 n->type = Node_val;
354 n->flags &= ~STRING;
355 n->stptr[0] = '0'; // STRCUR is still set
356 n->stlen = 1;
357
358 return n;
359 }
360
361 if ((n->flags & NUMCUR) != 0)
362 return n;
363 n->flags |= NUMCUR;
364
365 /* Trim leading white space, bailing out if there's nothing else */
366 for (cp = n->stptr, cpend = cp + n->stlen;
367 cp < cpend && isspace((unsigned char) *cp); cp++)
368 continue;
369
370 if (cp == cpend)
371 goto badnum;
372
373 /* At this point, we know the string is not entirely white space */
374 /* Trim trailing white space */
375 while (isspace((unsigned char) cpend[-1]))
376 cpend--;
377
378 /*
379 * 2/2007:
380 * POSIX, by way of severe language lawyering, seems to
381 * allow things like "inf" and "nan" to mean something.
382 * So if do_posix, the user gets what he deserves.
383 * This also allows hexadecimal floating point. Ugh.
384 */
385 if (! do_posix) {
386 if (is_alpha((unsigned char) *cp))
387 goto badnum;
388 else if (is_ieee_magic_val(cp)) {
389 if (cpend != cp + 4)
390 goto badnum;
391 /* else
392 fall through */
393 }
394 /* else
395 fall through */
396 }
397 /* else POSIX, so
398 fall through */
399
400 if (force_mpnum(n, (do_non_decimal_data && ! do_traditional), true)) {
401 if ((n->flags & USER_INPUT) != 0) {
402 /* leave USER_INPUT set to indicate a strnum */
403 n->flags &= ~STRING;
404 n->flags |= NUMBER;
405 }
406 } else
407 n->flags &= ~USER_INPUT;
408 return n;
409 badnum:
410 mpg_zero(n);
411 n->flags &= ~USER_INPUT;
412 return n;
413 }
414
415 /* mpg_format_val --- format a numeric value based on format */
416
417 static NODE *
418 mpg_format_val(const char *format, int index, NODE *s)
419 {
420 NODE *dummy[2], *r;
421 unsigned int oflags;
422
423 if (out_of_range(s)) {
424 const char *result = format_nan_inf(s, 'g');
425 return make_string(result, strlen(result));
426 }
427
428 /* create dummy node for a sole use of format_tree */
429 dummy[1] = s;
430 oflags = s->flags;
431
432 if (is_mpg_integer(s) || mpfr_integer_p(s->mpg_numbr)) {
433 /* integral value, use %d */
434 r = format_tree("%d", 2, dummy, 2);
435 s->stfmt = STFMT_UNUSED;
436 } else {
437 r = format_tree(format, fmt_list[index]->stlen, dummy, 2);
438 assert(r != NULL);
439 s->stfmt = index;
440 }
441 s->flags = oflags;
442 s->stlen = r->stlen;
443 if ((s->flags & (MALLOC|STRCUR)) == (MALLOC|STRCUR))
444 efree(s->stptr);
445 s->stptr = r->stptr;
446 s->flags |= STRCUR;
447 s->strndmode = MPFR_round_mode;
448 freenode(r); /* Do not unref(r)! We want to keep s->stptr == r->stpr. */
449 free_wstr(s);
450 return s;
451 }
452
453 /* mpg_cmp --- compare two numbers */
454
455 int
456 mpg_cmp(const NODE *t1, const NODE *t2)
457 {
458 /*
459 * For the purposes of sorting, NaN is considered greater than
460 * any other value, and all NaN values are considered equivalent and equal.
461 */
462
463 if (is_mpg_float(t1)) {
464 if (is_mpg_float(t2)) {
465 if (mpfr_nan_p(t1->mpg_numbr))
466 return ! mpfr_nan_p(t2->mpg_numbr);
467 if (mpfr_nan_p(t2->mpg_numbr))
468 return -1;
469 return mpfr_cmp(t1->mpg_numbr, t2->mpg_numbr);
470 }
471 if (mpfr_nan_p(t1->mpg_numbr))
472 return 1;
473 return mpfr_cmp_z(t1->mpg_numbr, t2->mpg_i);
474 } else if (is_mpg_float(t2)) {
475 int ret;
476 if (mpfr_nan_p(t2->mpg_numbr))
477 return -1;
478 ret = mpfr_cmp_z(t2->mpg_numbr, t1->mpg_i);
479 return ret > 0 ? -1 : (ret < 0);
480 } else if (is_mpg_integer(t1)) {
481 return mpz_cmp(t1->mpg_i, t2->mpg_i);
482 }
483
484 /* t1 and t2 are AWKNUMs */
485 return cmp_awknums(t1, t2);
486 }
487
488 /* mpg_cmp_as_numbers --- compare two numbers, similar to doubles */
489
490 bool
491 mpg_cmp_as_numbers(const NODE *t1, const NODE *t2, scalar_cmp_t comparison_type)
492 {
493 /*
494 * This routine provides numeric comparisons that should work
495 * the same as in C. It should NOT be used for sorting.
496 */
497
498 bool t1_nan = mpfr_nan_p(t1->mpg_numbr);
499 bool t2_nan = mpfr_nan_p(t2->mpg_numbr);
500 bool ret = false;
501
502 // MPFR is different than native doubles...
503 if (t1_nan || t2_nan)
504 return comparison_type == SCALAR_NEQ;
505
506 int di = mpg_cmp(t1, t2);
507
508 switch (comparison_type) {
509 case SCALAR_EQ:
510 ret = (di == 0);
511 break;
512 case SCALAR_NEQ:
513 ret = (di != 0);
514 break;
515 case SCALAR_LT:
516 ret = (di < 0);
517 break;
518 case SCALAR_LE:
519 ret = (di <= 0);
520 break;
521 case SCALAR_GT:
522 ret = (di > 0);
523 break;
524 case SCALAR_GE:
525 ret = (di >= 0);
526 break;
527 default:
528 cant_happen("invalid comparison type %d", (int) comparison_type);
529 break;
530 }
531
532 return ret;
533 }
534
535
536 /*
537 * mpg_update_var --- update NR or FNR.
538 * NR_node->var_value(mpz_t) = MNR(mpz_t) * LONG_MAX + NR(long)
539 */
540
541 NODE *
542 mpg_update_var(NODE *n)
543 {
544 NODE *val = n->var_value;
545 long nr = 0;
546 mpz_ptr nq = 0;
547
548 if (n == NR_node) {
549 nr = NR;
550 nq = MNR;
551 } else if (n == FNR_node) {
552 nr = FNR;
553 nq = MFNR;
554 } else
555 cant_happen("invalid node for mpg_update_var%s", "");
556
557 if (mpz_sgn(nq) == 0) {
558 /* Efficiency hack similar to that for AWKNUM */
559 if (is_mpg_float(val) || mpz_get_si(val->mpg_i) != nr) {
560 unref(n->var_value);
561 val = n->var_value = mpg_integer();
562 mpz_set_si(val->mpg_i, nr);
563 }
564 } else {
565 unref(n->var_value);
566 val = n->var_value = mpg_integer();
567 mpz_set_si(val->mpg_i, nr);
568 mpz_addmul_ui(val->mpg_i, nq, LONG_MAX); /* val->mpg_i += nq * LONG_MAX */
569 }
570 return val;
571 }
572
573 /* mpg_set_var --- set NR or FNR */
574
575 long
576 mpg_set_var(NODE *n)
577 {
578 long nr = 0;
579 mpz_ptr nq = 0, r;
580 NODE *val = n->var_value;
581
582 if (n == NR_node)
583 nq = MNR;
584 else if (n == FNR_node)
585 nq = MFNR;
586 else
587 cant_happen("invalid node for mpg_set_var%s", "");
588
589 if (is_mpg_integer(val))
590 r = val->mpg_i;
591 else {
592 /* convert float to integer */
593 mpfr_get_z(mpzval, val->mpg_numbr, MPFR_RNDZ);
594 r = mpzval;
595 }
596 nr = mpz_fdiv_q_ui(nq, r, LONG_MAX); /* nq (MNR or MFNR) is quotient */
597 return nr; /* remainder (NR or FNR) */
598 }
599
600 /* set_PREC --- update MPFR PRECISION related variables when PREC assigned to */
601
602 void
603 set_PREC()
604 {
605 long prec = 0;
606 NODE *val;
607 static const struct ieee_fmt {
608 const char *name;
609 mpfr_prec_t precision;
610 mpfr_exp_t emax;
611 mpfr_exp_t emin;
612 } ieee_fmts[] = {
613 { "half", 11, 16, -23 }, /* binary16 */
614 { "single", 24, 128, -148 }, /* binary32 */
615 { "double", 53, 1024, -1073 }, /* binary64 */
616 { "quad", 113, 16384, -16493 }, /* binary128 */
617 { "oct", 237, 262144, -262377 }, /* binary256, not in the IEEE 754-2008 standard */
618
619 /*
620 * For any bitwidth = 32 * k ( k >= 4),
621 * precision = 13 + bitwidth - int(4 * log2(bitwidth))
622 * emax = 1 << bitwidth - precision - 1
623 * emin = 4 - emax - precision
624 */
625 };
626
627 if (! do_mpfr)
628 return;
629
630 val = fixtype(PREC_node->var_value);
631
632 if ((val->flags & STRING) != 0) {
633 int i, j;
634
635 /* emulate IEEE-754 binary format */
636
637 for (i = 0, j = sizeof(ieee_fmts)/sizeof(ieee_fmts[0]); i < j; i++) {
638 if (strcasecmp(ieee_fmts[i].name, val->stptr) == 0)
639 break;
640 }
641
642 if (i < j) {
643 prec = ieee_fmts[i].precision;
644
645 /*
646 * We *DO NOT* change the MPFR exponent range using
647 * mpfr_set_{emin, emax} here. See format_ieee() for details.
648 */
649 max_exp = ieee_fmts[i].emax;
650 min_exp = ieee_fmts[i].emin;
651
652 do_ieee_fmt = true;
653 }
654 }
655
656 if (prec <= 0) {
657 force_number(val);
658 prec = get_number_si(val);
659 if (prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX) {
660 force_string(val);
661 warning(_("PREC value `%.*s' is invalid"), (int) val->stlen, val->stptr);
662 prec = 0;
663 } else
664 do_ieee_fmt = false;
665 }
666
667 if (prec > 0)
668 mpfr_set_default_prec(default_prec = prec);
669 }
670
671
672 /* get_rnd_mode --- convert string to MPFR rounding mode */
673
674 static mpfr_rnd_t
675 get_rnd_mode(const char rmode)
676 {
677 switch (rmode) {
678 case 'N':
679 case 'n':
680 return MPFR_RNDN; /* round to nearest (IEEE-754 roundTiesToEven) */
681 case 'Z':
682 case 'z':
683 return MPFR_RNDZ; /* round toward zero (IEEE-754 roundTowardZero) */
684 case 'U':
685 case 'u':
686 return MPFR_RNDU; /* round toward plus infinity (IEEE-754 roundTowardPositive) */
687 case 'D':
688 case 'd':
689 return MPFR_RNDD; /* round toward minus infinity (IEEE-754 roundTowardNegative) */
690 #if defined(MPFR_VERSION_MAJOR) && MPFR_VERSION_MAJOR > 2
691 case 'A':
692 case 'a':
693 return MPFR_RNDA; /* round away from zero */
694 #endif
695 default:
696 break;
697 }
698 return (mpfr_rnd_t) -1;
699 }
700
701 /*
702 * set_ROUNDMODE --- update MPFR rounding mode related variables
703 * when ROUNDMODE assigned to
704 */
705
706 void
707 set_ROUNDMODE()
708 {
709 if (do_mpfr) {
710 mpfr_rnd_t rndm = (mpfr_rnd_t) -1;
711 NODE *n;
712 n = force_string(ROUNDMODE_node->var_value);
713 if (n->stlen == 1)
714 rndm = get_rnd_mode(n->stptr[0]);
715 if (rndm != -1) {
716 mpfr_set_default_rounding_mode(rndm);
717 ROUND_MODE = rndm;
718 MPFR_round_mode = n->stptr[0];
719 } else
720 warning(_("ROUNDMODE value `%.*s' is invalid"), (int) n->stlen, n->stptr);
721 }
722 }
723
724
725 /* format_ieee --- make sure a number follows IEEE-754 floating-point standard */
726
727 int
728 format_ieee(mpfr_ptr x, int tval)
729 {
730 /*
731 * The MPFR doc says that it's our responsibility to make sure all numbers
732 * including those previously created are in range after we've changed the
733 * exponent range. Most MPFR operations and functions require
734 * the input arguments to have exponents within the current exponent range.
735 * Any argument outside the range results in a MPFR assertion failure
736 * like this:
737 *
738 * $ gawk -M 'BEGIN { x=1.0e-10000; print x+0; PREC="double"; print x+0}'
739 * 1e-10000
740 * init2.c:52: MPFR assertion failed ....
741 *
742 * A "naive" approach would be to keep track of the ternary state and
743 * the rounding mode for each number, and make sure it is in the current
744 * exponent range (using mpfr_check_range) before using it in an
745 * operation or function. Instead, we adopt the following strategy.
746 *
747 * When gawk starts, the exponent range is the MPFR default
748 * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT]. Any number that gawk
749 * creates must have exponent in this range (excluding infinities, NaNs and zeros).
750 * Each MPFR operation or function is performed with this default exponent
751 * range.
752 *
753 * When emulating IEEE-754 format, the exponents are *temporarily* changed,
754 * mpfr_check_range is called to make sure the number is in the new range,
755 * and mpfr_subnormalize is used to round following the rules of subnormal
756 * arithmetic. The exponent range is then *restored* to the original value
757 * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT].
758 */
759
760 (void) mpfr_set_emin(min_exp);
761 (void) mpfr_set_emax(max_exp);
762 tval = mpfr_check_range(x, tval, ROUND_MODE);
763 tval = mpfr_subnormalize(x, tval, ROUND_MODE);
764 (void) mpfr_set_emin(MPFR_EMIN_DEFAULT);
765 (void) mpfr_set_emax(MPFR_EMAX_DEFAULT);
766 return tval;
767 }
768
769
770 /* do_mpfr_atan2 --- do the atan2 function */
771
772 NODE *
773 do_mpfr_atan2(int nargs)
774 {
775 NODE *t1, *t2, *res;
776 mpfr_ptr p1, p2;
777 int tval;
778
779 check_exact_args(nargs, "atan2", 2);
780
781 t2 = POP_SCALAR();
782 t1 = POP_SCALAR();
783
784 if (do_lint) {
785 if ((fixtype(t1)->flags & NUMBER) == 0)
786 lintwarn(_("atan2: received non-numeric first argument"));
787 if ((fixtype(t2)->flags & NUMBER) == 0)
788 lintwarn(_("atan2: received non-numeric second argument"));
789 }
790 force_number(t1);
791 force_number(t2);
792
793 p1 = MP_FLOAT(t1);
794 p2 = MP_FLOAT(t2);
795 res = mpg_float();
796 /* See MPFR documentation for handling of special values like +inf as an argument */
797 tval = mpfr_atan2(res->mpg_numbr, p1, p2, ROUND_MODE);
798 IEEE_FMT(res->mpg_numbr, tval);
799
800 DEREF(t1);
801 DEREF(t2);
802 return res;
803 }
804
805 /* do_mpfr_func --- run an MPFR function - not inline, for debugging */
806
807 static inline NODE *
808 do_mpfr_func(const char *name,
809 int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t),
810 int nargs, bool warn_negative)
811 {
812 NODE *t1, *res;
813 mpfr_ptr p1;
814 int tval;
815 mpfr_prec_t argprec;
816
817 check_exact_args(nargs, name, 1);
818
819 t1 = POP_SCALAR();
820 if (do_lint && (fixtype(t1)->flags & NUMBER) == 0)
821 lintwarn(_("%s: received non-numeric argument"), name);
822
823 force_number(t1);
824 p1 = MP_FLOAT(t1);
825 if (warn_negative && mpfr_sgn(p1) < 0) {
826 force_string(t1);
827 warning(_("%s: received negative argument %.*s"), name, (int) t1->stlen, t1->stptr);
828 }
829 res = mpg_float();
830 if ((argprec = mpfr_get_prec(p1)) > default_prec)
831 mpfr_set_prec(res->mpg_numbr, argprec); /* needed at least for sqrt() */
832 tval = mpfr_func(res->mpg_numbr, p1, ROUND_MODE);
833 IEEE_FMT(res->mpg_numbr, tval);
834 DEREF(t1);
835 return res;
836 }
837
838 #define SPEC_MATH(X, WN) \
839 NODE *result; \
840 result = do_mpfr_func(#X, mpfr_##X, nargs, WN); \
841 return result
842
843 /* do_mpfr_sin --- do the sin function */
844
845 NODE *
846 do_mpfr_sin(int nargs)
847 {
848 SPEC_MATH(sin, false);
849 }
850
851 /* do_mpfr_cos --- do the cos function */
852
853 NODE *
854 do_mpfr_cos(int nargs)
855 {
856 SPEC_MATH(cos, false);
857 }
858
859 /* do_mpfr_exp --- exponential function */
860
861 NODE *
862 do_mpfr_exp(int nargs)
863 {
864 SPEC_MATH(exp, false);
865 }
866
867 /* do_mpfr_log --- the log function */
868
869 NODE *
870 do_mpfr_log(int nargs)
871 {
872 SPEC_MATH(log, true);
873 }
874
875 /* do_mpfr_sqrt --- do the sqrt function */
876
877 NODE *
878 do_mpfr_sqrt(int nargs)
879 {
880 SPEC_MATH(sqrt, true);
881 }
882
883 /* do_mpfr_int --- convert double to int for awk */
884
885 NODE *
886 do_mpfr_int(int nargs)
887 {
888 NODE *tmp, *r;
889
890 check_exact_args(nargs, "int", 1);
891
892 tmp = POP_SCALAR();
893 if (do_lint && (fixtype(tmp)->flags & NUMBER) == 0)
894 lintwarn(_("int: received non-numeric argument"));
895 force_number(tmp);
896
897 if (is_mpg_integer(tmp)) {
898 r = mpg_integer();
899 mpz_set(r->mpg_i, tmp->mpg_i);
900 } else {
901 if (! mpfr_number_p(tmp->mpg_numbr)) {
902 /* [+-]inf or NaN */
903 return tmp;
904 }
905
906 r = mpg_integer();
907 mpfr_get_z(r->mpg_i, tmp->mpg_numbr, MPFR_RNDZ);
908 }
909
910 DEREF(tmp);
911 return r;
912 }
913
914 /* do_mpfr_compl --- perform a ~ operation */
915
916 NODE *
917 do_mpfr_compl(int nargs)
918 {
919 NODE *tmp, *r;
920 mpz_ptr zptr;
921
922 check_exact_args(nargs, "compl", 1);
923
924 tmp = POP_SCALAR();
925 if (do_lint && (fixtype(tmp)->flags & NUMBER) == 0)
926 lintwarn(_("compl: received non-numeric argument"));
927
928 force_number(tmp);
929 if (is_mpg_float(tmp)) {
930 mpfr_ptr p = tmp->mpg_numbr;
931
932 if (! mpfr_number_p(p)) {
933 /* [+-]inf or NaN */
934 return tmp;
935 }
936 if (mpfr_sgn(p) < 0)
937 fatal("%s",
938 mpg_fmt(_("compl(%Rg): negative value is not allowed"), p)
939 );
940 if (do_lint) {
941 if (! mpfr_integer_p(p))
942 lintwarn("%s",
943 mpg_fmt(_("comp(%Rg): fractional value will be truncated"), p)
944 );
945 }
946
947 mpfr_get_z(mpzval, p, MPFR_RNDZ); /* float to integer conversion */
948 zptr = mpzval;
949 } else {
950 /* (tmp->flags & MPZN) != 0 */
951 zptr = tmp->mpg_i;
952 if (mpz_sgn(zptr) < 0)
953 fatal("%s",
954 mpg_fmt(_("compl(%Zd): negative values are not allowed"), zptr)
955 );
956 }
957
958 r = mpg_integer();
959 mpz_com(r->mpg_i, zptr);
960 DEREF(tmp);
961 return r;
962 }
963
964 /* get_intval --- get the (converted) integral operand of a binary function. */
965
966 static mpz_ptr
967 get_intval(NODE *t1, int argnum, const char *op)
968 {
969 mpz_ptr pz;
970
971 if (do_lint && (fixtype(t1)->flags & NUMBER) == 0)
972 lintwarn(_("%s: received non-numeric argument #%d"), op, argnum);
973
974 (void) force_number(t1);
975
976 if (is_mpg_float(t1)) {
977 mpfr_ptr left = t1->mpg_numbr;
978 if (! mpfr_number_p(left)) {
979 /* inf or NaN */
980 if (do_lint)
981 lintwarn("%s",
982 mpg_fmt(_("%s: argument #%d has invalid value %Rg, using 0"),
983 op, argnum, left)
984 );
985
986 emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
987 mpz_init(pz);
988 return pz; /* should be freed */
989 }
990
991 if (mpfr_sgn(left) < 0)
992 fatal("%s",
993 mpg_fmt(_("%s: argument #%d negative value %Rg is not allowed"),
994 op, argnum, left)
995 );
996
997 if (do_lint) {
998 if (! mpfr_integer_p(left))
999 lintwarn("%s",
1000 mpg_fmt(_("%s: argument #%d fractional value %Rg will be truncated"),
1001 op, argnum, left)
1002 );
1003 }
1004
1005 emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
1006 mpz_init(pz);
1007 mpfr_get_z(pz, left, MPFR_RNDZ); /* float to integer conversion */
1008 return pz; /* should be freed */
1009 }
1010 /* (t1->flags & MPZN) != 0 */
1011 pz = t1->mpg_i;
1012 if (mpz_sgn(pz) < 0)
1013 fatal("%s",
1014 mpg_fmt(_("%s: argument #%d negative value %Zd is not allowed"),
1015 op, argnum, pz)
1016 );
1017
1018 return pz; /* must not be freed */
1019 }
1020
1021
1022 /* free_intval --- free the converted integer value returned by get_intval() */
1023
1024 static inline void
1025 free_intval(NODE *t, mpz_ptr pz)
1026 {
1027 if ((t->flags & MPZN) == 0) {
1028 mpz_clear(pz);
1029 efree(pz);
1030 }
1031 }
1032
1033
1034 /* do_mpfr_lshift --- perform a << operation */
1035
1036 NODE *
1037 do_mpfr_lshift(int nargs)
1038 {
1039 NODE *t1, *t2, *res;
1040 unsigned long shift;
1041 mpz_ptr pz1, pz2;
1042
1043 check_exact_args(nargs, "lshift", 2);
1044
1045 t2 = POP_SCALAR();
1046 t1 = POP_SCALAR();
1047
1048 pz1 = get_intval(t1, 1, "lshift");
1049 pz2 = get_intval(t2, 2, "lshift");
1050
1051 /*
1052 * mpz_get_ui: If op is too big to fit an unsigned long then just
1053 * the least significant bits that do fit are returned.
1054 * The sign of op is ignored, only the absolute value is used.
1055 */
1056
1057 shift = mpz_get_ui(pz2); /* GMP integer => unsigned long conversion */
1058 res = mpg_integer();
1059 mpz_mul_2exp(res->mpg_i, pz1, shift); /* res = pz1 * 2^shift */
1060
1061 free_intval(t1, pz1);
1062 free_intval(t2, pz2);
1063 DEREF(t2);
1064 DEREF(t1);
1065 return res;
1066 }
1067
1068 /* do_mpfr_rshift --- perform a >> operation */
1069
1070 NODE *
1071 do_mpfr_rshift(int nargs)
1072 {
1073 NODE *t1, *t2, *res;
1074 unsigned long shift;
1075 mpz_ptr pz1, pz2;
1076
1077 check_exact_args(nargs, "rshift", 2);
1078
1079 t2 = POP_SCALAR();
1080 t1 = POP_SCALAR();
1081
1082 pz1 = get_intval(t1, 1, "rshift");
1083 pz2 = get_intval(t2, 2, "rshift");
1084
1085 /* N.B: See do_mpfp_lshift. */
1086 shift = mpz_get_ui(pz2); /* GMP integer => unsigned long conversion */
1087 res = mpg_integer();
1088 mpz_fdiv_q_2exp(res->mpg_i, pz1, shift); /* res = pz1 / 2^shift, round towards -inf */
1089
1090 free_intval(t1, pz1);
1091 free_intval(t2, pz2);
1092 DEREF(t2);
1093 DEREF(t1);
1094 return res;
1095 }
1096
1097
1098 /* do_mpfr_and --- perform an & operation */
1099
1100 NODE *
1101 do_mpfr_and(int nargs)
1102 {
1103 NODE *t1, *t2, *res;
1104 mpz_ptr pz1, pz2;
1105 int i;
1106
1107 if (nargs < 2)
1108 fatal(_("and: called with less than two arguments"));
1109
1110 t2 = POP_SCALAR();
1111 pz2 = get_intval(t2, nargs, "and");
1112
1113 res = mpg_integer();
1114 for (i = 1; i < nargs; i++) {
1115 t1 = POP_SCALAR();
1116 pz1 = get_intval(t1, nargs - i, "and");
1117 mpz_and(res->mpg_i, pz1, pz2);
1118 free_intval(t1, pz1);
1119 DEREF(t1);
1120 if (i == 1) {
1121 free_intval(t2, pz2);
1122 DEREF(t2);
1123 }
1124 pz2 = res->mpg_i;
1125 }
1126 return res;
1127 }
1128
1129
1130 /* do_mpfr_or --- perform an | operation */
1131
1132 NODE *
1133 do_mpfr_or(int nargs)
1134 {
1135 NODE *t1, *t2, *res;
1136 mpz_ptr pz1, pz2;
1137 int i;
1138
1139 if (nargs < 2)
1140 fatal(_("or: called with less than two arguments"));
1141
1142 t2 = POP_SCALAR();
1143 pz2 = get_intval(t2, nargs, "or");
1144
1145 res = mpg_integer();
1146 for (i = 1; i < nargs; i++) {
1147 t1 = POP_SCALAR();
1148 pz1 = get_intval(t1, nargs - i, "or");
1149 mpz_ior(res->mpg_i, pz1, pz2);
1150 free_intval(t1, pz1);
1151 DEREF(t1);
1152 if (i == 1) {
1153 free_intval(t2, pz2);
1154 DEREF(t2);
1155 }
1156 pz2 = res->mpg_i;
1157 }
1158 return res;
1159 }
1160
1161 /* do_mpfr_xor --- perform an ^ operation */
1162
1163 NODE *
1164 do_mpfr_xor(int nargs)
1165 {
1166 NODE *t1, *t2, *res;
1167 mpz_ptr pz1, pz2;
1168 int i;
1169
1170 if (nargs < 2)
1171 fatal(_("xor: called with less than two arguments"));
1172
1173 t2 = POP_SCALAR();
1174 pz2 = get_intval(t2, nargs, "xor");
1175
1176 res = mpg_integer();
1177 for (i = 1; i < nargs; i++) {
1178 t1 = POP_SCALAR();
1179 pz1 = get_intval(t1, nargs - i, "xor");
1180 mpz_xor(res->mpg_i, pz1, pz2);
1181 free_intval(t1, pz1);
1182 DEREF(t1);
1183 if (i == 1) {
1184 free_intval(t2, pz2);
1185 DEREF(t2);
1186 }
1187 pz2 = res->mpg_i;
1188 }
1189 return res;
1190 }
1191
1192 /* do_mpfr_strtonum --- the strtonum function */
1193
1194 NODE *
1195 do_mpfr_strtonum(int nargs)
1196 {
1197 NODE *tmp, *r;
1198
1199 check_exact_args(nargs, "strtonum", 1);
1200
1201 tmp = fixtype(POP_SCALAR());
1202 if ((tmp->flags & NUMBER) == 0) {
1203 r = mpg_integer(); /* will be changed to MPFR float if necessary in force_mpnum() */
1204 r->stptr = tmp->stptr;
1205 r->stlen = tmp->stlen;
1206 force_mpnum(r, true, use_lc_numeric);
1207 r->stptr = NULL;
1208 r->stlen = 0;
1209 r->wstptr = NULL;
1210 r->wstlen = 0;
1211 } else if (is_mpg_float(tmp)) {
1212 int tval;
1213 r = mpg_float();
1214 tval = mpfr_set(r->mpg_numbr, tmp->mpg_numbr, ROUND_MODE);
1215 IEEE_FMT(r->mpg_numbr, tval);
1216 } else {
1217 r = mpg_integer();
1218 mpz_set(r->mpg_i, tmp->mpg_i);
1219 }
1220
1221 DEREF(tmp);
1222 return r;
1223 }
1224
1225
1226 static bool firstrand = true;
1227 static gmp_randstate_t state;
1228 static mpz_t seed; /* current seed */
1229
1230 /* do_mpfr_rand --- do the rand function */
1231
1232 NODE *
1233 do_mpfr_rand(int nargs ATTRIBUTE_UNUSED)
1234 {
1235 NODE *res;
1236 int tval;
1237
1238 check_exact_args(nargs, "rand", 0);
1239
1240 if (firstrand) {
1241 #if 0
1242 /* Choose the default algorithm */
1243 gmp_randinit_default(state);
1244 #endif
1245 /*
1246 * Choose a specific (Mersenne Twister) algorithm in case the default
1247 * changes in the future.
1248 */
1249
1250 gmp_randinit_mt(state);
1251
1252 mpz_init(seed);
1253 mpz_set_ui(seed, 1);
1254 /* seed state */
1255 gmp_randseed(state, seed);
1256 firstrand = false;
1257 }
1258 res = mpg_float();
1259 tval = mpfr_urandomb(res->mpg_numbr, state);
1260 IEEE_FMT(res->mpg_numbr, tval);
1261 return res;
1262 }
1263
1264
1265 /* do_mpfr_srand --- seed the random number generator */
1266
1267 NODE *
1268 do_mpfr_srand(int nargs)
1269 {
1270 NODE *res;
1271
1272 if (firstrand) {
1273 #if 0
1274 /* Choose the default algorithm */
1275 gmp_randinit_default(state);
1276 #endif
1277 /*
1278 * Choose a specific algorithm (Mersenne Twister) in case default
1279 * changes in the future.
1280 */
1281
1282 gmp_randinit_mt(state);
1283
1284 mpz_init(seed);
1285 mpz_set_ui(seed, 1);
1286 /* No need to seed state, will change it below */
1287 firstrand = false;
1288 }
1289
1290 check_args_min_max(nargs, "srand", 0, 1);
1291
1292 res = mpg_integer();
1293 mpz_set(res->mpg_i, seed); /* previous seed */
1294
1295 if (nargs == 0)
1296 mpz_set_ui(seed, (unsigned long) time((time_t *) 0));
1297 else {
1298 NODE *tmp;
1299 tmp = POP_SCALAR();
1300 if (do_lint && (fixtype(tmp)->flags & NUMBER) == 0)
1301 lintwarn(_("srand: received non-numeric argument"));
1302 force_number(tmp);
1303 if (is_mpg_float(tmp))
1304 mpfr_get_z(seed, tmp->mpg_numbr, MPFR_RNDZ);
1305 else /* MP integer */
1306 mpz_set(seed, tmp->mpg_i);
1307 DEREF(tmp);
1308 }
1309
1310 gmp_randseed(state, seed);
1311 return res;
1312 }
1313
1314 #ifdef SUPPLY_INTDIV
1315 /* do_mpfr_intdiv --- do integer division, return quotient and remainder in dest array */
1316
1317 /*
1318 * We define the semantics as:
1319 * numerator = int(numerator)
1320 * denominator = int(denonmator)
1321 * quotient = int(numerator / denomator)
1322 * remainder = int(numerator % denomator)
1323 */
1324
1325 NODE *
1326 do_mpfr_intdiv(int nargs)
1327 {
1328 NODE *numerator, *denominator, *result;
1329 NODE *num, *denom;
1330 NODE *quotient, *remainder;
1331 NODE *sub, **lhs;
1332
1333 check_exact_args(nargs, "intdiv", 3);
1334
1335 result = POP_PARAM();
1336 if (result->type != Node_var_array)
1337 fatal(_("intdiv: third argument is not an array"));
1338 assoc_clear(result);
1339
1340 denominator = POP_SCALAR();
1341 numerator = POP_SCALAR();
1342
1343 if (do_lint) {
1344 if ((fixtype(numerator)->flags & NUMBER) == 0)
1345 lintwarn(_("intdiv: received non-numeric first argument"));
1346 if ((fixtype(denominator)->flags & NUMBER) == 0)
1347 lintwarn(_("intdiv: received non-numeric second argument"));
1348 }
1349
1350 (void) force_number(numerator);
1351 (void) force_number(denominator);
1352
1353 /* convert numerator and denominator to integer */
1354 if (is_mpg_integer(numerator)) {
1355 num = mpg_integer();
1356 mpz_set(num->mpg_i, numerator->mpg_i);
1357 } else {
1358 if (! mpfr_number_p(numerator->mpg_numbr)) {
1359 /* [+-]inf or NaN */
1360 unref(numerator);
1361 unref(denominator);
1362 return make_number((AWKNUM) -1);
1363 }
1364
1365 num = mpg_integer();
1366 mpfr_get_z(num->mpg_i, numerator->mpg_numbr, MPFR_RNDZ);
1367 }
1368
1369 if (is_mpg_integer(denominator)) {
1370 denom = mpg_integer();
1371 mpz_set(denom->mpg_i, denominator->mpg_i);
1372 } else {
1373 if (! mpfr_number_p(denominator->mpg_numbr)) {
1374 /* [+-]inf or NaN */
1375 unref(numerator);
1376 unref(denominator);
1377 unref(num);
1378 return make_number((AWKNUM) -1);
1379 }
1380
1381 denom = mpg_integer();
1382 mpfr_get_z(denom->mpg_i, denominator->mpg_numbr, MPFR_RNDZ);
1383 }
1384
1385 if (mpz_sgn(denom->mpg_i) == 0)
1386 fatal(_("intdiv: division by zero attempted"));
1387
1388 quotient = mpg_integer();
1389 remainder = mpg_integer();
1390
1391 /* do the division */
1392 mpz_tdiv_qr(quotient->mpg_i, remainder->mpg_i, num->mpg_i, denom->mpg_i);
1393 unref(num);
1394 unref(denom);
1395 unref(numerator);
1396 unref(denominator);
1397
1398 sub = make_string("quotient", 8);
1399 assoc_set(result, sub, quotient);
1400
1401 sub = make_string("remainder", 9);
1402 assoc_set(result, sub, remainder);
1403
1404 return make_number((AWKNUM) 0.0);
1405 }
1406 #endif /* SUPPLY_INTDIV */
1407
1408 /*
1409 * mpg_tofloat --- convert an arbitrary-precision integer operand to
1410 * a float without loss of precision. It is assumed that the
1411 * MPFR variable has already been initialized.
1412 */
1413
1414 static inline mpfr_ptr
1415 mpg_tofloat(mpfr_ptr mf, mpz_ptr mz)
1416 {
1417 size_t prec;
1418
1419 /*
1420 * When implicitely converting a GMP integer operand to a MPFR float, use
1421 * a precision sufficiently large to hold the converted value exactly.
1422 *
1423 * $ ./gawk -M 'BEGIN { print 13 % 2 }'
1424 * 1
1425 * If the user-specified precision is used to convert the integer 13 to a
1426 * float, one will get:
1427 * $ ./gawk -M 'BEGIN { PREC=2; print 13 % 2.0 }'
1428 * 0
1429 */
1430
1431 prec = mpz_sizeinbase(mz, 2); /* most significant 1 bit position starting at 1 */
1432 if (prec > PRECISION_MIN) {
1433 prec -= (size_t) mpz_scan1(mz, 0); /* least significant 1 bit index starting at 0 */
1434 if (prec > MPFR_PREC_MAX)
1435 prec = MPFR_PREC_MAX;
1436 else if (prec < PRECISION_MIN)
1437 prec = PRECISION_MIN;
1438 }
1439 else
1440 prec = PRECISION_MIN;
1441 /*
1442 * Always set the precision to avoid hysteresis, since do_mpfr_func
1443 * may copy our precision.
1444 */
1445 if (prec != mpfr_get_prec(mf))
1446 mpfr_set_prec(mf, prec);
1447
1448 mpfr_set_z(mf, mz, ROUND_MODE);
1449 return mf;
1450 }
1451
1452
1453 /* mpg_add --- add arbitrary-precision numbers */
1454
1455 static NODE *
1456 mpg_add(NODE *t1, NODE *t2)
1457 {
1458 NODE *r;
1459 int tval;
1460
1461 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1462 r = mpg_integer();
1463 mpz_add(r->mpg_i, t1->mpg_i, t2->mpg_i);
1464 } else {
1465 r = mpg_float();
1466 if (is_mpg_integer(t2))
1467 tval = mpfr_add_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1468 else if (is_mpg_integer(t1))
1469 tval = mpfr_add_z(r->mpg_numbr, t2->mpg_numbr, t1->mpg_i, ROUND_MODE);
1470 else
1471 tval = mpfr_add(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1472 IEEE_FMT(r->mpg_numbr, tval);
1473 }
1474 return r;
1475 }
1476
1477 /* mpg_sub --- subtract arbitrary-precision numbers */
1478
1479 static NODE *
1480 mpg_sub(NODE *t1, NODE *t2)
1481 {
1482 NODE *r;
1483 int tval;
1484
1485 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1486 r = mpg_integer();
1487 mpz_sub(r->mpg_i, t1->mpg_i, t2->mpg_i);
1488 } else {
1489 r = mpg_float();
1490 if (is_mpg_integer(t2))
1491 tval = mpfr_sub_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1492 else if (is_mpg_integer(t1)) {
1493 #if (!defined(MPFR_VERSION) || (MPFR_VERSION < MPFR_VERSION_NUM(3,1,0)))
1494 NODE *tmp = t1;
1495 t1 = t2;
1496 t2 = tmp;
1497 tval = mpfr_sub_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1498 tval = mpfr_neg(r->mpg_numbr, r->mpg_numbr, ROUND_MODE);
1499 t2 = t1;
1500 t1 = tmp;
1501 #else
1502 tval = mpfr_z_sub(r->mpg_numbr, t1->mpg_i, t2->mpg_numbr, ROUND_MODE);
1503 #endif
1504 } else
1505 tval = mpfr_sub(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1506 IEEE_FMT(r->mpg_numbr, tval);
1507 }
1508 return r;
1509 }
1510
1511 /* mpg_mul --- multiply arbitrary-precision numbers */
1512
1513 static NODE *
1514 mpg_mul(NODE *t1, NODE *t2)
1515 {
1516 NODE *r;
1517 int tval;
1518
1519 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1520 r = mpg_integer();
1521 mpz_mul(r->mpg_i, t1->mpg_i, t2->mpg_i);
1522 } else {
1523 r = mpg_float();
1524 if (is_mpg_integer(t2))
1525 tval = mpfr_mul_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1526 else if (is_mpg_integer(t1))
1527 tval = mpfr_mul_z(r->mpg_numbr, t2->mpg_numbr, t1->mpg_i, ROUND_MODE);
1528 else
1529 tval = mpfr_mul(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1530 IEEE_FMT(r->mpg_numbr, tval);
1531 }
1532 return r;
1533 }
1534
1535
1536 /* mpg_pow --- exponentiation involving arbitrary-precision numbers */
1537
1538 static NODE *
1539 mpg_pow(NODE *t1, NODE *t2)
1540 {
1541 NODE *r;
1542 int tval;
1543
1544 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1545 if (mpz_sgn(t2->mpg_i) >= 0 && mpz_fits_ulong_p(t2->mpg_i)) {
1546 r = mpg_integer();
1547 mpz_pow_ui(r->mpg_i, t1->mpg_i, mpz_get_ui(t2->mpg_i));
1548 } else {
1549 mpfr_ptr p1, p2;
1550 p1 = MP_FLOAT(t1);
1551 p2 = MP_FLOAT(t2);
1552 r = mpg_float();
1553 tval = mpfr_pow(r->mpg_numbr, p1, p2, ROUND_MODE);
1554 IEEE_FMT(r->mpg_numbr, tval);
1555 }
1556 } else {
1557 r = mpg_float();
1558 if (is_mpg_integer(t2))
1559 tval = mpfr_pow_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1560 else {
1561 mpfr_ptr p1;
1562 p1 = MP_FLOAT(t1);
1563 tval = mpfr_pow(r->mpg_numbr, p1, t2->mpg_numbr, ROUND_MODE);
1564 }
1565 IEEE_FMT(r->mpg_numbr, tval);
1566 }
1567 return r;
1568 }
1569
1570 /* mpg_div --- arbitrary-precision division */
1571
1572 static NODE *
1573 mpg_div(NODE *t1, NODE *t2)
1574 {
1575 NODE *r;
1576 int tval;
1577
1578 if (is_mpg_integer(t1) && is_mpg_integer(t2)
1579 && (mpz_sgn(t2->mpg_i) != 0) /* not dividing by 0 */
1580 && mpz_divisible_p(t1->mpg_i, t2->mpg_i)
1581 ) {
1582 r = mpg_integer();
1583 mpz_divexact(r->mpg_i, t1->mpg_i, t2->mpg_i);
1584 } else {
1585 mpfr_ptr p1, p2;
1586 p1 = MP_FLOAT(t1);
1587 p2 = MP_FLOAT(t2);
1588 if (mpfr_zero_p(p2))
1589 fatal(_("division by zero attempted"));
1590 r = mpg_float();
1591 tval = mpfr_div(r->mpg_numbr, p1, p2, ROUND_MODE);
1592 IEEE_FMT(r->mpg_numbr, tval);
1593 }
1594 return r;
1595 }
1596
1597 /* mpg_mod --- modulus operation with arbitrary-precision numbers */
1598
1599 static NODE *
1600 mpg_mod(NODE *t1, NODE *t2)
1601 {
1602 NODE *r;
1603 int tval;
1604
1605 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1606 /*
1607 * 8/2014: Originally, this was just
1608 *
1609 * r = mpg_integer();
1610 * mpz_mod(r->mpg_i, t1->mpg_i, t2->mpg_i);
1611 *
1612 * But that gave very strange results with negative numerator:
1613 *
1614 * $ ./gawk -M 'BEGIN { print -15 % 7 }'
1615 * 6
1616 *
1617 * So instead we use mpz_tdiv_qr() to get the correct result
1618 * and just throw away the quotient. We could not find any
1619 * reason why mpz_mod() wasn't working correctly.
1620 */
1621 NODE *dummy_quotient;
1622
1623 if (mpz_sgn(t2->mpg_i) == 0)
1624 fatal(_("division by zero attempted"));
1625 r = mpg_integer();
1626 dummy_quotient = mpg_integer();
1627 mpz_tdiv_qr(dummy_quotient->mpg_i, r->mpg_i, t1->mpg_i, t2->mpg_i);
1628 unref(dummy_quotient);
1629 } else {
1630 mpfr_ptr p1, p2;
1631 p1 = MP_FLOAT(t1);
1632 p2 = MP_FLOAT(t2);
1633 if (mpfr_zero_p(p2))
1634 fatal(_("division by zero attempted in `%%'"));
1635 r = mpg_float();
1636 tval = mpfr_fmod(r->mpg_numbr, p1, p2, ROUND_MODE);
1637 IEEE_FMT(r->mpg_numbr, tval);
1638 }
1639 return r;
1640 }
1641
1642 /*
1643 * mpg_interpret --- pre-exec hook in the interpreter. Handles
1644 * arithmetic operations with MPFR/GMP numbers.
1645 */
1646
1647 static int
1648 mpg_interpret(INSTRUCTION **cp)
1649 {
1650 INSTRUCTION *pc = *cp; /* current instruction */
1651 OPCODE op; /* current opcode */
1652 NODE *r = NULL;
1653 NODE *t1, *t2;
1654 NODE **lhs;
1655 int tval; /* the ternary value returned by a MPFR function */
1656
1657 op = pc->opcode;
1658 if (do_itrace) {
1659 switch (op) {
1660 case Op_plus_i:
1661 case Op_plus:
1662 case Op_minus_i:
1663 case Op_minus:
1664 case Op_times_i:
1665 case Op_times:
1666 case Op_exp_i:
1667 case Op_exp:
1668 case Op_quotient_i:
1669 case Op_quotient:
1670 case Op_mod_i:
1671 case Op_mod:
1672 case Op_preincrement:
1673 case Op_predecrement:
1674 case Op_postincrement:
1675 case Op_postdecrement:
1676 case Op_unary_minus:
1677 case Op_unary_plus:
1678 case Op_assign_plus:
1679 case Op_assign_minus:
1680 case Op_assign_times:
1681 case Op_assign_quotient:
1682 case Op_assign_mod:
1683 case Op_assign_exp:
1684 fprintf(stderr, "++ %s: mpg_interpret\n", opcode2str(op));
1685 fflush(stderr);
1686 break;
1687 default:
1688 return true; /* unhandled */
1689 }
1690 }
1691
1692 switch (op) {
1693 case Op_plus_i:
1694 t2 = force_number(pc->memory);
1695 goto plus;
1696 case Op_plus:
1697 t2 = POP_NUMBER();
1698 plus:
1699 t1 = TOP_NUMBER();
1700 r = mpg_add(t1, t2);
1701 DEREF(t1);
1702 if (op == Op_plus)
1703 DEREF(t2);
1704 REPLACE(r);
1705 break;
1706
1707 case Op_minus_i:
1708 t2 = force_number(pc->memory);
1709 goto minus;
1710 case Op_minus:
1711 t2 = POP_NUMBER();
1712 minus:
1713 t1 = TOP_NUMBER();
1714 r = mpg_sub(t1, t2);
1715 DEREF(t1);
1716 if (op == Op_minus)
1717 DEREF(t2);
1718 REPLACE(r);
1719 break;
1720
1721 case Op_times_i:
1722 t2 = force_number(pc->memory);
1723 goto times;
1724 case Op_times:
1725 t2 = POP_NUMBER();
1726 times:
1727 t1 = TOP_NUMBER();
1728 r = mpg_mul(t1, t2);
1729 DEREF(t1);
1730 if (op == Op_times)
1731 DEREF(t2);
1732 REPLACE(r);
1733 break;
1734
1735 case Op_exp_i:
1736 t2 = force_number(pc->memory);
1737 goto exp;
1738 case Op_exp:
1739 t2 = POP_NUMBER();
1740 exp:
1741 t1 = TOP_NUMBER();
1742 r = mpg_pow(t1, t2);
1743 DEREF(t1);
1744 if (op == Op_exp)
1745 DEREF(t2);
1746 REPLACE(r);
1747 break;
1748
1749 case Op_quotient_i:
1750 t2 = force_number(pc->memory);
1751 goto quotient;
1752 case Op_quotient:
1753 t2 = POP_NUMBER();
1754 quotient:
1755 t1 = TOP_NUMBER();
1756 r = mpg_div(t1, t2);
1757 DEREF(t1);
1758 if (op == Op_quotient)
1759 DEREF(t2);
1760 REPLACE(r);
1761 break;
1762
1763 case Op_mod_i:
1764 t2 = force_number(pc->memory);
1765 goto mod;
1766 case Op_mod:
1767 t2 = POP_NUMBER();
1768 mod:
1769 t1 = TOP_NUMBER();
1770 r = mpg_mod(t1, t2);
1771 DEREF(t1);
1772 if (op == Op_mod)
1773 DEREF(t2);
1774 REPLACE(r);
1775 break;
1776
1777 case Op_preincrement:
1778 case Op_predecrement:
1779 lhs = TOP_ADDRESS();
1780 t1 = *lhs;
1781 force_number(t1);
1782
1783 if (is_mpg_integer(t1)) {
1784 if (t1->valref == 1 && t1->flags == (MALLOC|MPZN|NUMCUR|NUMBER))
1785 /* Efficiency hack. Big speed-up (> 30%) in a tight loop */
1786 r = t1;
1787 else
1788 r = *lhs = mpg_integer();
1789 if (op == Op_preincrement)
1790 mpz_add_ui(r->mpg_i, t1->mpg_i, 1);
1791 else
1792 mpz_sub_ui(r->mpg_i, t1->mpg_i, 1);
1793 } else {
1794
1795 /*
1796 * An optimization like the one above is not going to work
1797 * for a floating-point number. With it,
1798 * gawk -M 'BEGIN { PREC=53; i=2^53+0.0; PREC=113; ++i; print i}'
1799 * will output 2^53 instead of 2^53+1.
1800 */
1801
1802 r = *lhs = mpg_float();
1803 tval = mpfr_add_si(r->mpg_numbr, t1->mpg_numbr,
1804 op == Op_preincrement ? 1 : -1,
1805 ROUND_MODE);
1806 IEEE_FMT(r->mpg_numbr, tval);
1807 }
1808 if (r != t1)
1809 unref(t1);
1810 UPREF(r);
1811 REPLACE(r);
1812 break;
1813
1814 case Op_postincrement:
1815 case Op_postdecrement:
1816 lhs = TOP_ADDRESS();
1817 t1 = *lhs;
1818 force_number(t1);
1819
1820 if (is_mpg_integer(t1)) {
1821 r = mpg_integer();
1822 mpz_set(r->mpg_i, t1->mpg_i);
1823 if (t1->valref == 1 && t1->flags == (MALLOC|MPZN|NUMCUR|NUMBER))
1824 /* Efficiency hack. Big speed-up (> 30%) in a tight loop */
1825 t2 = t1;
1826 else
1827 t2 = *lhs = mpg_integer();
1828 if (op == Op_postincrement)
1829 mpz_add_ui(t2->mpg_i, t1->mpg_i, 1);
1830 else
1831 mpz_sub_ui(t2->mpg_i, t1->mpg_i, 1);
1832 } else {
1833 r = mpg_float();
1834 tval = mpfr_set(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1835 IEEE_FMT(r->mpg_numbr, tval);
1836 t2 = *lhs = mpg_float();
1837 tval = mpfr_add_si(t2->mpg_numbr, t1->mpg_numbr,
1838 op == Op_postincrement ? 1 : -1,
1839 ROUND_MODE);
1840 IEEE_FMT(t2->mpg_numbr, tval);
1841 }
1842 if (t2 != t1)
1843 unref(t1);
1844 REPLACE(r);
1845 break;
1846
1847 case Op_unary_minus:
1848 t1 = TOP_NUMBER();
1849 if (is_mpg_float(t1)) {
1850 r = mpg_float();
1851 tval = mpfr_neg(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1852 IEEE_FMT(r->mpg_numbr, tval);
1853 } else {
1854 if (! is_zero(t1)) {
1855 r = mpg_integer();
1856 mpz_neg(r->mpg_i, t1->mpg_i);
1857 } else {
1858 // have to convert to MPFR for -0.0. sigh
1859 r = mpg_float();
1860 tval = mpfr_set_d(r->mpg_numbr, 0.0, ROUND_MODE);
1861 IEEE_FMT(r->mpg_numbr, tval);
1862 tval = mpfr_neg(r->mpg_numbr, r->mpg_numbr, ROUND_MODE);
1863 IEEE_FMT(r->mpg_numbr, tval);
1864 }
1865 }
1866 DEREF(t1);
1867 REPLACE(r);
1868 break;
1869
1870 case Op_unary_plus:
1871 t1 = TOP_NUMBER();
1872 if (is_mpg_float(t1)) {
1873 r = mpg_float();
1874 tval = mpfr_set(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1875 IEEE_FMT(r->mpg_numbr, tval);
1876 } else {
1877 r = mpg_integer();
1878 mpz_set(r->mpg_i, t1->mpg_i);
1879 }
1880 DEREF(t1);
1881 REPLACE(r);
1882 break;
1883
1884 case Op_assign_plus:
1885 case Op_assign_minus:
1886 case Op_assign_times:
1887 case Op_assign_quotient:
1888 case Op_assign_mod:
1889 case Op_assign_exp:
1890 lhs = POP_ADDRESS();
1891 t1 = *lhs;
1892 force_number(t1);
1893 t2 = TOP_NUMBER();
1894
1895 switch (op) {
1896 case Op_assign_plus:
1897 r = mpg_add(t1, t2);
1898 break;
1899 case Op_assign_minus:
1900 r = mpg_sub(t1, t2);
1901 break;
1902 case Op_assign_times:
1903 r = mpg_mul(t1, t2);
1904 break;
1905 case Op_assign_quotient:
1906 r = mpg_div(t1, t2);
1907 break;
1908 case Op_assign_mod:
1909 r = mpg_mod(t1, t2);
1910 break;
1911 case Op_assign_exp:
1912 r = mpg_pow(t1, t2);
1913 break;
1914 default:
1915 cant_happen("unexpected opcode %s", opcode2str(op));
1916 }
1917
1918 DEREF(t2);
1919 unref(*lhs);
1920 *lhs = r;
1921 UPREF(r);
1922 REPLACE(r);
1923 break;
1924
1925 default:
1926 return true; /* unhandled */
1927 }
1928
1929 *cp = pc->nexti; /* next instruction to execute */
1930 return false;
1931 }
1932
1933
1934 /* mpg_fmt --- output formatted string with special MPFR/GMP conversion specifiers */
1935
1936 const char *
1937 mpg_fmt(const char *mesg, ...)
1938 {
1939 static char *tmp = NULL;
1940 int ret;
1941 va_list args;
1942
1943 if (tmp != NULL) {
1944 mpfr_free_str(tmp);
1945 tmp = NULL;
1946 }
1947 va_start(args, mesg);
1948 ret = mpfr_vasprintf(& tmp, mesg, args);
1949 va_end(args);
1950 if (ret >= 0 && tmp != NULL)
1951 return tmp;
1952 return mesg;
1953 }
1954
1955 /* mpfr_unset --- clear out the MPFR values */
1956
1957 void
1958 mpfr_unset(NODE *n)
1959 {
1960 if (is_mpg_float(n))
1961 mpfr_clear(n->mpg_numbr);
1962 else if (is_mpg_integer(n))
1963 mpz_clear(n->mpg_i);
1964 }
1965
1966 /*
1967 * Custom memory allocation functions for GMP / MPFR. We need these so that the
1968 * persistent memory feature will also work with the -M option.
1969 *
1970 * These just call malloc/realloc/free; if we are using PMA then those are
1971 * redefined as macros to point at the pma functions, so all should "just work."
1972 */
1973
1974 /* mpfr_mem_alloc --- allocate memory */
1975
1976 void *
1977 mpfr_mem_alloc(size_t alloc_size)
1978 {
1979 return malloc(alloc_size);
1980 }
1981
1982 /* mpfr_mem_realloc --- reallocate memory */
1983
1984 void *
1985 mpfr_mem_realloc(void *ptr, size_t old_size, size_t new_size)
1986 {
1987 return realloc(ptr, new_size);
1988 }
1989
1990 /* mpfr_mem_free --- free memory */
1991
1992 void
1993 mpfr_mem_free(void *ptr, size_t size)
1994 {
1995 free(ptr);
1996 }
1997
1998 #else
1999
2000 void
2001 set_PREC()
2002 {
2003 /* dummy function */
2004 }
2005
2006 void
2007 set_ROUNDMODE()
2008 {
2009 /* dummy function */
2010 }
2011
2012 void
2013 mpfr_unset(NODE *n)
2014 {
2015 /* dummy function */
2016 }
2017 #endif