~vcs-imports/gawk/master

/*

 * mpfr.c - routines for arbitrary-precision number support in gawk.

 */

/*

 * Copyright (C) 2012, 2013, 2015, 2017, 2018,

 * the Free Software Foundation, Inc.

 *

 * This file is part of GAWK, the GNU implementation of the

 * AWK Programming Language.

 *

 * GAWK is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or

 * (at your option) any later version.

 *

 * GAWK is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA

 */

#include "awk.h"

#ifdef HAVE_MPFR

int MPFR_round_mode = 'N';	// default value

#if !defined(MPFR_VERSION_MAJOR) || MPFR_VERSION_MAJOR < 3

typedef mp_exp_t mpfr_exp_t;

#endif

extern NODE **fmt_list;          /* declared in eval.c */

mpz_t mpzval;	/* GMP integer type, used as temporary in few places */

mpz_t MNR;

mpz_t MFNR;

bool do_ieee_fmt;	/* IEEE-754 floating-point emulation */

mpfr_rnd_t ROUND_MODE;

static mpfr_prec_t default_prec;

static mpfr_rnd_t get_rnd_mode(const char rmode);

static NODE *mpg_force_number(NODE *n);

static NODE *mpg_make_number(double);

static NODE *mpg_format_val(const char *format, int index, NODE *s);

static int mpg_interpret(INSTRUCTION **cp);

static mpfr_exp_t min_exp = MPFR_EMIN_DEFAULT;

static mpfr_exp_t max_exp = MPFR_EMAX_DEFAULT;

/* temporary MPFR floats used to hold converted GMP integer operands */

static mpfr_t _mpf_t1;

static mpfr_t _mpf_t2;

/*

 * PRECISION_MIN is the precision used to initialize _mpf_t1 and _mpf_t2.

 * 64 bits should be enough for exact conversion of most integers to floats.

 */

#define PRECISION_MIN	64

/* mf = { _mpf_t1, _mpf_t2 } */

static inline mpfr_ptr mpg_tofloat(mpfr_ptr mf, mpz_ptr mz);

/* T = {t1, t2} */

#define MP_FLOAT(T) is_mpg_integer(T) ? mpg_tofloat(_mpf_##T, (T)->mpg_i) : (T)->mpg_numbr

/* init_mpfr --- set up MPFR related variables */

void

init_mpfr(mpfr_prec_t prec, const char *rmode)

{

	mpfr_set_default_prec(default_prec = prec);

	ROUND_MODE = get_rnd_mode(rmode[0]);

	mpfr_set_default_rounding_mode(ROUND_MODE);

	make_number = mpg_make_number;

	str2number = mpg_force_number;

	format_val = mpg_format_val;

	cmp_numbers = mpg_cmp;

	mpz_init(MNR);

	mpz_init(MFNR);

	do_ieee_fmt = false;

	mpfr_init2(_mpf_t1, PRECISION_MIN);

	mpfr_init2(_mpf_t2, PRECISION_MIN);

	mpz_init(mpzval);

	register_exec_hook(mpg_interpret, 0);

}

/* cleanup_mpfr --- clean stuff up, mainly for valgrind */

void

cleanup_mpfr(void)

{

	mpfr_clear(_mpf_t1);

	mpfr_clear(_mpf_t2);

}

/* mpg_node --- allocate a node to store MPFR float or GMP integer */

NODE *

mpg_node(unsigned int flags)

{

	NODE *r = make_number_node(flags);

	if (flags == MPFN)

		/* Initialize, set precision to the default precision, and value to NaN */

		mpfr_init(r->mpg_numbr);

	else

		/* Initialize and set value to 0 */

		mpz_init(r->mpg_i);

	return r;

}

/*

 * mpg_make_number --- make a arbitrary-precision number node

 *	and initialize with a C double

 */

static NODE *

mpg_make_number(double x)

{

	NODE *r;

	double ival;

	if ((ival = double_to_int(x)) != x) {

		int tval;

		r = mpg_float();

		tval = mpfr_set_d(r->mpg_numbr, x, ROUND_MODE);

		IEEE_FMT(r->mpg_numbr, tval);

	} else {

		r = mpg_integer();

		mpz_set_d(r->mpg_i, ival);

	}

	return r;

}

