use std use "fpmath"use "util" /*    See [Mul16] (6.2.1), [Tan89], and [Tan92]. While the exp and    expm1 functions are similar enough to be in the same file (Tang's    algorithms use the same S constants), they are not quite similar    enough to be in the same function. */pkg math =	pkglocal const exp32 : (x : flt32 -> flt32)	pkglocal const exp64 : (x : flt64 -> flt64) 	pkglocal const expm132 : (x : flt32 -> flt32)	pkglocal const expm164 : (x : flt64 -> flt64);; extern const horner_polyu32 : (x : flt32, a : uint32[:] -> flt32)extern const horner_polyu64 : (x : flt64, a : uint64[:] -> flt64) type fltdesc(@f, @u, @i) = struct	explode : (f : @f -> (bool, @i, @u))	assem : (n : bool, e : @i, s : @u -> @f)	horner : (f : @f, a : @u[:] -> @f)	sgnmask : @u	tobits : (f : @f -> @u)	frombits : (u : @u -> @f)	inf : @u	nan : @u 	/* For exp */	thresh_1_min : @u	thresh_1_max : @u	thresh_2 : @u	Ai : @u[:] 	/* For expm1 */	thresh_tiny : @u	thresh_huge_neg : @u	Log3_4 : @u	Log5_4 : @u	Bi : @u[:]	precision : @u 	/* For both */	nabs : @u	inv_L : @u	L1 : @u	L2 : @u	S : (@u, @u)[32];; const desc32 : fltdesc(flt32, uint32, int32) = [	.explode = std.flt32explode,	.assem = std.flt32assem,	.horner = horner_polyu32,	.sgnmask = (1 << 31),	.tobits = std.flt32bits,	.frombits = std.flt32frombits,	.inf = 0x7f800000,	.nan = 0x7fc00000,	.thresh_1_min = 0xc2cff1b4, /* ln(2^(-127 - 23)) ~= -103.9720770839917 */	.thresh_1_max = 0x42b17218, /* ln([2 - 2^(-24)] * 2^127) ~= 88.72283905206 */	.inv_L = 0x4238aa3b, /* L = log(2) / 32, this is 1/L in working precision */	.L1 = 0x3cb17200, /* L1 and L2 sum to give L in high precision, */	.L2 = 0x333fbe8e, /* and L1 has some trailing zeroes. */	.nabs = 9, /* L1 has 9 trailing zeroes in binary */	.Ai = [ /* Coefficients for approximating expm1 in a tiny range */		0x3f000044,		0x3e2aaaec,	][:],	.Log3_4 = 0xbe934b11, /* log(1 - 1/4) < x < log(1 + 1/4) */	.Log5_4 = 0x3e647fbf, /* triggers special expm1 case */	.thresh_tiny = 0x33000000, /* similar to thresh_1_{min,max}, but for expm1 */	.thresh_huge_neg = 0xc18aa122,	.Bi = [ /* Coefficients for approximating expm1 between log(3/4) and log(5/4) */		0x3e2aaaaa,		0x3d2aaaa0,		0x3c0889ff,		0x3ab64de5,		0x394ab327,	][:],	.S = [ /* Double-precise expansions of 2^(J/32) */		(0x3f800000, 0x00000000),		(0x3f82cd80, 0x35531585),		(0x3f85aac0, 0x34d9f312),		(0x3f889800, 0x35e8092e),		(0x3f8b95c0, 0x3471f546),		(0x3f8ea400, 0x36e62d17),		(0x3f91c3c0, 0x361b9d59),		(0x3f94f4c0, 0x36bea3fc),		(0x3f9837c0, 0x36c14637),		(0x3f9b8d00, 0x36e6e755),		(0x3f9ef500, 0x36c98247),		(0x3fa27040, 0x34c0c312),		(0x3fa5fec0, 0x36354d8b),		(0x3fa9a140, 0x3655a754),		(0x3fad5800, 0x36fba90b),		(0x3fb123c0, 0x36d6074b),		(0x3fb504c0, 0x36cccfe7),		(0x3fb8fb80, 0x36bd1d8c),		(0x3fbd0880, 0x368e7d60),		(0x3fc12c40, 0x35cca667),		(0x3fc56700, 0x36a84554),		(0x3fc9b980, 0x36f619b9),		(0x3fce2480, 0x35c151f8),		(0x3fd2a800, 0x366c8f89),		(0x3fd744c0, 0x36f32b5a),		(0x3fdbfb80, 0x36de5f6c),		(0x3fe0ccc0, 0x36776155),		(0x3fe5b900, 0x355cef90),		(0x3feac0c0, 0x355cfba5),		(0x3fefe480, 0x36e66f73),		(0x3ff52540, 0x36f45492),		(0x3ffa8380, 0x36cb6dc9),	],	.precision = 24, /* threshold to prevent underflow in a S[high] + 2^-m */] const desc64 : fltdesc(flt64, uint64, int64) = [	.explode = std.flt64explode,	.assem = std.flt64assem,	.horner = horner_polyu64,	.sgnmask = (1 << 63),	.tobits = std.flt64bits,	.frombits = std.flt64frombits,	.inf = 0x7ff0000000000000,	.nan = 0x7ff8000000000000,	.thresh_1_min = 0xc0874910d52d3052, /* ln(2^(-1023 - 52)) ~= -745.1332191019 */	.thresh_1_max = 0x40862e42fefa39ef, /* ln([2 - 2^(-53)] * 2^1023) ~= 709.78271289 */	.inv_L = 0x40471547652b82fe, /* below as in exp32 */	.L1 = 0x3f962e42fef00000,	.L2 = 0x3d8473de6af278ed,	.nabs = 20,	.Ai = [		0x3fe0000000000000,		0x3fc5555555548f7c,		0x3fa5555555545d4e,		0x3f811115b7aa905e,		0x3f56c1728d739765,	][:],	.Log3_4 = 0xbfd269621134db93,	.Log5_4 = 0x3fcc8ff7c79a9a22,	.thresh_tiny = 0x3c90000000000000,	.thresh_huge_neg = 0xc042b708872320e1,	.Bi = [		0x3fc5555555555549,		0x3fa55555555554b6,		0x3f8111111111a9f3,		0x3f56c16c16ce14c6,		0x3f2a01a01159dd2d,		0x3efa019f635825c4,		0x3ec71e14bfe3db59,		0x3e928295484734ea,		0x3e5a2836aa646b96,	][:],	.S = [		(0x3ff0000000000000, 0x0000000000000000),		(0x3ff059b0d3158540, 0x3d0a1d73e2a475b4),		(0x3ff0b5586cf98900, 0x3ceec5317256e308),		(0x3ff11301d0125b40, 0x3cf0a4ebbf1aed93),		(0x3ff172b83c7d5140, 0x3d0d6e6fbe462876),		(0x3ff1d4873168b980, 0x3d053c02dc0144c8),		(0x3ff2387a6e756200, 0x3d0c3360fd6d8e0b),		(0x3ff29e9df51fdec0, 0x3d009612e8afad12),		(0x3ff306fe0a31b700, 0x3cf52de8d5a46306),		(0x3ff371a7373aa9c0, 0x3ce54e28aa05e8a9),		(0x3ff3dea64c123400, 0x3d011ada0911f09f),		(0x3ff44e0860618900, 0x3d068189b7a04ef8),		(0x3ff4bfdad5362a00, 0x3d038ea1cbd7f621),		(0x3ff5342b569d4f80, 0x3cbdf0a83c49d86a),		(0x3ff5ab07dd485400, 0x3d04ac64980a8c8f),		(0x3ff6247eb03a5580, 0x3cd2c7c3e81bf4b7),		(0x3ff6a09e667f3bc0, 0x3ce921165f626cdd),		(0x3ff71f75e8ec5f40, 0x3d09ee91b8797785),		(0x3ff7a11473eb0180, 0x3cdb5f54408fdb37),		(0x3ff82589994cce00, 0x3cf28acf88afab35),		(0x3ff8ace5422aa0c0, 0x3cfb5ba7c55a192d),		(0x3ff93737b0cdc5c0, 0x3d027a280e1f92a0),		(0x3ff9c49182a3f080, 0x3cf01c7c46b071f3),		(0x3ffa5503b23e2540, 0x3cfc8b424491caf8),		(0x3ffae89f995ad380, 0x3d06af439a68bb99),		(0x3ffb7f76f2fb5e40, 0x3cdbaa9ec206ad4f),		(0x3ffc199bdd855280, 0x3cfc2220cb12a092),		(0x3ffcb720dcef9040, 0x3d048a81e5e8f4a5),		(0x3ffd5818dcfba480, 0x3cdc976816bad9b8),		(0x3ffdfc97337b9b40, 0x3cfeb968cac39ed3),		(0x3ffea4afa2a490c0, 0x3cf9858f73a18f5e),		(0x3fff50765b6e4540, 0x3c99d3e12dd8a18b),	],	.precision = 53,] const exp32 = {x : flt32	-> expgen(x, desc32)} const exp64 = {x : flt64	-> expgen(x, desc64)} generic expgen = {x : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i	var b = d.tobits(x)	var n, e, s	(n, e, s) = d.explode(x) 	/*	   Detect if exp(f) would round to outside representability.	   We don't do bias adjustment, so Tang's thresholds can	   be tightened.	 */	if !n && b >= d.thresh_1_max		-> d.frombits(d.inf)	elif n && b > d.thresh_1_min		-> (0.0 : @f)	;; 	/* Detect if exp(f) is close to 1. */	if (b & ~d.sgnmask) <= d.thresh_2		-> (1.0 : @f)	;; 	/* Argument reduction to [ -ln(2)/64, ln(2)/64 ] */	var inv_L = d.frombits(d.inv_L) 	var N = rn(x * inv_L)	var N2  = N % (32 : @i)	if N2 < 0		N2 += (32 : @i)	;;	var N1 = N - N2 	var R1 : @f, R2 : @f	var Nf = (N : @f) 	/*	   There are few enough significant bits that these are all	   exact, without need for an FMA. R1 + R2 will approximate	   (very well) f reduced into [ -ln(2)/64, ln(2)/64 ]	 */	if std.abs(N) >= (1 << d.nabs)		R1 = (x - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))	else		R1 = x - (N : @f) * d.frombits(d.L1)	;;	R2 = -1.0 * (N : @f) * d.frombits(d.L2) 	var M = (N1 : @i) / 32	var J = (N2 : std.size) 	/* Compute a polynomial approximation of exp1m */	var R = R1 + R2	var Q = R * R * d.horner(R, d.Ai)	var P = R1 + (R2 + Q) 	/* Expand out from reduced range */	var Su_hi, Su_lo	(Su_hi, Su_lo) = d.S[J]	var S_hi = d.frombits(Su_hi)	var S_lo = d.frombits(Su_lo) 	var S = S_hi + S_lo	var unscaled = S_hi + (S_lo + (P * S))	var nu, eu, su	(nu, eu, su) = d.explode(unscaled)	var exp = d.assem(nu, eu + M, su)	if (d.tobits(exp) == 0)		/* Make sure we don't quite return 0 */		-> d.frombits(1)	;; 	-> exp} const expm132 = {x : flt32	-> expm1gen(x, desc32)} const expm164 = {x : flt64	-> expm1gen(x, desc64)} generic expm1gen = {x : @f, d : fltdesc(@f, @u, @i) :: \    numeric,floating,std.equatable @f,    numeric,integral @u,    numeric,integral @i,    roundable @f -> @i	var b = d.tobits(x)	var n, e, s	(n, e, s) = d.explode(x) 	/* Special cases: +/- 0, inf, NaN, tiny, and huge */	if (b & ~d.sgnmask == 0)		-> x	elif n && (b & ~d.sgnmask == d.inf)		-> (-1.0 : @f)	elif (b & ~d.sgnmask == d.inf)		-> x	elif std.isnan(x)		-> d.frombits(d.nan)	elif (b & ~d.sgnmask) <= d.thresh_tiny		var two_to_large = d.assem(false, 100, 0)		var two_to_small = d.assem(false, -100, 0)		var abs_x = d.assem(false, e, s)		-> (two_to_large * x + abs_x) * two_to_small	elif !n && b >= d.thresh_1_max /* exp(x) = oo <=> expm1(x) = oo, as it turns out */		-> d.frombits(d.inf)	elif n && b >= d.thresh_huge_neg		-> (-1.0 : @f)	;; 	if ((n && b < d.Log3_4) || (!n && b < d.Log5_4))		/* Procedure 2 */ 		/* compute x^2  / 2 with extra precision */		var u = round(x, d)		var v = x - u		var y = u * u * (0.5 : @f)		var z = v * (x + u) * (0.5 : @f)		var q = x * x * x * d.horner(x, d.Bi) 		var yn, ye, ys		(yn, ye, ys) = d.explode(y)		if (ye >= -7)			-> (u + y) + (q + (v  + z))		else			-> x + (y + (q + z))		;;	;; 	/* Procedure 1 */	var inv_L = d.frombits(d.inv_L) 	var N = rn(x * inv_L)	var N2 = N % (32 : @i)	if N2 < 0		N2 += (32 : @i)	;;	var N1 = N - N2 	var R1 : @f, R2 : @f 	/*	   As in the exp case, R1 + R2 very well approximates f	   reduced into [ -ln(2)/64, ln(2)/64 ]	 */	if std.abs(N) >= (1 << d.nabs)		R1 = (x - (N1 : @f) * d.frombits(d.L1)) - ((N2 : @f) * d.frombits(d.L1))	else		R1 = x - (N : @f) * d.frombits(d.L1)	;;	R2 = -1.0 * (N : @f) * d.frombits(d.L2) 	var M = (N1 : @i) / 32	var J = (N2 : std.size) 	/* Compute a polynomial approximation of exp1m */	var R = R1 + R2	var Q = R * R * d.horner(R, d.Ai)	var P = R1 + (R2 + Q) 	/* Expand out */	var Su_hi, Su_lo	(Su_hi, Su_lo) = d.S[J]	var S_hi = d.frombits(Su_hi)	var S_lo = d.frombits(Su_lo)	var S = S_hi + S_lo 	/*	   Note that [Tan92] includes cases to avoid tripping the	   underflow flag when not technically required. We currently	   do not care about such flags, so that logic is omitted.	 */	if M >= d.precision		-> scale2(S_hi + (S * P + (S_lo - scale2(1.0, -M))), M)	elif M <= -8		-> scale2(S_hi + (S * P + S_lo), M) - (1.0 : @f)	;; 	-> scale2((S_hi - scale2(1.0, -M)) + (S_hi * P + S_lo * (P + (1.0 : @f))), M)} generic round = {f : @f, d : fltdesc(@f, @u, @i) :: numeric,floating,std.equatable @f, numeric,integral @u, numeric,integral @i, roundable @f -> @i	var b = d.tobits(f)	var n, e, s	(n, e, s) = d.explode(f)	var mask = ~((1 << d.nabs) - 1)	if need_round_away(0, ((s & mask) : uint64), (d.nabs : int64) + 1)		-> d.assem(n, e, 1 + s & mask)	;;	-> d.assem(n, e, s & mask)}