mpfr.c

   1 /*
   2  * mpfr.c - routines for arbitrary-precision number support in gawk.
   3  */
   4
   5 /*
   6  * Copyright (C) 2012, 2013 the Free Software Foundation, Inc.
   7  *
   8  * This file is part of GAWK, the GNU implementation of the
   9  * AWK Programming Language.
  10  *
  11  * GAWK is free software; you can redistribute it and/or modify
  12  * it under the terms of the GNU General Public License as published by
  13  * the Free Software Foundation; either version 3 of the License, or
  14  * (at your option) any later version.
  15  *
  16  * GAWK is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19  * GNU General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU General Public License
  22  * along with this program; if not, write to the Free Software
  23  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  24  */
  25
  26 #include "awk.h"
  27
  28 #ifdef HAVE_MPFR
  29
  30 #if !defined(MPFR_VERSION_MAJOR) || MPFR_VERSION_MAJOR < 3
  31 typedef mp_exp_t mpfr_exp_t;
  32 #endif
  33
  34 extern NODE **fmt_list;          /* declared in eval.c */
  35
  36 mpz_t mpzval;   /* GMP integer type, used as temporary in few places */
  37 mpz_t MNR;
  38 mpz_t MFNR;
  39 bool do_ieee_fmt;       /* IEEE-754 floating-point emulation */
  40 mpfr_rnd_t ROUND_MODE;
  41
  42 static mpfr_rnd_t get_rnd_mode(const char rmode);
  43 static NODE *mpg_force_number(NODE *n);
  44 static NODE *mpg_make_number(double);
  45 static NODE *mpg_format_val(const char *format, int index, NODE *s);
  46 static int mpg_interpret(INSTRUCTION **cp);
  47
  48 static mpfr_exp_t min_exp = MPFR_EMIN_DEFAULT;
  49 static mpfr_exp_t max_exp = MPFR_EMAX_DEFAULT;
  50
  51 /* temporary MPFR floats used to hold converted GMP integer operands */
  52 static mpfr_t _mpf_t1;
  53 static mpfr_t _mpf_t2;
  54
  55 /*
  56  * PRECISION_MIN is the precision used to initialize _mpf_t1 and _mpf_t2.
  57  * 64 bits should be enough for exact conversion of most integers to floats.
  58  */
  59
  60 #define PRECISION_MIN   64
  61
  62 /* mf = { _mpf_t1, _mpf_t2 } */
  63 static inline mpfr_ptr mpg_tofloat(mpfr_ptr mf, mpz_ptr mz);
  64 /* T = {t1, t2} */
  65 #define MP_FLOAT(T) is_mpg_integer(T) ? mpg_tofloat(_mpf_##T, (T)->mpg_i) : (T)->mpg_numbr
  66
  67
  68 /* init_mpfr --- set up MPFR related variables */
  69
  70 void
  71 init_mpfr(mpfr_prec_t prec, const char *rmode)
  72 {
  73         mpfr_set_default_prec(prec);
  74         ROUND_MODE = get_rnd_mode(rmode[0]);
  75         mpfr_set_default_rounding_mode(ROUND_MODE);
  76         make_number = mpg_make_number;
  77         str2number = mpg_force_number;
  78         format_val = mpg_format_val;
  79         cmp_numbers = mpg_cmp;
  80
  81         mpz_init(MNR);
  82         mpz_init(MFNR);
  83         do_ieee_fmt = false;
  84
  85         mpfr_init2(_mpf_t1, PRECISION_MIN);
  86         mpfr_init2(_mpf_t2, PRECISION_MIN);
  87         mpz_init(mpzval);
  88
  89         register_exec_hook(mpg_interpret, 0);
  90 }
  91
  92 /* mpg_node --- allocate a node to store MPFR float or GMP integer */
  93
  94 NODE *
  95 mpg_node(unsigned int tp)
  96 {
  97         NODE *r;
  98         getnode(r);
  99         r->type = Node_val;
 100
 101         if (tp == MPFN) {
 102                 /* Initialize, set precision to the default precision, and value to NaN */
 103                 mpfr_init(r->mpg_numbr);
 104                 r->flags = MPFN;
 105         } else {
 106                 /* Initialize and set value to 0 */
 107                 mpz_init(r->mpg_i);
 108                 r->flags = MPZN;
 109         }
 110
 111         r->valref = 1;
 112         r->flags |= MALLOC|NUMBER|NUMCUR;
 113         r->stptr = NULL;
 114         r->stlen = 0;
 115 #if MBS_SUPPORT
 116         r->wstptr = NULL;
 117         r->wstlen = 0;
 118 #endif /* defined MBS_SUPPORT */
 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                                 goto done;
 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 static inline 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) {
 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         if (do_nondec)
 297                 base = get_numbase(cp1, use_locale);
 298
 299         if (! mpg_maybe_float(cp1, use_locale)) {
 300                 mpg_zero(n);
 301                 errno = 0;
 302                 mpg_strtoui(n->mpg_i, cp1, cpend - cp1, & ptr, base);
 303                 if (*cp == '-')
 304                         mpz_neg(n->mpg_i, n->mpg_i);
 305                 goto done;
 306         }
 307
 308         if (is_mpg_integer(n)) {
 309                 mpz_clear(n->mpg_i);
 310                 n->flags &= ~MPZN;
 311         }
 312
 313         if (! is_mpg_float(n)) {
 314                 mpfr_init(n->mpg_numbr);
 315                 n->flags |= MPFN;
 316         }
 317
 318         errno = 0;
 319         tval = mpfr_strtofr(n->mpg_numbr, cp, & ptr, base, ROUND_MODE);
 320         IEEE_FMT(n->mpg_numbr, tval);
 321 done:
 322         /* trailing space is OK for NUMBER */
 323         while (isspace((unsigned char) *ptr))
 324                 ptr++;
 325         *cpend = save;
 326         if (errno == 0 && ptr == cpend)
 327                 return true;
 328         errno = 0;
 329         return false;
 330 }
 331
 332 /* mpg_force_number --- force a value to be a multiple-precision number */
 333
 334 static NODE *
 335 mpg_force_number(NODE *n)
 336 {
 337         unsigned int newflags = 0;
 338
 339         if (is_mpg_number(n) && (n->flags & NUMCUR) != 0)
 340                 return n;
 341
 342         if ((n->flags & MAYBE_NUM) != 0) {
 343                 n->flags &= ~MAYBE_NUM;
 344                 newflags = NUMBER;
 345         }
 346
 347         if (force_mpnum(n, (do_non_decimal_data && ! do_traditional), true)) {
 348                 n->flags |= newflags;
 349                 n->flags |= NUMCUR;
 350         }
 351         return n;
 352 }
 353
 354 /* mpg_format_val --- format a numeric value based on format */
 355
 356 static NODE *
 357 mpg_format_val(const char *format, int index, NODE *s)
 358 {
 359         NODE *dummy[2], *r;
 360         unsigned int oflags;
 361
 362         /* create dummy node for a sole use of format_tree */
 363         dummy[1] = s;
 364         oflags = s->flags;
 365
 366         if (is_mpg_integer(s) || mpfr_integer_p(s->mpg_numbr)) {
 367                 /* integral value, use %d */
 368                 r = format_tree("%d", 2, dummy, 2);
 369                 s->stfmt = -1;
 370         } else {
 371                 r = format_tree(format, fmt_list[index]->stlen, dummy, 2);
 372                 assert(r != NULL);
 373                 s->stfmt = (char) index;
 374         }
 375         s->flags = oflags;
 376         s->stlen = r->stlen;
 377         if ((s->flags & STRCUR) != 0)
 378                 efree(s->stptr);
 379         s->stptr = r->stptr;
 380         freenode(r);    /* Do not unref(r)! We want to keep s->stptr == r->stpr.  */
 381
 382         s->flags |= STRCUR;
 383         free_wstr(s);
 384         return s;
 385 }
 386
 387 /* mpg_cmp --- compare two numbers */
 388
 389 int
 390 mpg_cmp(const NODE *t1, const NODE *t2)
 391 {
 392         /*
 393          * For the purposes of sorting, NaN is considered greater than
 394          * any other value, and all NaN values are considered equivalent and equal.
 395          */
 396
 397         if (is_mpg_float(t1)) {
 398                 if (is_mpg_float(t2)) {
 399                         if (mpfr_nan_p(t1->mpg_numbr))
 400                                 return ! mpfr_nan_p(t2->mpg_numbr);
 401                         if (mpfr_nan_p(t2->mpg_numbr))
 402                                 return -1;
 403                         return mpfr_cmp(t1->mpg_numbr, t2->mpg_numbr);
 404                 }
 405                 if (mpfr_nan_p(t1->mpg_numbr))
 406                         return 1;
 407                 return mpfr_cmp_z(t1->mpg_numbr, t2->mpg_i);
 408         } else if (is_mpg_float(t2)) {
 409                 int ret;
 410                 if (mpfr_nan_p(t2->mpg_numbr))
 411                         return -1;
 412                 ret = mpfr_cmp_z(t2->mpg_numbr, t1->mpg_i);
 413                 return ret > 0 ? -1 : (ret < 0);
 414         } else if (is_mpg_integer(t1)) {
 415                 return mpz_cmp(t1->mpg_i, t2->mpg_i);
 416         }
 417
 418         /* t1 and t2 are AWKNUMs */
 419         return cmp_awknums(t1, t2);
 420 }
 421
 422
 423 /*
 424  * mpg_update_var --- update NR or FNR.
 425  *      NR_node->var_value(mpz_t) = MNR(mpz_t) * LONG_MAX + NR(long)
 426  */
 427
 428 NODE *
 429 mpg_update_var(NODE *n)
 430 {
 431         NODE *val = n->var_value;
 432         long nr = 0;
 433         mpz_ptr nq = 0;
 434
 435         if (n == NR_node) {
 436                 nr = NR;
 437                 nq = MNR;
 438         } else if (n == FNR_node) {
 439                 nr = FNR;
 440                 nq = MFNR;
 441         } else
 442                 cant_happen();
 443
 444         if (mpz_sgn(nq) == 0) {
 445                 /* Efficiency hack similar to that for AWKNUM */
 446                 if (is_mpg_float(val) || mpz_get_si(val->mpg_i) != nr) {
 447                         unref(n->var_value);
 448                         val = n->var_value = mpg_integer();
 449                         mpz_set_si(val->mpg_i, nr);
 450                 }
 451         } else {
 452                 unref(n->var_value);
 453                 val = n->var_value = mpg_integer();
 454                 mpz_set_si(val->mpg_i, nr);
 455                 mpz_addmul_ui(val->mpg_i, nq, LONG_MAX);        /* val->mpg_i += nq * LONG_MAX */
 456         }
 457         return val;
 458 }
 459
 460 /* mpg_set_var --- set NR or FNR */
 461
 462 long
 463 mpg_set_var(NODE *n)
 464 {
 465         long nr = 0;
 466         mpz_ptr nq = 0, r;
 467         NODE *val = n->var_value;
 468
 469         if (n == NR_node)
 470                 nq = MNR;
 471         else if (n == FNR_node)
 472                 nq = MFNR;
 473         else
 474                 cant_happen();
 475
 476         if (is_mpg_integer(val))
 477                 r = val->mpg_i;
 478         else {
 479                 /* convert float to integer */
 480                 mpfr_get_z(mpzval, val->mpg_numbr, MPFR_RNDZ);
 481                 r = mpzval;
 482         }
 483         nr = mpz_fdiv_q_ui(nq, r, LONG_MAX);    /* nq (MNR or MFNR) is quotient */
 484         return nr;      /* remainder (NR or FNR) */
 485 }
 486
 487 /* set_PREC --- update MPFR PRECISION related variables when PREC assigned to */
 488
 489 void
 490 set_PREC()
 491 {
 492         long prec = 0;
 493         NODE *val;
 494         static const struct ieee_fmt {
 495                 const char *name;
 496                 mpfr_prec_t precision;
 497                 mpfr_exp_t emax;
 498                 mpfr_exp_t emin;
 499         } ieee_fmts[] = {
 500 { "half",       11,     16,     -23     },      /* binary16 */
 501 { "single",     24,     128,    -148    },      /* binary32 */
 502 { "double",     53,     1024,   -1073   },      /* binary64 */
 503 { "quad",       113,    16384,  -16493  },      /* binary128 */
 504 { "oct",        237,    262144, -262377 },      /* binary256, not in the IEEE 754-2008 standard */
 505
 506                 /*
 507                  * For any bitwidth = 32 * k ( k >= 4),
 508                  * precision = 13 + bitwidth - int(4 * log2(bitwidth))
 509                  * emax = 1 << bitwidth - precision - 1
 510                  * emin = 4 - emax - precision
 511                  */
 512         };
 513
 514         if (! do_mpfr)
 515                 return;
 516
 517         val = PREC_node->var_value;
 518         if ((val->flags & MAYBE_NUM) != 0)
 519                 force_number(val);
 520
 521         if ((val->flags & (STRING|NUMBER)) == STRING) {
 522                 int i, j;
 523
 524                 /* emulate IEEE-754 binary format */
 525
 526                 for (i = 0, j = sizeof(ieee_fmts)/sizeof(ieee_fmts[0]); i < j; i++) {
 527                         if (strcasecmp(ieee_fmts[i].name, val->stptr) == 0)
 528                                 break;
 529                 }
 530
 531                 if (i < j) {
 532                         prec = ieee_fmts[i].precision;
 533
 534                         /*
 535                          * We *DO NOT* change the MPFR exponent range using
 536                          * mpfr_set_{emin, emax} here. See format_ieee() for details.
 537                          */
 538                         max_exp = ieee_fmts[i].emax;
 539                         min_exp = ieee_fmts[i].emin;
 540
 541                         do_ieee_fmt = true;
 542                 }
 543         }
 544
 545         if (prec <= 0) {
 546                 force_number(val);
 547                 prec = get_number_si(val);
 548                 if (prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX) {
 549                         force_string(val);
 550                         warning(_("PREC value `%.*s' is invalid"), (int) val->stlen, val->stptr);
 551                         prec = 0;
 552                 } else
 553                         do_ieee_fmt = false;
 554         }
 555
 556         if (prec > 0)
 557                 mpfr_set_default_prec(prec);
 558 }
 559
 560
 561 /* get_rnd_mode --- convert string to MPFR rounding mode */
 562
 563 static mpfr_rnd_t
 564 get_rnd_mode(const char rmode)
 565 {
 566         switch (rmode) {
 567         case 'N':
 568         case 'n':
 569                 return MPFR_RNDN;       /* round to nearest (IEEE-754 roundTiesToEven) */
 570         case 'Z':
 571         case 'z':
 572                 return MPFR_RNDZ;       /* round toward zero (IEEE-754 roundTowardZero) */
 573         case 'U':
 574         case 'u':
 575                 return MPFR_RNDU;       /* round toward plus infinity (IEEE-754 roundTowardPositive) */
 576         case 'D':
 577         case 'd':
 578                 return MPFR_RNDD;       /* round toward minus infinity (IEEE-754 roundTowardNegative) */
 579 #if defined(MPFR_VERSION_MAJOR) && MPFR_VERSION_MAJOR > 2
 580         case 'A':
 581         case 'a':
 582                 return MPFR_RNDA;       /* round away from zero (IEEE-754 roundTiesToAway) */
 583 #endif
 584         default:
 585                 break;
 586         }
 587         return -1;
 588 }
 589
 590 /*
 591  * set_ROUNDMODE --- update MPFR rounding mode related variables
 592  *      when ROUNDMODE assigned to
 593  */
 594
 595 void
 596 set_ROUNDMODE()
 597 {
 598         if (do_mpfr) {
 599                 mpfr_rnd_t rndm = -1;
 600                 NODE *n;
 601                 n = force_string(ROUNDMODE_node->var_value);
 602                 if (n->stlen == 1)
 603                         rndm = get_rnd_mode(n->stptr[0]);
 604                 if (rndm != -1) {
 605                         mpfr_set_default_rounding_mode(rndm);
 606                         ROUND_MODE = rndm;
 607                 } else
 608                         warning(_("RNDMODE value `%.*s' is invalid"), (int) n->stlen, n->stptr);
 609         }
 610 }
 611
 612
 613 /* format_ieee --- make sure a number follows IEEE-754 floating-point standard */
 614
 615 int
 616 format_ieee(mpfr_ptr x, int tval)
 617 {
 618         /*
 619          * The MPFR doc says that it's our responsibility to make sure all numbers
 620          * including those previously created are in range after we've changed the
 621          * exponent range. Most MPFR operations and functions require
 622          * the input arguments to have exponents within the current exponent range.
 623          * Any argument outside the range results in a MPFR assertion failure
 624          * like this:
 625          *
 626          *   $ gawk -M 'BEGIN { x=1.0e-10000; print x+0; PREC="double"; print x+0}'
 627          *   1e-10000
 628          *   init2.c:52: MPFR assertion failed ....
 629          *
 630          * A "naive" approach would be to keep track of the ternary state and
 631          * the rounding mode for each number, and make sure it is in the current
 632          * exponent range (using mpfr_check_range) before using it in an
 633          * operation or function. Instead, we adopt the following strategy.
 634          *
 635          * When gawk starts, the exponent range is the MPFR default
 636          * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT]. Any number that gawk
 637          * creates must have exponent in this range (excluding infinities, NaNs and zeros).
 638          * Each MPFR operation or function is performed with this default exponent
 639          * range.
 640          *
 641          * When emulating IEEE-754 format, the exponents are *temporarily* changed,
 642          * mpfr_check_range is called to make sure the number is in the new range,
 643          * and mpfr_subnormalize is used to round following the rules of subnormal
 644          * arithmetic. The exponent range is then *restored* to the original value
 645          * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT].
 646          */
 647
 648         (void) mpfr_set_emin(min_exp);
 649         (void) mpfr_set_emax(max_exp);
 650         tval = mpfr_check_range(x, tval, ROUND_MODE);
 651         tval = mpfr_subnormalize(x, tval, ROUND_MODE);
 652         (void) mpfr_set_emin(MPFR_EMIN_DEFAULT);
 653         (void) mpfr_set_emax(MPFR_EMAX_DEFAULT);
 654         return tval;
 655 }
 656
 657
 658 /* do_mpfr_atan2 --- do the atan2 function */
 659
 660 NODE *
 661 do_mpfr_atan2(int nargs)
 662 {
 663         NODE *t1, *t2, *res;
 664         mpfr_ptr p1, p2;
 665         int tval;
 666
 667         t2 = POP_SCALAR();
 668         t1 = POP_SCALAR();
 669
 670         if (do_lint) {
 671                 if ((t1->flags & (NUMCUR|NUMBER)) == 0)
 672                         lintwarn(_("atan2: received non-numeric first argument"));
 673                 if ((t2->flags & (NUMCUR|NUMBER)) == 0)
 674                         lintwarn(_("atan2: received non-numeric second argument"));
 675         }
 676         force_number(t1);
 677         force_number(t2);
 678
 679         p1 = MP_FLOAT(t1);
 680         p2 = MP_FLOAT(t2);
 681         res = mpg_float();
 682         /* See MPFR documentation for handling of special values like +inf as an argument */
 683         tval = mpfr_atan2(res->mpg_numbr, p1, p2, ROUND_MODE);
 684         IEEE_FMT(res->mpg_numbr, tval);
 685
 686         DEREF(t1);
 687         DEREF(t2);
 688         return res;
 689 }
 690
 691
 692 #define SPEC_MATH(X)                                            \
 693 NODE *t1, *res;                                                 \
 694 mpfr_ptr p1;                                                    \
 695 int tval;                                                       \
 696 t1 = POP_SCALAR();                                              \
 697 if (do_lint && (t1->flags & (NUMCUR|NUMBER)) == 0)              \
 698         lintwarn(_("%s: received non-numeric argument"), #X);   \
 699 force_number(t1);                                               \
 700 p1 = MP_FLOAT(t1);                                              \
 701 res = mpg_float();                                              \
 702 tval = mpfr_##X(res->mpg_numbr, p1, ROUND_MODE);                        \
 703 IEEE_FMT(res->mpg_numbr, tval);                                 \
 704 DEREF(t1);                                                      \
 705 return res
 706
 707
 708 /* do_mpfr_sin --- do the sin function */
 709
 710 NODE *
 711 do_mpfr_sin(int nargs)
 712 {
 713         SPEC_MATH(sin);
 714 }
 715
 716 /* do_mpfr_cos --- do the cos function */
 717
 718 NODE *
 719 do_mpfr_cos(int nargs)
 720 {
 721         SPEC_MATH(cos);
 722 }
 723
 724 /* do_mpfr_exp --- exponential function */
 725
 726 NODE *
 727 do_mpfr_exp(int nargs)
 728 {
 729         SPEC_MATH(exp);
 730 }
 731
 732 /* do_mpfr_log --- the log function */
 733
 734 NODE *
 735 do_mpfr_log(int nargs)
 736 {
 737         SPEC_MATH(log);
 738 }
 739
 740 /* do_mpfr_sqrt --- do the sqrt function */
 741
 742 NODE *
 743 do_mpfr_sqrt(int nargs)
 744 {
 745         SPEC_MATH(sqrt);
 746 }
 747
 748 /* do_mpfr_int --- convert double to int for awk */
 749
 750 NODE *
 751 do_mpfr_int(int nargs)
 752 {
 753         NODE *tmp, *r;
 754
 755         tmp = POP_SCALAR();
 756         if (do_lint && (tmp->flags & (NUMCUR|NUMBER)) == 0)
 757                 lintwarn(_("int: received non-numeric argument"));
 758         force_number(tmp);
 759
 760         if (is_mpg_integer(tmp)) {
 761                 r = mpg_integer();
 762                 mpz_set(r->mpg_i, tmp->mpg_i);
 763         } else {
 764                 if (! mpfr_number_p(tmp->mpg_numbr)) {
 765                         /* [+-]inf or NaN */
 766                         return tmp;
 767                 }
 768
 769                 r = mpg_integer();
 770                 mpfr_get_z(r->mpg_i, tmp->mpg_numbr, MPFR_RNDZ);
 771         }
 772
 773         DEREF(tmp);
 774         return r;
 775 }
 776
 777 /* do_mpfr_compl --- perform a ~ operation */
 778
 779 NODE *
 780 do_mpfr_compl(int nargs)
 781 {
 782         NODE *tmp, *r;
 783         mpz_ptr zptr;
 784
 785         tmp = POP_SCALAR();
 786         if (do_lint && (tmp->flags & (NUMCUR|NUMBER)) == 0)
 787                 lintwarn(_("compl: received non-numeric argument"));
 788
 789         force_number(tmp);
 790         if (is_mpg_float(tmp)) {
 791                 mpfr_ptr p = tmp->mpg_numbr;
 792
 793                 if (! mpfr_number_p(p)) {
 794                         /* [+-]inf or NaN */
 795                         return tmp;
 796                 }
 797                 if (do_lint) {
 798                         if (mpfr_sgn(p) < 0)
 799                                 lintwarn("%s",
 800                         mpg_fmt(_("compl(%Rg): negative value will give strange results"), p)
 801                                         );
 802                         if (! mpfr_integer_p(p))
 803                                 lintwarn("%s",
 804                         mpg_fmt(_("comp(%Rg): fractional value will be truncated"), p)
 805                                         );
 806                 }
 807
 808                 mpfr_get_z(mpzval, p, MPFR_RNDZ);       /* float to integer conversion */
 809                 zptr = mpzval;
 810         } else {
 811                 /* (tmp->flags & MPZN) != 0 */
 812                 zptr = tmp->mpg_i;
 813                 if (do_lint) {
 814                         if (mpz_sgn(zptr) < 0)
 815                                 lintwarn("%s",
 816                         mpg_fmt(_("cmpl(%Zd): negative values will give strange results"), zptr)
 817                                         );
 818                 }
 819         }
 820
 821         r = mpg_integer();
 822         mpz_com(r->mpg_i, zptr);
 823         DEREF(tmp);
 824         return r;
 825 }
 826
 827 /* get_intval --- get the (converted) integral operand of a binary function. */
 828
 829 static mpz_ptr
 830 get_intval(NODE *t1, int argnum, const char *op)
 831 {
 832         mpz_ptr pz;
 833
 834         if (do_lint && (t1->flags & (NUMCUR|NUMBER)) == 0)
 835                 lintwarn(_("%s: received non-numeric argument #%d"), op, argnum);
 836
 837         (void) force_number(t1);
 838
 839         if (is_mpg_float(t1)) {
 840                 mpfr_ptr left = t1->mpg_numbr;
 841                 if (! mpfr_number_p(left)) {
 842                         /* inf or NaN */
 843                         if (do_lint)
 844                                 lintwarn("%s",
 845                 mpg_fmt(_("%s: argument #%d has invalid value %Rg, using 0"),
 846                                         op, argnum, left)
 847                                 );
 848
 849                         emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
 850                         mpz_init(pz);
 851                         return pz;      /* should be freed */
 852                 }
 853
 854                 if (do_lint) {
 855                         if (mpfr_sgn(left) < 0)
 856                                 lintwarn("%s",
 857                 mpg_fmt(_("%s: argument #%d negative value %Rg will give strange results"),
 858                                         op, argnum, left)
 859                                 );
 860
 861                         if (! mpfr_integer_p(left))
 862                                 lintwarn("%s",
 863                 mpg_fmt(_("%s: argument #%d fractional value %Rg will be truncated"),
 864                                         op, argnum, left)
 865                                 );
 866                 }
 867
 868                 emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
 869                 mpz_init(pz);
 870                 mpfr_get_z(pz, left, MPFR_RNDZ);        /* float to integer conversion */
 871                 return pz;      /* should be freed */
 872         }
 873         /* (t1->flags & MPZN) != 0 */
 874         pz = t1->mpg_i;
 875         if (do_lint) {
 876                 if (mpz_sgn(pz) < 0)
 877                         lintwarn("%s",
 878         mpg_fmt(_("%s: argument #%d negative value %Zd will give strange results"),
 879                                         op, argnum, pz)
 880                                 );
 881         }
 882         return pz;      /* must not be freed */
 883 }
 884
 885
 886 /* free_intval --- free the converted integer value returned by get_intval() */
 887
 888 static inline void
 889 free_intval(NODE *t, mpz_ptr pz)
 890 {
 891         if ((t->flags & MPZN) == 0) {
 892                 mpz_clear(pz);
 893                 efree(pz);
 894         }
 895 }
 896
 897
 898 /* do_mpfr_lshift --- perform a << operation */
 899
 900 NODE *
 901 do_mpfr_lshift(int nargs)
 902 {
 903         NODE *t1, *t2, *res;
 904         unsigned long shift;
 905         mpz_ptr pz1, pz2;
 906
 907         t2 = POP_SCALAR();
 908         t1 = POP_SCALAR();
 909
 910         pz1 = get_intval(t1, 1, "lshift");
 911         pz2 = get_intval(t2, 2, "lshift");
 912
 913         /*
 914          * mpz_get_ui: If op is too big to fit an unsigned long then just
 915          * the least significant bits that do fit are returned.
 916          * The sign of op is ignored, only the absolute value is used.
 917          */
 918
 919         shift = mpz_get_ui(pz2);        /* GMP integer => unsigned long conversion */
 920         res = mpg_integer();
 921         mpz_mul_2exp(res->mpg_i, pz1, shift);           /* res = pz1 * 2^shift */
 922
 923         free_intval(t1, pz1);
 924         free_intval(t2, pz2);
 925         DEREF(t2);
 926         DEREF(t1);
 927         return res;
 928 }
 929
 930 /* do_mpfr_rshift --- perform a >> operation */
 931
 932 NODE *
 933 do_mpfr_rshift(int nargs)
 934 {
 935         NODE *t1, *t2, *res;
 936         unsigned long shift;
 937         mpz_ptr pz1, pz2;
 938
 939         t2 = POP_SCALAR();
 940         t1 = POP_SCALAR();
 941
 942         pz1 = get_intval(t1, 1, "rshift");
 943         pz2 = get_intval(t2, 2, "rshift");
 944
 945         /* N.B: See do_mpfp_lshift. */
 946         shift = mpz_get_ui(pz2);        /* GMP integer => unsigned long conversion */
 947         res = mpg_integer();
 948         mpz_fdiv_q_2exp(res->mpg_i, pz1, shift);        /* res = pz1 / 2^shift, round towards −inf */
 949
 950         free_intval(t1, pz1);
 951         free_intval(t2, pz2);
 952         DEREF(t2);
 953         DEREF(t1);
 954         return res;
 955 }
 956
 957
 958 /* do_mpfr_and --- perform an & operation */
 959
 960 NODE *
 961 do_mpfr_and(int nargs)
 962 {
 963         NODE *t1, *t2, *res;
 964         mpz_ptr pz1, pz2;
 965         int i;
 966
 967         if (nargs < 2)
 968                 fatal(_("and: called with less than two arguments"));
 969
 970         t2 = POP_SCALAR();
 971         pz2 = get_intval(t2, nargs, "and");
 972
 973         res = mpg_integer();
 974         for (i = 1; i < nargs; i++) {
 975                 t1 = POP_SCALAR();
 976                 pz1 = get_intval(t1, nargs - i, "and");
 977                 mpz_and(res->mpg_i, pz1, pz2);
 978                 free_intval(t1, pz1);
 979                 DEREF(t1);
 980                 if (i == 1) {
 981                         free_intval(t2, pz2);
 982                         DEREF(t2);
 983                 }
 984                 pz2 = res->mpg_i;
 985         }
 986         return res;
 987 }
 988
 989
 990 /* do_mpfr_or --- perform an | operation */
 991
 992 NODE *
 993 do_mpfr_or(int nargs)
 994 {
 995         NODE *t1, *t2, *res;
 996         mpz_ptr pz1, pz2;
 997         int i;
 998
 999         if (nargs < 2)
1000                 fatal(_("or: called with less than two arguments"));
1001
1002         t2 = POP_SCALAR();
1003         pz2 = get_intval(t2, nargs, "or");
1004
1005         res = mpg_integer();
1006         for (i = 1; i < nargs; i++) {
1007                 t1 = POP_SCALAR();
1008                 pz1 = get_intval(t1, nargs - i, "or");
1009                 mpz_ior(res->mpg_i, pz1, pz2);
1010                 free_intval(t1, pz1);
1011                 DEREF(t1);
1012                 if (i == 1) {
1013                         free_intval(t2, pz2);
1014                         DEREF(t2);
1015                 }
1016                 pz2 = res->mpg_i;
1017         }
1018         return res;
1019 }
1020
1021 /* do_mpfr_xor --- perform an ^ operation */
1022
1023 NODE *
1024 do_mpfr_xor(int nargs)
1025 {
1026         NODE *t1, *t2, *res;
1027         mpz_ptr pz1, pz2;
1028         int i;
1029
1030         if (nargs < 2)
1031                 fatal(_("xor: called with less than two arguments"));
1032
1033         t2 = POP_SCALAR();
1034         pz2 = get_intval(t2, nargs, "xor");
1035
1036         res = mpg_integer();
1037         for (i = 1; i < nargs; i++) {
1038                 t1 = POP_SCALAR();
1039                 pz1 = get_intval(t1, nargs - i, "xor");
1040                 mpz_xor(res->mpg_i, pz1, pz2);
1041                 free_intval(t1, pz1);
1042                 DEREF(t1);
1043                 if (i == 1) {
1044                         free_intval(t2, pz2);
1045                         DEREF(t2);
1046                 }
1047                 pz2 = res->mpg_i;
1048         }
1049         return res;
1050 }
1051
1052 /* do_mpfr_strtonum --- the strtonum function */
1053
1054 NODE *
1055 do_mpfr_strtonum(int nargs)
1056 {
1057         NODE *tmp, *r;
1058
1059         tmp = POP_SCALAR();
1060         if ((tmp->flags & (NUMBER|NUMCUR)) == 0) {
1061                 r = mpg_integer();      /* will be changed to MPFR float if necessary in force_mpnum() */
1062                 r->stptr = tmp->stptr;
1063                 r->stlen = tmp->stlen;
1064                 force_mpnum(r, true, use_lc_numeric);
1065                 r->stptr = NULL;
1066                 r->stlen = 0;
1067         } else {
1068                 (void) force_number(tmp);
1069                 if (is_mpg_float(tmp)) {
1070                         int tval;
1071                         r = mpg_float();
1072                         tval = mpfr_set(r->mpg_numbr, tmp->mpg_numbr, ROUND_MODE);
1073                         IEEE_FMT(r->mpg_numbr, tval);
1074                 } else {
1075                         r = mpg_integer();
1076                         mpz_set(r->mpg_i, tmp->mpg_i);
1077                 }
1078         }
1079
1080         DEREF(tmp);
1081         return r;
1082 }
1083
1084
1085 static bool firstrand = true;
1086 static gmp_randstate_t state;
1087 static mpz_t seed;      /* current seed */
1088
1089 /* do_mpfr_rand --- do the rand function */
1090
1091 NODE *
1092 do_mpfr_rand(int nargs ATTRIBUTE_UNUSED)
1093 {
1094         NODE *res;
1095         int tval;
1096
1097         if (firstrand) {
1098 #if 0
1099                 /* Choose the default algorithm */
1100                 gmp_randinit_default(state);
1101 #endif
1102                 /*
1103                  * Choose a specific (Mersenne Twister) algorithm in case the default
1104                  * changes in the future.
1105                  */
1106
1107                 gmp_randinit_mt(state);
1108
1109                 mpz_init(seed);
1110                 mpz_set_ui(seed, 1);
1111                 /* seed state */
1112                 gmp_randseed(state, seed);
1113                 firstrand = false;
1114         }
1115         res = mpg_float();
1116         tval = mpfr_urandomb(res->mpg_numbr, state);
1117         IEEE_FMT(res->mpg_numbr, tval);
1118         return res;
1119 }
1120
1121
1122 /* do_mpfr_srand --- seed the random number generator */
1123
1124 NODE *
1125 do_mpfr_srand(int nargs)
1126 {
1127         NODE *res;
1128
1129         if (firstrand) {
1130 #if 0
1131                 /* Choose the default algorithm */
1132                 gmp_randinit_default(state);
1133 #endif
1134                 /*
1135                  * Choose a specific algorithm (Mersenne Twister) in case default
1136                  * changes in the future.
1137                  */
1138
1139                 gmp_randinit_mt(state);
1140
1141                 mpz_init(seed);
1142                 mpz_set_ui(seed, 1);
1143                 /* No need to seed state, will change it below */
1144                 firstrand = false;
1145         }
1146
1147         res = mpg_integer();
1148         mpz_set(res->mpg_i, seed);      /* previous seed */
1149
1150         if (nargs == 0)
1151                 mpz_set_ui(seed, (unsigned long) time((time_t *) 0));
1152         else {
1153                 NODE *tmp;
1154                 tmp = POP_SCALAR();
1155                 if (do_lint && (tmp->flags & (NUMCUR|NUMBER)) == 0)
1156                         lintwarn(_("srand: received non-numeric argument"));
1157                 force_number(tmp);
1158                 if (is_mpg_float(tmp))
1159                         mpfr_get_z(seed, tmp->mpg_numbr, MPFR_RNDZ);
1160                 else /* MP integer */
1161                         mpz_set(seed, tmp->mpg_i);
1162                 DEREF(tmp);
1163         }
1164
1165         gmp_randseed(state, seed);
1166         return res;
1167 }
1168
1169 /*
1170  * mpg_tofloat --- convert an arbitrary-precision integer operand to
1171  *      a float without loss of precision. It is assumed that the
1172  *      MPFR variable has already been initialized.
1173  */
1174
1175 static inline mpfr_ptr
1176 mpg_tofloat(mpfr_ptr mf, mpz_ptr mz)
1177 {
1178         size_t prec;
1179
1180         /*
1181          * When implicitely converting a GMP integer operand to a MPFR float, use
1182          * a precision sufficiently large to hold the converted value exactly.
1183          *
1184          *      $ ./gawk -M 'BEGIN { print 13 % 2 }'
1185          *      1
1186          * If the user-specified precision is used to convert the integer 13 to a
1187          * float, one will get:
1188          *      $ ./gawk -M 'BEGIN { PREC=2; print 13 % 2.0 }'
1189          *      0
1190          */
1191
1192         prec = mpz_sizeinbase(mz, 2);   /* most significant 1 bit position starting at 1 */
1193         if (prec > PRECISION_MIN) {
1194                 prec -= (size_t) mpz_scan1(mz, 0);      /* least significant 1 bit index starting at 0 */
1195                 if (prec > MPFR_PREC_MAX)
1196                         prec = MPFR_PREC_MAX;
1197                 if (prec > PRECISION_MIN)
1198                         mpfr_set_prec(mf, prec);
1199         }
1200
1201         mpfr_set_z(mf, mz, ROUND_MODE);
1202         return mf;
1203 }
1204
1205
1206 /* mpg_add --- add arbitrary-precision numbers */
1207
1208 static NODE *
1209 mpg_add(NODE *t1, NODE *t2)
1210 {
1211         NODE *r;
1212         int tval;
1213
1214         if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1215                 r = mpg_integer();
1216                 mpz_add(r->mpg_i, t1->mpg_i, t2->mpg_i);
1217         } else {
1218                 r = mpg_float();
1219                 if (is_mpg_integer(t2))
1220                         tval = mpfr_add_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1221                 else if (is_mpg_integer(t1))
1222                         tval = mpfr_add_z(r->mpg_numbr, t2->mpg_numbr, t1->mpg_i, ROUND_MODE);
1223                 else
1224                         tval = mpfr_add(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1225                 IEEE_FMT(r->mpg_numbr, tval);
1226         }
1227         return r;
1228 }
1229
1230 /* mpg_sub --- subtract arbitrary-precision numbers */
1231
1232 static NODE *
1233 mpg_sub(NODE *t1, NODE *t2)
1234 {
1235         NODE *r;
1236         int tval;
1237
1238         if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1239                 r = mpg_integer();
1240                 mpz_sub(r->mpg_i, t1->mpg_i, t2->mpg_i);
1241         } else {
1242                 r = mpg_float();
1243                 if (is_mpg_integer(t2))
1244                         tval = mpfr_sub_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1245                 else if (is_mpg_integer(t1)) {
1246 #if (!defined(MPFR_VERSION) || (MPFR_VERSION < MPFR_VERSION_NUM(3,1,0)))
1247                         NODE *tmp = t1;
1248                         t1 = t2;
1249                         t2 = tmp;
1250                         tval = mpfr_sub_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1251                         tval = mpfr_neg(r->mpg_numbr, r->mpg_numbr, ROUND_MODE);
1252                         t2 = t1;
1253                         t1 = tmp;
1254 #else
1255                         tval = mpfr_z_sub(r->mpg_numbr, t1->mpg_i, t2->mpg_numbr, ROUND_MODE);
1256 #endif
1257                 } else
1258                         tval = mpfr_sub(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1259                 IEEE_FMT(r->mpg_numbr, tval);
1260         }
1261         return r;
1262 }
1263
1264 /* mpg_mul --- multiply arbitrary-precision numbers */
1265
1266 static NODE *
1267 mpg_mul(NODE *t1, NODE *t2)
1268 {
1269         NODE *r;
1270         int tval;
1271
1272         if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1273                 r = mpg_integer();
1274                 mpz_mul(r->mpg_i, t1->mpg_i, t2->mpg_i);
1275         } else {
1276                 r = mpg_float();
1277                 if (is_mpg_integer(t2))
1278                         tval = mpfr_mul_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1279                 else if (is_mpg_integer(t1))
1280                         tval = mpfr_mul_z(r->mpg_numbr, t2->mpg_numbr, t1->mpg_i, ROUND_MODE);
1281                 else
1282                         tval = mpfr_mul(r->mpg_numbr, t1->mpg_numbr, t2->mpg_numbr, ROUND_MODE);
1283                 IEEE_FMT(r->mpg_numbr, tval);
1284         }
1285         return r;
1286 }
1287
1288
1289 /* mpg_pow --- exponentiation involving arbitrary-precision numbers */
1290
1291 static NODE *
1292 mpg_pow(NODE *t1, NODE *t2)
1293 {
1294         NODE *r;
1295         int tval;
1296
1297         if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1298                 if (mpz_sgn(t2->mpg_i) >= 0 && mpz_fits_ulong_p(t2->mpg_i)) {
1299                         r = mpg_integer();
1300                         mpz_pow_ui(r->mpg_i, t1->mpg_i, mpz_get_ui(t2->mpg_i));
1301                 } else {
1302                         mpfr_ptr p1, p2;
1303                         p1 = MP_FLOAT(t1);
1304                         p2 = MP_FLOAT(t2);
1305                         r = mpg_float();
1306                         tval = mpfr_pow(r->mpg_numbr, p1, p2, ROUND_MODE);
1307                         IEEE_FMT(r->mpg_numbr, tval);
1308                 }
1309         } else {
1310                 r = mpg_float();
1311                 if (is_mpg_integer(t2))
1312                         tval = mpfr_pow_z(r->mpg_numbr, t1->mpg_numbr, t2->mpg_i, ROUND_MODE);
1313                 else {
1314                         mpfr_ptr p1;
1315                         p1 = MP_FLOAT(t1);
1316                         tval = mpfr_pow(r->mpg_numbr, p1, t2->mpg_numbr, ROUND_MODE);
1317                 }
1318                 IEEE_FMT(r->mpg_numbr, tval);
1319         }
1320         return r;
1321 }
1322
1323 /* mpg_div --- arbitrary-precision division */
1324
1325 static NODE *
1326 mpg_div(NODE *t1, NODE *t2)
1327 {
1328         NODE *r;
1329         int tval;
1330
1331         if (is_mpg_integer(t1) && is_mpg_integer(t2)
1332                         && (mpz_sgn(t2->mpg_i) != 0)    /* not dividing by 0 */
1333                         && mpz_divisible_p(t1->mpg_i, t2->mpg_i)
1334         ) {
1335                 r = mpg_integer();
1336                 mpz_divexact(r->mpg_i, t1->mpg_i, t2->mpg_i);
1337         } else {
1338                 mpfr_ptr p1, p2;
1339                 p1 = MP_FLOAT(t1);
1340                 p2 = MP_FLOAT(t2);
1341                 r = mpg_float();
1342                 tval = mpfr_div(r->mpg_numbr, p1, p2, ROUND_MODE);
1343                 IEEE_FMT(r->mpg_numbr, tval);
1344         }
1345         return r;
1346 }
1347
1348 /* mpg_mod --- modulus operation with arbitrary-precision numbers */
1349
1350 static NODE *
1351 mpg_mod(NODE *t1, NODE *t2)
1352 {
1353         NODE *r;
1354         int tval;
1355
1356         if (is_mpg_integer(t1) && is_mpg_integer(t2)) {
1357                 r = mpg_integer();
1358                 mpz_mod(r->mpg_i, t1->mpg_i, t2->mpg_i);
1359         } else {
1360                 mpfr_ptr p1, p2;
1361                 p1 = MP_FLOAT(t1);
1362                 p2 = MP_FLOAT(t2);
1363                 r = mpg_float();
1364                 tval = mpfr_fmod(r->mpg_numbr, p1, p2, ROUND_MODE);
1365                 IEEE_FMT(r->mpg_numbr, tval);
1366         }
1367         return r;
1368 }
1369
1370 /*
1371  * mpg_interpret --- pre-exec hook in the interpreter. Handles
1372  *      arithmetic operations with MPFR/GMP numbers.
1373  */
1374
1375 static int
1376 mpg_interpret(INSTRUCTION **cp)
1377 {
1378         INSTRUCTION *pc = *cp;  /* current instruction */
1379         OPCODE op;      /* current opcode */
1380         NODE *r = NULL;
1381         NODE *t1, *t2;
1382         NODE **lhs;
1383         int tval;       /* the ternary value returned by a MPFR function */
1384
1385         switch ((op = pc->opcode)) {
1386         case Op_plus_i:
1387                 t2 = force_number(pc->memory);
1388                 goto plus;
1389         case Op_plus:
1390                 t2 = POP_NUMBER();
1391 plus:
1392                 t1 = TOP_NUMBER();
1393                 r = mpg_add(t1, t2);
1394                 DEREF(t1);
1395                 if (op == Op_plus)
1396                         DEREF(t2);
1397                 REPLACE(r);
1398                 break;
1399
1400         case Op_minus_i:
1401                 t2 = force_number(pc->memory);
1402                 goto minus;
1403         case Op_minus:
1404                 t2 = POP_NUMBER();
1405 minus:
1406                 t1 = TOP_NUMBER();
1407                 r = mpg_sub(t1, t2);
1408                 DEREF(t1);
1409                 if (op == Op_minus)
1410                         DEREF(t2);
1411                 REPLACE(r);
1412                 break;
1413
1414         case Op_times_i:
1415                 t2 = force_number(pc->memory);
1416                 goto times;
1417         case Op_times:
1418                 t2 = POP_NUMBER();
1419 times:
1420                 t1 = TOP_NUMBER();
1421                 r = mpg_mul(t1, t2);
1422                 DEREF(t1);
1423                 if (op == Op_times)
1424                         DEREF(t2);
1425                 REPLACE(r);
1426                 break;
1427
1428         case Op_exp_i:
1429                 t2 = force_number(pc->memory);
1430                 goto exp;
1431         case Op_exp:
1432                 t2 = POP_NUMBER();
1433 exp:
1434                 t1 = TOP_NUMBER();
1435                 r = mpg_pow(t1, t2);
1436                 DEREF(t1);
1437                 if (op == Op_exp)
1438                         DEREF(t2);
1439                 REPLACE(r);
1440                 break;
1441
1442         case Op_quotient_i:
1443                 t2 = force_number(pc->memory);
1444                 goto quotient;
1445         case Op_quotient:
1446                 t2 = POP_NUMBER();
1447 quotient:
1448                 t1 = TOP_NUMBER();
1449                 r = mpg_div(t1, t2);
1450                 DEREF(t1);
1451                 if (op == Op_quotient)
1452                         DEREF(t2);
1453                 REPLACE(r);
1454                 break;
1455
1456         case Op_mod_i:
1457                 t2 = force_number(pc->memory);
1458                 goto mod;
1459         case Op_mod:
1460                 t2 = POP_NUMBER();
1461 mod:
1462                 t1 = TOP_NUMBER();
1463                 r = mpg_mod(t1, t2);
1464                 DEREF(t1);
1465                 if (op == Op_mod)
1466                         DEREF(t2);
1467                 REPLACE(r);
1468                 break;
1469
1470         case Op_preincrement:
1471         case Op_predecrement:
1472                 lhs = TOP_ADDRESS();
1473                 t1 = *lhs;
1474                 force_number(t1);
1475
1476                 if (is_mpg_integer(t1)) {
1477                         if (t1->valref == 1 && t1->flags == (MALLOC|MPZN|NUMCUR|NUMBER))
1478                         /* Efficiency hack. Big speed-up (> 30%) in a tight loop */
1479                                 r = t1;
1480                         else
1481                                 r = *lhs = mpg_integer();
1482                         if (op == Op_preincrement)
1483                                 mpz_add_ui(r->mpg_i, t1->mpg_i, 1);
1484                         else
1485                                 mpz_sub_ui(r->mpg_i, t1->mpg_i, 1);
1486                 } else {
1487
1488                         /*
1489                          * An optimization like the one above is not going to work
1490                          * for a floating-point number. With it,
1491                          *   gawk -M 'BEGIN { PREC=53; i=2^53+0.0; PREC=113; ++i; print i}'
1492                          * will output 2^53 instead of 2^53+1.
1493                          */
1494
1495                         r = *lhs = mpg_float();
1496                         tval = mpfr_add_si(r->mpg_numbr, t1->mpg_numbr,
1497                                         op == Op_preincrement ? 1 : -1,
1498                                         ROUND_MODE);
1499                         IEEE_FMT(r->mpg_numbr, tval);
1500                 }
1501                 if (r != t1)
1502                         unref(t1);
1503                 UPREF(r);
1504                 REPLACE(r);
1505                 break;
1506
1507         case Op_postincrement:
1508         case Op_postdecrement:
1509                 lhs = TOP_ADDRESS();
1510                 t1 = *lhs;
1511                 force_number(t1);
1512
1513                 if (is_mpg_integer(t1)) {
1514                         r = mpg_integer();
1515                         mpz_set(r->mpg_i, t1->mpg_i);
1516                         if (t1->valref == 1 && t1->flags == (MALLOC|MPZN|NUMCUR|NUMBER))
1517                         /* Efficiency hack. Big speed-up (> 30%) in a tight loop */
1518                                 t2 = t1;
1519                         else
1520                                 t2 = *lhs = mpg_integer();
1521                         if (op == Op_postincrement)
1522                                 mpz_add_ui(t2->mpg_i, t1->mpg_i, 1);
1523                         else
1524                                 mpz_sub_ui(t2->mpg_i, t1->mpg_i, 1);
1525                 } else {
1526                         r = mpg_float();
1527                         tval = mpfr_set(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1528                         IEEE_FMT(r->mpg_numbr, tval);
1529                         t2 = *lhs = mpg_float();
1530                         tval = mpfr_add_si(t2->mpg_numbr, t1->mpg_numbr,
1531                                         op == Op_postincrement ? 1 : -1,
1532                                         ROUND_MODE);
1533                         IEEE_FMT(t2->mpg_numbr, tval);
1534                 }
1535                 if (t2 != t1)
1536                         unref(t1);
1537                 REPLACE(r);
1538                 break;
1539
1540         case Op_unary_minus:
1541                 t1 = TOP_NUMBER();
1542                 if (is_mpg_float(t1)) {
1543                         r = mpg_float();
1544                         tval = mpfr_neg(r->mpg_numbr, t1->mpg_numbr, ROUND_MODE);
1545                         IEEE_FMT(r->mpg_numbr, tval);
1546                 } else {
1547                         r = mpg_integer();
1548                         mpz_neg(r->mpg_i, t1->mpg_i);
1549                 }
1550                 DEREF(t1);
1551                 REPLACE(r);
1552                 break;
1553
1554         case Op_assign_plus:
1555         case Op_assign_minus:
1556         case Op_assign_times:
1557         case Op_assign_quotient:
1558         case Op_assign_mod:
1559         case Op_assign_exp:
1560                 lhs = POP_ADDRESS();
1561                 t1 = *lhs;
1562                 force_number(t1);
1563                 t2 = TOP_NUMBER();
1564
1565                 switch (op) {
1566                 case Op_assign_plus:
1567                         r = mpg_add(t1, t2);
1568                         break;
1569                 case Op_assign_minus:
1570                         r = mpg_sub(t1, t2);
1571                         break;
1572                 case Op_assign_times:
1573                         r = mpg_mul(t1, t2);
1574                         break;
1575                 case Op_assign_quotient:
1576                         r = mpg_div(t1, t2);
1577                         break;
1578                 case Op_assign_mod:
1579                         r = mpg_mod(t1, t2);
1580                         break;
1581                 case Op_assign_exp:
1582                         r = mpg_pow(t1, t2);
1583                         break;
1584                 default:
1585                         cant_happen();
1586                 }
1587
1588                 DEREF(t2);
1589                 unref(*lhs);
1590                 *lhs = r;
1591                 UPREF(r);
1592                 REPLACE(r);
1593                 break;
1594
1595         default:
1596                 return true;    /* unhandled */
1597         }
1598
1599         *cp = pc->nexti;        /* next instruction to execute */
1600         return false;
1601 }
1602
1603
1604 /* mpg_fmt --- output formatted string with special MPFR/GMP conversion specifiers */
1605
1606 const char *
1607 mpg_fmt(const char *mesg, ...)
1608 {
1609         static char *tmp = NULL;
1610         int ret;
1611         va_list args;
1612
1613         if (tmp != NULL) {
1614                 mpfr_free_str(tmp);
1615                 tmp = NULL;
1616         }
1617         va_start(args, mesg);
1618         ret = mpfr_vasprintf(& tmp, mesg, args);
1619         va_end(args);
1620         if (ret >= 0 && tmp != NULL)
1621                 return tmp;
1622         return mesg;
1623 }
1624
1625 /* mpfr_unset --- clear out the MPFR values */
1626
1627 void
1628 mpfr_unset(NODE *n)
1629 {
1630         if (is_mpg_float(n))
1631                 mpfr_clear(n->mpg_numbr);
1632         else if (is_mpg_integer(n))
1633                 mpz_clear(n->mpg_i);
1634 }
1635
1636 #else
1637
1638 void
1639 set_PREC()
1640 {
1641         /* dummy function */
1642 }
1643
1644 void
1645 set_ROUNDMODE()
1646 {
1647         /* dummy function */
1648 }
1649
1650 void
1651 mpfr_unset(NODE *n)
1652 {
1653         /* dummy function */
1654 }
1655 #endif