/* mpg_strtoui --- assign arbitrary-precision integral value from a string */

int

mpg_strtoui(mpz_ptr zi, char *str, size_t len, char **end, int base)

{

	char *s = str;

	char *start;

	int ret = -1;

	/*

	 * mpz_set_str does not like leading 0x or 0X for hex (or 0 for octal)

	 * with a non-zero base argument.

	*/

	if (base == 16 && len >= 2 && *s == '0' && (s[1] == 'x' || s[1] == 'X')) {

		s += 2; len -= 2;

	} else if (base == 8 && len >= 1 && *s == '0') {

		s++; len--;

	}

	start = s;

	while (len > 0) {

		switch (*s) {

		case '0':

		case '1':

		case '2':

		case '3':

		case '4':

		case '5':

		case '6':

		case '7':

			break;

		case '8':

		case '9':

			if (base == 8)

				goto done;

			break;

		case 'a':

		case 'b':

		case 'c':

		case 'd':

		case 'e':

		case 'f':

		case 'A':

		case 'B':

		case 'C':

		case 'D':

		case 'E':

		case 'F':

			if (base == 16)

				break;

		default:

			goto done;

		}

		s++; len--;

	}

done:

	if (s > start) {

		char save = *s;

		*s = '\0';

		ret = mpz_set_str(zi, start, base);

		*s = save;

	}

	if (end != NULL)

		*end = s;

	return ret;

}

/* mpg_maybe_float --- test if a string may contain arbitrary-precision float */

static int

mpg_maybe_float(const char *str, int use_locale)

{

	int dec_point = '.';

	const char *s = str;

#if defined(HAVE_LOCALE_H)

	/*

	 * loc.decimal_point may not have been initialized yet,

	 * so double check it before using it.

	 */

	if (use_locale && loc.decimal_point != NULL && loc.decimal_point[0] != '\0')

		dec_point = loc.decimal_point[0];	/* XXX --- assumes one char */

#endif

	if (strlen(s) >= 3

		&& (   (   (s[0] == 'i' || s[0] == 'I')

			&& (s[1] == 'n' || s[1] == 'N')

			&& (s[2] == 'f' || s[2] == 'F'))

		    || (   (s[0] == 'n' || s[0] == 'N')

			&& (s[1] == 'a' || s[1] == 'A')

			&& (s[2] == 'n' || s[2] == 'N'))))

		return true;

	for (; *s != '\0'; s++) {

		if (*s == dec_point || *s == 'e' || *s == 'E')

			return true;

	}

	return false;

}

/* mpg_zero --- initialize with arbitrary-precision integer(GMP) and set value to zero */

static inline void

mpg_zero(NODE *n)

{

	if (is_mpg_float(n)) {

		mpfr_clear(n->mpg_numbr);

		n->flags &= ~MPFN;

	}

	if (! is_mpg_integer(n)) {

		mpz_init(n->mpg_i);	/* this also sets its value to 0 */

		n->flags |= MPZN;

	} else

		mpz_set_si(n->mpg_i, 0);

}

/* force_mpnum --- force a value to be a GMP integer or MPFR float */

static int

force_mpnum(NODE *n, int do_nondec, int use_locale)

{

	char *cp, *cpend, *ptr, *cp1;

	char save;

	int tval, base = 10;

	if (n->stlen == 0) {

		mpg_zero(n);

		return false;

	}

	cp = n->stptr;

	cpend = n->stptr + n->stlen;

	while (cp < cpend && isspace((unsigned char) *cp))

		cp++;

	if (cp == cpend) {	/* only spaces */

		mpg_zero(n);

		return false;

	}

	save = *cpend;

	*cpend = '\0';

	if (*cp == '+' || *cp == '-')

		cp1 = cp + 1;

	else

		cp1 = cp;

	if (do_nondec)

		base = get_numbase(cp1, cpend - cp1, use_locale);

	if (! mpg_maybe_float(cp1, use_locale)) {

		mpg_zero(n);

		errno = 0;

		mpg_strtoui(n->mpg_i, cp1, cpend - cp1, & ptr, base);

		if (*cp == '-')

			mpz_neg(n->mpg_i, n->mpg_i);

		goto done;

	}

	if (is_mpg_integer(n)) {

		mpz_clear(n->mpg_i);

		n->flags &= ~MPZN;

	}

	if (! is_mpg_float(n)) {

		mpfr_init(n->mpg_numbr);

		n->flags |= MPFN;

	}

	errno = 0;

	tval = mpfr_strtofr(n->mpg_numbr, cp, & ptr, base, ROUND_MODE);

	IEEE_FMT(n->mpg_numbr, tval);

done:

	/* trailing space is OK for NUMBER */

	while (ptr < cpend && isspace((unsigned char) *ptr))

		ptr++;

	*cpend = save;

	if (errno == 0 && ptr == cpend)

		return true;

	errno = 0;

	return false;

}

