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
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
init_mpfr(mpfr_prec_t prec,const char * rmode)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
cleanup_mpfr(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 *
mpg_node(unsigned int flags)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 *
mpg_make_number(double x)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
mpg_strtoui(mpz_ptr zi,char * str,size_t len,char ** end,int base)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
mpg_maybe_float(const char * str,int use_locale)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
mpg_zero(NODE * n)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
force_mpnum(NODE * n,int do_nondec,int use_locale)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 *
mpg_force_number(NODE * n)348 mpg_force_number(NODE *n)
349 {
350 char *cp, *cpend;
351
352 if ((n->flags & NUMCUR) != 0)
353 return n;
354 n->flags |= NUMCUR;
355
356 /* Trim leading white space, bailing out if there's nothing else */
357 for (cp = n->stptr, cpend = cp + n->stlen;
358 cp < cpend && isspace((unsigned char) *cp); cp++)
359 continue;
360
361 if (cp == cpend)
362 goto badnum;
363
364 /* At this point, we know the string is not entirely white space */
365 /* Trim trailing white space */
366 while (isspace((unsigned char) cpend[-1]))
367 cpend--;
368
369 /*
370 * 2/2007:
371 * POSIX, by way of severe language lawyering, seems to
372 * allow things like "inf" and "nan" to mean something.
373 * So if do_posix, the user gets what he deserves.
374 * This also allows hexadecimal floating point. Ugh.
375 */
376 if (! do_posix) {
377 if (is_alpha((unsigned char) *cp))
378 goto badnum;
379 else if (is_ieee_magic_val(cp)) {
380 if (cpend != cp + 4)
381 goto badnum;
382 /* else
383 fall through */
384 }
385 /* else
386 fall through */
387 }
388 /* else POSIX, so
389 fall through */
390
391 if (force_mpnum(n, (do_non_decimal_data && ! do_traditional), true)) {
392 if ((n->flags & USER_INPUT) != 0) {
393 /* leave USER_INPUT set to indicate a strnum */
394 n->flags &= ~STRING;
395 n->flags |= NUMBER;
396 }
397 } else
398 n->flags &= ~USER_INPUT;
399 return n;
400 badnum:
401 mpg_zero(n);
402 n->flags &= ~USER_INPUT;
403 return n;
404 }
405
406 /* mpg_format_val --- format a numeric value based on format */
407
408 static NODE *
mpg_format_val(const char * format,int index,NODE * s)409 mpg_format_val(const char *format, int index, NODE *s)
410 {
411 NODE *dummy[2], *r;
412 unsigned int oflags;
413
414 if (out_of_range(s)) {
415 const char *result = format_nan_inf(s, 'g');
416 return make_string(result, strlen(result));
417 }
418
419 /* create dummy node for a sole use of format_tree */
420 dummy[1] = s;
421 oflags = s->flags;
422
423 if (is_mpg_integer(s) || mpfr_integer_p(s->mpg_numbr)) {
424 /* integral value, use %d */
425 r = format_tree("%d", 2, dummy, 2);
426 s->stfmt = STFMT_UNUSED;
427 } else {
428 r = format_tree(format, fmt_list[index]->stlen, dummy, 2);
429 assert(r != NULL);
430 s->stfmt = index;
431 }
432 s->flags = oflags;
433 s->stlen = r->stlen;
434 if ((s->flags & (MALLOC|STRCUR)) == (MALLOC|STRCUR))
435 efree(s->stptr);
436 s->stptr = r->stptr;
437 s->flags |= STRCUR;
438 s->strndmode = MPFR_round_mode;
439 freenode(r); /* Do not unref(r)! We want to keep s->stptr == r->stpr. */
440 free_wstr(s);
441 return s;
442 }
443
444 /* mpg_cmp --- compare two numbers */
445
446 int
mpg_cmp(const NODE * t1,const NODE * t2)447 mpg_cmp(const NODE *t1, const NODE *t2)
448 {
449 /*
450 * For the purposes of sorting, NaN is considered greater than
451 * any other value, and all NaN values are considered equivalent and equal.
452 */
453
454 if (is_mpg_float(t1)) {
455 if (is_mpg_float(t2)) {
456 if (mpfr_nan_p(t1->mpg_numbr))
457 return ! mpfr_nan_p(t2->mpg_numbr);
458 if (mpfr_nan_p(t2->mpg_numbr))
459 return -1;
460 return mpfr_cmp(t1->mpg_numbr, t2->mpg_numbr);
461 }
462 if (mpfr_nan_p(t1->mpg_numbr))
463 return 1;
464 return mpfr_cmp_z(t1->mpg_numbr, t2->mpg_i);
465 } else if (is_mpg_float(t2)) {
466 int ret;
467 if (mpfr_nan_p(t2->mpg_numbr))
468 return -1;
469 ret = mpfr_cmp_z(t2->mpg_numbr, t1->mpg_i);
470 return ret > 0 ? -1 : (ret < 0);
471 } else if (is_mpg_integer(t1)) {
472 return mpz_cmp(t1->mpg_i, t2->mpg_i);
473 }
474
475 /* t1 and t2 are AWKNUMs */
476 return cmp_awknums(t1, t2);
477 }
478
479
480 /*
481 * mpg_update_var --- update NR or FNR.
482 * NR_node->var_value(mpz_t) = MNR(mpz_t) * LONG_MAX + NR(long)
483 */
484
485 NODE *
mpg_update_var(NODE * n)486 mpg_update_var(NODE *n)
487 {
488 NODE *val = n->var_value;
489 long nr = 0;
490 mpz_ptr nq = 0;
491
492 if (n == NR_node) {
493 nr = NR;
494 nq = MNR;
495 } else if (n == FNR_node) {
496 nr = FNR;
497 nq = MFNR;
498 } else
499 cant_happen();
500
501 if (mpz_sgn(nq) == 0) {
502 /* Efficiency hack similar to that for AWKNUM */
503 if (is_mpg_float(val) || mpz_get_si(val->mpg_i) != nr) {
504 unref(n->var_value);
505 val = n->var_value = mpg_integer();
506 mpz_set_si(val->mpg_i, nr);
507 }
508 } else {
509 unref(n->var_value);
510 val = n->var_value = mpg_integer();
511 mpz_set_si(val->mpg_i, nr);
512 mpz_addmul_ui(val->mpg_i, nq, LONG_MAX); /* val->mpg_i += nq * LONG_MAX */
513 }
514 return val;
515 }
516
517 /* mpg_set_var --- set NR or FNR */
518
519 long
mpg_set_var(NODE * n)520 mpg_set_var(NODE *n)
521 {
522 long nr = 0;
523 mpz_ptr nq = 0, r;
524 NODE *val = n->var_value;
525
526 if (n == NR_node)
527 nq = MNR;
528 else if (n == FNR_node)
529 nq = MFNR;
530 else
531 cant_happen();
532
533 if (is_mpg_integer(val))
534 r = val->mpg_i;
535 else {
536 /* convert float to integer */
537 mpfr_get_z(mpzval, val->mpg_numbr, MPFR_RNDZ);
538 r = mpzval;
539 }
540 nr = mpz_fdiv_q_ui(nq, r, LONG_MAX); /* nq (MNR or MFNR) is quotient */
541 return nr; /* remainder (NR or FNR) */
542 }
543
544 /* set_PREC --- update MPFR PRECISION related variables when PREC assigned to */
545
546 void
set_PREC()547 set_PREC()
548 {
549 long prec = 0;
550 NODE *val;
551 static const struct ieee_fmt {
552 const char *name;
553 mpfr_prec_t precision;
554 mpfr_exp_t emax;
555 mpfr_exp_t emin;
556 } ieee_fmts[] = {
557 { "half", 11, 16, -23 }, /* binary16 */
558 { "single", 24, 128, -148 }, /* binary32 */
559 { "double", 53, 1024, -1073 }, /* binary64 */
560 { "quad", 113, 16384, -16493 }, /* binary128 */
561 { "oct", 237, 262144, -262377 }, /* binary256, not in the IEEE 754-2008 standard */
562
563 /*
564 * For any bitwidth = 32 * k ( k >= 4),
565 * precision = 13 + bitwidth - int(4 * log2(bitwidth))
566 * emax = 1 << bitwidth - precision - 1
567 * emin = 4 - emax - precision
568 */
569 };
570
571 if (! do_mpfr)
572 return;
573
574 val = fixtype(PREC_node->var_value);
575
576 if ((val->flags & STRING) != 0) {
577 int i, j;
578
579 /* emulate IEEE-754 binary format */
580
581 for (i = 0, j = sizeof(ieee_fmts)/sizeof(ieee_fmts[0]); i < j; i++) {
582 if (strcasecmp(ieee_fmts[i].name, val->stptr) == 0)
583 break;
584 }
585
586 if (i < j) {
587 prec = ieee_fmts[i].precision;
588
589 /*
590 * We *DO NOT* change the MPFR exponent range using
591 * mpfr_set_{emin, emax} here. See format_ieee() for details.
592 */
593 max_exp = ieee_fmts[i].emax;
594 min_exp = ieee_fmts[i].emin;
595
596 do_ieee_fmt = true;
597 }
598 }
599
600 if (prec <= 0) {
601 force_number(val);
602 prec = get_number_si(val);
603 if (prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX) {
604 force_string(val);
605 warning(_("PREC value `%.*s' is invalid"), (int) val->stlen, val->stptr);
606 prec = 0;
607 } else
608 do_ieee_fmt = false;
609 }
610
611 if (prec > 0)
612 mpfr_set_default_prec(default_prec = prec);
613 }
614
615
616 /* get_rnd_mode --- convert string to MPFR rounding mode */
617
618 static mpfr_rnd_t
get_rnd_mode(const char rmode)619 get_rnd_mode(const char rmode)
620 {
621 switch (rmode) {
622 case 'N':
623 case 'n':
624 return MPFR_RNDN; /* round to nearest (IEEE-754 roundTiesToEven) */
625 case 'Z':
626 case 'z':
627 return MPFR_RNDZ; /* round toward zero (IEEE-754 roundTowardZero) */
628 case 'U':
629 case 'u':
630 return MPFR_RNDU; /* round toward plus infinity (IEEE-754 roundTowardPositive) */
631 case 'D':
632 case 'd':
633 return MPFR_RNDD; /* round toward minus infinity (IEEE-754 roundTowardNegative) */
634 #if defined(MPFR_VERSION_MAJOR) && MPFR_VERSION_MAJOR > 2
635 case 'A':
636 case 'a':
637 return MPFR_RNDA; /* round away from zero */
638 #endif
639 default:
640 break;
641 }
642 return (mpfr_rnd_t) -1;
643 }
644
645 /*
646 * set_ROUNDMODE --- update MPFR rounding mode related variables
647 * when ROUNDMODE assigned to
648 */
649
650 void
set_ROUNDMODE()651 set_ROUNDMODE()
652 {
653 if (do_mpfr) {
654 mpfr_rnd_t rndm = (mpfr_rnd_t) -1;
655 NODE *n;
656 n = force_string(ROUNDMODE_node->var_value);
657 if (n->stlen == 1)
658 rndm = get_rnd_mode(n->stptr[0]);
659 if (rndm != -1) {
660 mpfr_set_default_rounding_mode(rndm);
661 ROUND_MODE = rndm;
662 MPFR_round_mode = n->stptr[0];
663 } else
664 warning(_("ROUNDMODE value `%.*s' is invalid"), (int) n->stlen, n->stptr);
665 }
666 }
667
668
669 /* format_ieee --- make sure a number follows IEEE-754 floating-point standard */
670
671 int
format_ieee(mpfr_ptr x,int tval)672 format_ieee(mpfr_ptr x, int tval)
673 {
674 /*
675 * The MPFR doc says that it's our responsibility to make sure all numbers
676 * including those previously created are in range after we've changed the
677 * exponent range. Most MPFR operations and functions require
678 * the input arguments to have exponents within the current exponent range.
679 * Any argument outside the range results in a MPFR assertion failure
680 * like this:
681 *
682 * $ gawk -M 'BEGIN { x=1.0e-10000; print x+0; PREC="double"; print x+0}'
683 * 1e-10000
684 * init2.c:52: MPFR assertion failed ....
685 *
686 * A "naive" approach would be to keep track of the ternary state and
687 * the rounding mode for each number, and make sure it is in the current
688 * exponent range (using mpfr_check_range) before using it in an
689 * operation or function. Instead, we adopt the following strategy.
690 *
691 * When gawk starts, the exponent range is the MPFR default
692 * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT]. Any number that gawk
693 * creates must have exponent in this range (excluding infinities, NaNs and zeros).
694 * Each MPFR operation or function is performed with this default exponent
695 * range.
696 *
697 * When emulating IEEE-754 format, the exponents are *temporarily* changed,
698 * mpfr_check_range is called to make sure the number is in the new range,
699 * and mpfr_subnormalize is used to round following the rules of subnormal
700 * arithmetic. The exponent range is then *restored* to the original value
701 * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT].
702 */
703
704 (void) mpfr_set_emin(min_exp);
705 (void) mpfr_set_emax(max_exp);
706 tval = mpfr_check_range(x, tval, ROUND_MODE);
707 tval = mpfr_subnormalize(x, tval, ROUND_MODE);
708 (void) mpfr_set_emin(MPFR_EMIN_DEFAULT);
709 (void) mpfr_set_emax(MPFR_EMAX_DEFAULT);
710 return tval;
711 }
712
713
714 /* do_mpfr_atan2 --- do the atan2 function */
715
716 NODE *
do_mpfr_atan2(int nargs)717 do_mpfr_atan2(int nargs)
718 {
719 NODE *t1, *t2, *res;
720 mpfr_ptr p1, p2;
721 int tval;
722
723 t2 = POP_SCALAR();
724 t1 = POP_SCALAR();
725
726 if (do_lint) {
727 if ((fixtype(t1)->flags & NUMBER) == 0)
728 lintwarn(_("atan2: received non-numeric first argument"));
729 if ((fixtype(t2)->flags & NUMBER) == 0)
730 lintwarn(_("atan2: received non-numeric second argument"));
731 }
732 force_number(t1);
733 force_number(t2);
734
735 p1 = MP_FLOAT(t1);
736 p2 = MP_FLOAT(t2);
737 res = mpg_float();
738 /* See MPFR documentation for handling of special values like +inf as an argument */
739 tval = mpfr_atan2(res->mpg_numbr, p1, p2, ROUND_MODE);
740 IEEE_FMT(res->mpg_numbr, tval);
741
742 DEREF(t1);
743 DEREF(t2);
744 return res;
745 }
746
747 /* do_mpfr_func --- run an MPFR function - not inline, for debugging */
748
749 static inline NODE *
do_mpfr_func(const char * name,int (* mpfr_func)(mpfr_ptr,mpfr_srcptr,mpfr_rnd_t),int nargs,bool warn_negative)750 do_mpfr_func(const char *name,
751 int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t),
752 int nargs, bool warn_negative)
753 {
754 NODE *t1, *res;
755 mpfr_ptr p1;
756 int tval;
757 mpfr_prec_t argprec;
758
759 t1 = POP_SCALAR();
760 if (do_lint && (fixtype(t1)->flags & NUMBER) == 0)
761 lintwarn(_("%s: received non-numeric argument"), name);
762
763 force_number(t1);
764 p1 = MP_FLOAT(t1);
765 if (warn_negative && mpfr_sgn(p1) < 0) {
766 force_string(t1);
767 warning(_("%s: received negative argument %.*s"), name, (int) t1->stlen, t1->stptr);
768 }
769 res = mpg_float();
770 if ((argprec = mpfr_get_prec(p1)) > default_prec)
771 mpfr_set_prec(res->mpg_numbr, argprec); /* needed at least for sqrt() */
772 tval = mpfr_func(res->mpg_numbr, p1, ROUND_MODE);
773 IEEE_FMT(res->mpg_numbr, tval);
774 DEREF(t1);
775 return res;
776 }
777
778 #define SPEC_MATH(X, WN) \
779 NODE *result; \
780 result = do_mpfr_func(#X, mpfr_##X, nargs, WN); \
781 return result
782
783 /* do_mpfr_sin --- do the sin function */
784
785 NODE *
do_mpfr_sin(int nargs)786 do_mpfr_sin(int nargs)
787 {
788 SPEC_MATH(sin, false);
789 }
790
791 /* do_mpfr_cos --- do the cos function */
792
793 NODE *
do_mpfr_cos(int nargs)794 do_mpfr_cos(int nargs)
795 {
796 SPEC_MATH(cos, false);
797 }
798
799 /* do_mpfr_exp --- exponential function */
800
801 NODE *
do_mpfr_exp(int nargs)802 do_mpfr_exp(int nargs)
803 {
804 SPEC_MATH(exp, false);
805 }
806
807 /* do_mpfr_log --- the log function */
808
809 NODE *
do_mpfr_log(int nargs)810 do_mpfr_log(int nargs)
811 {
812 SPEC_MATH(log, true);
813 }
814
815 /* do_mpfr_sqrt --- do the sqrt function */
816
817 NODE *
do_mpfr_sqrt(int nargs)818 do_mpfr_sqrt(int nargs)
819 {
820 SPEC_MATH(sqrt, true);
821 }
822
823 /* do_mpfr_int --- convert double to int for awk */
824
825 NODE *
do_mpfr_int(int nargs)826 do_mpfr_int(int nargs)
827 {
828 NODE *tmp, *r;
829
830 tmp = POP_SCALAR();
831 if (do_lint && (fixtype(tmp)->flags & NUMBER) == 0)
832 lintwarn(_("int: received non-numeric argument"));
833 force_number(tmp);
834
835 if (is_mpg_integer(tmp)) {
836 r = mpg_integer();
837 mpz_set(r->mpg_i, tmp->mpg_i);
838 } else {
839 if (! mpfr_number_p(tmp->mpg_numbr)) {
840 /* [+-]inf or NaN */
841 return tmp;
842 }
843
844 r = mpg_integer();
845 mpfr_get_z(r->mpg_i, tmp->mpg_numbr, MPFR_RNDZ);
846 }
847
848 DEREF(tmp);
849 return r;
850 }
851
852 /* do_mpfr_compl --- perform a ~ operation */
853
854 NODE *
do_mpfr_compl(int nargs)855 do_mpfr_compl(int nargs)
856 {
857 NODE *tmp, *r;
858 mpz_ptr zptr;
859
860 tmp = POP_SCALAR();
861 if (do_lint && (fixtype(tmp)->flags & NUMBER) == 0)
862 lintwarn(_("compl: received non-numeric argument"));
863
864 force_number(tmp);
865 if (is_mpg_float(tmp)) {
866 mpfr_ptr p = tmp->mpg_numbr;
867
868 if (! mpfr_number_p(p)) {
869 /* [+-]inf or NaN */
870 return tmp;
871 }
872 if (mpfr_sgn(p) < 0)
873 fatal("%s",
874 mpg_fmt(_("compl(%Rg): negative value is not allowed"), p)
875 );
876 if (do_lint) {
877 if (! mpfr_integer_p(p))
878 lintwarn("%s",
879 mpg_fmt(_("comp(%Rg): fractional value will be truncated"), p)
880 );
881 }
882
883 mpfr_get_z(mpzval, p, MPFR_RNDZ); /* float to integer conversion */
884 zptr = mpzval;
885 } else {
886 /* (tmp->flags & MPZN) != 0 */
887 zptr = tmp->mpg_i;
888 if (mpz_sgn(zptr) < 0)
889 fatal("%s",
890 mpg_fmt(_("compl(%Zd): negative values are not allowed"), zptr)
891 );
892 }
893
894 r = mpg_integer();
895 mpz_com(r->mpg_i, zptr);
896 DEREF(tmp);
897 return r;
898 }
899
900 /* get_intval --- get the (converted) integral operand of a binary function. */
901
902 static mpz_ptr
get_intval(NODE * t1,int argnum,const char * op)903 get_intval(NODE *t1, int argnum, const char *op)
904 {
905 mpz_ptr pz;
906
907 if (do_lint && (fixtype(t1)->flags & NUMBER) == 0)
908 lintwarn(_("%s: received non-numeric argument #%d"), op, argnum);
909
910 (void) force_number(t1);
911
912 if (is_mpg_float(t1)) {
913 mpfr_ptr left = t1->mpg_numbr;
914 if (! mpfr_number_p(left)) {
915 /* inf or NaN */
916 if (do_lint)
917 lintwarn("%s",
918 mpg_fmt(_("%s: argument #%d has invalid value %Rg, using 0"),
919 op, argnum, left)
920 );
921
922 emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
923 mpz_init(pz);
924 return pz; /* should be freed */
925 }
926
927 if (mpfr_sgn(left) < 0)
928 fatal("%s",
929 mpg_fmt(_("%s: argument #%d negative value %Rg is not allowed"),
930 op, argnum, left)
931 );
932
933 if (do_lint) {
934 if (! mpfr_integer_p(left))
935 lintwarn("%s",
936 mpg_fmt(_("%s: argument #%d fractional value %Rg will be truncated"),
937 op, argnum, left)
938 );
939 }
940
941 emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
942 mpz_init(pz);
943 mpfr_get_z(pz, left, MPFR_RNDZ); /* float to integer conversion */
944 return pz; /* should be freed */
945 }
946 /* (t1->flags & MPZN) != 0 */
947 pz = t1->mpg_i;
948 if (mpz_sgn(pz) < 0)
949 fatal("%s",
950 mpg_fmt(_("%s: argument #%d negative value %Zd is not allowed"),
951 op, argnum, pz)
952 );
953
954 return pz; /* must not be freed */
955 }
956
957
958 /* free_intval --- free the converted integer value returned by get_intval() */
959
960 static inline void
free_intval(NODE * t,mpz_ptr pz)961 free_intval(NODE *t, mpz_ptr pz)
962 {
963 if ((t->flags & MPZN) == 0) {
964 mpz_clear(pz);
965 efree(pz);
966 }
967 }
968
969
970 /* do_mpfr_lshift --- perform a << operation */
971
972 NODE *
do_mpfr_lshift(int nargs)973 do_mpfr_lshift(int nargs)
974 {
975 NODE *t1, *t2, *res;
976 unsigned long shift;
977 mpz_ptr pz1, pz2;
978
979 t2 = POP_SCALAR();
980 t1 = POP_SCALAR();
981
982 pz1 = get_intval(t1, 1, "lshift");
983 pz2 = get_intval(t2, 2, "lshift");
984
985 /*
986 * mpz_get_ui: If op is too big to fit an unsigned long then just
987 * the least significant bits that do fit are returned.
988 * The sign of op is ignored, only the absolute value is used.
989 */
990
991 shift = mpz_get_ui(pz2); /* GMP integer => unsigned long conversion */
992 res = mpg_integer();
993 mpz_mul_2exp(res->mpg_i, pz1, shift); /* res = pz1 * 2^shift */
994
995 free_intval(t1, pz1);
996 free_intval(t2, pz2);
997 DEREF(t2);
998 DEREF(t1);
999 return res;
1000 }
1001
1002 /* do_mpfr_rshift --- perform a >> operation */
1003
1004 NODE *
do_mpfr_rshift(int nargs)1005 do_mpfr_rshift(int nargs)
1006 {
1007 NODE *t1, *t2, *res;
1008 unsigned long shift;
1009 mpz_ptr pz1, pz2;
1010
1011 t2 = POP_SCALAR();
1012 t1 = POP_SCALAR();
1013
1014 pz1 = get_intval(t1, 1, "rshift");
1015 pz2 = get_intval(t2, 2, "rshift");
1016
1017 /* N.B: See do_mpfp_lshift. */
1018 shift = mpz_get_ui(pz2); /* GMP integer => unsigned long conversion */
1019 res = mpg_integer();
1020 mpz_fdiv_q_2exp(res->mpg_i, pz1, shift); /* res = pz1 / 2^shift, round towards -inf */
1021
1022 free_intval(t1, pz1);
1023 free_intval(t2, pz2);
1024 DEREF(t2);
1025 DEREF(t1);
1026 return res;
1027 }
1028
1029
1030 /* do_mpfr_and --- perform an & operation */
1031
1032 NODE *
do_mpfr_and(int nargs)1033 do_mpfr_and(int nargs)
1034 {
1035 NODE *t1, *t2, *res;
1036 mpz_ptr pz1, pz2;
1037 int i;
1038
1039 if (nargs < 2)
1040 fatal(_("and: called with less than two arguments"));
1041
1042 t2 = POP_SCALAR();
1043 pz2 = get_intval(t2, nargs, "and");
1044
1045 res = mpg_integer();
1046 for (i = 1; i < nargs; i++) {
1047 t1 = POP_SCALAR();
1048 pz1 = get_intval(t1, nargs - i, "and");
1049 mpz_and(res->mpg_i, pz1, pz2);
1050 free_intval(t1, pz1);
1051 DEREF(t1);
1052 if (i == 1) {
1053 free_intval(t2, pz2);
1054 DEREF(t2);
1055 }
1056 pz2 = res->mpg_i;
1057 }
1058 return res;
1059 }
1060
1061
1062 /* do_mpfr_or --- perform an | operation */
1063
1064 NODE *
do_mpfr_or(int nargs)1065 do_mpfr_or(int nargs)
1066 {
1067 NODE *t1, *t2, *res;
1068 mpz_ptr pz1, pz2;
1069 int i;
1070
1071 if (nargs < 2)
1072 fatal(_("or: called with less than two arguments"));
1073
1074 t2 = POP_SCALAR();
1075 pz2 = get_intval(t2, nargs, "or");
1076
1077 res = mpg_integer();
1078 for (i = 1; i < nargs; i++) {
1079 t1 = POP_SCALAR();
1080 pz1 = get_intval(t1, nargs - i, "or");
1081 mpz_ior(res->mpg_i, pz1, pz2);
1082 free_intval(t1, pz1);
1083 DEREF(t1);
1084 if (i == 1) {
1085 free_intval(t2, pz2);
1086 DEREF(t2);
1087 }
1088 pz2 = res->mpg_i;
1089 }
1090 return res;
1091 }
1092
1093 /* do_mpfr_xor --- perform an ^ operation */
1094
1095 NODE *
do_mpfr_xor(int nargs)1096 do_mpfr_xor(int nargs)
1097 {
1098 NODE *t1, *t2, *res;
1099 mpz_ptr pz1, pz2;
1100 int i;
1101
1102 if (nargs < 2)
1103 fatal(_("xor: called with less than two arguments"));
1104
1105 t2 = POP_SCALAR();
1106 pz2 = get_intval(t2, nargs, "xor");
1107
1108 res = mpg_integer();
1109 for (i = 1; i < nargs; i++) {
1110 t1 = POP_SCALAR();
1111 pz1 = get_intval(t1, nargs - i, "xor");
1112 mpz_xor(res->mpg_i, pz1, pz2);
1113 free_intval(t1, pz1);
1114 DEREF(t1);
1115 if (i == 1) {
1116 free_intval(t2, pz2);
1117 DEREF(t2);
1118 }
1119 pz2 = res->mpg_i;
1120 }
1121 return res;
1122 }
1123
1124 /* do_mpfr_strtonum --- the strtonum function */
1125
1126 NODE *
do_mpfr_strtonum(int nargs)1127 do_mpfr_strtonum(int nargs)
1128 {
1129 NODE *tmp, *r;
1130
1131 tmp = fixtype(POP_SCALAR());
1132 if ((tmp->flags & NUMBER) == 0) {
1133 r = mpg_integer(); /* will be changed to MPFR float if necessary in force_mpnum() */
1134 r->stptr = tmp->stptr;
1135 r->stlen = tmp->stlen;
1136 force_mpnum(r, true, use_lc_numeric);
1137 r->stptr = NULL;
1138 r->stlen = 0;
1139 r->wstptr = NULL;
1140 r->wstlen = 0;
1141 } else if (is_mpg_float(tmp)) {
1142 int tval;
1143 r = mpg_float();
1144 tval = mpfr_set(r->mpg_numbr, tmp->mpg_numbr, ROUND_MODE);
1145 IEEE_FMT(r->mpg_numbr, tval);
1146 } else {
1147 r = mpg_integer();
1148 mpz_set(r->mpg_i, tmp->mpg_i);
1149 }
1150
1151 DEREF(tmp);
1152 return r;
1153 }
1154
1155
1156 static bool firstrand = true;
1157 static gmp_randstate_t state;
1158 static mpz_t seed; /* current seed */
1159
1160 /* do_mpfr_rand --- do the rand function */
1161
1162 NODE *
do_mpfr_rand(int nargs ATTRIBUTE_UNUSED)1163 do_mpfr_rand(int nargs ATTRIBUTE_UNUSED)
1164 {
1165 NODE *res;
1166 int tval;
1167
1168 if (firstrand) {
1169 #if 0
1170 /* Choose the default algorithm */
1171 gmp_randinit_default(state);
1172 #endif
1173 /*
1174 * Choose a specific (Mersenne Twister) algorithm in case the default
1175 * changes in the future.
1176 */
1177
1178 gmp_randinit_mt(state);
1179
1180 mpz_init(seed);
1181 mpz_set_ui(seed, 1);
1182 /* seed state */
1183 gmp_randseed(state, seed);
1184 firstrand = false;
1185 }
1186 res = mpg_float();
1187 tval = mpfr_urandomb(res->mpg_numbr, state);
1188 IEEE_FMT(res->mpg_numbr, tval);
1189 return res;
1190 }
1191
1192
1193 /* do_mpfr_srand --- seed the random number generator */
1194
1195 NODE *
do_mpfr_srand(int nargs)1196 do_mpfr_srand(int nargs)
1197 {
1198 NODE *res;
1199
1200 if (firstrand) {
1201 #if 0
1202 /* Choose the default algorithm */
1203 gmp_randinit_default(state);
1204 #endif
1205 /*
1206 * Choose a specific algorithm (Mersenne Twister) in case default
1207 * changes in the future.
1208 */
1209
1210 gmp_randinit_mt(state);
1211
1212 mpz_init(seed);
1213 mpz_set_ui(seed, 1);
1214 /* No need to seed state, will change it below */
1215 firstrand = false;
1216 }
1217
1218 res = mpg_integer();
1219 mpz_set(res->mpg_i, seed); /* previous seed */
1220
1221 if (nargs == 0)
1222 mpz_set_ui(seed, (unsigned long) time((time_t *) 0));
1223 else {
1224 NODE *tmp;
1225 tmp = POP_SCALAR();
1226 if (do_lint && (fixtype(tmp)->flags & NUMBER) == 0)
1227 lintwarn(_("srand: received non-numeric argument"));
1228 force_number(tmp);
1229 if (is_mpg_float(tmp))
1230 mpfr_get_z(seed, tmp->mpg_numbr, MPFR_RNDZ);
1231 else /* MP integer */
1232 mpz_set(seed, tmp->mpg_i);
1233 DEREF(tmp);
1234 }
1235
1236 gmp_randseed(state, seed);
1237 return res;
1238 }
1239
1240 #ifdef SUPPLY_INTDIV
1241 /* do_mpfr_intdiv --- do integer division, return quotient and remainder in dest array */
1242
1243 /*
1244 * We define the semantics as:
1245 * numerator = int(numerator)
1246 * denominator = int(denonmator)
1247 * quotient = int(numerator / denomator)
1248 * remainder = int(numerator % denomator)
1249 */
1250
1251 NODE *
do_mpfr_intdiv(int nargs)1252 do_mpfr_intdiv(int nargs)
1253 {
1254 NODE *numerator, *denominator, *result;
1255 NODE *num, *denom;
1256 NODE *quotient, *remainder;
1257 NODE *sub, **lhs;
1258
1259 result = POP_PARAM();
1260 if (result->type != Node_var_array)
1261 fatal(_("intdiv: third argument is not an array"));
1262 assoc_clear(result);
1263
1264 denominator = POP_SCALAR();
1265 numerator = POP_SCALAR();
1266
1267 if (do_lint) {
1268 if ((fixtype(numerator)->flags & NUMBER) == 0)
1269 lintwarn(_("intdiv: received non-numeric first argument"));
1270 if ((fixtype(denominator)->flags & NUMBER) == 0)
1271 lintwarn(_("intdiv: received non-numeric second argument"));
1272 }
1273
1274 (void) force_number(numerator);
1275 (void) force_number(denominator);
1276
1277 /* convert numerator and denominator to integer */
1278 if (is_mpg_integer(numerator)) {
1279 num = mpg_integer();
1280 mpz_set(num->mpg_i, numerator->mpg_i);
1281 } else {
1282 if (! mpfr_number_p(numerator->mpg_numbr)) {
1283 /* [+-]inf or NaN */
1284 unref(numerator);
1285 unref(denominator);
1286 return make_number((AWKNUM) -1);
1287 }
1288
1289 num = mpg_integer();
1290 mpfr_get_z(num->mpg_i, numerator->mpg_numbr, MPFR_RNDZ);
1291 }
1292
1293 if (is_mpg_integer(denominator)) {
1294 denom = mpg_integer();
1295 mpz_set(denom->mpg_i, denominator->mpg_i);
1296 } else {
1297 if (! mpfr_number_p(denominator->mpg_numbr)) {
1298 /* [+-]inf or NaN */
1299 unref(numerator);
1300 unref(denominator);
1301 unref(num);
1302 return make_number((AWKNUM) -1);
1303 }
1304
1305 denom = mpg_integer();
1306 mpfr_get_z(denom->mpg_i, denominator->mpg_numbr, MPFR_RNDZ);
1307 }
1308
1309 if (mpz_sgn(denom->mpg_i) == 0)
1310 fatal(_("intdiv: division by zero attempted"));
1311
1312 quotient = mpg_integer();
1313 remainder = mpg_integer();
1314
1315 /* do the division */
1316 mpz_tdiv_qr(quotient->mpg_i, remainder->mpg_i, num->mpg_i, denom->mpg_i);
1317 unref(num);
1318 unref(denom);
1319 unref(numerator);
1320 unref(denominator);
1321
1322 sub = make_string("quotient", 8);
1323 assoc_set(result, sub, quotient);
1324
1325 sub = make_string("remainder", 9);
1326 assoc_set(result, sub, remainder);
1327
1328 return make_number((AWKNUM) 0.0);
1329 }
1330 #endif /* SUPPLY_INTDIV */
1331
1332 /*
1333 * mpg_tofloat --- convert an arbitrary-precision integer operand to
1334 * a float without loss of precision. It is assumed that the
1335 * MPFR variable has already been initialized.
1336 */
1337
1338 static inline mpfr_ptr
mpg_tofloat(mpfr_ptr mf,mpz_ptr mz)1339 mpg_tofloat(mpfr_ptr mf, mpz_ptr mz)
1340 {
1341 size_t prec;
1342
1343 /*
1344 * When implicitely converting a GMP integer operand to a MPFR float, use
1345 * a precision sufficiently large to hold the converted value exactly.
1346 *
1347 * $ ./gawk -M 'BEGIN { print 13 % 2 }'
1348 * 1
1349 * If the user-specified precision is used to convert the integer 13 to a
1350 * float, one will get:
1351 * $ ./gawk -M 'BEGIN { PREC=2; print 13 % 2.0 }'
1352 * 0
1353 */
1354
1355 prec = mpz_sizeinbase(mz, 2); /* most significant 1 bit position starting at 1 */
1356 if (prec > PRECISION_MIN) {
1357 prec -= (size_t) mpz_scan1(mz, 0); /* least significant 1 bit index starting at 0 */
1358 if (prec > MPFR_PREC_MAX)
1359 prec = MPFR_PREC_MAX;
1360 else if (prec < PRECISION_MIN)
1361 prec = PRECISION_MIN;
1362 }
1363 else
1364 prec = PRECISION_MIN;
1365 /*
1366 * Always set the precision to avoid hysteresis, since do_mpfr_func
1367 * may copy our precision.
1368 */
1369 if (prec != mpfr_get_prec(mf))
1370 mpfr_set_prec(mf, prec);
1371
1372 mpfr_set_z(mf, mz, ROUND_MODE);
1373 return mf;
1374 }
1375
1376
1377 /* mpg_add --- add arbitrary-precision numbers */
1378
1379 static NODE *
mpg_add(NODE * t1,NODE * t2)1380 mpg_add(NODE *t1, NODE *t2)
1381 {
1382 NODE *r;
1383 int tval;
1384
1385 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1386 r = mpg_integer();
1387 mpz_add(r->mpg_i, t1->mpg_i, t2->mpg_i);
1388 } else {
1389 r = mpg_float();
1390 if (is_mpg_integer(t2))
1391 tval = mpfr_add_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1392 else if (is_mpg_integer(t1))
1393 tval = mpfr_add_z(r->mpg_numbr, t2->mpg_numbr, t1->mpg_i, ROUND_MODE);
1394 else
1395 tval = mpfr_add(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1396 IEEE_FMT(r->mpg_numbr, tval);
1397 }
1398 return r;
1399 }
1400
1401 /* mpg_sub --- subtract arbitrary-precision numbers */
1402
1403 static NODE *
mpg_sub(NODE * t1,NODE * t2)1404 mpg_sub(NODE *t1, NODE *t2)
1405 {
1406 NODE *r;
1407 int tval;
1408
1409 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1410 r = mpg_integer();
1411 mpz_sub(r->mpg_i, t1->mpg_i, t2->mpg_i);
1412 } else {
1413 r = mpg_float();
1414 if (is_mpg_integer(t2))
1415 tval = mpfr_sub_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1416 else if (is_mpg_integer(t1)) {
1417 #if (!defined(MPFR_VERSION) || (MPFR_VERSION < MPFR_VERSION_NUM(3,1,0)))
1418 NODE *tmp = t1;
1419 t1 = t2;
1420 t2 = tmp;
1421 tval = mpfr_sub_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1422 tval = mpfr_neg(r->mpg_numbr, r->mpg_numbr, ROUND_MODE);
1423 t2 = t1;
1424 t1 = tmp;
1425 #else
1426 tval = mpfr_z_sub(r->mpg_numbr, t1->mpg_i, t2->mpg_numbr, ROUND_MODE);
1427 #endif
1428 } else
1429 tval = mpfr_sub(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1430 IEEE_FMT(r->mpg_numbr, tval);
1431 }
1432 return r;
1433 }
1434
1435 /* mpg_mul --- multiply arbitrary-precision numbers */
1436
1437 static NODE *
mpg_mul(NODE * t1,NODE * t2)1438 mpg_mul(NODE *t1, NODE *t2)
1439 {
1440 NODE *r;
1441 int tval;
1442
1443 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1444 r = mpg_integer();
1445 mpz_mul(r->mpg_i, t1->mpg_i, t2->mpg_i);
1446 } else {
1447 r = mpg_float();
1448 if (is_mpg_integer(t2))
1449 tval = mpfr_mul_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1450 else if (is_mpg_integer(t1))
1451 tval = mpfr_mul_z(r->mpg_numbr, t2->mpg_numbr, t1->mpg_i, ROUND_MODE);
1452 else
1453 tval = mpfr_mul(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1454 IEEE_FMT(r->mpg_numbr, tval);
1455 }
1456 return r;
1457 }
1458
1459
1460 /* mpg_pow --- exponentiation involving arbitrary-precision numbers */
1461
1462 static NODE *
mpg_pow(NODE * t1,NODE * t2)1463 mpg_pow(NODE *t1, NODE *t2)
1464 {
1465 NODE *r;
1466 int tval;
1467
1468 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1469 if (mpz_sgn(t2->mpg_i) >= 0 && mpz_fits_ulong_p(t2->mpg_i)) {
1470 r = mpg_integer();
1471 mpz_pow_ui(r->mpg_i, t1->mpg_i, mpz_get_ui(t2->mpg_i));
1472 } else {
1473 mpfr_ptr p1, p2;
1474 p1 = MP_FLOAT(t1);
1475 p2 = MP_FLOAT(t2);
1476 r = mpg_float();
1477 tval = mpfr_pow(r->mpg_numbr, p1, p2, ROUND_MODE);
1478 IEEE_FMT(r->mpg_numbr, tval);
1479 }
1480 } else {
1481 r = mpg_float();
1482 if (is_mpg_integer(t2))
1483 tval = mpfr_pow_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1484 else {
1485 mpfr_ptr p1;
1486 p1 = MP_FLOAT(t1);
1487 tval = mpfr_pow(r->mpg_numbr, p1, t2->mpg_numbr, ROUND_MODE);
1488 }
1489 IEEE_FMT(r->mpg_numbr, tval);
1490 }
1491 return r;
1492 }
1493
1494 /* mpg_div --- arbitrary-precision division */
1495
1496 static NODE *
mpg_div(NODE * t1,NODE * t2)1497 mpg_div(NODE *t1, NODE *t2)
1498 {
1499 NODE *r;
1500 int tval;
1501
1502 if (is_mpg_integer(t1) && is_mpg_integer(t2)
1503 && (mpz_sgn(t2->mpg_i) != 0) /* not dividing by 0 */
1504 && mpz_divisible_p(t1->mpg_i, t2->mpg_i)
1505 ) {
1506 r = mpg_integer();
1507 mpz_divexact(r->mpg_i, t1->mpg_i, t2->mpg_i);
1508 } else {
1509 mpfr_ptr p1, p2;
1510 p1 = MP_FLOAT(t1);
1511 p2 = MP_FLOAT(t2);
1512 if (mpfr_zero_p(p2))
1513 fatal(_("division by zero attempted"));
1514 r = mpg_float();
1515 tval = mpfr_div(r->mpg_numbr, p1, p2, ROUND_MODE);
1516 IEEE_FMT(r->mpg_numbr, tval);
1517 }
1518 return r;
1519 }
1520
1521 /* mpg_mod --- modulus operation with arbitrary-precision numbers */
1522
1523 static NODE *
mpg_mod(NODE * t1,NODE * t2)1524 mpg_mod(NODE *t1, NODE *t2)
1525 {
1526 NODE *r;
1527 int tval;
1528
1529 if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1530 /*
1531 * 8/2014: Originally, this was just
1532 *
1533 * r = mpg_integer();
1534 * mpz_mod(r->mpg_i, t1->mpg_i, t2->mpg_i);
1535 *
1536 * But that gave very strange results with negative numerator:
1537 *
1538 * $ ./gawk -M 'BEGIN { print -15 % 7 }'
1539 * 6
1540 *
1541 * So instead we use mpz_tdiv_qr() to get the correct result
1542 * and just throw away the quotient. We could not find any
1543 * reason why mpz_mod() wasn't working correctly.
1544 */
1545 NODE *dummy_quotient;
1546
1547 if (mpz_sgn(t2->mpg_i) == 0)
1548 fatal(_("division by zero attempted"));
1549 r = mpg_integer();
1550 dummy_quotient = mpg_integer();
1551 mpz_tdiv_qr(dummy_quotient->mpg_i, r->mpg_i, t1->mpg_i, t2->mpg_i);
1552 unref(dummy_quotient);
1553 } else {
1554 mpfr_ptr p1, p2;
1555 p1 = MP_FLOAT(t1);
1556 p2 = MP_FLOAT(t2);
1557 if (mpfr_zero_p(p2))
1558 fatal(_("division by zero attempted in `%%'"));
1559 r = mpg_float();
1560 tval = mpfr_fmod(r->mpg_numbr, p1, p2, ROUND_MODE);
1561 IEEE_FMT(r->mpg_numbr, tval);
1562 }
1563 return r;
1564 }
1565
1566 /*
1567 * mpg_interpret --- pre-exec hook in the interpreter. Handles
1568 * arithmetic operations with MPFR/GMP numbers.
1569 */
1570
1571 static int
mpg_interpret(INSTRUCTION ** cp)1572 mpg_interpret(INSTRUCTION **cp)
1573 {
1574 INSTRUCTION *pc = *cp; /* current instruction */
1575 OPCODE op; /* current opcode */
1576 NODE *r = NULL;
1577 NODE *t1, *t2;
1578 NODE **lhs;
1579 int tval; /* the ternary value returned by a MPFR function */
1580
1581 op = pc->opcode;
1582 if (do_itrace) {
1583 switch (op) {
1584 case Op_plus_i:
1585 case Op_plus:
1586 case Op_minus_i:
1587 case Op_minus:
1588 case Op_times_i:
1589 case Op_times:
1590 case Op_exp_i:
1591 case Op_exp:
1592 case Op_quotient_i:
1593 case Op_quotient:
1594 case Op_mod_i:
1595 case Op_mod:
1596 case Op_preincrement:
1597 case Op_predecrement:
1598 case Op_postincrement:
1599 case Op_postdecrement:
1600 case Op_unary_minus:
1601 case Op_unary_plus:
1602 case Op_assign_plus:
1603 case Op_assign_minus:
1604 case Op_assign_times:
1605 case Op_assign_quotient:
1606 case Op_assign_mod:
1607 case Op_assign_exp:
1608 fprintf(stderr, "++ %s: mpg_interpret\n", opcode2str(op));
1609 fflush(stderr);
1610 break;
1611 default:
1612 return true; /* unhandled */
1613 }
1614 }
1615
1616 switch (op) {
1617 case Op_plus_i:
1618 t2 = force_number(pc->memory);
1619 goto plus;
1620 case Op_plus:
1621 t2 = POP_NUMBER();
1622 plus:
1623 t1 = TOP_NUMBER();
1624 r = mpg_add(t1, t2);
1625 DEREF(t1);
1626 if (op == Op_plus)
1627 DEREF(t2);
1628 REPLACE(r);
1629 break;
1630
1631 case Op_minus_i:
1632 t2 = force_number(pc->memory);
1633 goto minus;
1634 case Op_minus:
1635 t2 = POP_NUMBER();
1636 minus:
1637 t1 = TOP_NUMBER();
1638 r = mpg_sub(t1, t2);
1639 DEREF(t1);
1640 if (op == Op_minus)
1641 DEREF(t2);
1642 REPLACE(r);
1643 break;
1644
1645 case Op_times_i:
1646 t2 = force_number(pc->memory);
1647 goto times;
1648 case Op_times:
1649 t2 = POP_NUMBER();
1650 times:
1651 t1 = TOP_NUMBER();
1652 r = mpg_mul(t1, t2);
1653 DEREF(t1);
1654 if (op == Op_times)
1655 DEREF(t2);
1656 REPLACE(r);
1657 break;
1658
1659 case Op_exp_i:
1660 t2 = force_number(pc->memory);
1661 goto exp;
1662 case Op_exp:
1663 t2 = POP_NUMBER();
1664 exp:
1665 t1 = TOP_NUMBER();
1666 r = mpg_pow(t1, t2);
1667 DEREF(t1);
1668 if (op == Op_exp)
1669 DEREF(t2);
1670 REPLACE(r);
1671 break;
1672
1673 case Op_quotient_i:
1674 t2 = force_number(pc->memory);
1675 goto quotient;
1676 case Op_quotient:
1677 t2 = POP_NUMBER();
1678 quotient:
1679 t1 = TOP_NUMBER();
1680 r = mpg_div(t1, t2);
1681 DEREF(t1);
1682 if (op == Op_quotient)
1683 DEREF(t2);
1684 REPLACE(r);
1685 break;
1686
1687 case Op_mod_i:
1688 t2 = force_number(pc->memory);
1689 goto mod;
1690 case Op_mod:
1691 t2 = POP_NUMBER();
1692 mod:
1693 t1 = TOP_NUMBER();
1694 r = mpg_mod(t1, t2);
1695 DEREF(t1);
1696 if (op == Op_mod)
1697 DEREF(t2);
1698 REPLACE(r);
1699 break;
1700
1701 case Op_preincrement:
1702 case Op_predecrement:
1703 lhs = TOP_ADDRESS();
1704 t1 = *lhs;
1705 force_number(t1);
1706
1707 if (is_mpg_integer(t1)) {
1708 if (t1->valref == 1 && t1->flags == (MALLOC|MPZN|NUMCUR|NUMBER))
1709 /* Efficiency hack. Big speed-up (> 30%) in a tight loop */
1710 r = t1;
1711 else
1712 r = *lhs = mpg_integer();
1713 if (op == Op_preincrement)
1714 mpz_add_ui(r->mpg_i, t1->mpg_i, 1);
1715 else
1716 mpz_sub_ui(r->mpg_i, t1->mpg_i, 1);
1717 } else {
1718
1719 /*
1720 * An optimization like the one above is not going to work
1721 * for a floating-point number. With it,
1722 * gawk -M 'BEGIN { PREC=53; i=2^53+0.0; PREC=113; ++i; print i}'
1723 * will output 2^53 instead of 2^53+1.
1724 */
1725
1726 r = *lhs = mpg_float();
1727 tval = mpfr_add_si(r->mpg_numbr, t1->mpg_numbr,
1728 op == Op_preincrement ? 1 : -1,
1729 ROUND_MODE);
1730 IEEE_FMT(r->mpg_numbr, tval);
1731 }
1732 if (r != t1)
1733 unref(t1);
1734 UPREF(r);
1735 REPLACE(r);
1736 break;
1737
1738 case Op_postincrement:
1739 case Op_postdecrement:
1740 lhs = TOP_ADDRESS();
1741 t1 = *lhs;
1742 force_number(t1);
1743
1744 if (is_mpg_integer(t1)) {
1745 r = mpg_integer();
1746 mpz_set(r->mpg_i, t1->mpg_i);
1747 if (t1->valref == 1 && t1->flags == (MALLOC|MPZN|NUMCUR|NUMBER))
1748 /* Efficiency hack. Big speed-up (> 30%) in a tight loop */
1749 t2 = t1;
1750 else
1751 t2 = *lhs = mpg_integer();
1752 if (op == Op_postincrement)
1753 mpz_add_ui(t2->mpg_i, t1->mpg_i, 1);
1754 else
1755 mpz_sub_ui(t2->mpg_i, t1->mpg_i, 1);
1756 } else {
1757 r = mpg_float();
1758 tval = mpfr_set(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1759 IEEE_FMT(r->mpg_numbr, tval);
1760 t2 = *lhs = mpg_float();
1761 tval = mpfr_add_si(t2->mpg_numbr, t1->mpg_numbr,
1762 op == Op_postincrement ? 1 : -1,
1763 ROUND_MODE);
1764 IEEE_FMT(t2->mpg_numbr, tval);
1765 }
1766 if (t2 != t1)
1767 unref(t1);
1768 REPLACE(r);
1769 break;
1770
1771 case Op_unary_minus:
1772 t1 = TOP_NUMBER();
1773 if (is_mpg_float(t1)) {
1774 r = mpg_float();
1775 tval = mpfr_neg(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1776 IEEE_FMT(r->mpg_numbr, tval);
1777 } else {
1778 r = mpg_integer();
1779 mpz_neg(r->mpg_i, t1->mpg_i);
1780 }
1781 DEREF(t1);
1782 REPLACE(r);
1783 break;
1784
1785 case Op_unary_plus:
1786 t1 = TOP_NUMBER();
1787 if (is_mpg_float(t1)) {
1788 r = mpg_float();
1789 tval = mpfr_set(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1790 IEEE_FMT(r->mpg_numbr, tval);
1791 } else {
1792 r = mpg_integer();
1793 mpz_set(r->mpg_i, t1->mpg_i);
1794 }
1795 DEREF(t1);
1796 REPLACE(r);
1797 break;
1798
1799 case Op_assign_plus:
1800 case Op_assign_minus:
1801 case Op_assign_times:
1802 case Op_assign_quotient:
1803 case Op_assign_mod:
1804 case Op_assign_exp:
1805 lhs = POP_ADDRESS();
1806 t1 = *lhs;
1807 force_number(t1);
1808 t2 = TOP_NUMBER();
1809
1810 switch (op) {
1811 case Op_assign_plus:
1812 r = mpg_add(t1, t2);
1813 break;
1814 case Op_assign_minus:
1815 r = mpg_sub(t1, t2);
1816 break;
1817 case Op_assign_times:
1818 r = mpg_mul(t1, t2);
1819 break;
1820 case Op_assign_quotient:
1821 r = mpg_div(t1, t2);
1822 break;
1823 case Op_assign_mod:
1824 r = mpg_mod(t1, t2);
1825 break;
1826 case Op_assign_exp:
1827 r = mpg_pow(t1, t2);
1828 break;
1829 default:
1830 cant_happen();
1831 }
1832
1833 DEREF(t2);
1834 unref(*lhs);
1835 *lhs = r;
1836 UPREF(r);
1837 REPLACE(r);
1838 break;
1839
1840 default:
1841 return true; /* unhandled */
1842 }
1843
1844 *cp = pc->nexti; /* next instruction to execute */
1845 return false;
1846 }
1847
1848
1849 /* mpg_fmt --- output formatted string with special MPFR/GMP conversion specifiers */
1850
1851 const char *
mpg_fmt(const char * mesg,...)1852 mpg_fmt(const char *mesg, ...)
1853 {
1854 static char *tmp = NULL;
1855 int ret;
1856 va_list args;
1857
1858 if (tmp != NULL) {
1859 mpfr_free_str(tmp);
1860 tmp = NULL;
1861 }
1862 va_start(args, mesg);
1863 ret = mpfr_vasprintf(& tmp, mesg, args);
1864 va_end(args);
1865 if (ret >= 0 && tmp != NULL)
1866 return tmp;
1867 return mesg;
1868 }
1869
1870 /* mpfr_unset --- clear out the MPFR values */
1871
1872 void
mpfr_unset(NODE * n)1873 mpfr_unset(NODE *n)
1874 {
1875 if (is_mpg_float(n))
1876 mpfr_clear(n->mpg_numbr);
1877 else if (is_mpg_integer(n))
1878 mpz_clear(n->mpg_i);
1879 }
1880
1881 #else
1882
1883 void
set_PREC()1884 set_PREC()
1885 {
1886 /* dummy function */
1887 }
1888
1889 void
set_ROUNDMODE()1890 set_ROUNDMODE()
1891 {
1892 /* dummy function */
1893 }
1894
1895 void
mpfr_unset(NODE * n)1896 mpfr_unset(NODE *n)
1897 {
1898 /* dummy function */
1899 }
1900 #endif
1901