Commit 36e2eaf2bd6d7c69ea6ffbf866c57ea96d6026f1 | zig

author	David Rubin <david@vortan.dev> 2026-04-09 17:38:17+02:00
committer	Andrew Kelley <andrew@ziglang.org> 2026-04-09 17:38:17+02:00
log	36e2eaf2bd6d7c69ea6ffbf866c57ea96d6026f1
tree	a6d10f4ed66c7f71aacce20a0595fc3ce38dd82f
parent	5ccfeb9268a5cb42efc7cd3350ae73bd960eafc3

revert zig libc implementation of pow (#31794)

Partially reverts commit b5fbed986b.

Reverts commit 7aae7dd3f4.

Closes #31207

Reviewed-on: https://codeberg.org/ziglang/zig/pulls/31794
Co-authored-by: David Rubin <david@vortan.dev>
Co-committed-by: David Rubin <david@vortan.dev>

8 files changed, 759 insertions(+), 5 deletions(-)

lib
- c
  - math.zig
- libc
  - musl
    src
    math
    exp_data.c
    exp_data.h
    pow.c
    pow_data.c
    pow_data.h
src
- libs
  - musl.zig
  - wasi_libc.zig

lib/c/math.zig-5

...	...	@@ -76,7 +76,6 @@ comptime {
76	76	symbol(&modf, "modf");
77	77	symbol(&nan, "nan");
78	78	symbol(&nanf, "nanf");
79		symbol(&pow, "pow");
80	79	symbol(&pow10, "pow10");
81	80	symbol(&pow10f, "pow10f");
82	81	symbol(&tanh, "tanh");
...	...	@@ -346,10 +345,6 @@ fn nanl(_: [*:0]const c_char) callconv(.c) c_longdouble {
346	345	return math.nan(c_longdouble);
347	346	}
348	347
349		fn pow(x: f64, y: f64) callconv(.c) f64 {
350		return math.pow(f64, x, y);
351		}
352
353	348	fn pow10(x: f64) callconv(.c) f64 {
354	349	return exp10(x);
355	350	}

lib/libc/musl/src/math/exp_data.c created+182

...	...	@@ -0,0 +1,182 @@
	1	/*
	2	* Shared data between exp, exp2 and pow.
	3	*
	4	* Copyright (c) 2018, Arm Limited.
	5	* SPDX-License-Identifier: MIT
	6	*/
	7
	8	#include "exp_data.h"
	9
	10	#define N (1 << EXP_TABLE_BITS)
	11
	12	const struct exp_data __exp_data = {
	13	// N/ln2
	14	.invln2N = 0x1.71547652b82fep0 * N,
	15	// -ln2/N
	16	.negln2hiN = -0x1.62e42fefa0000p-8,
	17	.negln2loN = -0x1.cf79abc9e3b3ap-47,
	18	// Used for rounding when !TOINT_INTRINSICS
	19	#if EXP_USE_TOINT_NARROW
	20	.shift = 0x1800000000.8p0,
	21	#else
	22	.shift = 0x1.8p52,
	23	#endif
	24	// exp polynomial coefficients.
	25	.poly = {
	26	// abs error: 1.555*2^-66
	27	// ulp error: 0.509 (0.511 without fma)
	28	// if \|x\| < ln2/256+eps
	29	// abs error if \|x\| < ln2/256+0x1p-15: 1.09*2^-65
	30	// abs error if \|x\| < ln2/128: 1.7145*2^-56
	31	0x1.ffffffffffdbdp-2,
	32	0x1.555555555543cp-3,
	33	0x1.55555cf172b91p-5,
	34	0x1.1111167a4d017p-7,
	35	},
	36	.exp2_shift = 0x1.8p52 / N,
	37	// exp2 polynomial coefficients.
	38	.exp2_poly = {
	39	// abs error: 1.2195*2^-65
	40	// ulp error: 0.507 (0.511 without fma)
	41	// if \|x\| < 1/256
	42	// abs error if \|x\| < 1/128: 1.9941*2^-56
	43	0x1.62e42fefa39efp-1,
	44	0x1.ebfbdff82c424p-3,
	45	0x1.c6b08d70cf4b5p-5,
	46	0x1.3b2abd24650ccp-7,
	47	0x1.5d7e09b4e3a84p-10,
	48	},
	49	// 2^(k/N) ~= H[k]*(1 + T[k]) for int k in [0,N)
	50	// tab[2*k] = asuint64(T[k])
	51	// tab[2*k+1] = asuint64(H[k]) - (k << 52)/N
	52	.tab = {
	53	0x0, 0x3ff0000000000000,
	54	0x3c9b3b4f1a88bf6e, 0x3feff63da9fb3335,
	55	0xbc7160139cd8dc5d, 0x3fefec9a3e778061,
	56	0xbc905e7a108766d1, 0x3fefe315e86e7f85,
	57	0x3c8cd2523567f613, 0x3fefd9b0d3158574,
	58	0xbc8bce8023f98efa, 0x3fefd06b29ddf6de,
	59	0x3c60f74e61e6c861, 0x3fefc74518759bc8,
	60	0x3c90a3e45b33d399, 0x3fefbe3ecac6f383,
	61	0x3c979aa65d837b6d, 0x3fefb5586cf9890f,
	62	0x3c8eb51a92fdeffc, 0x3fefac922b7247f7,
	63	0x3c3ebe3d702f9cd1, 0x3fefa3ec32d3d1a2,
	64	0xbc6a033489906e0b, 0x3fef9b66affed31b,
	65	0xbc9556522a2fbd0e, 0x3fef9301d0125b51,
	66	0xbc5080ef8c4eea55, 0x3fef8abdc06c31cc,
	67	0xbc91c923b9d5f416, 0x3fef829aaea92de0,
	68	0x3c80d3e3e95c55af, 0x3fef7a98c8a58e51,
	69	0xbc801b15eaa59348, 0x3fef72b83c7d517b,
	70	0xbc8f1ff055de323d, 0x3fef6af9388c8dea,
	71	0x3c8b898c3f1353bf, 0x3fef635beb6fcb75,
	72	0xbc96d99c7611eb26, 0x3fef5be084045cd4,
	73	0x3c9aecf73e3a2f60, 0x3fef54873168b9aa,
	74	0xbc8fe782cb86389d, 0x3fef4d5022fcd91d,
	75	0x3c8a6f4144a6c38d, 0x3fef463b88628cd6,
	76	0x3c807a05b0e4047d, 0x3fef3f49917ddc96,
	77	0x3c968efde3a8a894, 0x3fef387a6e756238,
	78	0x3c875e18f274487d, 0x3fef31ce4fb2a63f,
	79	0x3c80472b981fe7f2, 0x3fef2b4565e27cdd,
	80	0xbc96b87b3f71085e, 0x3fef24dfe1f56381,
	81	0x3c82f7e16d09ab31, 0x3fef1e9df51fdee1,
	82	0xbc3d219b1a6fbffa, 0x3fef187fd0dad990,
	83	0x3c8b3782720c0ab4, 0x3fef1285a6e4030b,
	84	0x3c6e149289cecb8f, 0x3fef0cafa93e2f56,
	85	0x3c834d754db0abb6, 0x3fef06fe0a31b715,
	86	0x3c864201e2ac744c, 0x3fef0170fc4cd831,
	87	0x3c8fdd395dd3f84a, 0x3feefc08b26416ff,
	88	0xbc86a3803b8e5b04, 0x3feef6c55f929ff1,
	89	0xbc924aedcc4b5068, 0x3feef1a7373aa9cb,
	90	0xbc9907f81b512d8e, 0x3feeecae6d05d866,
	91	0xbc71d1e83e9436d2, 0x3feee7db34e59ff7,
	92	0xbc991919b3ce1b15, 0x3feee32dc313a8e5,
	93	0x3c859f48a72a4c6d, 0x3feedea64c123422,
	94	0xbc9312607a28698a, 0x3feeda4504ac801c,
	95	0xbc58a78f4817895b, 0x3feed60a21f72e2a,
	96	0xbc7c2c9b67499a1b, 0x3feed1f5d950a897,
	97	0x3c4363ed60c2ac11, 0x3feece086061892d,
	98	0x3c9666093b0664ef, 0x3feeca41ed1d0057,
	99	0x3c6ecce1daa10379, 0x3feec6a2b5c13cd0,
	100	0x3c93ff8e3f0f1230, 0x3feec32af0d7d3de,
	101	0x3c7690cebb7aafb0, 0x3feebfdad5362a27,
	102	0x3c931dbdeb54e077, 0x3feebcb299fddd0d,
	103	0xbc8f94340071a38e, 0x3feeb9b2769d2ca7,
	104	0xbc87deccdc93a349, 0x3feeb6daa2cf6642,
	105	0xbc78dec6bd0f385f, 0x3feeb42b569d4f82,
	106	0xbc861246ec7b5cf6, 0x3feeb1a4ca5d920f,
	107	0x3c93350518fdd78e, 0x3feeaf4736b527da,
	108	0x3c7b98b72f8a9b05, 0x3feead12d497c7fd,
	109	0x3c9063e1e21c5409, 0x3feeab07dd485429,
	110	0x3c34c7855019c6ea, 0x3feea9268a5946b7,
	111	0x3c9432e62b64c035, 0x3feea76f15ad2148,
	112	0xbc8ce44a6199769f, 0x3feea5e1b976dc09,
	113	0xbc8c33c53bef4da8, 0x3feea47eb03a5585,
	114	0xbc845378892be9ae, 0x3feea34634ccc320,
	115	0xbc93cedd78565858, 0x3feea23882552225,
	116	0x3c5710aa807e1964, 0x3feea155d44ca973,
	117	0xbc93b3efbf5e2228, 0x3feea09e667f3bcd,
	118	0xbc6a12ad8734b982, 0x3feea012750bdabf,
	119	0xbc6367efb86da9ee, 0x3fee9fb23c651a2f,
	120	0xbc80dc3d54e08851, 0x3fee9f7df9519484,
	121	0xbc781f647e5a3ecf, 0x3fee9f75e8ec5f74,
	122	0xbc86ee4ac08b7db0, 0x3fee9f9a48a58174,
	123	0xbc8619321e55e68a, 0x3fee9feb564267c9,
	124	0x3c909ccb5e09d4d3, 0x3feea0694fde5d3f,
	125	0xbc7b32dcb94da51d, 0x3feea11473eb0187,
	126	0x3c94ecfd5467c06b, 0x3feea1ed0130c132,
	127	0x3c65ebe1abd66c55, 0x3feea2f336cf4e62,
	128	0xbc88a1c52fb3cf42, 0x3feea427543e1a12,
	129	0xbc9369b6f13b3734, 0x3feea589994cce13,
	130	0xbc805e843a19ff1e, 0x3feea71a4623c7ad,
	131	0xbc94d450d872576e, 0x3feea8d99b4492ed,
	132	0x3c90ad675b0e8a00, 0x3feeaac7d98a6699,
	133	0x3c8db72fc1f0eab4, 0x3feeace5422aa0db,
	134	0xbc65b6609cc5e7ff, 0x3feeaf3216b5448c,
	135	0x3c7bf68359f35f44, 0x3feeb1ae99157736,
	136	0xbc93091fa71e3d83, 0x3feeb45b0b91ffc6,
	137	0xbc5da9b88b6c1e29, 0x3feeb737b0cdc5e5,
	138	0xbc6c23f97c90b959, 0x3feeba44cbc8520f,
	139	0xbc92434322f4f9aa, 0x3feebd829fde4e50,
	140	0xbc85ca6cd7668e4b, 0x3feec0f170ca07ba,
	141	0x3c71affc2b91ce27, 0x3feec49182a3f090,
	142	0x3c6dd235e10a73bb, 0x3feec86319e32323,
	143	0xbc87c50422622263, 0x3feecc667b5de565,
	144	0x3c8b1c86e3e231d5, 0x3feed09bec4a2d33,
	145	0xbc91bbd1d3bcbb15, 0x3feed503b23e255d,
	146	0x3c90cc319cee31d2, 0x3feed99e1330b358,
	147	0x3c8469846e735ab3, 0x3feede6b5579fdbf,
	148	0xbc82dfcd978e9db4, 0x3feee36bbfd3f37a,
	149	0x3c8c1a7792cb3387, 0x3feee89f995ad3ad,
	150	0xbc907b8f4ad1d9fa, 0x3feeee07298db666,
	151	0xbc55c3d956dcaeba, 0x3feef3a2b84f15fb,
	152	0xbc90a40e3da6f640, 0x3feef9728de5593a,
	153	0xbc68d6f438ad9334, 0x3feeff76f2fb5e47,
	154	0xbc91eee26b588a35, 0x3fef05b030a1064a,
	155	0x3c74ffd70a5fddcd, 0x3fef0c1e904bc1d2,
	156	0xbc91bdfbfa9298ac, 0x3fef12c25bd71e09,
	157	0x3c736eae30af0cb3, 0x3fef199bdd85529c,
	158	0x3c8ee3325c9ffd94, 0x3fef20ab5fffd07a,
	159	0x3c84e08fd10959ac, 0x3fef27f12e57d14b,
	160	0x3c63cdaf384e1a67, 0x3fef2f6d9406e7b5,
	161	0x3c676b2c6c921968, 0x3fef3720dcef9069,
	162	0xbc808a1883ccb5d2, 0x3fef3f0b555dc3fa,
	163	0xbc8fad5d3ffffa6f, 0x3fef472d4a07897c,
	164	0xbc900dae3875a949, 0x3fef4f87080d89f2,
	165	0x3c74a385a63d07a7, 0x3fef5818dcfba487,
	166	0xbc82919e2040220f, 0x3fef60e316c98398,
	167	0x3c8e5a50d5c192ac, 0x3fef69e603db3285,
	168	0x3c843a59ac016b4b, 0x3fef7321f301b460,
	169	0xbc82d52107b43e1f, 0x3fef7c97337b9b5f,
	170	0xbc892ab93b470dc9, 0x3fef864614f5a129,
	171	0x3c74b604603a88d3, 0x3fef902ee78b3ff6,
	172	0x3c83c5ec519d7271, 0x3fef9a51fbc74c83,
	173	0xbc8ff7128fd391f0, 0x3fefa4afa2a490da,
	174	0xbc8dae98e223747d, 0x3fefaf482d8e67f1,
	175	0x3c8ec3bc41aa2008, 0x3fefba1bee615a27,
	176	0x3c842b94c3a9eb32, 0x3fefc52b376bba97,
	177	0x3c8a64a931d185ee, 0x3fefd0765b6e4540,
	178	0xbc8e37bae43be3ed, 0x3fefdbfdad9cbe14,
	179	0x3c77893b4d91cd9d, 0x3fefe7c1819e90d8,
	180	0x3c5305c14160cc89, 0x3feff3c22b8f71f1,
	181	},
	182	};

lib/libc/musl/src/math/exp_data.h created+26

...	...	@@ -0,0 +1,26 @@
	1	/*
	2	* Copyright (c) 2018, Arm Limited.
	3	* SPDX-License-Identifier: MIT
	4	*/
	5	#ifndef _EXP_DATA_H
	6	#define _EXP_DATA_H
	7
	8	#include <features.h>
	9	#include <stdint.h>
	10
	11	#define EXP_TABLE_BITS 7
	12	#define EXP_POLY_ORDER 5
	13	#define EXP_USE_TOINT_NARROW 0
	14	#define EXP2_POLY_ORDER 5
	15	extern hidden const struct exp_data {
	16	double invln2N;
	17	double shift;
	18	double negln2hiN;
	19	double negln2loN;
	20	double poly[4]; /* Last four coefficients. */
	21	double exp2_shift;
	22	double exp2_poly[EXP2_POLY_ORDER];
	23	uint64_t tab[2*(1 << EXP_TABLE_BITS)];
	24	} __exp_data;
	25
	26	#endif

lib/libc/musl/src/math/pow.c created+343

...	...	@@ -0,0 +1,343 @@
	1	/*
	2	* Double-precision x^y function.
	3	*
	4	* Copyright (c) 2018, Arm Limited.
	5	* SPDX-License-Identifier: MIT
	6	*/
	7
	8	#include <math.h>
	9	#include <stdint.h>
	10	#include "libm.h"
	11	#include "exp_data.h"
	12	#include "pow_data.h"
	13
	14	/*
	15	Worst-case error: 0.54 ULP (~= ulperr_exp + 1024Ln2relerr_log*2^53)
	16	relerr_log: 1.3 * 2^-68 (Relative error of log, 1.5 * 2^-68 without fma)
	17	ulperr_exp: 0.509 ULP (ULP error of exp, 0.511 ULP without fma)
	18	*/
	19
	20	#define T __pow_log_data.tab
	21	#define A __pow_log_data.poly
	22	#define Ln2hi __pow_log_data.ln2hi
	23	#define Ln2lo __pow_log_data.ln2lo
	24	#define N (1 << POW_LOG_TABLE_BITS)
	25	#define OFF 0x3fe6955500000000
	26
	27	/* Top 12 bits of a double (sign and exponent bits). */
	28	static inline uint32_t top12(double x)
	29	{
	30	return asuint64(x) >> 52;
	31	}
	32
	33	/* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
	34	additional 15 bits precision. IX is the bit representation of x, but
	35	normalized in the subnormal range using the sign bit for the exponent. */
	36	static inline double_t log_inline(uint64_t ix, double_t *tail)
	37	{
	38	/* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
	39	double_t z, r, y, invc, logc, logctail, kd, hi, t1, t2, lo, lo1, lo2, p;
	40	uint64_t iz, tmp;
	41	int k, i;
	42
	43	/* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
	44	The range is split into N subintervals.
	45	The ith subinterval contains z and c is near its center. */
	46	tmp = ix - OFF;
	47	i = (tmp >> (52 - POW_LOG_TABLE_BITS)) % N;
	48	k = (int64_t)tmp >> 52; /* arithmetic shift */
	49	iz = ix - (tmp & 0xfffULL << 52);
	50	z = asdouble(iz);
	51	kd = (double_t)k;
	52
	53	/* log(x) = kLn2 + log(c) + log1p(z/c-1). /
	54	invc = T[i].invc;
	55	logc = T[i].logc;
	56	logctail = T[i].logctail;
	57
	58	/* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
	59	\|z/c - 1\| < 1/N, so r = z/c - 1 is exactly representible. */
	60	#if __FP_FAST_FMA
	61	r = __builtin_fma(z, invc, -1.0);
	62	#else
	63	/* Split z such that rhi, rlo and rhirhi are exact and \|rlo\| <= \|r\|. /
	64	double_t zhi = asdouble((iz + (1ULL << 31)) & (-1ULL << 32));
	65	double_t zlo = z - zhi;
	66	double_t rhi = zhi * invc - 1.0;
	67	double_t rlo = zlo * invc;
	68	r = rhi + rlo;
	69	#endif
	70
	71	/* kLn2 + log(c) + r. /
	72	t1 = kd * Ln2hi + logc;
	73	t2 = t1 + r;
	74	lo1 = kd * Ln2lo + logctail;
	75	lo2 = t1 - t2 + r;
	76
	77	/* Evaluation is optimized assuming superscalar pipelined execution. */
	78	double_t ar, ar2, ar3, lo3, lo4;
	79	ar = A[0] * r; /* A[0] = -0.5. */
	80	ar2 = r * ar;
	81	ar3 = r * ar2;
	82	/* kLn2 + log(c) + r + A[0]rr. /
	83	#if __FP_FAST_FMA
	84	hi = t2 + ar2;
	85	lo3 = __builtin_fma(ar, r, -ar2);
	86	lo4 = t2 - hi + ar2;
	87	#else
	88	double_t arhi = A[0] * rhi;
	89	double_t arhi2 = rhi * arhi;
	90	hi = t2 + arhi2;
	91	lo3 = rlo * (ar + arhi);
	92	lo4 = t2 - hi + arhi2;
	93	#endif
	94	/* p = log1p(r) - r - A[0]rr. */
	95	p = (ar3 * (A[1] + r * A[2] +
	96	ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r * A[6]))));
	97	lo = lo1 + lo2 + lo3 + lo4 + p;
	98	y = hi + lo;
	99	*tail = hi - y + lo;
	100	return y;
	101	}
	102
	103	#undef N
	104	#undef T
	105	#define N (1 << EXP_TABLE_BITS)
	106	#define InvLn2N __exp_data.invln2N
	107	#define NegLn2hiN __exp_data.negln2hiN
	108	#define NegLn2loN __exp_data.negln2loN
	109	#define Shift __exp_data.shift
	110	#define T __exp_data.tab
	111	#define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
	112	#define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
	113	#define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
	114	#define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
	115	#define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
	116
	117	/* Handle cases that may overflow or underflow when computing the result that
	118	is scale*(1+TMP) without intermediate rounding. The bit representation of
	119	scale is in SBITS, however it has a computed exponent that may have
	120	overflown into the sign bit so that needs to be adjusted before using it as
	121	a double. (int32_t)KI is the k used in the argument reduction and exponent
	122	adjustment of scale, positive k here means the result may overflow and
	123	negative k means the result may underflow. */
	124	static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki)
	125	{
	126	double_t scale, y;
	127
	128	if ((ki & 0x80000000) == 0) {
	129	/* k > 0, the exponent of scale might have overflowed by <= 460. */
	130	sbits -= 1009ull << 52;
	131	scale = asdouble(sbits);
	132	y = 0x1p1009 * (scale + scale * tmp);
	133	return eval_as_double(y);
	134	}
	135	/* k < 0, need special care in the subnormal range. */
	136	sbits += 1022ull << 52;
	137	/* Note: sbits is signed scale. */
	138	scale = asdouble(sbits);
	139	y = scale + scale * tmp;
	140	if (fabs(y) < 1.0) {
	141	/* Round y to the right precision before scaling it into the subnormal
	142	range to avoid double rounding that can cause 0.5+E/2 ulp error where
	143	E is the worst-case ulp error outside the subnormal range. So this
	144	is only useful if the goal is better than 1 ulp worst-case error. */
	145	double_t hi, lo, one = 1.0;
	146	if (y < 0.0)
	147	one = -1.0;
	148	lo = scale - y + scale * tmp;
	149	hi = one + y;
	150	lo = one - hi + y + lo;
	151	y = eval_as_double(hi + lo) - one;
	152	/* Fix the sign of 0. */
	153	if (y == 0.0)
	154	y = asdouble(sbits & 0x8000000000000000);
	155	/* The underflow exception needs to be signaled explicitly. */
	156	fp_force_eval(fp_barrier(0x1p-1022) * 0x1p-1022);
	157	}
	158	y = 0x1p-1022 * y;
	159	return eval_as_double(y);
	160	}
	161
	162	#define SIGN_BIAS (0x800 << EXP_TABLE_BITS)
	163
	164	/* Computes sign*exp(x+xtail) where \|xtail\| < 2^-8/N and \|xtail\| <= \|x\|.
	165	The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1. */
	166	static inline double exp_inline(double_t x, double_t xtail, uint32_t sign_bias)
	167	{
	168	uint32_t abstop;
	169	uint64_t ki, idx, top, sbits;
	170	/* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
	171	double_t kd, z, r, r2, scale, tail, tmp;
	172
	173	abstop = top12(x) & 0x7ff;
	174	if (predict_false(abstop - top12(0x1p-54) >=
	175	top12(512.0) - top12(0x1p-54))) {
	176	if (abstop - top12(0x1p-54) >= 0x80000000) {
	177	/* Avoid spurious underflow for tiny x. */
	178	/* Note: 0 is common input. */
	179	double_t one = WANT_ROUNDING ? 1.0 + x : 1.0;
	180	return sign_bias ? -one : one;
	181	}
	182	if (abstop >= top12(1024.0)) {
	183	/* Note: inf and nan are already handled. */
	184	if (asuint64(x) >> 63)
	185	return __math_uflow(sign_bias);
	186	else
	187	return __math_oflow(sign_bias);
	188	}
	189	/* Large x is special cased below. */
	190	abstop = 0;
	191	}
	192
	193	/* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
	194	/* x = ln2/Nk + r, with int k and r in [-ln2/2N, ln2/2N]. /
	195	z = InvLn2N * x;
	196	#if TOINT_INTRINSICS
	197	kd = roundtoint(z);
	198	ki = converttoint(z);
	199	#elif EXP_USE_TOINT_NARROW
	200	/* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes. */
	201	kd = eval_as_double(z + Shift);
	202	ki = asuint64(kd) >> 16;
	203	kd = (double_t)(int32_t)ki;
	204	#else
	205	/* z - kd is in [-1, 1] in non-nearest rounding modes. */
	206	kd = eval_as_double(z + Shift);
	207	ki = asuint64(kd);
	208	kd -= Shift;
	209	#endif
	210	r = x + kd * NegLn2hiN + kd * NegLn2loN;
	211	/* The code assumes 2^-200 < \|xtail\| < 2^-8/N. */
	212	r += xtail;
	213	/* 2^(k/N) ~= scale * (1 + tail). */
	214	idx = 2 * (ki % N);
	215	top = (ki + sign_bias) << (52 - EXP_TABLE_BITS);
	216	tail = asdouble(T[idx]);
	217	/* This is only a valid scale when -1023N < k < 1024N. */
	218	sbits = T[idx + 1] + top;
	219	/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */
	220	/* Evaluation is optimized assuming superscalar pipelined execution. */
	221	r2 = r * r;
	222	/* Without fma the worst case error is 0.25/N ulp larger. */
	223	/* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp. */
	224	tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
	225	if (predict_false(abstop == 0))
	226	return specialcase(tmp, sbits, ki);
	227	scale = asdouble(sbits);
	228	/* Note: tmp == 0 or \|tmp\| > 2^-200 and scale > 2^-739, so there
	229	is no spurious underflow here even without fma. */
	230	return eval_as_double(scale + scale * tmp);
	231	}
	232
	233	/* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is
	234	the bit representation of a non-zero finite floating-point value. */
	235	static inline int checkint(uint64_t iy)
	236	{
	237	int e = iy >> 52 & 0x7ff;
	238	if (e < 0x3ff)
	239	return 0;
	240	if (e > 0x3ff + 52)
	241	return 2;
	242	if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
	243	return 0;
	244	if (iy & (1ULL << (0x3ff + 52 - e)))
	245	return 1;
	246	return 2;
	247	}
	248
	249	/* Returns 1 if input is the bit representation of 0, infinity or nan. */
	250	static inline int zeroinfnan(uint64_t i)
	251	{
	252	return 2 * i - 1 >= 2 * asuint64(INFINITY) - 1;
	253	}
	254
	255	double pow(double x, double y)
	256	{
	257	uint32_t sign_bias = 0;
	258	uint64_t ix, iy;
	259	uint32_t topx, topy;
	260
	261	ix = asuint64(x);
	262	iy = asuint64(y);
	263	topx = top12(x);
	264	topy = top12(y);
	265	if (predict_false(topx - 0x001 >= 0x7ff - 0x001 \|\|
	266	(topy & 0x7ff) - 0x3be >= 0x43e - 0x3be)) {
	267	/* Note: if \|y\| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
	268	and if \|y\| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */
	269	/* Special cases: (x < 0x1p-126 or inf or nan) or
	270	(\|y\| < 0x1p-65 or \|y\| >= 0x1p63 or nan). */
	271	if (predict_false(zeroinfnan(iy))) {
	272	if (2 * iy == 0)
	273	return issignaling_inline(x) ? x + y : 1.0;
	274	if (ix == asuint64(1.0))
	275	return issignaling_inline(y) ? x + y : 1.0;
	276	if (2 * ix > 2 * asuint64(INFINITY) \|\|
	277	2 * iy > 2 * asuint64(INFINITY))
	278	return x + y;
	279	if (2 * ix == 2 * asuint64(1.0))
	280	return 1.0;
	281	if ((2 * ix < 2 * asuint64(1.0)) == !(iy >> 63))
	282	return 0.0; /* \|x\|<1 && y==inf or \|x\|>1 && y==-inf. */
	283	return y * y;
	284	}
	285	if (predict_false(zeroinfnan(ix))) {
	286	double_t x2 = x * x;
	287	if (ix >> 63 && checkint(iy) == 1)
	288	x2 = -x2;
	289	/* Without the barrier some versions of clang hoist the 1/x2 and
	290	thus division by zero exception can be signaled spuriously. */
	291	return iy >> 63 ? fp_barrier(1 / x2) : x2;
	292	}
	293	/* Here x and y are non-zero finite. */
	294	if (ix >> 63) {
	295	/* Finite x < 0. */
	296	int yint = checkint(iy);
	297	if (yint == 0)
	298	return __math_invalid(x);
	299	if (yint == 1)
	300	sign_bias = SIGN_BIAS;
	301	ix &= 0x7fffffffffffffff;
	302	topx &= 0x7ff;
	303	}
	304	if ((topy & 0x7ff) - 0x3be >= 0x43e - 0x3be) {
	305	/* Note: sign_bias == 0 here because y is not odd. */
	306	if (ix == asuint64(1.0))
	307	return 1.0;
	308	if ((topy & 0x7ff) < 0x3be) {
	309	/* \|y\| < 2^-65, x^y ~= 1 + ylog(x). /
	310	if (WANT_ROUNDING)
	311	return ix > asuint64(1.0) ? 1.0 + y :
	312	1.0 - y;
	313	else
	314	return 1.0;
	315	}
	316	return (ix > asuint64(1.0)) == (topy < 0x800) ?
	317	__math_oflow(0) :
	318	__math_uflow(0);
	319	}
	320	if (topx == 0) {
	321	/* Normalize subnormal x so exponent becomes negative. */
	322	ix = asuint64(x * 0x1p52);
	323	ix &= 0x7fffffffffffffff;
	324	ix -= 52ULL << 52;
	325	}
	326	}
	327
	328	double_t lo;
	329	double_t hi = log_inline(ix, &lo);
	330	double_t ehi, elo;
	331	#if __FP_FAST_FMA
	332	ehi = y * hi;
	333	elo = y * lo + __builtin_fma(y, hi, -ehi);
	334	#else
	335	double_t yhi = asdouble(iy & -1ULL << 27);
	336	double_t ylo = y - yhi;
	337	double_t lhi = asdouble(asuint64(hi) & -1ULL << 27);
	338	double_t llo = hi - lhi + lo;
	339	ehi = yhi * lhi;
	340	elo = ylo * lhi + y * llo; /* \|elo\| < \|ehi\| * 2^-25. */
	341	#endif
	342	return exp_inline(ehi, elo, sign_bias);
	343	}

lib/libc/musl/src/math/pow_data.c created+180

...	...	@@ -0,0 +1,180 @@
	1	/*
	2	* Data for the log part of pow.
	3	*
	4	* Copyright (c) 2018, Arm Limited.
	5	* SPDX-License-Identifier: MIT
	6	*/
	7
	8	#include "pow_data.h"
	9
	10	#define N (1 << POW_LOG_TABLE_BITS)
	11
	12	const struct pow_log_data __pow_log_data = {
	13	.ln2hi = 0x1.62e42fefa3800p-1,
	14	.ln2lo = 0x1.ef35793c76730p-45,
	15	.poly = {
	16	// relative error: 0x1.11922ap-70
	17	// in -0x1.6bp-8 0x1.6bp-8
	18	// Coefficients are scaled to match the scaling during evaluation.
	19	-0x1p-1,
	20	0x1.555555555556p-2 * -2,
	21	-0x1.0000000000006p-2 * -2,
	22	0x1.999999959554ep-3 * 4,
	23	-0x1.555555529a47ap-3 * 4,
	24	0x1.2495b9b4845e9p-3 * -8,
	25	-0x1.0002b8b263fc3p-3 * -8,
	26	},
	27	/* Algorithm:
	28
	29	x = 2^k z
	30	log(x) = k ln2 + log(c) + log(z/c)
	31	log(z/c) = poly(z/c - 1)
	32
	33	where z is in [0x1.69555p-1; 0x1.69555p0] which is split into N subintervals
	34	and z falls into the ith one, then table entries are computed as
	35
	36	tab[i].invc = 1/c
	37	tab[i].logc = round(0x1p43*log(c))/0x1p43
	38	tab[i].logctail = (double)(log(c) - logc)
	39
	40	where c is chosen near the center of the subinterval such that 1/c has only a
	41	few precision bits so z/c - 1 is exactly representible as double:
	42
	43	1/c = center < 1 ? round(N/center)/N : round(2*N/center)/N/2
	44
	45	Note: \|z/c - 1\| < 1/N for the chosen c, \|log(c) - logc - logctail\| < 0x1p-97,
	46	the last few bits of logc are rounded away so k*ln2hi + logc has no rounding
	47	error and the interval for z is selected such that near x == 1, where log(x)
	48	is tiny, large cancellation error is avoided in logc + poly(z/c - 1). */
	49	.tab = {
	50	#define A(a, b, c) {a, 0, b, c},
	51	A(0x1.6a00000000000p+0, -0x1.62c82f2b9c800p-2, 0x1.ab42428375680p-48)
	52	A(0x1.6800000000000p+0, -0x1.5d1bdbf580800p-2, -0x1.ca508d8e0f720p-46)
	53	A(0x1.6600000000000p+0, -0x1.5767717455800p-2, -0x1.362a4d5b6506dp-45)
	54	A(0x1.6400000000000p+0, -0x1.51aad872df800p-2, -0x1.684e49eb067d5p-49)
	55	A(0x1.6200000000000p+0, -0x1.4be5f95777800p-2, -0x1.41b6993293ee0p-47)
	56	A(0x1.6000000000000p+0, -0x1.4618bc21c6000p-2, 0x1.3d82f484c84ccp-46)
	57	A(0x1.5e00000000000p+0, -0x1.404308686a800p-2, 0x1.c42f3ed820b3ap-50)
	58	A(0x1.5c00000000000p+0, -0x1.3a64c55694800p-2, 0x1.0b1c686519460p-45)
	59	A(0x1.5a00000000000p+0, -0x1.347dd9a988000p-2, 0x1.5594dd4c58092p-45)
	60	A(0x1.5800000000000p+0, -0x1.2e8e2bae12000p-2, 0x1.67b1e99b72bd8p-45)
	61	A(0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46)
	62	A(0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46)
	63	A(0x1.5400000000000p+0, -0x1.22941fbcf7800p-2, -0x1.65a242853da76p-46)
	64	A(0x1.5200000000000p+0, -0x1.1c898c1699800p-2, -0x1.fafbc68e75404p-46)
	65	A(0x1.5000000000000p+0, -0x1.1675cababa800p-2, 0x1.f1fc63382a8f0p-46)
	66	A(0x1.4e00000000000p+0, -0x1.1058bf9ae4800p-2, -0x1.6a8c4fd055a66p-45)
	67	A(0x1.4c00000000000p+0, -0x1.0a324e2739000p-2, -0x1.c6bee7ef4030ep-47)
	68	A(0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48)
	69	A(0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48)
	70	A(0x1.4800000000000p+0, -0x1.fb9186d5e4000p-3, 0x1.d572aab993c87p-47)
	71	A(0x1.4600000000000p+0, -0x1.ef0adcbdc6000p-3, 0x1.b26b79c86af24p-45)
	72	A(0x1.4400000000000p+0, -0x1.e27076e2af000p-3, -0x1.72f4f543fff10p-46)
	73	A(0x1.4200000000000p+0, -0x1.d5c216b4fc000p-3, 0x1.1ba91bbca681bp-45)
	74	A(0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45)
	75	A(0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45)
	76	A(0x1.3e00000000000p+0, -0x1.bc286742d9000p-3, 0x1.94eb0318bb78fp-46)
	77	A(0x1.3c00000000000p+0, -0x1.af3c94e80c000p-3, 0x1.a4e633fcd9066p-52)
	78	A(0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45)
	79	A(0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45)
	80	A(0x1.3800000000000p+0, -0x1.9525a9cf45000p-3, -0x1.ad1d904c1d4e3p-45)
	81	A(0x1.3600000000000p+0, -0x1.87fa06520d000p-3, 0x1.bbdbf7fdbfa09p-45)
	82	A(0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45)
	83	A(0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45)
	84	A(0x1.3200000000000p+0, -0x1.6d60fe719d000p-3, -0x1.0e46aa3b2e266p-46)
	85	A(0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46)
	86	A(0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46)
	87	A(0x1.2e00000000000p+0, -0x1.526e5e3a1b000p-3, -0x1.0de8b90075b8fp-45)
	88	A(0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46)
	89	A(0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46)
	90	A(0x1.2a00000000000p+0, -0x1.371fc201e9000p-3, 0x1.178864d27543ap-48)
	91	A(0x1.2800000000000p+0, -0x1.29552f81ff000p-3, -0x1.48d301771c408p-45)
	92	A(0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45)
	93	A(0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45)
	94	A(0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47)
	95	A(0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47)
	96	A(0x1.2200000000000p+0, -0x1.fec9131dbe000p-4, -0x1.575545ca333f2p-45)
	97	A(0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45)
	98	A(0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45)
	99	A(0x1.1e00000000000p+0, -0x1.c5e548f5bc000p-4, -0x1.d0c57585fbe06p-46)
	100	A(0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45)
	101	A(0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45)
	102	A(0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46)
	103	A(0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46)
	104	A(0x1.1800000000000p+0, -0x1.6f0d28ae56000p-4, -0x1.69737c93373dap-45)
	105	A(0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46)
	106	A(0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46)
	107	A(0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45)
	108	A(0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45)
	109	A(0x1.1200000000000p+0, -0x1.16536eea38000p-4, 0x1.47c5e768fa309p-46)
	110	A(0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45)
	111	A(0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45)
	112	A(0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46)
	113	A(0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46)
	114	A(0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45)
	115	A(0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45)
	116	A(0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48)
	117	A(0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48)
	118	A(0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45)
	119	A(0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45)
	120	A(0x1.0600000000000p+0, -0x1.7b91b07d58000p-6, -0x1.88d5493faa639p-45)
	121	A(0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50)
	122	A(0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50)
	123	A(0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46)
	124	A(0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46)
	125	A(0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0)
	126	A(0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0)
	127	A(0x1.fc00000000000p-1, 0x1.0101575890000p-7, -0x1.0c76b999d2be8p-46)
	128	A(0x1.f800000000000p-1, 0x1.0205658938000p-6, -0x1.3dc5b06e2f7d2p-45)
	129	A(0x1.f400000000000p-1, 0x1.8492528c90000p-6, -0x1.aa0ba325a0c34p-45)
	130	A(0x1.f000000000000p-1, 0x1.0415d89e74000p-5, 0x1.111c05cf1d753p-47)
	131	A(0x1.ec00000000000p-1, 0x1.466aed42e0000p-5, -0x1.c167375bdfd28p-45)
	132	A(0x1.e800000000000p-1, 0x1.894aa149fc000p-5, -0x1.97995d05a267dp-46)
	133	A(0x1.e400000000000p-1, 0x1.ccb73cdddc000p-5, -0x1.a68f247d82807p-46)
	134	A(0x1.e200000000000p-1, 0x1.eea31c006c000p-5, -0x1.e113e4fc93b7bp-47)
	135	A(0x1.de00000000000p-1, 0x1.1973bd1466000p-4, -0x1.5325d560d9e9bp-45)
	136	A(0x1.da00000000000p-1, 0x1.3bdf5a7d1e000p-4, 0x1.cc85ea5db4ed7p-45)
	137	A(0x1.d600000000000p-1, 0x1.5e95a4d97a000p-4, -0x1.c69063c5d1d1ep-45)
	138	A(0x1.d400000000000p-1, 0x1.700d30aeac000p-4, 0x1.c1e8da99ded32p-49)
	139	A(0x1.d000000000000p-1, 0x1.9335e5d594000p-4, 0x1.3115c3abd47dap-45)
	140	A(0x1.cc00000000000p-1, 0x1.b6ac88dad6000p-4, -0x1.390802bf768e5p-46)
	141	A(0x1.ca00000000000p-1, 0x1.c885801bc4000p-4, 0x1.646d1c65aacd3p-45)
	142	A(0x1.c600000000000p-1, 0x1.ec739830a2000p-4, -0x1.dc068afe645e0p-45)
	143	A(0x1.c400000000000p-1, 0x1.fe89139dbe000p-4, -0x1.534d64fa10afdp-45)
	144	A(0x1.c000000000000p-1, 0x1.1178e8227e000p-3, 0x1.1ef78ce2d07f2p-45)
	145	A(0x1.be00000000000p-1, 0x1.1aa2b7e23f000p-3, 0x1.ca78e44389934p-45)
	146	A(0x1.ba00000000000p-1, 0x1.2d1610c868000p-3, 0x1.39d6ccb81b4a1p-47)
	147	A(0x1.b800000000000p-1, 0x1.365fcb0159000p-3, 0x1.62fa8234b7289p-51)
	148	A(0x1.b400000000000p-1, 0x1.4913d8333b000p-3, 0x1.5837954fdb678p-45)
	149	A(0x1.b200000000000p-1, 0x1.527e5e4a1b000p-3, 0x1.633e8e5697dc7p-45)
	150	A(0x1.ae00000000000p-1, 0x1.6574ebe8c1000p-3, 0x1.9cf8b2c3c2e78p-46)
	151	A(0x1.ac00000000000p-1, 0x1.6f0128b757000p-3, -0x1.5118de59c21e1p-45)
	152	A(0x1.aa00000000000p-1, 0x1.7898d85445000p-3, -0x1.c661070914305p-46)
	153	A(0x1.a600000000000p-1, 0x1.8beafeb390000p-3, -0x1.73d54aae92cd1p-47)
	154	A(0x1.a400000000000p-1, 0x1.95a5adcf70000p-3, 0x1.7f22858a0ff6fp-47)
	155	A(0x1.a000000000000p-1, 0x1.a93ed3c8ae000p-3, -0x1.8724350562169p-45)
	156	A(0x1.9e00000000000p-1, 0x1.b31d8575bd000p-3, -0x1.c358d4eace1aap-47)
	157	A(0x1.9c00000000000p-1, 0x1.bd087383be000p-3, -0x1.d4bc4595412b6p-45)
	158	A(0x1.9a00000000000p-1, 0x1.c6ffbc6f01000p-3, -0x1.1ec72c5962bd2p-48)
	159	A(0x1.9600000000000p-1, 0x1.db13db0d49000p-3, -0x1.aff2af715b035p-45)
	160	A(0x1.9400000000000p-1, 0x1.e530effe71000p-3, 0x1.212276041f430p-51)
	161	A(0x1.9200000000000p-1, 0x1.ef5ade4dd0000p-3, -0x1.a211565bb8e11p-51)
	162	A(0x1.9000000000000p-1, 0x1.f991c6cb3b000p-3, 0x1.bcbecca0cdf30p-46)
	163	A(0x1.8c00000000000p-1, 0x1.07138604d5800p-2, 0x1.89cdb16ed4e91p-48)
	164	A(0x1.8a00000000000p-1, 0x1.0c42d67616000p-2, 0x1.7188b163ceae9p-45)
	165	A(0x1.8800000000000p-1, 0x1.1178e8227e800p-2, -0x1.c210e63a5f01cp-45)
	166	A(0x1.8600000000000p-1, 0x1.16b5ccbacf800p-2, 0x1.b9acdf7a51681p-45)
	167	A(0x1.8400000000000p-1, 0x1.1bf99635a6800p-2, 0x1.ca6ed5147bdb7p-45)
	168	A(0x1.8200000000000p-1, 0x1.214456d0eb800p-2, 0x1.a87deba46baeap-47)
	169	A(0x1.7e00000000000p-1, 0x1.2bef07cdc9000p-2, 0x1.a9cfa4a5004f4p-45)
	170	A(0x1.7c00000000000p-1, 0x1.314f1e1d36000p-2, -0x1.8e27ad3213cb8p-45)
	171	A(0x1.7a00000000000p-1, 0x1.36b6776be1000p-2, 0x1.16ecdb0f177c8p-46)
	172	A(0x1.7800000000000p-1, 0x1.3c25277333000p-2, 0x1.83b54b606bd5cp-46)
	173	A(0x1.7600000000000p-1, 0x1.419b423d5e800p-2, 0x1.8e436ec90e09dp-47)
	174	A(0x1.7400000000000p-1, 0x1.4718dc271c800p-2, -0x1.f27ce0967d675p-45)
	175	A(0x1.7200000000000p-1, 0x1.4c9e09e173000p-2, -0x1.e20891b0ad8a4p-45)
	176	A(0x1.7000000000000p-1, 0x1.522ae0738a000p-2, 0x1.ebe708164c759p-45)
	177	A(0x1.6e00000000000p-1, 0x1.57bf753c8d000p-2, 0x1.fadedee5d40efp-46)
	178	A(0x1.6c00000000000p-1, 0x1.5d5bddf596000p-2, -0x1.a0b2a08a465dcp-47)
	179	},
	180	};

lib/libc/musl/src/math/pow_data.h created+22

...	...	@@ -0,0 +1,22 @@
	1	/*
	2	* Copyright (c) 2018, Arm Limited.
	3	* SPDX-License-Identifier: MIT
	4	*/
	5	#ifndef _POW_DATA_H
	6	#define _POW_DATA_H
	7
	8	#include <features.h>
	9
	10	#define POW_LOG_TABLE_BITS 7
	11	#define POW_LOG_POLY_ORDER 8
	12	extern hidden const struct pow_log_data {
	13	double ln2hi;
	14	double ln2lo;
	15	double poly[POW_LOG_POLY_ORDER - 1]; /* First coefficient is 1. */
	16	/* Note: the pad field is unused, but allows slightly faster indexing. */
	17	struct {
	18	double invc, pad, logc, logctail;
	19	} tab[1 << POW_LOG_TABLE_BITS];
	20	} __pow_log_data;
	21
	22	#endif

src/libs/musl.zig+3

...	...	@@ -816,6 +816,7 @@ const src_files = [_][]const u8{
816	816	"musl/src/math/exp10l.c",
817	817	"musl/src/math/exp2f_data.c",
818	818	"musl/src/math/exp2l.c",
	819	"musl/src/math/exp_data.c",
819	820	"musl/src/math/expl.c",
820	821	"musl/src/math/expm1.c",
821	822	"musl/src/math/expm1f.c",
...	...	@@ -925,6 +926,8 @@ const src_files = [_][]const u8{
925	926	"musl/src/math/nexttowardf.c",
926	927	"musl/src/math/nexttowardl.c",
927	928	"musl/src/math/__polevll.c",
	929	"musl/src/math/pow.c",
	930	"musl/src/math/pow_data.c",
928	931	"musl/src/math/powerpc64/fma.c",
929	932	"musl/src/math/powerpc64/fmaf.c",
930	933	"musl/src/math/powerpc64/lround.c",

src/libs/wasi_libc.zig+3

...	...	@@ -687,6 +687,7 @@ const libc_top_half_src_files = [_][]const u8{
687	687	"musl/src/math/exp10l.c",
688	688	"musl/src/math/exp2f_data.c",
689	689	"musl/src/math/exp2l.c",
	690	"musl/src/math/exp_data.c",
690	691	"musl/src/math/expl.c",
691	692	"musl/src/math/expm1.c",
692	693	"musl/src/math/expm1f.c",
...	...	@@ -751,6 +752,8 @@ const libc_top_half_src_files = [_][]const u8{
751	752	"musl/src/math/nexttowardf.c",
752	753	"musl/src/math/nexttowardl.c",
753	754	"musl/src/math/__polevll.c",
	755	"musl/src/math/pow.c",
	756	"musl/src/math/pow_data.c",
754	757	"musl/src/math/powf.c",
755	758	"musl/src/math/powf_data.c",
756	759	"musl/src/math/remainder.c",