/* mpg_force_number --- force a value to be a multiple-precision number */

static NODE *

mpg_force_number(NODE *n)

{

	if ((n->flags & NUMCUR) != 0)

		return n;

	n->flags |= NUMCUR;

	if (force_mpnum(n, (do_non_decimal_data && ! do_traditional), true)) {

		if ((n->flags & USER_INPUT) != 0) {

			/* leave USER_INPUT set to indicate a strnum */

			n->flags &= ~STRING;

			n->flags |= NUMBER;

		}

	} else

		n->flags &= ~USER_INPUT;

	return n;

}

/* mpg_format_val --- format a numeric value based on format */

static NODE *

mpg_format_val(const char *format, int index, NODE *s)

{

	NODE *dummy[2], *r;

	unsigned int oflags;

	/* create dummy node for a sole use of format_tree */

	dummy[1] = s;

	oflags = s->flags;

	if (is_mpg_integer(s) || mpfr_integer_p(s->mpg_numbr)) {

		/* integral value, use %d */

		r = format_tree("%d", 2, dummy, 2);

		s->stfmt = STFMT_UNUSED;

	} else {

		r = format_tree(format, fmt_list[index]->stlen, dummy, 2);

		assert(r != NULL);

		s->stfmt = index;

	}

	s->flags = oflags;

	s->stlen = r->stlen;

	if ((s->flags & (MALLOC|STRCUR)) == (MALLOC|STRCUR))

		efree(s->stptr);

	s->stptr = r->stptr;

	s->flags |= STRCUR;

	s->strndmode = MPFR_round_mode;

	freenode(r);	/* Do not unref(r)! We want to keep s->stptr == r->stpr.  */

	free_wstr(s);

	return s;

}

/* mpg_cmp --- compare two numbers */

int

mpg_cmp(const NODE *t1, const NODE *t2)

{

	/*

	 * For the purposes of sorting, NaN is considered greater than

	 * any other value, and all NaN values are considered equivalent and equal.

	 */

	if (is_mpg_float(t1)) {

		if (is_mpg_float(t2)) {

			if (mpfr_nan_p(t1->mpg_numbr))

				return ! mpfr_nan_p(t2->mpg_numbr);

			if (mpfr_nan_p(t2->mpg_numbr))

				return -1;

			return mpfr_cmp(t1->mpg_numbr, t2->mpg_numbr);

		}

		if (mpfr_nan_p(t1->mpg_numbr))

			return 1;

		return mpfr_cmp_z(t1->mpg_numbr, t2->mpg_i);

	} else if (is_mpg_float(t2)) {

		int ret;

		if (mpfr_nan_p(t2->mpg_numbr))

			return -1;

		ret = mpfr_cmp_z(t2->mpg_numbr, t1->mpg_i);

		return ret > 0 ? -1 : (ret < 0);

	} else if (is_mpg_integer(t1)) {

		return mpz_cmp(t1->mpg_i, t2->mpg_i);

	}

	/* t1 and t2 are AWKNUMs */

	return cmp_awknums(t1, t2);

}

/*

 * mpg_update_var --- update NR or FNR.

 *	NR_node->var_value(mpz_t) = MNR(mpz_t) * LONG_MAX + NR(long)

 */

NODE *

mpg_update_var(NODE *n)

