aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754')
-rw-r--r--sysdeps/ieee754/dbl-64/s_erfc.c1
-rw-r--r--sysdeps/ieee754/float128/s_erfcf128.c1
-rw-r--r--sysdeps/ieee754/flt-32/e_gammaf_r.c2
-rw-r--r--sysdeps/ieee754/flt-32/e_lgammaf_r.c575
-rw-r--r--sysdeps/ieee754/flt-32/k_tanf.c102
-rw-r--r--sysdeps/ieee754/flt-32/lgamma_negf.c283
-rw-r--r--sysdeps/ieee754/flt-32/math_config.h27
-rw-r--r--sysdeps/ieee754/flt-32/s_cbrtf.c136
-rw-r--r--sysdeps/ieee754/flt-32/s_erfcf.c187
-rw-r--r--sysdeps/ieee754/flt-32/s_erff.c470
-rw-r--r--sysdeps/ieee754/flt-32/s_expm1f.c2
-rw-r--r--sysdeps/ieee754/flt-32/s_tanf.c224
-rw-r--r--sysdeps/ieee754/ldbl-128/s_erfcl.c1
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/s_erfcl.c1
-rw-r--r--sysdeps/ieee754/ldbl-96/s_erfcl.c1
15 files changed, 1068 insertions, 945 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_erfc.c b/sysdeps/ieee754/dbl-64/s_erfc.c
new file mode 100644
index 0000000..95d17c8
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/s_erfc.c
@@ -0,0 +1 @@
+/* Not required. */
diff --git a/sysdeps/ieee754/float128/s_erfcf128.c b/sysdeps/ieee754/float128/s_erfcf128.c
new file mode 100644
index 0000000..95d17c8
--- /dev/null
+++ b/sysdeps/ieee754/float128/s_erfcf128.c
@@ -0,0 +1 @@
+/* Not required. */
diff --git a/sysdeps/ieee754/flt-32/e_gammaf_r.c b/sysdeps/ieee754/flt-32/e_gammaf_r.c
index 6b1f95d..66e8cae 100644
--- a/sysdeps/ieee754/flt-32/e_gammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_gammaf_r.c
@@ -140,7 +140,7 @@ __ieee754_gammaf_r (float x, int *signgamp)
};
double m = z - 0x1.7p+1;
- double i = roundeven (m);
+ double i = roundeven_finite (m);
double step = copysign (1.0, i);
double d = m - i, d2 = d * d, d4 = d2 * d2, d8 = d4 * d4;
double f = (c[0] + d * c[1]) + d2 * (c[2] + d * c[3])
diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
index a1a3a60..75ec25f 100644
--- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
@@ -1,247 +1,366 @@
-/* e_lgammaf_r.c -- float version of e_lgamma_r.c.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+/* Correctly-rounded logarithm of the absolute value of the gamma function
+ for binary32 value.
+Copyright (c) 2023, 2024 Alexei Sibidanov.
+
+This file is part of the CORE-MATH project
+project (file src/binary32/lgamma/lgammaf.c, revision bc385c2).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+/* Changes with respect to the original CORE-MATH code:
+ - removed the dealing with errno
+ (this is done in the wrapper math/w_lgammaf_compat2.c).
+ - usage of math_narrow_eval to deal with underflow/overflow.
+ - deal with signamp. */
+
+#include <array_length.h>
+#include <stdint.h>
#include <math.h>
-#include <math-narrow-eval.h>
-#include <math_private.h>
-#include <libc-diag.h>
#include <libm-alias-finite.h>
+#include <limits.h>
+#include <math-narrow-eval.h>
+#include "math_config.h"
-static const float
-two23= 8.3886080000e+06, /* 0x4b000000 */
-half= 5.0000000000e-01, /* 0x3f000000 */
-one = 1.0000000000e+00, /* 0x3f800000 */
-pi = 3.1415927410e+00, /* 0x40490fdb */
-a0 = 7.7215664089e-02, /* 0x3d9e233f */
-a1 = 3.2246702909e-01, /* 0x3ea51a66 */
-a2 = 6.7352302372e-02, /* 0x3d89f001 */
-a3 = 2.0580807701e-02, /* 0x3ca89915 */
-a4 = 7.3855509982e-03, /* 0x3bf2027e */
-a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */
-a6 = 1.1927076848e-03, /* 0x3a9c54a1 */
-a7 = 5.1006977446e-04, /* 0x3a05b634 */
-a8 = 2.2086278477e-04, /* 0x39679767 */
-a9 = 1.0801156895e-04, /* 0x38e28445 */
-a10 = 2.5214456400e-05, /* 0x37d383a2 */
-a11 = 4.4864096708e-05, /* 0x383c2c75 */
-tc = 1.4616321325e+00, /* 0x3fbb16c3 */
-tf = -1.2148628384e-01, /* 0xbdf8cdcd */
-/* tt = -(tail of tf) */
-tt = 6.6971006518e-09, /* 0x31e61c52 */
-t0 = 4.8383611441e-01, /* 0x3ef7b95e */
-t1 = -1.4758771658e-01, /* 0xbe17213c */
-t2 = 6.4624942839e-02, /* 0x3d845a15 */
-t3 = -3.2788541168e-02, /* 0xbd064d47 */
-t4 = 1.7970675603e-02, /* 0x3c93373d */
-t5 = -1.0314224288e-02, /* 0xbc28fcfe */
-t6 = 6.1005386524e-03, /* 0x3bc7e707 */
-t7 = -3.6845202558e-03, /* 0xbb7177fe */
-t8 = 2.2596477065e-03, /* 0x3b141699 */
-t9 = -1.4034647029e-03, /* 0xbab7f476 */
-t10 = 8.8108185446e-04, /* 0x3a66f867 */
-t11 = -5.3859531181e-04, /* 0xba0d3085 */
-t12 = 3.1563205994e-04, /* 0x39a57b6b */
-t13 = -3.1275415677e-04, /* 0xb9a3f927 */
-t14 = 3.3552918467e-04, /* 0x39afe9f7 */
-u0 = -7.7215664089e-02, /* 0xbd9e233f */
-u1 = 6.3282704353e-01, /* 0x3f2200f4 */
-u2 = 1.4549225569e+00, /* 0x3fba3ae7 */
-u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */
-u4 = 2.2896373272e-01, /* 0x3e6a7578 */
-u5 = 1.3381091878e-02, /* 0x3c5b3c5e */
-v1 = 2.4559779167e+00, /* 0x401d2ebe */
-v2 = 2.1284897327e+00, /* 0x4008392d */
-v3 = 7.6928514242e-01, /* 0x3f44efdf */
-v4 = 1.0422264785e-01, /* 0x3dd572af */
-v5 = 3.2170924824e-03, /* 0x3b52d5db */
-s0 = -7.7215664089e-02, /* 0xbd9e233f */
-s1 = 2.1498242021e-01, /* 0x3e5c245a */
-s2 = 3.2577878237e-01, /* 0x3ea6cc7a */
-s3 = 1.4635047317e-01, /* 0x3e15dce6 */
-s4 = 2.6642270386e-02, /* 0x3cda40e4 */
-s5 = 1.8402845599e-03, /* 0x3af135b4 */
-s6 = 3.1947532989e-05, /* 0x3805ff67 */
-r1 = 1.3920053244e+00, /* 0x3fb22d3b */
-r2 = 7.2193557024e-01, /* 0x3f38d0c5 */
-r3 = 1.7193385959e-01, /* 0x3e300f6e */
-r4 = 1.8645919859e-02, /* 0x3c98bf54 */
-r5 = 7.7794247773e-04, /* 0x3a4beed6 */
-r6 = 7.3266842264e-06, /* 0x36f5d7bd */
-w0 = 4.1893854737e-01, /* 0x3ed67f1d */
-w1 = 8.3333335817e-02, /* 0x3daaaaab */
-w2 = -2.7777778450e-03, /* 0xbb360b61 */
-w3 = 7.9365057172e-04, /* 0x3a500cfd */
-w4 = -5.9518753551e-04, /* 0xba1c065c */
-w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
-w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
-
-static const float zero= 0.0000000000e+00;
-
-static float
-sin_pif(float x)
+static double
+as_r7 (double x, const double *c)
{
- float y,z;
- int n,ix;
-
- GET_FLOAT_WORD(ix,x);
- ix &= 0x7fffffff;
-
- if(ix<0x3e800000) return __sinf (pi*x);
- y = -x; /* x is assume negative */
-
- /*
- * argument reduction, make sure inexact flag not raised if input
- * is an integer
- */
- z = floorf(y);
- if(z!=y) { /* inexact anyway */
- y *= (float)0.5;
- y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */
- n = (int) (y*(float)4.0);
- } else {
- if(ix>=0x4b800000) {
- y = zero; n = 0; /* y must be even */
- } else {
- if(ix<0x4b000000) z = y+two23; /* exact */
- GET_FLOAT_WORD(n,z);
- n &= 1;
- y = n;
- n<<= 2;
- }
- }
- switch (n) {
- case 0: y = __sinf (pi*y); break;
- case 1:
- case 2: y = __cosf (pi*((float)0.5-y)); break;
- case 3:
- case 4: y = __sinf (pi*(one-y)); break;
- case 5:
- case 6: y = -__cosf (pi*(y-(float)1.5)); break;
- default: y = __sinf (pi*(y-(float)2.0)); break;
- }
- return -y;
+ return (((x - c[0]) * (x - c[1])) * ((x - c[2]) * (x - c[3])))
+ * (((x - c[4]) * (x - c[5])) * ((x - c[6])));
}
+static double
+as_r8 (double x, const double *c)
+{
+ return (((x - c[0]) * (x - c[1])) * ((x - c[2]) * (x - c[3])))
+ * (((x - c[4]) * (x - c[5])) * ((x - c[6]) * (x - c[7])));
+}
+
+static double
+as_sinpi (double x)
+{
+ static const double c[] =
+ {
+ 0x1p+2, -0x1.de9e64df22ea4p+1, 0x1.472be122401f8p+0,
+ -0x1.d4fcd82df91bp-3, 0x1.9f05c97e0aab2p-6, -0x1.f3091c427b611p-10,
+ 0x1.b22c9bfdca547p-14, -0x1.15484325ef569p-18
+ };
+ x -= 0.5;
+ double x2 = x * x, x4 = x2 * x2, x8 = x4 * x4;
+ return (0.25 - x2)
+ * ((c[0] + x2 * c[1]) + x4 * (c[2] + x2 * c[3])
+ + x8 * ((c[4] + x2 * c[5]) + x4 * (c[6] + x2 * c[7])));
+}
+
+static double
+as_ln (double x)
+{
+ uint64_t t = asuint64 (x);
+ int e = (t >> 52) - 0x3ff;
+ static const double c[] =
+ {
+ 0x1.fffffffffff24p-1, -0x1.ffffffffd1d67p-2, 0x1.55555537802dep-2,
+ -0x1.ffffeca81b866p-3, 0x1.999611761d772p-3, -0x1.54f3e581b61bfp-3,
+ 0x1.1e642b4cb5143p-3, -0x1.9115a5af1e1edp-4
+ };
+ static const double il[] =
+ {
+ 0x1.59caeec280116p-57, 0x1.f0a30c01162aap-5, 0x1.e27076e2af2ebp-4,
+ 0x1.5ff3070a793d6p-3, 0x1.c8ff7c79a9a2p-3, 0x1.1675cababa60fp-2,
+ 0x1.4618bc21c5ec2p-2, 0x1.739d7f6bbd007p-2, 0x1.9f323ecbf984dp-2,
+ 0x1.c8ff7c79a9a21p-2, 0x1.f128f5faf06ecp-2, 0x1.0be72e4252a83p-1,
+ 0x1.1e85f5e7040d1p-1, 0x1.307d7334f10bep-1, 0x1.41d8fe84672afp-1,
+ 0x1.52a2d265bc5abp-1
+ };
+ static const double ix[] =
+ {
+ 0x1p+0, 0x1.e1e1e1e1e1e1ep-1, 0x1.c71c71c71c71cp-1,
+ 0x1.af286bca1af28p-1, 0x1.999999999999ap-1, 0x1.8618618618618p-1,
+ 0x1.745d1745d1746p-1, 0x1.642c8590b2164p-1, 0x1.5555555555555p-1,
+ 0x1.47ae147ae147bp-1, 0x1.3b13b13b13b14p-1, 0x1.2f684bda12f68p-1,
+ 0x1.2492492492492p-1, 0x1.1a7b9611a7b96p-1, 0x1.1111111111111p-1,
+ 0x1.0842108421084p-1
+ };
+ int i = (t >> 48) & 0xf;
+ t = (t & (~UINT64_C(0) >> 12)) | (INT64_C(0x3ff) << 52);
+ double z = ix[i] * asdouble (t) - 1;
+ double z2 = z * z, z4 = z2 * z2;
+ return e * 0x1.62e42fefa39efp-1 + il[i]
+ + z * ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3])
+ + z4 * ((c[4] + z * c[5]) + z2 * (c[6] + z * c[7])));
+}
float
-__ieee754_lgammaf_r(float x, int *signgamp)
+__ieee754_lgammaf_r (float x, int *signgamp)
{
- float t,y,z,nadj,p,p1,p2,p3,q,r,w;
- int i,hx,ix;
-
- GET_FLOAT_WORD(hx,x);
-
- /* purge off +-inf, NaN, +-0, and negative arguments */
- *signgamp = 1;
- ix = hx&0x7fffffff;
- if(__builtin_expect(ix>=0x7f800000, 0)) return x*x;
- if(__builtin_expect(ix==0, 0))
- {
- if (hx < 0)
- *signgamp = -1;
- return one/fabsf(x);
- }
- if(__builtin_expect(ix<0x30800000, 0)) {
- /* |x|<2**-30, return -log(|x|) */
- if(hx<0) {
- *signgamp = -1;
- return -__ieee754_logf(-x);
- } else return -__ieee754_logf(x);
+ static const struct
+ {
+ float x;
+ float f;
+ float df;
+ } tb[] = {
+ { -0x1.efc2a2p+14, -0x1.222dbcp+18, -0x1p-7 },
+ { -0x1.627346p+7, -0x1.73235ep+9, -0x1p-16 },
+ { -0x1.08b14p+4, -0x1.f0cbe6p+4, -0x1p-21 },
+ { -0x1.69d628p+3, -0x1.0eac2ap+4, -0x1p-21 },
+ { -0x1.904902p+2, -0x1.65532cp+2, 0x1p-23 },
+ { -0x1.9272d2p+1, -0x1.170b98p-8, 0x1p-33 },
+ { -0x1.625edap+1, 0x1.6a6c4ap-5, -0x1p-30 },
+ { -0x1.5fc2aep+1, 0x1.c0a484p-11, -0x1p-36 },
+ { -0x1.5fb43ep+1, 0x1.5b697p-17, 0x1p-42 },
+ { -0x1.5fa20cp+1, -0x1.132f7ap-10, 0x1p-35 },
+ { -0x1.580c1ep+1, -0x1.5787c6p-4, 0x1p-29 },
+ { -0x1.3a7fcap+1, -0x1.e4cf24p-24, -0x1p-49 },
+ { -0x1.c2f04p-30, 0x1.43a6f6p+4, 0x1p-21 },
+ { -0x1.ade594p-30, 0x1.446ab2p+4, -0x1p-21 },
+ { -0x1.437e74p-40, 0x1.b7dec2p+4, -0x1p-21 },
+ { -0x1.d85bfep-43, 0x1.d31592p+4, -0x1p-21 },
+ { -0x1.f51c8ep-49, 0x1.0a572ap+5, -0x1p-20 },
+ { -0x1.108a5ap-66, 0x1.6d7b18p+5, -0x1p-20 },
+ { -0x1.ecf3fep-73, 0x1.8f8e5ap+5, -0x1p-20 },
+ { -0x1.25cb66p-123, 0x1.547a44p+6, -0x1p-19 },
+ { 0x1.ecf3fep-73, 0x1.8f8e5ap+5, -0x1p-20 },
+ { 0x1.108a5ap-66, 0x1.6d7b18p+5, -0x1p-20 },
+ { 0x1.a68bbcp-42, 0x1.c9c6e8p+4, 0x1p-21 },
+ { 0x1.ddfd06p-12, 0x1.ec5ba8p+2, -0x1p-23 },
+ { 0x1.f8a754p-9, 0x1.63acc2p+2, 0x1p-23 },
+ { 0x1.8d16b2p+5, 0x1.1e4b4ep+7, 0x1p-18 },
+ { 0x1.359e0ep+10, 0x1.d9ad02p+12, -0x1p-13 },
+ { 0x1.a82a2cp+13, 0x1.c38036p+16, 0x1p-9 },
+ { 0x1.62c646p+14, 0x1.9075bep+17, -0x1p-8 },
+ { 0x1.7f298p+31, 0x1.f44946p+35, -0x1p+10 },
+ { 0x1.a45ea4p+33, 0x1.25dcbcp+38, -0x1p+13 },
+ { 0x1.f9413ep+76, 0x1.9d5ab4p+82, -0x1p+57 },
+ { 0x1.dcbbaap+99, 0x1.fc5772p+105, 0x1p+80 },
+ { 0x1.58ace8p+112, 0x1.9e4f66p+118, -0x1p+93 },
+ { 0x1.87bdfp+115, 0x1.e465aep+121, 0x1p+96 },
+ };
+
+ float fx = floor (x);
+ float ax = fabsf (x);
+ uint32_t t = asuint (ax);
+ if (__glibc_unlikely (t >= (0xffu << 23)))
+ {
+ *signgamp = 1;
+ if (t == (0xffu << 23))
+ return INFINITY;
+ return x + x; /* nan */
+ }
+ if (__glibc_unlikely (fx == x))
+ {
+ if (x <= 0.0f)
+ {
+ *signgamp = asuint (x) >> 31 ? -1 : 1;
+ return 1.0f / 0.0f;
}
- if(hx<0) {
- if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */
- return fabsf (x)/zero;
- if (ix > 0x40000000 /* X < 2.0f. */
- && ix < 0x41700000 /* X > -15.0f. */)
- return __lgamma_negf (x, signgamp);
- t = sin_pif(x);
- if(t==zero) return one/fabsf(t); /* -integer */
- nadj = __ieee754_logf(pi/fabsf(t*x));
- if(t<zero) *signgamp = -1;
- x = -x;
+ if (x == 1.0f || x == 2.0f)
+ {
+ *signgamp = 1;
+ return 0.0f;
}
+ }
+
+ /* Check the value of fx to avoid a spurious invalid exception.
+ Note that for a binary32 |x| >= 2^23, x is necessarily an integer,
+ and we already dealed with negative integers, thus now:
+ -2^23 < x < +Inf and x is not a negative integer nor 0, 1, 2. */
+ if (__glibc_likely (fx >= 0))
+ *signgamp = 1;
+ else
+ /* gamma(x) is negative in (-2n-1,-2n), thus when fx is odd. */
+ *signgamp = 1 - ((((int) fx) & 1) << 1);
- /* purge off 1 and 2 */
- if (ix==0x3f800000||ix==0x40000000) r = 0;
- /* for x < 2.0 */
- else if(ix<0x40000000) {
- if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
- r = -__ieee754_logf(x);
- if(ix>=0x3f3b4a20) {y = one-x; i= 0;}
- else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;}
- else {y = x; i=2;}
- } else {
- r = zero;
- if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */
- else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */
- else {y=x-one;i=2;}
+ double z = ax, f;
+ if (__glibc_unlikely (ax < 0x1.52p-1f))
+ {
+ static const double rn[] =
+ {
+ -0x1.505bdf4b65acp+4, -0x1.51c80eb47e068p+2,
+ 0x1.0000000007cb8p+0, -0x1.4ac529250a1fcp+1,
+ -0x1.a8c99dbe1621ap+0, -0x1.4abdcc74115eap+0,
+ -0x1.1b87fe5a5b923p+0, -0x1.05b8a4d47ff64p+0
+ };
+ const double c0 = 0x1.0fc0fad268c4dp+2;
+ static const double rd[] =
+ {
+ -0x1.4db2cfe9a5265p+5, -0x1.062e99d1c4f27p+3,
+ -0x1.c81bc2ecf25f6p+1, -0x1.108e55c10091bp+1,
+ -0x1.7dd25af0b83d4p+0, -0x1.36bf1880125fcp+0,
+ -0x1.1379fc8023d9cp+0, -0x1.03712e41525d2p+0
+ };
+ double s = x;
+ f = (c0 * s) * as_r8 (s, rn) / as_r8 (s, rd) - as_ln (z);
+ }
+ else
+ {
+ if (ax > 0x1.afc1ap+1f)
+ {
+ if (__glibc_unlikely (x > 0x1.895f1cp+121f))
+ return math_narrow_eval (0x1p127f * 0x1p127f);
+ /* |x|>=2**23, must be -integer */
+ if (__glibc_unlikely (x < 0.0f && ax > 0x1p+23))
+ return ax / 0.0f;
+ double lz = as_ln (z);
+ f = (z - 0.5) * (lz - 1) + 0x1.acfe390c97d69p-2;
+ if (ax < 0x1.0p+20f)
+ {
+ double iz = 1.0 / z, iz2 = iz * iz;
+ if (ax > 1198.0f)
+ f += iz * (1. / 12.);
+ else if (ax > 0x1.279a7p+6f)
+ {
+ static const double c[] =
+ {
+ 0x1.555555547fbadp-4, -0x1.6c0fd270c465p-9
+ };
+ f += iz * (c[0] + iz2 * c[1]);
+ }
+ else if (ax > 0x1.555556p+3f)
+ {
+ static const double c[] =
+ {
+ 0x1.555555554de0bp-4, -0x1.6c16bdc45944fp-9,
+ 0x1.a0077f300ecb3p-11, -0x1.2e9cfff3b29c2p-11
+ };
+ double iz4 = iz2 * iz2;
+ f += iz * ((c[0] + iz2 * c[1]) + iz4 * (c[2] + iz2 * c[3]));
+ }
+ else
+ {
+ static const double c[] =
+ {
+ 0x1.5555555551286p-4, -0x1.6c16c0e7c4cf4p-9,
+ 0x1.a0193267fe6f2p-11, -0x1.37e87ec19cb45p-11,
+ 0x1.b40011dfff081p-11, -0x1.c16c8946b19b6p-10,
+ 0x1.e9f47ace150d8p-9, -0x1.4f5843a71a338p-8
+ };
+ double iz4 = iz2 * iz2, iz8 = iz4 * iz4;
+ double p = ((c[0] + iz2 * c[1]) + iz4 * (c[2] + iz2 * c[3]))
+ + iz8 * ((c[4] + iz2 * c[5])
+ + iz4 * (c[6] + iz2 * c[7]));
+ f += iz * p;
+ }
}
- switch(i) {
- case 0:
- z = y*y;
- p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
- p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
- p = y*p1+p2;
- r += (p-(float)0.5*y); break;
- case 1:
- z = y*y;
- w = z*y;
- p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
- p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
- p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
- p = z*p1-(tt-w*(p2+y*p3));
- r += (tf + p); break;
- case 2:
- p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
- p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
- r += (-(float)0.5*y + p1/p2);
+ if (x < 0.0f)
+ {
+ f = 0x1.250d048e7a1bdp+0 - f - lz;
+ double lp = as_ln (as_sinpi (x - fx));
+ f -= lp;
}
}
- else if(ix<0x41000000) { /* x < 8.0 */
- i = (int)x;
- t = zero;
- y = x-(float)i;
- p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
- q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
- r = half*y+p/q;
- z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
- switch(i) {
- case 7: z *= (y+(float)6.0); /* FALLTHRU */
- case 6: z *= (y+(float)5.0); /* FALLTHRU */
- case 5: z *= (y+(float)4.0); /* FALLTHRU */
- case 4: z *= (y+(float)3.0); /* FALLTHRU */
- case 3: z *= (y+(float)2.0); /* FALLTHRU */
- r += __ieee754_logf(z); break;
+ else
+ {
+ static const double rn[] =
+ {
+ -0x1.667923ff14df7p+5, -0x1.2d35f25ad8f64p+3,
+ -0x1.b8c9eab9d5bd3p+1, -0x1.7a4a97f494127p+0,
+ -0x1.3a6c8295b4445p-1, -0x1.da44e8b810024p-3,
+ -0x1.9061e81c77e4ap-5
+ };
+ if (x < 0.0f)
+ {
+ int ni = floorf (-2 * x);
+ if ((ni & 1) == 0 && ni == -2 * x)
+ return 1.0f / 0.0f;
+ }
+ const double c0 = 0x1.3cc0e6a0106b3p+2;
+ static const double rd[] =
+ {
+ -0x1.491a899e84c52p+6, -0x1.d202961b9e098p+3,
+ -0x1.4ced68c631ed6p+2, -0x1.2589eedf40738p+1,
+ -0x1.1302e3337271p+0, -0x1.c36b802f26dffp-2,
+ -0x1.3ded448acc39dp-3, -0x1.bffc491078eafp-6
+ };
+ f = (z - 1) * (z - 2) * c0 * as_r7 (z, rn) / as_r8 (z, rd);
+ if (x < 0.0f)
+ {
+ if (__glibc_unlikely (t < 0x40301b93u && t > 0x402f95c2u))
+ {
+ double h = (x + 0x1.5fb410a1bd901p+1)
+ - 0x1.a19a96d2e6f85p-54;
+ double h2 = h * h;
+ double h4 = h2 * h2;
+ static const double c[] =
+ {
+ -0x1.ea12da904b18cp+0, 0x1.3267f3c265a54p+3,
+ -0x1.4185ac30cadb3p+4, 0x1.f504accc3f2e4p+5,
+ -0x1.8588444c679b4p+7, 0x1.43740491dc22p+9,
+ -0x1.12400ea23f9e6p+11, 0x1.dac829f365795p+12
+ };
+ f = h * ((c[0] + h * c[1]) + h2 * (c[2] + h * c[3])
+ + h4 * ((c[4] + h * c[5]) + h2 * (c[6] + h * c[7])));
+ }
+ else if (__glibc_unlikely (t > 0x401ceccbu && t < 0x401d95cau))
+ {
+ double h = (x + 0x1.3a7fc9600f86cp+1)
+ + 0x1.55f64f98af8dp-55;
+ double h2 = h * h;
+ double h4 = h2 * h2;
+ static const double c[] =
+ {
+ 0x1.83fe966af535fp+0, 0x1.36eebb002f61ap+2,
+ 0x1.694a60589a0b3p+0, 0x1.1718d7aedb0b5p+3,
+ 0x1.733a045eca0d3p+2, 0x1.8d4297421205bp+4,
+ 0x1.7feea5fb29965p+4
+ };
+ f = h
+ * ((c[0] + h * c[1]) + h2 * (c[2] + h * c[3])
+ + h4 * ((c[4] + h * c[5]) + h2 * (c[6])));
+ }
+ else if (__glibc_unlikely (t > 0x40492009u && t < 0x404940efu))
+ {
+ double h = (x + 0x1.9260dbc9e59afp+1)
+ + 0x1.f717cd335a7b3p-53;
+ double h2 = h * h;
+ double h4 = h2 * h2;
+ static const double c[] =
+ {
+ 0x1.f20a65f2fac55p+2, 0x1.9d4d297715105p+4,
+ 0x1.c1137124d5b21p+6, 0x1.267203d24de38p+9,
+ 0x1.99a63399a0b44p+11, 0x1.2941214faaf0cp+14,
+ 0x1.bb912c0c9cdd1p+16
+ };
+ f = h * ((c[0] + h * c[1]) + h2 * (c[2] + h * c[3])
+ + h4 * ((c[4] + h * c[5]) + h2 * (c[6])));
+ }
+ else
+ {
+ f = 0x1.250d048e7a1bdp+0 - f;
+ double lp = as_ln (as_sinpi (x - fx) * z);
+ f -= lp;
+ }
}
- /* 8.0 <= x < 2**26 */
- } else if (ix < 0x4c800000) {
- t = __ieee754_logf(x);
- z = one/x;
- y = z*z;
- w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
- r = (x-half)*(t-one)+w;
- } else
- /* 2**26 <= x <= inf */
- r = math_narrow_eval (x*(__ieee754_logf(x)-one));
- /* NADJ is set for negative arguments but not otherwise,
- resulting in warnings that it may be used uninitialized
- although in the cases where it is used it has always been
- set. */
- DIAG_PUSH_NEEDS_COMMENT;
- DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized");
- if(hx<0) r = nadj - r;
- DIAG_POP_NEEDS_COMMENT;
- return r;
+ }
+ }
+
+ uint64_t tl = (asuint64 (f) + 5) & 0xfffffff;
+ float r = f;
+ if (__glibc_unlikely (tl <= 31u))
+ {
+ t = asuint (x);
+ for (unsigned i = 0; i < array_length (tb); i++)
+ {
+ if (t == asuint (tb[i].x))
+ return tb[i].f + tb[i].df;
+ }
+ }
+ return r;
}
libm_alias_finite (__ieee754_lgammaf_r, __lgammaf_r)
diff --git a/sysdeps/ieee754/flt-32/k_tanf.c b/sysdeps/ieee754/flt-32/k_tanf.c
index e1c9d14..1cc8931 100644
--- a/sysdeps/ieee754/flt-32/k_tanf.c
+++ b/sysdeps/ieee754/flt-32/k_tanf.c
@@ -1,101 +1 @@
-/* k_tanf.c -- float version of k_tan.c
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: k_tanf.c,v 1.4 1995/05/10 20:46:39 jtc Exp $";
-#endif
-
-#include <float.h>
-#include <math.h>
-#include <math_private.h>
-#include <math-underflow.h>
-static const float
-one = 1.0000000000e+00, /* 0x3f800000 */
-pio4 = 7.8539812565e-01, /* 0x3f490fda */
-pio4lo= 3.7748947079e-08, /* 0x33222168 */
-T[] = {
- 3.3333334327e-01, /* 0x3eaaaaab */
- 1.3333334029e-01, /* 0x3e088889 */
- 5.3968254477e-02, /* 0x3d5d0dd1 */
- 2.1869488060e-02, /* 0x3cb327a4 */
- 8.8632395491e-03, /* 0x3c11371f */
- 3.5920790397e-03, /* 0x3b6b6916 */
- 1.4562094584e-03, /* 0x3abede48 */
- 5.8804126456e-04, /* 0x3a1a26c8 */
- 2.4646313977e-04, /* 0x398137b9 */
- 7.8179444245e-05, /* 0x38a3f445 */
- 7.1407252108e-05, /* 0x3895c07a */
- -1.8558637748e-05, /* 0xb79bae5f */
- 2.5907305826e-05, /* 0x37d95384 */
-};
-
-float __kernel_tanf(float x, float y, int iy)
-{
- float z,r,v,w,s;
- int32_t ix,hx;
- GET_FLOAT_WORD(hx,x);
- ix = hx&0x7fffffff; /* high word of |x| */
- if(ix<0x39000000) /* x < 2**-13 */
- {if((int)x==0) { /* generate inexact */
- if((ix|(iy+1))==0) return one/fabsf(x);
- else if (iy == 1)
- {
- math_check_force_underflow (x);
- return x;
- }
- else
- return -one / x;
- }
- }
- if(ix>=0x3f2ca140) { /* |x|>=0.6744 */
- if(hx<0) {x = -x; y = -y;}
- z = pio4-x;
- w = pio4lo-y;
- x = z+w; y = 0.0;
- if (fabsf (x) < 0x1p-13f)
- return (1 - ((hx >> 30) & 2)) * iy * (1.0f - 2 * iy * x);
- }
- z = x*x;
- w = z*z;
- /* Break x^5*(T[1]+x^2*T[2]+...) into
- * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
- * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
- */
- r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11]))));
- v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12])))));
- s = z*x;
- r = y + z*(s*(r+v)+y);
- r += T[0]*s;
- w = x+r;
- if(ix>=0x3f2ca140) {
- v = (float)iy;
- return (float)(1-((hx>>30)&2))*(v-(float)2.0*(x-(w*w/(w+v)-r)));
- }
- if(iy==1) return w;
- else { /* if allow error up to 2 ulp,
- simply return -1.0/(x+r) here */
- /* compute -1.0/(x+r) accurately */
- float a,t;
- int32_t i;
- z = w;
- GET_FLOAT_WORD(i,z);
- SET_FLOAT_WORD(z,i&0xfffff000);
- v = r-(z - x); /* z+v = r+x */
- t = a = -(float)1.0/w; /* a = -1.0/w */
- GET_FLOAT_WORD(i,t);
- SET_FLOAT_WORD(t,i&0xfffff000);
- s = (float)1.0+t*z;
- return t+a*(s+t*v);
- }
-}
+/* Not needed. */
diff --git a/sysdeps/ieee754/flt-32/lgamma_negf.c b/sysdeps/ieee754/flt-32/lgamma_negf.c
index a8aa74e..1cc8931 100644
--- a/sysdeps/ieee754/flt-32/lgamma_negf.c
+++ b/sysdeps/ieee754/flt-32/lgamma_negf.c
@@ -1,282 +1 @@
-/* lgammaf expanding around zeros.
- Copyright (C) 2015-2024 Free Software Foundation, Inc.
- This file is part of the GNU C Library.
-
- The GNU C Library is free software; you can redistribute it and/or
- modify it under the terms of the GNU Lesser General Public
- License as published by the Free Software Foundation; either
- version 2.1 of the License, or (at your option) any later version.
-
- The GNU C Library 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
- Lesser General Public License for more details.
-
- You should have received a copy of the GNU Lesser General Public
- License along with the GNU C Library; if not, see
- <https://www.gnu.org/licenses/>. */
-
-#include <float.h>
-#include <math.h>
-#include <math-narrow-eval.h>
-#include <math_private.h>
-#include <fenv_private.h>
-
-static const float lgamma_zeros[][2] =
- {
- { -0x2.74ff94p+0f, 0x1.3fe0f2p-24f },
- { -0x2.bf682p+0f, -0x1.437b2p-24f },
- { -0x3.24c1b8p+0f, 0x6.c34cap-28f },
- { -0x3.f48e2cp+0f, 0x1.707a04p-24f },
- { -0x4.0a13ap+0f, 0x1.e99aap-24f },
- { -0x4.fdd5ep+0f, 0x1.64454p-24f },
- { -0x5.021a98p+0f, 0x2.03d248p-24f },
- { -0x5.ffa4cp+0f, 0x2.9b82fcp-24f },
- { -0x6.005ac8p+0f, -0x1.625f24p-24f },
- { -0x6.fff3p+0f, 0x2.251e44p-24f },
- { -0x7.000dp+0f, 0x8.48078p-28f },
- { -0x7.fffe6p+0f, 0x1.fa98c4p-28f },
- { -0x8.0001ap+0f, -0x1.459fcap-28f },
- { -0x8.ffffdp+0f, -0x1.c425e8p-24f },
- { -0x9.00003p+0f, 0x1.c44b82p-24f },
- { -0xap+0f, 0x4.9f942p-24f },
- { -0xap+0f, -0x4.9f93b8p-24f },
- { -0xbp+0f, 0x6.b9916p-28f },
- { -0xbp+0f, -0x6.b9915p-28f },
- { -0xcp+0f, 0x8.f76c8p-32f },
- { -0xcp+0f, -0x8.f76c7p-32f },
- { -0xdp+0f, 0xb.09231p-36f },
- { -0xdp+0f, -0xb.09231p-36f },
- { -0xep+0f, 0xc.9cba5p-40f },
- { -0xep+0f, -0xc.9cba5p-40f },
- { -0xfp+0f, 0xd.73f9fp-44f },
- };
-
-static const float e_hi = 0x2.b7e15p+0f, e_lo = 0x1.628aeep-24f;
-
-/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
- approximation to lgamma function. */
-
-static const float lgamma_coeff[] =
- {
- 0x1.555556p-4f,
- -0xb.60b61p-12f,
- 0x3.403404p-12f,
- };
-
-#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
-
-/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
- the integer end-point of the half-integer interval containing x and
- x0 is the zero of lgamma in that half-integer interval. Each
- polynomial is expressed in terms of x-xm, where xm is the midpoint
- of the interval for which the polynomial applies. */
-
-static const float poly_coeff[] =
- {
- /* Interval [-2.125, -2] (polynomial degree 5). */
- -0x1.0b71c6p+0f,
- -0xc.73a1ep-4f,
- -0x1.ec8462p-4f,
- -0xe.37b93p-4f,
- -0x1.02ed36p-4f,
- -0xe.cbe26p-4f,
- /* Interval [-2.25, -2.125] (polynomial degree 5). */
- -0xf.29309p-4f,
- -0xc.a5cfep-4f,
- 0x3.9c93fcp-4f,
- -0x1.02a2fp+0f,
- 0x9.896bep-4f,
- -0x1.519704p+0f,
- /* Interval [-2.375, -2.25] (polynomial degree 5). */
- -0xd.7d28dp-4f,
- -0xe.6964cp-4f,
- 0xb.0d4f1p-4f,
- -0x1.9240aep+0f,
- 0x1.dadabap+0f,
- -0x3.1778c4p+0f,
- /* Interval [-2.5, -2.375] (polynomial degree 6). */
- -0xb.74ea2p-4f,
- -0x1.2a82cp+0f,
- 0x1.880234p+0f,
- -0x3.320c4p+0f,
- 0x5.572a38p+0f,
- -0x9.f92bap+0f,
- 0x1.1c347ep+4f,
- /* Interval [-2.625, -2.5] (polynomial degree 6). */
- -0x3.d10108p-4f,
- 0x1.cd5584p+0f,
- 0x3.819c24p+0f,
- 0x6.84cbb8p+0f,
- 0xb.bf269p+0f,
- 0x1.57fb12p+4f,
- 0x2.7b9854p+4f,
- /* Interval [-2.75, -2.625] (polynomial degree 6). */
- -0x6.b5d25p-4f,
- 0x1.28d604p+0f,
- 0x1.db6526p+0f,
- 0x2.e20b38p+0f,
- 0x4.44c378p+0f,
- 0x6.62a08p+0f,
- 0x9.6db3ap+0f,
- /* Interval [-2.875, -2.75] (polynomial degree 5). */
- -0x8.a41b2p-4f,
- 0xc.da87fp-4f,
- 0x1.147312p+0f,
- 0x1.7617dap+0f,
- 0x1.d6c13p+0f,
- 0x2.57a358p+0f,
- /* Interval [-3, -2.875] (polynomial degree 5). */
- -0xa.046d6p-4f,
- 0x9.70b89p-4f,
- 0xa.a89a6p-4f,
- 0xd.2f2d8p-4f,
- 0xd.e32b4p-4f,
- 0xf.fb741p-4f,
- };
-
-static const size_t poly_deg[] =
- {
- 5,
- 5,
- 5,
- 6,
- 6,
- 6,
- 5,
- 5,
- };
-
-static const size_t poly_end[] =
- {
- 5,
- 11,
- 17,
- 24,
- 31,
- 38,
- 44,
- 50,
- };
-
-/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
-
-static float
-lg_sinpi (float x)
-{
- if (x <= 0.25f)
- return __sinf (M_PIf * x);
- else
- return __cosf (M_PIf * (0.5f - x));
-}
-
-/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
-
-static float
-lg_cospi (float x)
-{
- if (x <= 0.25f)
- return __cosf (M_PIf * x);
- else
- return __sinf (M_PIf * (0.5f - x));
-}
-
-/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
-
-static float
-lg_cotpi (float x)
-{
- return lg_cospi (x) / lg_sinpi (x);
-}
-
-/* Compute lgamma of a negative argument -15 < X < -2, setting
- *SIGNGAMP accordingly. */
-
-float
-__lgamma_negf (float x, int *signgamp)
-{
- /* Determine the half-integer region X lies in, handle exact
- integers and determine the sign of the result. */
- int i = floorf (-2 * x);
- if ((i & 1) == 0 && i == -2 * x)
- return 1.0f / 0.0f;
- float xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
- i -= 4;
- *signgamp = ((i & 2) == 0 ? -1 : 1);
-
- SET_RESTORE_ROUNDF (FE_TONEAREST);
-
- /* Expand around the zero X0 = X0_HI + X0_LO. */
- float x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1];
- float xdiff = x - x0_hi - x0_lo;
-
- /* For arguments in the range -3 to -2, use polynomial
- approximations to an adjusted version of the gamma function. */
- if (i < 2)
- {
- int j = floorf (-8 * x) - 16;
- float xm = (-33 - 2 * j) * 0.0625f;
- float x_adj = x - xm;
- size_t deg = poly_deg[j];
- size_t end = poly_end[j];
- float g = poly_coeff[end];
- for (size_t j = 1; j <= deg; j++)
- g = g * x_adj + poly_coeff[end - j];
- return __log1pf (g * xdiff / (x - xn));
- }
-
- /* The result we want is log (sinpi (X0) / sinpi (X))
- + log (gamma (1 - X0) / gamma (1 - X)). */
- float x_idiff = fabsf (xn - x), x0_idiff = fabsf (xn - x0_hi - x0_lo);
- float log_sinpi_ratio;
- if (x0_idiff < x_idiff * 0.5f)
- /* Use log not log1p to avoid inaccuracy from log1p of arguments
- close to -1. */
- log_sinpi_ratio = __ieee754_logf (lg_sinpi (x0_idiff)
- / lg_sinpi (x_idiff));
- else
- {
- /* Use log1p not log to avoid inaccuracy from log of arguments
- close to 1. X0DIFF2 has positive sign if X0 is further from
- XN than X is from XN, negative sign otherwise. */
- float x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5f;
- float sx0d2 = lg_sinpi (x0diff2);
- float cx0d2 = lg_cospi (x0diff2);
- log_sinpi_ratio = __log1pf (2 * sx0d2
- * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff)));
- }
-
- float log_gamma_ratio;
- float y0 = math_narrow_eval (1 - x0_hi);
- float y0_eps = -x0_hi + (1 - y0) - x0_lo;
- float y = math_narrow_eval (1 - x);
- float y_eps = -x + (1 - y);
- /* We now wish to compute LOG_GAMMA_RATIO
- = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
- accurately approximates the difference Y0 + Y0_EPS - Y -
- Y_EPS. Use Stirling's approximation. */
- float log_gamma_high
- = (xdiff * __log1pf ((y0 - e_hi - e_lo + y0_eps) / e_hi)
- + (y - 0.5f + y_eps) * __log1pf (xdiff / y));
- /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
- float y0r = 1 / y0, yr = 1 / y;
- float y0r2 = y0r * y0r, yr2 = yr * yr;
- float rdiff = -xdiff / (y * y0);
- float bterm[NCOEFF];
- float dlast = rdiff, elast = rdiff * yr * (yr + y0r);
- bterm[0] = dlast * lgamma_coeff[0];
- for (size_t j = 1; j < NCOEFF; j++)
- {
- float dnext = dlast * y0r2 + elast;
- float enext = elast * yr2;
- bterm[j] = dnext * lgamma_coeff[j];
- dlast = dnext;
- elast = enext;
- }
- float log_gamma_low = 0;
- for (size_t j = 0; j < NCOEFF; j++)
- log_gamma_low += bterm[NCOEFF - 1 - j];
- log_gamma_ratio = log_gamma_high + log_gamma_low;
-
- return log_sinpi_ratio + log_gamma_ratio;
-}
+/* Not needed. */
diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
index dc07ebd..b30a03e 100644
--- a/sysdeps/ieee754/flt-32/math_config.h
+++ b/sysdeps/ieee754/flt-32/math_config.h
@@ -57,6 +57,33 @@ static inline int32_t
converttoint (double_t x);
#endif
+#ifndef ROUNDEVEN_INTRINSICS
+/* When set, roundeven_finite will route to the internal roundeven function. */
+# define ROUNDEVEN_INTRINSICS 1
+#endif
+
+/* Round x to nearest integer value in floating-point format, rounding halfway
+ cases to even. If the input is non finite the result is unspecified. */
+static inline double
+roundeven_finite (double x)
+{
+ if (!isfinite (x))
+ __builtin_unreachable ();
+#if ROUNDEVEN_INTRINSICS
+ return roundeven (x);
+#else
+ double y = round (x);
+ if (fabs (x - y) == 0.5)
+ {
+ union { double f; uint64_t i; } u = {y};
+ union { double f; uint64_t i; } v = {y - copysign (1.0, x)};
+ if (__builtin_ctzll (v.i) > __builtin_ctzll (u.i))
+ y = v.f;
+ }
+ return y;
+#endif
+}
+
static inline uint32_t
asuint (float f)
{
diff --git a/sysdeps/ieee754/flt-32/s_cbrtf.c b/sysdeps/ieee754/flt-32/s_cbrtf.c
index 68b8b0e..5a7a9a9 100644
--- a/sysdeps/ieee754/flt-32/s_cbrtf.c
+++ b/sysdeps/ieee754/flt-32/s_cbrtf.c
@@ -1,61 +1,99 @@
-/* Compute cubic root of float value.
- Copyright (C) 1997-2024 Free Software Foundation, Inc.
- This file is part of the GNU C Library.
+/* Correctly-rounded cubic root of binary32 value.
- The GNU C Library is free software; you can redistribute it and/or
- modify it under the terms of the GNU Lesser General Public
- License as published by the Free Software Foundation; either
- version 2.1 of the License, or (at your option) any later version.
+Copyright (c) 2023, 2024 Alexei Sibidanov.
- The GNU C Library 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
- Lesser General Public License for more details.
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/cbrt/cbrtf.c, revision bc385c2).
- You should have received a copy of the GNU Lesser General Public
- License along with the GNU C Library; if not, see
- <https://www.gnu.org/licenses/>. */
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
-#include <math.h>
-#include <libm-alias-float.h>
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
-#define CBRT2 1.2599210498948731648 /* 2^(1/3) */
-#define SQR_CBRT2 1.5874010519681994748 /* 2^(2/3) */
-
-static const double factor[5] =
-{
- 1.0 / SQR_CBRT2,
- 1.0 / CBRT2,
- 1.0,
- CBRT2,
- SQR_CBRT2
-};
-
+#include <fenv.h>
+#include <libm-alias-float.h>
+#include <math.h>
+#include <stdint.h>
+#include "math_config.h"
float
__cbrtf (float x)
{
- float xm, ym, u, t2;
- int xe;
-
- /* Reduce X. XM now is an range 1.0 to 0.5. */
- xm = __frexpf (fabsf (x), &xe);
-
- /* If X is not finite or is null return it (with raising exceptions
- if necessary.
- Note: *Our* version of `frexp' sets XE to zero if the argument is
- Inf or NaN. This is not portable but faster. */
- if (xe == 0 && fpclassify (x) <= FP_ZERO)
- return x + x;
-
- u = (0.492659620528969547 + (0.697570460207922770
- - 0.191502161678719066 * xm) * xm);
-
- t2 = u * u * u;
-
- ym = u * (t2 + 2.0 * xm) / (2.0 * t2 + xm) * factor[2 + xe % 3];
-
- return __ldexpf (x > 0.0 ? ym : -ym, xe / 3);
+ static const union
+ {
+ double d;
+ uint64_t u;
+ } escale[3] =
+ {
+ { .d = 1.0 },
+ { .d = 0x1.428a2f98d728bp+0 }, /* 2^(1/3) */
+ { .d = 0x1.965fea53d6e3dp+0 }, /* 2^(2/3) */
+ };
+ uint32_t u = asuint (x);
+ uint32_t au = u << 1;
+ uint32_t sgn = u >> 31;
+ uint32_t e = au >> 24;
+ if (__glibc_unlikely (au < 1u << 24 || au >= 0xffu << 24))
+ {
+ if (au >= 0xffu << 24)
+ return x + x; /* inf, nan */
+ if (au == 0)
+ return x; /* +-0 */
+ int nz = __builtin_clz (au) - 7; /* subnormal */
+ au <<= nz;
+ e -= nz - 1;
+ }
+ uint32_t mant = au & 0xffffff;
+ e += 899;
+ uint32_t et = e / 3, it = e % 3;
+ uint64_t isc = escale[it].u;
+ isc += (int64_t) (et - 342) << 52;
+ isc |= (int64_t) sgn << 63;
+ double cvt2 = asdouble (isc);
+ static const double c[] =
+ {
+ 0x1.2319d352ea5d5p-1, 0x1.67ad8ee258d1ap-1, -0x1.9342edf9cbad9p-2,
+ 0x1.b6388fc510a75p-3, -0x1.6002455599e2fp-4, 0x1.7b096936192c4p-6,
+ -0x1.e5577187e8bf8p-9, 0x1.169ef81d6c34ep-12
+ };
+ double z = asdouble ((uint64_t) mant << 28 | UINT64_C(0x3ff) << 52);
+ double r0 = -0x1.9931c6c2d19d1p-6 / z;
+ double z2 = z * z;
+ double z4 = z2 * z2;
+ double f = ((c[0] + z * c[1]) + z2 * (c[2] + z * c[3]))
+ + z4 * ((c[4] + z * c[5]) + z2 * (c[6] + z * c[7])) + r0;
+ double r = f * cvt2;
+ float ub = r;
+ float lb = r - cvt2 * 1.4182e-9;
+ if (__glibc_likely (ub == lb))
+ return ub;
+ const double u0 = -0x1.ab16ec65d138fp+3;
+ double h = f * f * f - z;
+ f -= (f * r0 * u0) * h;
+ r = f * cvt2;
+ uint64_t cvt1 = asuint64 (r);
+ ub = r;
+ int64_t m0 = cvt1 << 19;
+ int64_t m1 = m0 >> 63;
+ if (__glibc_unlikely ((m0 ^ m1) < (UINT64_C(1) << 31)))
+ {
+ cvt1 = (cvt1 + (UINT64_C(1) << 31)) & UINT64_C(0xffffffff00000000);
+ ub = asdouble (cvt1);
+ }
+ return ub;
}
libm_alias_float (__cbrt, cbrt)
diff --git a/sysdeps/ieee754/flt-32/s_erfcf.c b/sysdeps/ieee754/flt-32/s_erfcf.c
new file mode 100644
index 0000000..3dae2a0
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/s_erfcf.c
@@ -0,0 +1,187 @@
+/* Correctly-rounded complementary error function for the binary32 format
+
+Copyright (c) 2023, 2024 Alexei Sibidanov.
+
+This file is part of the CORE-MATH project
+project (file src/binary32/erfc/erfcf.c revision bc385c2).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+#include <errno.h>
+#include <math.h>
+#include <stdint.h>
+#include <libm-alias-float.h>
+#include "math_config.h"
+
+static const double E[] =
+ {
+ 0x1p+0, 0x1.0163da9fb3335p+0, 0x1.02c9a3e778061p+0,
+ 0x1.04315e86e7f85p+0, 0x1.059b0d3158574p+0, 0x1.0706b29ddf6dep+0,
+ 0x1.0874518759bc8p+0, 0x1.09e3ecac6f383p+0, 0x1.0b5586cf9890fp+0,
+ 0x1.0cc922b7247f7p+0, 0x1.0e3ec32d3d1a2p+0, 0x1.0fb66affed31bp+0,
+ 0x1.11301d0125b51p+0, 0x1.12abdc06c31ccp+0, 0x1.1429aaea92dep+0,
+ 0x1.15a98c8a58e51p+0, 0x1.172b83c7d517bp+0, 0x1.18af9388c8deap+0,
+ 0x1.1a35beb6fcb75p+0, 0x1.1bbe084045cd4p+0, 0x1.1d4873168b9aap+0,
+ 0x1.1ed5022fcd91dp+0, 0x1.2063b88628cd6p+0, 0x1.21f49917ddc96p+0,
+ 0x1.2387a6e756238p+0, 0x1.251ce4fb2a63fp+0, 0x1.26b4565e27cddp+0,
+ 0x1.284dfe1f56381p+0, 0x1.29e9df51fdee1p+0, 0x1.2b87fd0dad99p+0,
+ 0x1.2d285a6e4030bp+0, 0x1.2ecafa93e2f56p+0, 0x1.306fe0a31b715p+0,
+ 0x1.32170fc4cd831p+0, 0x1.33c08b26416ffp+0, 0x1.356c55f929ff1p+0,
+ 0x1.371a7373aa9cbp+0, 0x1.38cae6d05d866p+0, 0x1.3a7db34e59ff7p+0,
+ 0x1.3c32dc313a8e5p+0, 0x1.3dea64c123422p+0, 0x1.3fa4504ac801cp+0,
+ 0x1.4160a21f72e2ap+0, 0x1.431f5d950a897p+0, 0x1.44e086061892dp+0,
+ 0x1.46a41ed1d0057p+0, 0x1.486a2b5c13cdp+0, 0x1.4a32af0d7d3dep+0,
+ 0x1.4bfdad5362a27p+0, 0x1.4dcb299fddd0dp+0, 0x1.4f9b2769d2ca7p+0,
+ 0x1.516daa2cf6642p+0, 0x1.5342b569d4f82p+0, 0x1.551a4ca5d920fp+0,
+ 0x1.56f4736b527dap+0, 0x1.58d12d497c7fdp+0, 0x1.5ab07dd485429p+0,
+ 0x1.5c9268a5946b7p+0, 0x1.5e76f15ad2148p+0, 0x1.605e1b976dc09p+0,
+ 0x1.6247eb03a5585p+0, 0x1.6434634ccc32p+0, 0x1.6623882552225p+0,
+ 0x1.68155d44ca973p+0, 0x1.6a09e667f3bcdp+0, 0x1.6c012750bdabfp+0,
+ 0x1.6dfb23c651a2fp+0, 0x1.6ff7df9519484p+0, 0x1.71f75e8ec5f74p+0,
+ 0x1.73f9a48a58174p+0, 0x1.75feb564267c9p+0, 0x1.780694fde5d3fp+0,
+ 0x1.7a11473eb0187p+0, 0x1.7c1ed0130c132p+0, 0x1.7e2f336cf4e62p+0,
+ 0x1.80427543e1a12p+0, 0x1.82589994cce13p+0, 0x1.8471a4623c7adp+0,
+ 0x1.868d99b4492edp+0, 0x1.88ac7d98a6699p+0, 0x1.8ace5422aa0dbp+0,
+ 0x1.8cf3216b5448cp+0, 0x1.8f1ae99157736p+0, 0x1.9145b0b91ffc6p+0,
+ 0x1.93737b0cdc5e5p+0, 0x1.95a44cbc8520fp+0, 0x1.97d829fde4e5p+0,
+ 0x1.9a0f170ca07bap+0, 0x1.9c49182a3f09p+0, 0x1.9e86319e32323p+0,
+ 0x1.a0c667b5de565p+0, 0x1.a309bec4a2d33p+0, 0x1.a5503b23e255dp+0,
+ 0x1.a799e1330b358p+0, 0x1.a9e6b5579fdbfp+0, 0x1.ac36bbfd3f37ap+0,
+ 0x1.ae89f995ad3adp+0, 0x1.b0e07298db666p+0, 0x1.b33a2b84f15fbp+0,
+ 0x1.b59728de5593ap+0, 0x1.b7f76f2fb5e47p+0, 0x1.ba5b030a1064ap+0,
+ 0x1.bcc1e904bc1d2p+0, 0x1.bf2c25bd71e09p+0, 0x1.c199bdd85529cp+0,
+ 0x1.c40ab5fffd07ap+0, 0x1.c67f12e57d14bp+0, 0x1.c8f6d9406e7b5p+0,
+ 0x1.cb720dcef9069p+0, 0x1.cdf0b555dc3fap+0, 0x1.d072d4a07897cp+0,
+ 0x1.d2f87080d89f2p+0, 0x1.d5818dcfba487p+0, 0x1.d80e316c98398p+0,
+ 0x1.da9e603db3285p+0, 0x1.dd321f301b46p+0, 0x1.dfc97337b9b5fp+0,
+ 0x1.e264614f5a129p+0, 0x1.e502ee78b3ff6p+0, 0x1.e7a51fbc74c83p+0,
+ 0x1.ea4afa2a490dap+0, 0x1.ecf482d8e67f1p+0, 0x1.efa1bee615a27p+0,
+ 0x1.f252b376bba97p+0, 0x1.f50765b6e454p+0, 0x1.f7bfdad9cbe14p+0,
+ 0x1.fa7c1819e90d8p+0, 0x1.fd3c22b8f71f1p+0
+ };
+
+float
+__erfcf (float xf)
+{
+ float axf = fabsf (xf);
+ double axd = axf;
+ double x2 = axd * axd;
+ uint32_t t = asuint (xf);
+ unsigned int at = t & (~0u >> 1);
+ unsigned int sgn = t >> 31;
+ int64_t i = at > 0x40051000;
+ /* for x < -0x1.ea8f94p+1, erfc(x) rounds to 2 (to nearest) */
+ if (__glibc_unlikely (t > 0xc07547ca))
+ { /* xf < -0x1.ea8f94p+1 */
+ if (__glibc_unlikely (t >= 0xff800000))
+ { /* -Inf or NaN */
+ if (t == 0xff800000)
+ return 2.0f; /* -Inf */
+ return xf + xf; /* NaN */
+ }
+ return 2.0f - 0x1p-25f; /* rounds to 2 or nextbelow(2) */
+ }
+ /* at is the absolute value of xf
+ for x >= 0x1.41bbf8p+3, erfc(x) < 2^-150, thus rounds to 0 or to 2^-149
+ depending on the rounding mode */
+ if (__glibc_unlikely (at >= 0x4120ddfc))
+ { /* |xf| >= 0x1.41bbf8p+3 */
+ if (__glibc_unlikely (at >= 0x7f800000))
+ { /* +Inf or NaN */
+ if (at == 0x7f800000)
+ return 0.0f; /* +Inf */
+ return xf + xf; /* NaN */
+ }
+ __set_errno (ERANGE);
+ /* 0x1p-149f * 0.25f rounds to 0 or 2^-149 depending on rounding */
+ return 0x1p-149f * 0.25f;
+ }
+ if (__glibc_unlikely (at <= 0x3db80000))
+ { /* |x| <= 0x1.7p-4 */
+ if (__glibc_unlikely (t == 0xb76c9f62))
+ return 0x1.00010ap+0f + 0x1p-25f; /* exceptional case */
+ /* for |x| <= 0x1.c5bf88p-26. erfc(x) rounds to 1 (to nearest) */
+ if (__glibc_unlikely (at <= 0x32e2dfc4))
+ { /* |x| <= 0x1.c5bf88p-26 */
+ if (__glibc_unlikely (at == 0))
+ return 1.0f;
+ static const float d[] = { -0x1p-26, 0x1p-25 };
+ return 1.0f + d[sgn];
+ }
+ /* around 0, erfc(x) behaves as 1 - (odd polynomial) */
+ static const double c[] =
+ {
+ 0x1.20dd750429b6dp+0, -0x1.812746b03610bp-2, 0x1.ce2f218831d2fp-4,
+ -0x1.b82c609607dcbp-6, 0x1.553af09b8008ep-8
+ };
+ double f0 = xf
+ * (c[0] + x2 * (c[1] + x2 * (c[2] + x2 * (c[3] + x2 * (c[4])))));
+ return 1.0 - f0;
+ }
+
+ /* now -0x1.ea8f94p+1 <= x <= 0x1.41bbf8p+3, with |x| > 0x1.7p-4 */
+ const double iln2 = 0x1.71547652b82fep+0;
+ const double ln2h = 0x1.62e42fefap-8;
+ const double ln2l = 0x1.cf79abd6f5dc8p-47;
+ uint64_t jt = asuint64 (fma (x2, iln2, -(1024 + 0x1p-8)));
+ int64_t j = (int64_t) (jt << 12) >> 48;
+ double S = asdouble (((j >> 7) + (0x3ff | sgn << 11)) << 52);
+ static const double ch[] =
+ {
+ -0x1.ffffffffff333p-2, 0x1.5555555556a14p-3, -0x1.55556666659b4p-5,
+ 0x1.1111074cc7b22p-7
+ };
+ double d = (x2 + ln2h * j) + ln2l * j;
+ double d2 = d * d;
+ double e0 = E[j & 127];
+ double f = d + d2 * ((ch[0] + d * ch[1]) + d2 * (ch[2] + d * ch[3]));
+ static const double ct[][16] =
+ {
+ {
+ 0x1.c162355429b28p-1, 0x1.d99999999999ap+1, 0x1.da951cece2b85p-2,
+ -0x1.70ef6cff4bcc4p+0, 0x1.3d7f7b3d617dep+1, -0x1.9d0aa47537c51p+1,
+ 0x1.9754ea9a3fcb1p+1, -0x1.27a5453fcc015p+1, 0x1.1ef2e0531aebap+0,
+ -0x1.eca090f5a1c06p-3, -0x1.7a3cd173a063cp-4, 0x1.30fa68a68fdddp-4,
+ 0x1.55ad9a326993ap-10, -0x1.07e7b0bb39fbfp-6, 0x1.2328706c0e95p-10,
+ 0x1.d6aa0b7b19cfep-9
+ },
+ {
+ 0x1.137c8983f8516p+2, 0x1.799999999999ap+1, 0x1.05b53aa241333p-3,
+ -0x1.a3f53872bf87p-3, 0x1.de4c30742c9d5p-4, -0x1.cb24bfa591986p-5,
+ 0x1.666aec059ca5fp-6, -0x1.a61250eb26b0bp-8, 0x1.2b28b7924b34dp-10,
+ 0x1.41b13a9d45013p-15, -0x1.6dd5e8a273613p-14, 0x1.09ce8ea5e8da5p-16,
+ 0x1.33923b4102981p-18, -0x1.1dfd161e3f984p-19, -0x1.c87618fcae3b3p-23,
+ 0x1.e8a6ffa0ba2c7p-23
+ }
+ };
+ double z = (axd - ct[i][0]) / (axd + ct[i][1]);
+ double z2 = z * z, z4 = z2 * z2;
+ double z8 = z4 * z4;
+ const double *c = ct[i] + 3;
+ double s = (((c[0] + z * c[1]) + z2 * (c[2] + z * c[3]))
+ + z4 * ((c[4] + z * c[5]) + z2 * (c[6] + z * c[7])))
+ + z8 * (((c[8] + z * c[9]) + z2 * (c[10] + z * c[11])) + z4 * (c[12]));
+ s = ct[i][2] + z * s;
+ static const double off[] = { 0, 2 };
+ double r = (S * (e0 - f * e0)) * s;
+ double y = off[sgn] + r;
+ return y;
+}
+libm_alias_float (__erfc, erfc)
diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c
index ba29734..025c207 100644
--- a/sysdeps/ieee754/flt-32/s_erff.c
+++ b/sysdeps/ieee754/flt-32/s_erff.c
@@ -1,232 +1,256 @@
-/* s_erff.c -- float version of s_erf.c.
- */
+/* Correctly-rounded error function for binary32 value.
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+Copyright (c) 2022-2024 Alexei Sibidanov.
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_erff.c,v 1.4 1995/05/10 20:47:07 jtc Exp $";
-#endif
+This file is part of the CORE-MATH project
+project (file src/binary32/erf/erff.c revision bc385c2).
-#include <errno.h>
-#include <float.h>
-#include <math.h>
-#include <math-narrow-eval.h>
-#include <math_private.h>
-#include <math-underflow.h>
-#include <libm-alias-float.h>
-#include <fix-int-fp-convert-zero.h>
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
-static const float
-tiny = 1e-30,
-half= 5.0000000000e-01, /* 0x3F000000 */
-one = 1.0000000000e+00, /* 0x3F800000 */
-two = 2.0000000000e+00, /* 0x40000000 */
- /* c = (subfloat)0.84506291151 */
-erx = 8.4506291151e-01, /* 0x3f58560b */
-/*
- * Coefficients for approximation to erf on [0,0.84375]
- */
-efx = 1.2837916613e-01, /* 0x3e0375d4 */
-pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
-pp1 = -3.2504209876e-01, /* 0xbea66beb */
-pp2 = -2.8481749818e-02, /* 0xbce9528f */
-pp3 = -5.7702702470e-03, /* 0xbbbd1489 */
-pp4 = -2.3763017452e-05, /* 0xb7c756b1 */
-qq1 = 3.9791721106e-01, /* 0x3ecbbbce */
-qq2 = 6.5022252500e-02, /* 0x3d852a63 */
-qq3 = 5.0813062117e-03, /* 0x3ba68116 */
-qq4 = 1.3249473704e-04, /* 0x390aee49 */
-qq5 = -3.9602282413e-06, /* 0xb684e21a */
-/*
- * Coefficients for approximation to erf in [0.84375,1.25]
- */
-pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */
-pa1 = 4.1485610604e-01, /* 0x3ed46805 */
-pa2 = -3.7220788002e-01, /* 0xbebe9208 */
-pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */
-pa4 = -1.1089469492e-01, /* 0xbde31cc2 */
-pa5 = 3.5478305072e-02, /* 0x3d1151b3 */
-pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */
-qa1 = 1.0642088205e-01, /* 0x3dd9f331 */
-qa2 = 5.4039794207e-01, /* 0x3f0a5785 */
-qa3 = 7.1828655899e-02, /* 0x3d931ae7 */
-qa4 = 1.2617121637e-01, /* 0x3e013307 */
-qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */
-qa6 = 1.1984500103e-02, /* 0x3c445aa3 */
-/*
- * Coefficients for approximation to erfc in [1.25,1/0.35]
- */
-ra0 = -9.8649440333e-03, /* 0xbc21a093 */
-ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */
-ra2 = -1.0558626175e+01, /* 0xc128f022 */
-ra3 = -6.2375331879e+01, /* 0xc2798057 */
-ra4 = -1.6239666748e+02, /* 0xc322658c */
-ra5 = -1.8460508728e+02, /* 0xc3389ae7 */
-ra6 = -8.1287437439e+01, /* 0xc2a2932b */
-ra7 = -9.8143291473e+00, /* 0xc11d077e */
-sa1 = 1.9651271820e+01, /* 0x419d35ce */
-sa2 = 1.3765776062e+02, /* 0x4309a863 */
-sa3 = 4.3456588745e+02, /* 0x43d9486f */
-sa4 = 6.4538726807e+02, /* 0x442158c9 */
-sa5 = 4.2900814819e+02, /* 0x43d6810b */
-sa6 = 1.0863500214e+02, /* 0x42d9451f */
-sa7 = 6.5702495575e+00, /* 0x40d23f7c */
-sa8 = -6.0424413532e-02, /* 0xbd777f97 */
-/*
- * Coefficients for approximation to erfc in [1/.35,28]
- */
-rb0 = -9.8649431020e-03, /* 0xbc21a092 */
-rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */
-rb2 = -1.7757955551e+01, /* 0xc18e104b */
-rb3 = -1.6063638306e+02, /* 0xc320a2ea */
-rb4 = -6.3756646729e+02, /* 0xc41f6441 */
-rb5 = -1.0250950928e+03, /* 0xc480230b */
-rb6 = -4.8351919556e+02, /* 0xc3f1c275 */
-sb1 = 3.0338060379e+01, /* 0x41f2b459 */
-sb2 = 3.2579251099e+02, /* 0x43a2e571 */
-sb3 = 1.5367296143e+03, /* 0x44c01759 */
-sb4 = 3.1998581543e+03, /* 0x4547fdbb */
-sb5 = 2.5530502930e+03, /* 0x451f90ce */
-sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
-sb7 = -2.2440952301e+01; /* 0xc1b38712 */
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
-float __erff(float x)
-{
- int32_t hx,ix,i;
- float R,S,P,Q,s,y,z,r;
- GET_FLOAT_WORD(hx,x);
- ix = hx&0x7fffffff;
- if(ix>=0x7f800000) { /* erf(nan)=nan */
- i = ((uint32_t)hx>>31)<<1;
- return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */
- }
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
- if(ix < 0x3f580000) { /* |x|<0.84375 */
- if(ix < 0x31800000) { /* |x|<2**-28 */
- if (ix < 0x04000000)
- {
- /* Avoid spurious underflow. */
- float ret = 0.0625f * (16.0f * x + (16.0f * efx) * x);
- math_check_force_underflow (ret);
- return ret;
- }
- return x + efx*x;
- }
- z = x*x;
- r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
- s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
- y = r/s;
- return x + x*y;
- }
- if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
- s = fabsf(x)-one;
- P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
- if(hx>=0) return erx + P/Q; else return -erx - P/Q;
- }
- if (ix >= 0x40c00000) { /* inf>|x|>=6 */
- if(hx>=0) return one-tiny; else return tiny-one;
- }
- x = fabsf(x);
- s = one/(x*x);
- if(ix< 0x4036DB6E) { /* |x| < 1/0.35 */
- R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
- ra5+s*(ra6+s*ra7))))));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
- sa5+s*(sa6+s*(sa7+s*sa8)))))));
- } else { /* |x| >= 1/0.35 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
- rb5+s*rb6)))));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
- sb5+s*(sb6+s*sb7))))));
- }
- GET_FLOAT_WORD(ix,x);
- SET_FLOAT_WORD(z,ix&0xfffff000);
- r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S);
- if(hx>=0) return one-r/x; else return r/x-one;
-}
-libm_alias_float (__erf, erf)
+#include <math.h>
+#include <stdint.h>
+#include <libm-alias-float.h>
+#include "math_config.h"
-float __erfcf(float x)
+float
+__erff (float x)
{
- int32_t hx,ix;
- float R,S,P,Q,s,y,z,r;
- GET_FLOAT_WORD(hx,x);
- ix = hx&0x7fffffff;
- if(ix>=0x7f800000) { /* erfc(nan)=nan */
- /* erfc(+-inf)=0,2 */
- float ret = (float)(((uint32_t)hx>>31)<<1)+one/x;
- if (FIX_INT_FP_CONVERT_ZERO && ret == 0.0f)
- return 0.0f;
- return ret;
- }
-
- if(ix < 0x3f580000) { /* |x|<0.84375 */
- if(ix < 0x32800000) /* |x|<2**-26 */
- return one-x;
- z = x*x;
- r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
- s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
- y = r/s;
- if(hx < 0x3e800000) { /* x<1/4 */
- return one-(x+x*y);
- } else {
- r = x*y;
- r += (x-half);
- return half - r ;
- }
- }
- if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
- s = fabsf(x)-one;
- P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
- if(hx>=0) {
- z = one-erx; return z - P/Q;
- } else {
- z = erx+P/Q; return one+z;
- }
- }
- if (ix < 0x41e00000) { /* |x|<28 */
- x = fabsf(x);
- s = one/(x*x);
- if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
- R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
- ra5+s*(ra6+s*ra7))))));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
- sa5+s*(sa6+s*(sa7+s*sa8)))))));
- } else { /* |x| >= 1/.35 ~ 2.857143 */
- if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
- rb5+s*rb6)))));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
- sb5+s*(sb6+s*sb7))))));
- }
- GET_FLOAT_WORD(ix,x);
- SET_FLOAT_WORD(z,ix&0xffffe000);
- r = __ieee754_expf(-z*z-(float)0.5625)*
- __ieee754_expf((z-x)*(z+x)+R/S);
- if(hx>0) {
- float ret = math_narrow_eval (r/x);
- if (ret == 0)
- __set_errno (ERANGE);
- return ret;
- } else
- return two-r/x;
- } else {
- if(hx>0) {
- __set_errno (ERANGE);
- return tiny*tiny;
- } else
- return two-tiny;
- }
+ /* for 7 <= i < 63, C[i-7] is a degree-7 polynomial approximation of
+ erf(i/16+1/32+x) for -1/32 <= x <= 1/32 */
+ static const double C[56][8] = {
+ { 0x1.f86faa9428f9cp-2, 0x1.cfc41e36c7dfap-1, -0x1.b2c7dc53508b9p-2,
+ -0x1.5a9de93fa556ep-3, 0x1.731793dbb01b5p-3, 0x1.133e06426cf18p-6,
+ -0x1.a12a6289cafd8p-5, 0x1.717d6f1d6f557p-9 },
+ { 0x1.1855a5fd3dd50p-1, 0x1.b3aafcc27502fp-1, -0x1.cee5ac8e92bb2p-2,
+ -0x1.fa02983ca2d79p-4, 0x1.77cd746cb1922p-3, -0x1.fa6f277886487p-10,
+ -0x1.8de75458db416p-5, 0x1.00899c98551c9p-7 },
+ { 0x1.32a54cb8db67ap-1, 0x1.96164fafd8de5p-1, -0x1.e23a7ea0c9ad3p-2,
+ -0x1.3f5ee15671cf4p-4, 0x1.70e468a3d72d9p-3, -0x1.3da68037cfc99p-6,
+ -0x1.69ed9ba1f9839p-5, 0x1.8cab9244a4ff4p-7 },
+ { 0x1.4b13713ad3513p-1, 0x1.7791b886e7405p-1, -0x1.ecef423109bf5p-2,
+ -0x1.15c3c5cec6847p-5, 0x1.5f688fc931ba6p-3, -0x1.1da63ed190037p-5,
+ -0x1.38427ca63cca4p-5, 0x1.fa00e52525e17p-7 },
+ { 0x1.61955607dd15dp-1, 0x1.58a445da7c74ep-1, -0x1.ef6c246a0f66cp-2,
+ 0x1.e83e0d9d61330p-8, 0x1.44cc65535bc9fp-3, -0x1.87d3c4860435dp-5,
+ -0x1.f90b10501169bp-6, 0x1.22295856d427ap-6 },
+ { 0x1.762870f720c6fp-1, 0x1.39ccc1b136d5cp-1, -0x1.ea4feea4e4744p-2,
+ 0x1.715e5952ebfbap-5, 0x1.22cdbd83c75c4p-3, -0x1.da50aa1d925b6p-5,
+ -0x1.754dc0a29b4ddp-6, 0x1.350b6bef9392cp-6 },
+ { 0x1.88d1cd474a2e0p-1, 0x1.1b7e98fe26219p-1, -0x1.de65a22ce1419p-2,
+ 0x1.40686a3f16400p-4, 0x1.f6b0cbb216b2bp-4, -0x1.09c7c903edd57p-4,
+ -0x1.da7529fde641p-7, 0x1.362a7a0588eabp-6 },
+ { 0x1.999d4192a5717p-1, 0x1.fc3ee5d1524b3p-2, -0x1.cc990045b55c8p-2,
+ 0x1.b37338e68b37dp-4, 0x1.a0d120c872ea7p-4, -0x1.19bb2b07ecff6p-4,
+ -0x1.a110f5f593aafp-8, 0x1.272c15a57720ep-6 },
+ { 0x1.a89c850b7d54dp-1, 0x1.c40b0729ed54ap-2, -0x1.b5eaaef0a2346p-2,
+ 0x1.0847c7dacbae1p-3, 0x1.47de0ba6d18fbp-4, -0x1.1d9de77a4b648p-4,
+ 0x1.30ffbe56f0726p-10, 0x1.0a9cb99feea01p-6 },
+ { 0x1.b5e62fce16096p-1, 0x1.8eed36b886d95p-2, -0x1.9b64a06e50705p-2,
+ 0x1.2bb6e2c744df5p-3, 0x1.dee3261ca61bcp-5, -0x1.16996004f7da5p-4,
+ 0x1.fdff37bae983ep-8, 0x1.c750083e65f9ap-7 },
+ { 0x1.c194b1d49a184p-1, 0x1.5d4fd33729015p-2, -0x1.7e0f4f045addbp-2,
+ 0x1.444bc66c31a1bp-3, 0x1.356dbf16ec8f1p-5, -0x1.0643de0906cd8p-4,
+ 0x1.b281af7bd3a2cp-7, 0x1.6b97eaa2c6abdp-7 },
+ { 0x1.cbc54b476248ep-1, 0x1.2f7cc3fe6f423p-2, -0x1.5ee8429e36de8p-2,
+ 0x1.52a8395f96177p-3, 0x1.313761ba257dcp-6, -0x1.dcf844d5fed8fp-5,
+ 0x1.1e1420f475fa9p-6, 0x1.091c7dc1e18b2p-7 },
+ { 0x1.d4970f9ce00d9p-1, 0x1.059f59af7a905p-2, -0x1.3eda354de36c3p-2,
+ 0x1.57b85ad439779p-3, 0x1.8e913b9778136p-10, -0x1.a2893bd3435f4p-5,
+ 0x1.4d3a90e37164ap-6, 0x1.4ce7f6e19a902p-8 },
+ { 0x1.dc29fb60715b0p-1, 0x1.bf8e1b1ca2277p-3, -0x1.1eb7095e5d6d2p-2,
+ 0x1.549ea6f7a64f4p-3, -0x1.b10f12f3877a3p-7, -0x1.61420c8f7156ap-5,
+ 0x1.674f1f92a8812p-6, 0x1.25543ffd74d52p-9 },
+ { 0x1.e29e22a89d767p-1, 0x1.7bd5c7df3fe99p-3, -0x1.fe674494077bfp-3,
+ 0x1.4a9feacf86578p-3, -0x1.a008269076644p-6, -0x1.1cf0e8fb4f1cbp-5,
+ 0x1.6e0d2ef105fb3p-6, -0x1.367205fbd7876p-12 },
+ { 0x1.e812fc64db36ap-1, 0x1.3fda6bc016991p-3, -0x1.c1cb278627920p-3,
+ 0x1.3b10512314f1ep-3, -0x1.1e6457bb1b9a9p-5, -0x1.b1f6474e2388cp-6,
+ 0x1.640a5345f7ec7p-6, -0x1.3dae5a997fdbp-9 },
+ { 0x1.eca6ccd709544p-1, 0x1.0b3f52ce8c380p-3, -0x1.8885019f63c6dp-3,
+ 0x1.274275fc91a05p-3, -0x1.57f73699a8372p-5, -0x1.3076a305fc7cep-6,
+ 0x1.4c6ae04843a41p-6, -0x1.0be5fcf5ecc91p-8 },
+ { 0x1.f0762fde45ee7p-1, 0x1.bb1c972f23e4ap-4, -0x1.5341e3c01b58dp-3,
+ 0x1.107929f6f0b60p-3, -0x1.7e1b34f976c02p-5, -0x1.73b62589c234ap-7,
+ 0x1.2a97ee1876486p-6, -0x1.595f40a3150fep-8 },
+ { 0x1.f39bc242e43e6p-1, 0x1.6c7e64e7281c5p-4, -0x1.2274b86835fd3p-3,
+ 0x1.efb890e5c770dp-4, -0x1.92c7db16847e0p-5, -0x1.45477db5e2dd4p-8,
+ 0x1.01fc6165fc866p-6, -0x1.8845509030c2cp-8 },
+ { 0x1.f62fe80272419p-1, 0x1.297db960e4f5dp-4, -0x1.ecb83b087c04fp-4,
+ 0x1.bce18363ca3d1p-4, -0x1.985aaf776482cp-5, 0x1.cd953efdae886p-12,
+ 0x1.ab9a0b89b54ffp-7, -0x1.9b5e576ccc31cp-8 },
+ { 0x1.f848acb544e95p-1, 0x1.e1d4cf1e24501p-5, -0x1.9e12e1fde5552p-4,
+ 0x1.8a27806df3d1bp-4, -0x1.91674e5eb3319p-5, 0x1.3bc75595b2db8p-8,
+ 0x1.51bc537ac61afp-7, -0x1.96b23b19ea04dp-8 },
+ { 0x1.f9f9ba8d3c733p-1, 0x1.83298d7172108p-5, -0x1.58d101f905a75p-4,
+ 0x1.58f1456f8639bp-4, -0x1.808d1850b8231p-5, 0x1.0c1bd99c348a7p-7,
+ 0x1.f61e9d7bc48cap-8, -0x1.7f07c13441774p-8 },
+ { 0x1.fb54641aebbc9p-1, 0x1.34ac36ad8dafap-5, -0x1.1c8ec267f9405p-4,
+ 0x1.2a52c5d841848p-4, -0x1.68541c02b3b6bp-5, 0x1.5afe400196379p-7,
+ 0x1.565b2d6eda3d6p-8, -0x1.596aaff29e739p-8 },
+ { 0x1.fc67bcf2d7b8fp-1, 0x1.e85c449e377efp-6, -0x1.d177f166c07c6p-5,
+ 0x1.fe23b7584b504p-5, -0x1.4b12109613313p-5, 0x1.8d9905c0acf7dp-7,
+ 0x1.9265032a669dap-9, -0x1.2ac4a6dbcbf3ep-8 },
+ { 0x1.fd40bd6d7a785p-1, 0x1.7f5188610ddc7p-6, -0x1.7954423f7c998p-5,
+ 0x1.af5baae33887fp-5, -0x1.2ad77c7cbc474p-5, 0x1.a7b8c47ec2a51p-7,
+ 0x1.46646ee094bccp-10, -0x1.ef19d8db8673p-9 },
+ { 0x1.fdea6e062d0c9p-1, 0x1.2a875b5ffab58p-6, -0x1.2f3178cd6dcd5p-5,
+ 0x1.68d1c45b94182p-5, -0x1.09648ed3aeaefp-5, 0x1.ad8b150d38164p-7,
+ -0x1.e9a6023d9429fp-13, -0x1.8722d19ee2e8ep-9 },
+ { 0x1.fe6e1742f7cf5p-1, 0x1.cd5ec93c1243ap-7, -0x1.e2ff3aaacb386p-6,
+ 0x1.2aa4e5823cc89p-5, -0x1.d049842dbe399p-6, 0x1.a34edb21ab302p-7,
+ -0x1.676e5996c7f9bp-10, -0x1.23b01a35140bfp-9 },
+ { 0x1.fed37386190fbp-1, 0x1.61beae53b72c2p-7, -0x1.7d6193f22c3c1p-6,
+ 0x1.e947279e3bb7dp-6, -0x1.906031b97ca97p-6, 0x1.8d14d62561755p-7,
+ -0x1.1f245e7178882p-9, -0x1.9257d4eb47685p-10 },
+ { 0x1.ff20e0a7ba8c2p-1, 0x1.0d1d69569b839p-7, -0x1.2a8ca0dc02752p-6,
+ 0x1.8cc071b709751p-6, -0x1.54a149f1b070cp-6, 0x1.6e9137b13412cp-7,
+ -0x1.6577ed3d8e83bp-9, -0x1.e9c1a5178a289p-11 },
+ { 0x1.ff5b8fb26f5f6p-1, 0x1.9646f35a7663cp-8, -0x1.cf68ed9311b0bp-7,
+ 0x1.3e8735b5a694fp-6, -0x1.1e1612d026fdfp-6, 0x1.4afd8e6ca636dp-7,
+ -0x1.8c375170ccb22p-9, -0x1.c799443c4fd3bp-12 },
+ { 0x1.ff87b1913e853p-1, 0x1.30499b5039596p-8, -0x1.64964201ec8bap-7,
+ 0x1.fa73d7eafba98p-7, -0x1.daa3022141fbbp-7, 0x1.2509444c063b7p-7,
+ -0x1.99482a2f8a0a1p-9, -0x1.403d1f76c9454p-15 },
+ { 0x1.ffa89fe5b3625p-1, 0x1.c4412bf4b8f35p-9, -0x1.100f347126cf0p-7,
+ 0x1.8ebda07671d40p-7, -0x1.850c6a31c98c1p-7, 0x1.fdac860c67d21p-8,
+ -0x1.927d03d2ba12cp-9, 0x1.0ff620b4190fep-12 },
+ { 0x1.ffc10194fcb64p-1, 0x1.4d78bba8ca621p-9, -0x1.9ba107a443e02p-8,
+ 0x1.36f273fbc04ccp-7, -0x1.3b38716ac7e6fp-7, 0x1.b3fe0181914acp-8,
+ -0x1.7d3fe7de98c5cp-9, 0x1.ea31f8e5317f7p-12 },
+ { 0x1.ffd2eae369a07p-1, 0x1.e7f232d9e266cp-10, -0x1.34c7442dd48d9p-8,
+ 0x1.e066bed070a0bp-8, -0x1.f914f3c42fc0dp-8, 0x1.6f4664ed2260fp-8,
+ -0x1.5e59910761d24p-9, 0x1.39cbb6e84c126p-11 },
+ { 0x1.ffdff92db56e5p-1, 0x1.6235fbd7a4373p-10, -0x1.cb5e029b9e56ap-9,
+ 0x1.6fa4c7ef274dap-8, -0x1.903a089a835f3p-8, 0x1.30f12e0ca1901p-8,
+ -0x1.39d21b6957f99p-9, 0x1.5d3f8495a703cp-11 },
+ { 0x1.ffe96a78a04a9p-1, 0x1.fe41cd9bb4f2cp-11, -0x1.52d7b28966c0cp-9,
+ 0x1.16c192d86a1a7p-8, -0x1.39bfce951100cp-8, 0x1.f376a7869f9e3p-9,
+ -0x1.12e6cef999c4fp-9, 0x1.66acd4d667b5p-11 },
+ { 0x1.fff0312b010b5p-1, 0x1.6caa0d3583018p-11, -0x1.efb729f4cf75bp-10,
+ 0x1.a2da7cebe12acp-9, -0x1.e6c27a24bc759p-9, 0x1.93b1f4d8ea65p-9,
+ -0x1.d82050aa94a08p-10, 0x1.5cd7dc75d6cbap-11 },
+ { 0x1.fff50456dab8cp-1, 0x1.0295ef6591865p-11, -0x1.679880e95a4dap-10,
+ 0x1.37d38e3a5c8ebp-9, -0x1.75b3708aebb8fp-9, 0x1.4231c4b4b0296p-9,
+ -0x1.8e26476489318p-10, 0x1.45c3b570dd924p-11 },
+ { 0x1.fff86cfd3e657p-1, 0x1.6be02102b353dp-12, -0x1.02b157780d6aep-10,
+ 0x1.cc1d886861133p-10, -0x1.1bff6f12ec9abp-9, 0x1.fc0f77bd9c736p-10,
+ -0x1.4a3320bd0959dp-10, 0x1.267f8b4f95d2p-11 },
+ { 0x1.fffad0b901755p-1, 0x1.fc0d55470cf5ep-13, -0x1.7121aff5e820ep-11,
+ 0x1.506d6992f7de5p-10, -0x1.ab595d3ecd0d6p-10, 0x1.8bdd79daaf754p-10,
+ -0x1.0d9b090f997c1p-10, 0x1.031ab9fd1c7dap-11 },
+ { 0x1.fffc7a37857d2p-1, 0x1.5feada379d8a5p-13, -0x1.05304df58f3aap-11,
+ 0x1.e79c081b8600fp-11, -0x1.3e5dbe33232e0p-10, 0x1.30eb208200729p-10,
+ -0x1.b1d493b147945p-11, 0x1.bd587bbc071bep-12 },
+ { 0x1.fffd9fdeabccep-1, 0x1.e3bcf436a1a49p-14, -0x1.6e953111ef0a1p-12,
+ 0x1.5e3edf6768654p-11, -0x1.d5be67c0547a4p-11, 0x1.d07d9ffa1d435p-11,
+ -0x1.58328f5f358cap-11, 0x1.76d42d95c42c4p-12 },
+ { 0x1.fffe68f4fa777p-1, 0x1.49e17724f4cddp-14, -0x1.fe48c44e229c1p-13,
+ 0x1.f2bd95d76f188p-12, -0x1.57388cb12d011p-11, 0x1.5decc25c5c079p-11,
+ -0x1.0d7499d1b0d2dp-11, 0x1.359332c94ecdcp-12 },
+ { 0x1.fffef1960d85dp-1, 0x1.be6abbb10a4cdp-15, -0x1.6040381a8c313p-13,
+ 0x1.5fff1dde9ee9dp-12, -0x1.f0c933efa9971p-12, 0x1.04cbf4a5cd760p-11,
+ -0x1.a07f150af6dadp-12, 0x1.f68dd183426bap-13 },
+ { 0x1.ffff4db27f146p-1, 0x1.2bb5cc22e5cd8p-15, -0x1.e25894899f526p-14,
+ 0x1.ec8a8e5a72757p-13, -0x1.64256ae0a3cf9p-12, 0x1.80a836c18c46cp-12,
+ -0x1.3dea401af6775p-12, 0x1.915ddff3fe0d1p-13 },
+ { 0x1.ffff8b500e77cp-1, 0x1.8f4ccca7fc769p-16, -0x1.478cffe305946p-14,
+ 0x1.559f04adde504p-13, -0x1.f9e1577d6961dp-13, 0x1.18bda53c14716p-12,
+ -0x1.df8634c35541cp-13, 0x1.3bb5c6b616337p-13 },
+ { 0x1.ffffb43555b5fp-1, 0x1.07ebd2a2d26c8p-16, -0x1.b93e442a37f2bp-15,
+ 0x1.d5cf15159ce28p-14, -0x1.63f5e1469c006p-13, 0x1.95a03acebac18p-13,
+ -0x1.656e5e2a1f8e2p-13, 0x1.e98c437189bdep-14 },
+ { 0x1.ffffcf23ff5fcp-1, 0x1.5a2adfa0b492cp-17, -0x1.26c88270759f0p-15,
+ 0x1.40473572b99a8p-14, -0x1.f057cbde578a5p-14, 0x1.22178d1c3c948p-13,
+ -0x1.0765b61a0d859p-13, 0x1.765b3ea03ddbep-14 },
+ { 0x1.ffffe0bd3e852p-1, 0x1.c282cd3957a72p-18, -0x1.86ad6dfa44faap-16,
+ 0x1.b0f313f03a029p-15, -0x1.56e44abecd255p-14, 0x1.9ad1ecfe34a89p-14,
+ -0x1.7fe4033478618p-14, 0x1.1a8184e049fbfp-14 },
+ { 0x1.ffffec2641a9ep-1, 0x1.22df29821407ep-18, -0x1.00c902a6cfd98p-16,
+ 0x1.22234eb88671fp-15, -0x1.d57a181c9e6e1p-15, 0x1.200c283b54a90p-14,
+ -0x1.14b4c3295a7d0p-14, 0x1.a4f966f713bdep-15 },
+ { 0x1.fffff37d63a36p-1, 0x1.74adc8f405eecp-19, -0x1.4ed4228e44858p-17,
+ 0x1.81918baea92bap-16, -0x1.3e81b17a0009cp-15, 0x1.9004a36116436p-15,
+ -0x1.8aa1ba400e076p-15, 0x1.35cd4e2340a9ep-15 },
+ { 0x1.fffff82cdcf1bp-1, 0x1.d9c73698fa87dp-20, -0x1.b11017ec67115p-18,
+ 0x1.fc0dfadf653f8p-17, -0x1.ac4e03cd2dfc2p-16, 0x1.131806b5abbc5p-15,
+ -0x1.1672ef66fcaafp-15, 0x1.c2882c7debed7p-16 },
+ { 0x1.fffffb248c39dp-1, 0x1.2acee2f5ec66ap-20, -0x1.15cc570408a36p-18,
+ 0x1.4be757bbb75a3p-17, -0x1.1d6aa5f8d2940p-16, 0x1.76c5937d5105ep-16,
+ -0x1.84dffc3ca9302p-16, 0x1.43c8315f2c30ap-16 },
+ { 0x1.fffffd01f36afp-1, 0x1.75fa8dbc840bap-21, -0x1.6186da0133f5ap-19,
+ 0x1.ae023231e1af5p-18, -0x1.790812f7ca394p-17, 0x1.f9c25656d0ef2p-17,
+ -0x1.0cc66682e304cp-16, 0x1.cc170a75d6f9cp-17 },
+ { 0x1.fffffe2ba0ea5p-1, 0x1.d06ad6ecde88ep-22, -0x1.be46aa8edc9a1p-20,
+ 0x1.143860c7840b8p-18, -0x1.edaba78fb1260p-18, 0x1.52138a96ecee2p-17,
+ -0x1.6fca538c4e2eep-17, 0x1.434040640bcefp-17 },
+ { 0x1.fffffee3cc32cp-1, 0x1.1e1e857adb8ddp-22, -0x1.1769ce5f2a6e8p-20,
+ 0x1.5fe5d479b0543p-19, -0x1.405d865c94c2ap-18, 0x1.bfc94feb96afcp-18,
+ -0x1.f245d5f3e8358p-18, 0x1.c142456acf443p-18 },
+ };
+ float ax = fabsf (x);
+ uint32_t ux = asuint (ax);
+ double s = x;
+ double z = ax;
+ /* 0x407ad444 corresponds to x = 0x1.f5a888p+1 = 3.91921..., which is the
+ largest float such that erf(x) does not round to 1 (to nearest). */
+ if (__glibc_unlikely (ux > 0x407ad444u))
+ {
+ float os = copysignf (1.0f, x);
+ if (ux > (0xffu << 23))
+ return x + x; /* nan */
+ if (ux == (0xffu << 23))
+ return os; /* +-inf */
+ return os - 0x1p-25f * os;
+ }
+ double v = floor (16.0 * z);
+ uint32_t i = 16.0f * ax;
+ /* 0x3ee00000 corresponds to x = 0.4375, for smaller x we have i < 7. */
+ if (__glibc_unlikely (ux < 0x3ee00000u))
+ {
+ static const double c[] =
+ {
+ 0x1.20dd750429b6dp+0, -0x1.812746b0375fbp-2,
+ 0x1.ce2f219fd6f45p-4, -0x1.b82ce2cbf0838p-6,
+ 0x1.565bb655adb85p-8, -0x1.c025bfc879c94p-11,
+ 0x1.f81718f61309cp-14, -0x1.cc67bd88f5867p-17
+ };
+ double z2 = s * s, z4 = z2 * z2, z8 = z4 * z4;
+ double c0 = c[0] + z2 * c[1];
+ double c2 = c[2] + z2 * c[3];
+ double c4 = c[4] + z2 * c[5];
+ double c6 = c[6] + z2 * c[7];
+ c0 += z4 * c2;
+ c4 += z4 * c6;
+ c0 += z8 * c4;
+ return s * c0;
+ }
+ z = (z - 0.03125) - 0.0625 * v;
+ const double *c = C[i - 7];
+ double z2 = z * z, z4 = z2 * z2;
+ double c0 = c[0] + z * c[1];
+ double c2 = c[2] + z * c[3];
+ double c4 = c[4] + z * c[5];
+ double c6 = c[6] + z * c[7];
+ c0 += z2 * c2;
+ c4 += z2 * c6;
+ c0 += z4 * c4;
+ return copysign (c0, s);
}
-libm_alias_float (__erfc, erfc)
+libm_alias_float (__erf, erf)
diff --git a/sysdeps/ieee754/flt-32/s_expm1f.c b/sysdeps/ieee754/flt-32/s_expm1f.c
index edd7c9a..a36e578 100644
--- a/sysdeps/ieee754/flt-32/s_expm1f.c
+++ b/sysdeps/ieee754/flt-32/s_expm1f.c
@@ -95,7 +95,7 @@ __expm1f (float x)
return __math_oflowf (0);
}
double a = iln2 * z;
- double ia = roundeven (a);
+ double ia = roundeven_finite (a);
double h = a - ia;
double h2 = h * h;
uint64_t u = asuint64 (ia + big);
diff --git a/sysdeps/ieee754/flt-32/s_tanf.c b/sysdeps/ieee754/flt-32/s_tanf.c
index ae6600b..dfe56fc 100644
--- a/sysdeps/ieee754/flt-32/s_tanf.c
+++ b/sysdeps/ieee754/flt-32/s_tanf.c
@@ -1,76 +1,180 @@
-/* s_tanf.c -- float version of s_tan.c.
- */
+/* Correctly-rounded tangent of binary32 value.
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+Copyright (c) 2022-2024 Alexei Sibidanov.
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_tanf.c,v 1.4 1995/05/10 20:48:20 jtc Exp $";
-#endif
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/tan/tanf.c, revision 59d21d7).
-#include <errno.h>
-#include <math.h>
-#include <math_private.h>
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+#include <array_length.h>
+#include <stdint.h>
#include <libm-alias-float.h>
-#include "s_sincosf.h"
+#include "math_config.h"
+#include <math_uint128.h>
-/* Reduce range of X to a multiple of PI/2. The modulo result is between
- -PI/4 and PI/4 and returned as a high part y[0] and a low part y[1].
- The low bit in the return value indicates the first or 2nd half of tanf. */
-static inline int32_t
-rem_pio2f (float x, float *y)
+/* argument reduction
+ for |z| < 2^28, return r such that 2/pi*x = q + r */
+static inline double
+rltl (float z, int *q)
{
- double dx = x;
- int n;
- const sincos_t *p = &__sincosf_table[0];
+ double x = z;
+ double idl = -0x1.b1bbead603d8bp-32 * x;
+ double idh = 0x1.45f306ep-1 * x;
+ double id = roundeven_finite (idh);
+ *q = (int64_t) id;
+ return (idh - id) + idl;
+}
- if (__glibc_likely (abstop12 (x) < abstop12 (120.0f)))
- dx = reduce_fast (dx, p, &n);
- else
+/* argument reduction
+ same as rltl, but for |x| >= 2^28 */
+static double __attribute__ ((noinline))
+rbig (uint32_t u, int *q)
+{
+ static const uint64_t ipi[] =
{
- uint32_t xi = asuint (x);
- int sign = xi >> 31;
-
- dx = reduce_large (xi, &n);
- dx = sign ? -dx : dx;
+ 0xfe5163abdebbc562, 0xdb6295993c439041,
+ 0xfc2757d1f534ddc0, 0xa2f9836e4e441529
+ };
+ int e = (u >> 23) & 0xff, i;
+ uint64_t m = (u & (~0u >> 9)) | 1 << 23;
+ u128 p0 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[0]));
+ u128 p1 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[1]));
+ p1 = u128_add (p1, u128_rshift (p0, 64));
+ u128 p2 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[2]));
+ p2 = u128_add (p2, u128_rshift (p1, 64));
+ u128 p3 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[3]));
+ p3 = u128_add (p3, u128_rshift (p2, 64));
+ uint64_t p3h = u128_high (p3);
+ uint64_t p3l = u128_low (p3);
+ uint64_t p2l = u128_low (p2);
+ uint64_t p1l = u128_low (p1);
+ int64_t a;
+ int k = e - 127, s = k - 23;
+ /* in ctanf(), rbig() is called in the case 127+28 <= e < 0xff
+ thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
+ if (s < 64)
+ {
+ i = p3h << s | p3l >> (64 - s);
+ a = p3l << s | p2l >> (64 - s);
}
-
- y[0] = dx;
- y[1] = dx - y[0];
- return n;
+ else if (s == 64)
+ {
+ i = p3l;
+ a = p2l;
+ }
+ else
+ { /* s > 64 */
+ i = p3l << (s - 64) | p2l >> (128 - s);
+ a = p2l << (s - 64) | p1l >> (128 - s);
+ }
+ int sgn = u;
+ sgn >>= 31;
+ int64_t sm = a >> 63;
+ i -= sm;
+ double z = (a ^ sgn) * 0x1p-64;
+ i = (i ^ sgn) - sgn;
+ *q = i;
+ return z;
}
-float __tanf(float x)
+float
+__tanf (float x)
{
- float y[2],z=0.0;
- int32_t n, ix;
-
- GET_FLOAT_WORD(ix,x);
-
- /* |x| ~< pi/4 */
- ix &= 0x7fffffff;
- if(ix <= 0x3f490fda) return __kernel_tanf(x,z,1);
-
- /* tan(Inf or NaN) is NaN */
- else if (ix>=0x7f800000) {
- if (ix==0x7f800000)
- __set_errno (EDOM);
- return x-x; /* NaN */
+ uint32_t t = asuint (x);
+ int e = (t >> 23) & 0xff;
+ int i;
+ double z;
+ if (__glibc_likely (e < 127 + 28)) /* |x| < 2^28 */
+ {
+ if (__glibc_unlikely (e < 115))
+ {
+ if (__glibc_unlikely (e < 102))
+ return fmaf (x, fabsf (x), x);
+ float x2 = x * x;
+ return fmaf (x, 0x1.555556p-2f * x2, x);
}
-
- /* argument reduction needed */
- else {
- n = rem_pio2f(x,y);
- return __kernel_tanf(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even
- -1 -- n odd */
+ z = rltl (x, &i);
+ }
+ else if (e < 0xff)
+ z = rbig (t, &i);
+ else
+ {
+ if (t << 9)
+ return x + x; /* nan */
+ return __math_invalidf (x);
+ }
+ double z2 = z * z;
+ double z4 = z2 * z2;
+ static const double cn[] =
+ {
+ 0x1.921fb54442d18p+0, -0x1.fd226e573289fp-2,
+ 0x1.b7a60c8dac9f6p-6, -0x1.725beb40f33e5p-13
+ };
+ static const double cd[] =
+ {
+ 0x1p+0, -0x1.2395347fb829dp+0,
+ 0x1.2313660f29c36p-3, -0x1.9a707ab98d1c1p-9
+ };
+ static const double s[] = { 0, 1 };
+ double n = cn[0] + z2 * cn[1];
+ double n2 = cn[2] + z2 * cn[3];
+ n += z4 * n2;
+ double d = cd[0] + z2 * cd[1];
+ double d2 = cd[2] + z2 * cd[3];
+ d += z4 * d2;
+ n *= z;
+ double s0 = s[i & 1];
+ double s1 = s[1 - (i & 1)];
+ double r1 = (n * s1 - d * s0) / (n * s0 + d * s1);
+ uint64_t tail = (asuint64 (r1) + 7) & (~UINT64_C(0) >> 35);
+ if (__glibc_unlikely (tail <= 14))
+ {
+ static const struct
+ {
+ float arg;
+ float rh;
+ float rl;
+ } st[] = {
+ { 0x1.143ec4p+0f, 0x1.ddf9f6p+0f, -0x1.891d24p-52f },
+ { 0x1.ada6aap+27f, 0x1.e80304p-3f, 0x1.419f46p-58f },
+ { 0x1.af61dap+48f, 0x1.60d1c8p-2f, -0x1.2d6c3ap-55f },
+ { 0x1.0088bcp+52f, 0x1.ca1edp+0f, 0x1.f6053p-53f },
+ { 0x1.f90dfcp+72f, 0x1.597f9cp-1f, 0x1.925978p-53f },
+ { 0x1.cc4e22p+85f, -0x1.f33584p+1f, 0x1.d7254ap-51f },
+ { 0x1.a6ce12p+86f, -0x1.c5612ep-1f, -0x1.26c33ep-53f },
+ { 0x1.6a0b76p+102f, -0x1.e42a1ep+0f, -0x1.1dc906p-52f },
+ };
+ uint32_t ax = t & (~0u >> 1);
+ uint32_t sgn = t >> 31;
+ for (int j = 0; j < array_length (st); j++)
+ {
+ if (__glibc_unlikely (asfloat (st[j].arg) == ax))
+ {
+ if (sgn)
+ return -st[j].rh - st[j].rl;
+ else
+ return st[j].rh + st[j].rl;
+ }
}
+ }
+ return r1;
}
libm_alias_float (__tan, tan)
diff --git a/sysdeps/ieee754/ldbl-128/s_erfcl.c b/sysdeps/ieee754/ldbl-128/s_erfcl.c
new file mode 100644
index 0000000..95d17c8
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_erfcl.c
@@ -0,0 +1 @@
+/* Not required. */
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfcl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfcl.c
new file mode 100644
index 0000000..95d17c8
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128ibm/s_erfcl.c
@@ -0,0 +1 @@
+/* Not required. */
diff --git a/sysdeps/ieee754/ldbl-96/s_erfcl.c b/sysdeps/ieee754/ldbl-96/s_erfcl.c
new file mode 100644
index 0000000..95d17c8
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-96/s_erfcl.c
@@ -0,0 +1 @@
+/* Not required. */
generated by cgit v1.1 at 2024-12-03 17:48:29 +0000