~vcs-imports/gawk/master

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