{

	NODE *val = n->var_value;

	long nr = 0;

	mpz_ptr nq = 0;

	if (n == NR_node) {

		nr = NR;

		nq = MNR;

	} else if (n == FNR_node) {

		nr = FNR;

		nq = MFNR;

	} else

		cant_happen();

	if (mpz_sgn(nq) == 0) {

		/* Efficiency hack similar to that for AWKNUM */

		if (is_mpg_float(val) || mpz_get_si(val->mpg_i) != nr) {

			unref(n->var_value);

			val = n->var_value = mpg_integer();

			mpz_set_si(val->mpg_i, nr);

		}

	} else {

		unref(n->var_value);

		val = n->var_value = mpg_integer();

		mpz_set_si(val->mpg_i, nr);

		mpz_addmul_ui(val->mpg_i, nq, LONG_MAX);	/* val->mpg_i += nq * LONG_MAX */

	}

	return val;

}

/* mpg_set_var --- set NR or FNR */

long

mpg_set_var(NODE *n)

{

	long nr = 0;

	mpz_ptr nq = 0, r;

	NODE *val = n->var_value;

	if (n == NR_node)

		nq = MNR;

	else if (n == FNR_node)

		nq = MFNR;

	else

		cant_happen();

	if (is_mpg_integer(val))

		r = val->mpg_i;

	else {

		/* convert float to integer */

		mpfr_get_z(mpzval, val->mpg_numbr, MPFR_RNDZ);

		r = mpzval;

	}

	nr = mpz_fdiv_q_ui(nq, r, LONG_MAX);	/* nq (MNR or MFNR) is quotient */

	return nr;	/* remainder (NR or FNR) */

}

/* set_PREC --- update MPFR PRECISION related variables when PREC assigned to */

void

set_PREC()

{

	long prec = 0;

	NODE *val;

	static const struct ieee_fmt {

		const char *name;

		mpfr_prec_t precision;

		mpfr_exp_t emax;

		mpfr_exp_t emin;

	} ieee_fmts[] = {

		{ "half",	11,	16,	-23	},	/* binary16 */

		{ "single",	24,	128,	-148	},	/* binary32 */

		{ "double",	53,	1024,	-1073	},	/* binary64 */

		{ "quad",	113,	16384,	-16493	},	/* binary128 */

		{ "oct",	237,	262144,	-262377	},	/* binary256, not in the IEEE 754-2008 standard */

		/*

 		 * For any bitwidth = 32 * k ( k >= 4),

 		 * precision = 13 + bitwidth - int(4 * log2(bitwidth))

		 * emax = 1 << bitwidth - precision - 1

		 * emin = 4 - emax - precision

 		 */

	};

	if (! do_mpfr)

		return;

	val = fixtype(PREC_node->var_value);

	if ((val->flags & STRING) != 0) {

		int i, j;

		/* emulate IEEE-754 binary format */

		for (i = 0, j = sizeof(ieee_fmts)/sizeof(ieee_fmts[0]); i < j; i++) {

			if (strcasecmp(ieee_fmts[i].name, val->stptr) == 0)

				break;

		}

		if (i < j) {

			prec = ieee_fmts[i].precision;

			/*

			 * We *DO NOT* change the MPFR exponent range using

			 * mpfr_set_{emin, emax} here. See format_ieee() for details.

			 */

			max_exp = ieee_fmts[i].emax;

			min_exp = ieee_fmts[i].emin;

			do_ieee_fmt = true;

		}

	}

	if (prec <= 0) {

		force_number(val);

		prec = get_number_si(val);

		if (prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX) {

			force_string(val);

			warning(_("PREC value `%.*s' is invalid"), (int) val->stlen, val->stptr);

			prec = 0;

		} else

			do_ieee_fmt = false;

	}

	if (prec > 0)

		mpfr_set_default_prec(default_prec = prec);

}

/* get_rnd_mode --- convert string to MPFR rounding mode */

static mpfr_rnd_t

get_rnd_mode(const char rmode)

