Sync upstream FreeBSD libm (real changes).

Test: treehugger
Change-Id: Icf591ba195a3f5080203b157aa7f43d518a9cc69
This commit is contained in:
Elliott Hughes 2023-07-19 14:11:58 -07:00
parent a2cdb247e1
commit 8810bd7585
7 changed files with 14 additions and 12 deletions

View file

@ -22,15 +22,15 @@ __FBSDID("$FreeBSD$");
* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
* Note 2. About jn(n,x), yn(n,x) * Note 2. About jn(n,x), yn(n,x)
* For n=0, j0(x) is called, * For n=0, j0(x) is called.
* for n=1, j1(x) is called, * For n=1, j1(x) is called.
* for n<x, forward recursion us used starting * For n<x, forward recursion is used starting
* from values of j0(x) and j1(x). * from values of j0(x) and j1(x).
* for n>x, a continued fraction approximation to * For n>x, a continued fraction approximation to
* j(n,x)/j(n-1,x) is evaluated and then backward * j(n,x)/j(n-1,x) is evaluated and then backward
* recursion is used starting from a supposed value * recursion is used starting from a supposed value
* for j(n,x). The resulting value of j(0,x) is * for j(n,x). The resulting values of j(0,x) or j(1,x) are
* compared with the actual value to correct the * compared with the actual values to correct the
* supposed value of j(n,x). * supposed value of j(n,x).
* *
* yn(n,x) is similar in all respects, except * yn(n,x) is similar in all respects, except

View file

@ -27,7 +27,7 @@ __FBSDID("$FreeBSD$");
* = log(6.3*5.3) + lgamma(5.3) * = log(6.3*5.3) + lgamma(5.3)
* = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
* 2. Polynomial approximation of lgamma around its * 2. Polynomial approximation of lgamma around its
* minimun ymin=1.461632144968362245 to maintain monotonicity. * minimum ymin=1.461632144968362245 to maintain monotonicity.
* On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
* Let z = x-ymin; * Let z = x-ymin;
* lgamma(x) = -1.214862905358496078218 + z^2*poly(z) * lgamma(x) = -1.214862905358496078218 + z^2*poly(z)

View file

@ -50,7 +50,7 @@ __kernel_tandf(double x, int iy)
* We add the small terms from lowest degree up for efficiency on * We add the small terms from lowest degree up for efficiency on
* non-sequential machines (the lowest degree terms tend to be ready * non-sequential machines (the lowest degree terms tend to be ready
* earlier). Apart from this, we don't care about order of * earlier). Apart from this, we don't care about order of
* operations, and don't need to to care since we have precision to * operations, and don't need to care since we have precision to
* spare. However, the chosen splitting is good for accuracy too, * spare. However, the chosen splitting is good for accuracy too,
* and would give results as accurate as Horner's method if the * and would give results as accurate as Horner's method if the
* small terms were added from highest degree down. * small terms were added from highest degree down.

View file

@ -108,7 +108,7 @@ cbrt(double x)
r=x/s; /* error <= 0.5 ulps; |r| < |t| */ r=x/s; /* error <= 0.5 ulps; |r| < |t| */
w=t+t; /* t+t is exact */ w=t+t; /* t+t is exact */
r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */ r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ t=t+t*r; /* error <= (0.5 + 0.5/3) * ulp */
return(t); return(t);
} }

View file

@ -136,7 +136,7 @@ cbrtl(long double x)
r=x/s; /* error <= 0.5 ulps; |r| < |t| */ r=x/s; /* error <= 0.5 ulps; |r| < |t| */
w=t+t; /* t+t is exact */ w=t+t; /* t+t is exact */
r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */ r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ t=t+t*r; /* error <= (0.5 + 0.5/3) * ulp */
t *= v.e; t *= v.e;
RETURNI(t); RETURNI(t);

View file

@ -138,7 +138,8 @@ cospi(double x)
return (j0 & 1 ? -c : c); return (j0 & 1 ? -c : c);
} }
if (ix >= 0x7f800000) /* x = +-inf or nan. */
if (ix >= 0x7ff00000)
return (vzero / vzero); return (vzero / vzero);
/* /*

View file

@ -155,7 +155,8 @@ sinpi(double x)
return ((hx & 0x80000000) ? -s : s); return ((hx & 0x80000000) ? -s : s);
} }
if (ix >= 0x7f800000) /* x = +-inf or nan. */
if (ix >= 0x7ff00000)
return (vzero / vzero); return (vzero / vzero);
/* /*