{

	switch (rmode) {

	case 'N':

	case 'n':

		return MPFR_RNDN;	/* round to nearest (IEEE-754 roundTiesToEven) */

	case 'Z':

	case 'z':

		return MPFR_RNDZ;	/* round toward zero (IEEE-754 roundTowardZero) */

	case 'U':

	case 'u':

		return MPFR_RNDU;	/* round toward plus infinity (IEEE-754 roundTowardPositive) */

	case 'D':

	case 'd':

		return MPFR_RNDD;	/* round toward minus infinity (IEEE-754 roundTowardNegative) */

#if defined(MPFR_VERSION_MAJOR) && MPFR_VERSION_MAJOR > 2

	case 'A':

	case 'a':

		return MPFR_RNDA;	/* round away from zero */

#endif

	default:

		break;

	}

	return -1;

}

/*

 * set_ROUNDMODE --- update MPFR rounding mode related variables

 *	when ROUNDMODE assigned to

 */

void

set_ROUNDMODE()

{

	if (do_mpfr) {

		mpfr_rnd_t rndm = -1;

		NODE *n;

		n = force_string(ROUNDMODE_node->var_value);

		if (n->stlen == 1)

			rndm = get_rnd_mode(n->stptr[0]);

		if (rndm != -1) {

			mpfr_set_default_rounding_mode(rndm);

			ROUND_MODE = rndm;

			MPFR_round_mode = n->stptr[0];

		} else

			warning(_("RNDMODE value `%.*s' is invalid"), (int) n->stlen, n->stptr);

	}

}

/* format_ieee --- make sure a number follows IEEE-754 floating-point standard */

int

format_ieee(mpfr_ptr x, int tval)

{

	/*

	 * The MPFR doc says that it's our responsibility to make sure all numbers

	 * including those previously created are in range after we've changed the

	 * exponent range. Most MPFR operations and functions require

	 * the input arguments to have exponents within the current exponent range.

	 * Any argument outside the range results in a MPFR assertion failure

	 * like this:

	 *

	 *   $ gawk -M 'BEGIN { x=1.0e-10000; print x+0; PREC="double"; print x+0}'

	 *   1e-10000

	 *   init2.c:52: MPFR assertion failed ....

	 *

	 * A "naive" approach would be to keep track of the ternary state and

	 * the rounding mode for each number, and make sure it is in the current

	 * exponent range (using mpfr_check_range) before using it in an

	 * operation or function. Instead, we adopt the following strategy.

	 *

	 * When gawk starts, the exponent range is the MPFR default

	 * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT]. Any number that gawk

	 * creates must have exponent in this range (excluding infinities, NaNs and zeros).

	 * Each MPFR operation or function is performed with this default exponent

	 * range.

	 *

	 * When emulating IEEE-754 format, the exponents are *temporarily* changed,

	 * mpfr_check_range is called to make sure the number is in the new range,

	 * and mpfr_subnormalize is used to round following the rules of subnormal

	 * arithmetic. The exponent range is then *restored* to the original value

	 * [MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT].

	 */

	(void) mpfr_set_emin(min_exp);

	(void) mpfr_set_emax(max_exp);

	tval = mpfr_check_range(x, tval, ROUND_MODE);

	tval = mpfr_subnormalize(x, tval, ROUND_MODE);

	(void) mpfr_set_emin(MPFR_EMIN_DEFAULT);

	(void) mpfr_set_emax(MPFR_EMAX_DEFAULT);

	return tval;

}

/* do_mpfr_atan2 --- do the atan2 function */

NODE *

do_mpfr_atan2(int nargs)

{

	NODE *t1, *t2, *res;

	mpfr_ptr p1, p2;

	int tval;

	t2 = POP_SCALAR();

	t1 = POP_SCALAR();

	if (do_lint) {

		if ((fixtype(t1)->flags & NUMBER) == 0)

			lintwarn(_("atan2: received non-numeric first argument"));

		if ((fixtype(t2)->flags & NUMBER) == 0)

			lintwarn(_("atan2: received non-numeric second argument"));

	}

	force_number(t1);

	force_number(t2);

	p1 = MP_FLOAT(t1);

	p2 = MP_FLOAT(t2);

	res = mpg_float();

	/* See MPFR documentation for handling of special values like +inf as an argument */

	tval = mpfr_atan2(res->mpg_numbr, p1, p2, ROUND_MODE);

	IEEE_FMT(res->mpg_numbr, tval);

	DEREF(t1);