tequilaOS/platform_bionic
2
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Merge "Sync upstream FreeBSD libm." into main

Browse source
This commit is contained in:
Elliott Hughes 2023-08-04 14:51:38 +00:00 committed by Gerrit Code Review
commit 9a80691e38
56 changed files with 245 additions and 385 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

28
libm/upstream-freebsd/lib/msun/ld128/s_expl.c
View file

@ -65,8 +65,6 @@ expl(long double x)
int k;
uint16_t hx, ix;
DOPRINT_START(&x);
/* Filter out exceptional cases. */
u.e = x;
hx = u.xbits.expsign;
@ -74,15 +72,15 @@ expl(long double x)
if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */
if (ix == BIAS + LDBL_MAX_EXP) {
if (hx & 0x8000) /* x is -Inf or -NaN */
RETURNP(-1 / x);
RETURNP(x + x); /* x is +Inf or +NaN */
RETURNF(-1 / x);
RETURNF(x + x); /* x is +Inf or +NaN */
}
if (x > o_threshold)
RETURNP(huge * huge);
RETURNF(huge * huge);
if (x < u_threshold)
RETURNP(tiny * tiny);
RETURNF(tiny * tiny);
} else if (ix < BIAS - 114) { /* |x| < 0x1p-114 */
RETURN2P(1, x); /* 1 with inexact iff x != 0 */
RETURNF(1 + x); /* 1 with inexact iff x != 0 */
}
ENTERI();
@ -210,8 +208,6 @@ expm1l(long double x)
int k, n, n2;
uint16_t hx, ix;
DOPRINT_START(&x);
/* Filter out exceptional cases. */
u.e = x;
hx = u.xbits.expsign;
@ -219,11 +215,11 @@ expm1l(long double x)
if (ix >= BIAS + 7) { /* |x| >= 128 or x is NaN */
if (ix == BIAS + LDBL_MAX_EXP) {
if (hx & 0x8000) /* x is -Inf or -NaN */
RETURNP(-1 / x - 1);
RETURNP(x + x); /* x is +Inf or +NaN */
RETURNF(-1 / x - 1);
RETURNF(x + x); /* x is +Inf or +NaN */
}
if (x > o_threshold)
RETURNP(huge * huge);
RETURNF(huge * huge);
/*
* expm1l() never underflows, but it must avoid
* unrepresentable large negative exponents. We used a
@ -232,7 +228,7 @@ expm1l(long double x)
* in the same way as large ones here.
*/
if (hx & 0x8000) /* x <= -128 */
RETURN2P(tiny, -1); /* good for x < -114ln2 - eps */
RETURNF(tiny - 1); /* good for x < -114ln2 - eps */
}
ENTERI();
@ -244,7 +240,7 @@ expm1l(long double x)
if (x < T3) {
if (ix < BIAS - 113) { /* |x| < 0x1p-113 */
/* x (rounded) with inexact if x != 0: */
RETURNPI(x == 0 ? x :
RETURNI(x == 0 ? x :
(0x1p200 * x + fabsl(x)) * 0x1p-200);
}
q = x * x2 * C3 + x2 * x2 * (C4 + x * (C5 + x * (C6 +
@ -265,9 +261,9 @@ expm1l(long double x)
hx2_hi = x_hi * x_hi / 2;
hx2_lo = x_lo * (x + x_hi) / 2;
if (ix >= BIAS - 7)
RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q);
RETURNI((hx2_hi + x_hi) + (hx2_lo + x_lo + q));
else
RETURN2PI(x, hx2_lo + q + hx2_hi);
RETURNI(x + (hx2_lo + q + hx2_hi));
}
/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */

24
libm/upstream-freebsd/lib/msun/ld128/s_logl.c
View file

@ -573,24 +573,23 @@ log1pl(long double x)
int i, k;
int16_t ax, hx;
DOPRINT_START(&x);
EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
if (hx < 0x3fff) { /* x < 1, or x neg NaN */
ax = hx & 0x7fff;
if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */
if (ax == 0x3fff && (lx | llx) == 0)
RETURNP(-1 / zero); /* log1p(-1) = -Inf */
RETURNF(-1 / zero); /* log1p(-1) = -Inf */
/* log1p(x < 1, or x NaN) = qNaN: */
RETURNP((x - x) / (x - x));
RETURNF((x - x) / (x - x));
}
if (ax <= 0x3f8d) { /* |x| < 2**-113 */
if ((int)x == 0)
RETURNP(x); /* x with inexact if x != 0 */
RETURNF(x); /* x with inexact if x != 0 */
}
f_hi = 1;
f_lo = x;
} else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */
RETURNP(x + x); /* log1p(Inf or NaN) = Inf or qNaN */
RETURNF(x + x); /* log1p(Inf or NaN) = Inf or qNaN */
} else if (hx < 0x40e1) { /* 1 <= x < 2**226 */
f_hi = x;
f_lo = 1;
@ -669,7 +668,7 @@ log1pl(long double x)
#endif
_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
RETURN2PI(val_hi, val_lo);
RETURNI(val_hi + val_lo);
}
#ifdef STRUCT_RETURN
@ -680,7 +679,6 @@ logl(long double x)
struct ld r;
ENTERI();
DOPRINT_START(&x);
k_logl(x, &r);
RETURNSPI(&r);
}
@ -708,15 +706,13 @@ log10l(long double x)
long double hi, lo;
ENTERI();
DOPRINT_START(&x);
k_logl(x, &r);
if (!r.lo_set)
RETURNPI(r.hi);
RETURNI(r.hi);
_2sumF(r.hi, r.lo);
hi = (float)r.hi;
lo = r.lo + (r.hi - hi);
RETURN2PI(invln10_hi * hi,
invln10_lo_plus_hi * lo + invln10_lo * hi);
RETURNI(invln10_hi * hi + (invln10_lo_plus_hi * lo + invln10_lo * hi));
}
long double
@ -726,15 +722,13 @@ log2l(long double x)
long double hi, lo;
ENTERI();
DOPRINT_START(&x);
k_logl(x, &r);
if (!r.lo_set)
RETURNPI(r.hi);
RETURNI(r.hi);
_2sumF(r.hi, r.lo);
hi = (float)r.hi;
lo = r.lo + (r.hi - hi);
RETURN2PI(invln2_hi * hi,
invln2_lo_plus_hi * lo + invln2_lo * hi);
RETURNI(invln2_hi * hi + (invln2_lo_plus_hi * lo + invln2_lo * hi));
}
#endif /* STRUCT_RETURN */

4
libm/upstream-freebsd/lib/msun/src/e_acos.c
View file

@ -14,7 +14,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_acos(x)
/* acos(x)
* Method :
* acos(x) = pi/2 - asin(x)
* acos(-x) = pi/2 + asin(x)
@ -62,7 +62,7 @@ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
double
__ieee754_acos(double x)
acos(double x)
{
double z,p,q,r,w,s,c,df;
int32_t hx,ix;

2
libm/upstream-freebsd/lib/msun/src/e_acosf.c
View file

@ -32,7 +32,7 @@ pS2 = -8.6563630030e-03,
qS1 = -7.0662963390e-01;
float
__ieee754_acosf(float x)
acosf(float x)
{
float z,p,q,r,w,s,c,df;
int32_t hx,ix;

8
libm/upstream-freebsd/lib/msun/src/e_acosh.c
View file

@ -15,7 +15,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_acosh(x)
/* acosh(x)
* Method :
* Based on
* acosh(x) = log [ x + sqrt(x*x-1) ]
@ -39,7 +39,7 @@ one = 1.0,
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double
__ieee754_acosh(double x)
acosh(double x)
{
double t;
int32_t hx;
@ -51,12 +51,12 @@ __ieee754_acosh(double x)
if(hx >=0x7ff00000) { /* x is inf of NaN */
return x+x;
} else
return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */
return log(x)+ln2; /* acosh(huge)=log(2x) */
} else if(((hx-0x3ff00000)|lx)==0) {
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t=x*x;
return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
return log(2.0*x-one/(x+sqrt(t-one)));
} else { /* 1<x<2 */
t = x-one;
return log1p(t+sqrt(2.0*t+t*t));

8
libm/upstream-freebsd/lib/msun/src/e_acoshf.c
View file

@ -24,7 +24,7 @@ one = 1.0,
ln2 = 6.9314718246e-01; /* 0x3f317218 */
float
__ieee754_acoshf(float x)
acoshf(float x)
{
float t;
int32_t hx;
@ -35,14 +35,14 @@ __ieee754_acoshf(float x)
if(hx >=0x7f800000) { /* x is inf of NaN */
return x+x;
} else
return __ieee754_logf(x)+ln2; /* acosh(huge)=log(2x) */
return logf(x)+ln2; /* acosh(huge)=log(2x) */
} else if (hx==0x3f800000) {
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t=x*x;
return __ieee754_logf((float)2.0*x-one/(x+__ieee754_sqrtf(t-one)));
return logf((float)2.0*x-one/(x+sqrtf(t-one)));
} else { /* 1<x<2 */
t = x-one;
return log1pf(t+__ieee754_sqrtf((float)2.0*t+t*t));
return log1pf(t+sqrtf((float)2.0*t+t*t));
}
}

4
libm/upstream-freebsd/lib/msun/src/e_asin.c
View file

@ -14,7 +14,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_asin(x)
/* asin(x)
* Method :
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
* we approximate asin(x) on [0,0.5] by
@ -68,7 +68,7 @@ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
double
__ieee754_asin(double x)
asin(double x)
{
double t=0.0,w,p,q,c,r,s;
int32_t hx,ix;

2
libm/upstream-freebsd/lib/msun/src/e_asinf.c
View file

@ -32,7 +32,7 @@ static const double
pio2 = 1.570796326794896558e+00;
float
__ieee754_asinf(float x)
asinf(float x)
{
double s;
float t,w,p,q;

4
libm/upstream-freebsd/lib/msun/src/e_atan2.c
View file

@ -15,7 +15,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_atan2(y,x)
/* atan2(y,x)
* Method :
* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
* 2. Reduce x to positive by (if x and y are unexceptional):
@ -58,7 +58,7 @@ static volatile double
pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
double
__ieee754_atan2(double y, double x)
atan2(double y, double x)
{
double z;
int32_t k,m,hx,hy,ix,iy;

2
libm/upstream-freebsd/lib/msun/src/e_atan2f.c
View file

@ -30,7 +30,7 @@ static volatile float
pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
float
__ieee754_atan2f(float y, float x)
atan2f(float y, float x)
{
float z;
int32_t k,m,hx,hy,ix,iy;

4
libm/upstream-freebsd/lib/msun/src/e_atanh.c
View file

@ -15,7 +15,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_atanh(x)
/* atanh(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
@ -42,7 +42,7 @@ static const double one = 1.0, huge = 1e300;
static const double zero = 0.0;
double
__ieee754_atanh(double x)
atanh(double x)
{
double t;
int32_t hx,ix;

2
libm/upstream-freebsd/lib/msun/src/e_atanhf.c
View file

@ -24,7 +24,7 @@ static const float one = 1.0, huge = 1e30;
static const float zero = 0.0;
float
__ieee754_atanhf(float x)
atanhf(float x)
{
float t;
int32_t hx,ix;

8
libm/upstream-freebsd/lib/msun/src/e_cosh.c
View file

@ -14,7 +14,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_cosh(x)
/* cosh(x)
* Method :
* mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
* 1. Replace x by |x| (cosh(x) = cosh(-x)).
@ -43,7 +43,7 @@ __FBSDID("$FreeBSD$");
static const double one = 1.0, half=0.5, huge = 1.0e300;
double
__ieee754_cosh(double x)
cosh(double x)
{
double t,w;
int32_t ix;
@ -65,12 +65,12 @@ __ieee754_cosh(double x)
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ix < 0x40360000) {
t = __ieee754_exp(fabs(x));
t = exp(fabs(x));
return half*t+half/t;
}
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x));
if (ix < 0x40862E42) return half*exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
if (ix<=0x408633CE)

6
libm/upstream-freebsd/lib/msun/src/e_coshf.c
View file

@ -22,7 +22,7 @@ __FBSDID("$FreeBSD$");
static const float one = 1.0, half=0.5, huge = 1.0e30;
float
__ieee754_coshf(float x)
coshf(float x)
{
float t,w;
int32_t ix;
@ -43,12 +43,12 @@ __ieee754_coshf(float x)
/* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x41100000) {
t = __ieee754_expf(fabsf(x));
t = expf(fabsf(x));
return half*t+half/t;
}
/* |x| in [9, log(maxfloat)] return half*exp(|x|) */
if (ix < 0x42b17217) return half*__ieee754_expf(fabsf(x));
if (ix < 0x42b17217) return half*expf(fabsf(x));
/* |x| in [log(maxfloat), overflowthresold] */
if (ix<=0x42b2d4fc)

4
libm/upstream-freebsd/lib/msun/src/e_exp.c
View file

@ -13,7 +13,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_exp(x)
/* exp(x)
* Returns the exponential of x.
*
* Method
@ -102,7 +102,7 @@ huge = 1.0e+300,
twom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/
double
__ieee754_exp(double x) /* default IEEE double exp */
exp(double x) /* default IEEE double exp */
{
double y,hi=0.0,lo=0.0,c,t,twopk;
int32_t k=0,xsb;

2
libm/upstream-freebsd/lib/msun/src/e_expf.c
View file

@ -43,7 +43,7 @@ huge = 1.0e+30,
twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */
float
__ieee754_expf(float x)
expf(float x)
{
float y,hi=0.0,lo=0.0,c,t,twopk;
int32_t k=0,xsb;

4
libm/upstream-freebsd/lib/msun/src/e_fmod.c
View file

@ -15,7 +15,7 @@
__FBSDID("$FreeBSD$");
/*
* __ieee754_fmod(x,y)
* fmod(x,y)
* Return x mod y in exact arithmetic
* Method: shift and subtract
*/
@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
static const double one = 1.0, Zero[] = {0.0, -0.0,};
double
__ieee754_fmod(double x, double y)
fmod(double x, double y)
{
int32_t n,hx,hy,hz,ix,iy,sx,i;
u_int32_t lx,ly,lz;

4
libm/upstream-freebsd/lib/msun/src/e_fmodf.c
View file

@ -17,7 +17,7 @@
__FBSDID("$FreeBSD$");
/*
* __ieee754_fmodf(x,y)
* fmodf(x,y)
* Return x mod y in exact arithmetic
* Method: shift and subtract
*/
@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
static const float one = 1.0, Zero[] = {0.0, -0.0,};
float
__ieee754_fmodf(float x, float y)
fmodf(float x, float y)
{
int32_t n,hx,hy,hz,ix,iy,sx,i;

8
libm/upstream-freebsd/lib/msun/src/e_gamma.c
View file

@ -15,10 +15,10 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_gamma(x)
/* gamma(x)
* Return the logarithm of the Gamma function of x.
*
* Method: call __ieee754_gamma_r
* Method: call gamma_r
*/
#include "math.h"
@ -27,7 +27,7 @@ __FBSDID("$FreeBSD$");
extern int signgam;
double
__ieee754_gamma(double x)
gamma(double x)
{
return __ieee754_gamma_r(x,&signgam);
return gamma_r(x,&signgam);
}

8
libm/upstream-freebsd/lib/msun/src/e_gamma_r.c
View file

@ -15,18 +15,18 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_gamma_r(x, signgamp)
/* gamma_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
* Method: See __ieee754_lgamma_r
* Method: See lgamma_r
*/
#include "math.h"
#include "math_private.h"
double
__ieee754_gamma_r(double x, int *signgamp)
gamma_r(double x, int *signgamp)
{
return __ieee754_lgamma_r(x,signgamp);
return lgamma_r(x,signgamp);
}

8
libm/upstream-freebsd/lib/msun/src/e_gammaf.c
View file

@ -16,10 +16,10 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_gammaf(x)
/* gammaf(x)
* Return the logarithm of the Gamma function of x.
*
* Method: call __ieee754_gammaf_r
* Method: call gammaf_r
*/
#include "math.h"
@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
extern int signgam;
float
__ieee754_gammaf(float x)
gammaf(float x)
{
return __ieee754_gammaf_r(x,&signgam);
return gammaf_r(x,&signgam);
}

8
libm/upstream-freebsd/lib/msun/src/e_gammaf_r.c
View file

@ -16,18 +16,18 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_gammaf_r(x, signgamp)
/* gammaf_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
* Method: See __ieee754_lgammaf_r
* Method: See lgammaf_r
*/
#include "math.h"
#include "math_private.h"
float
__ieee754_gammaf_r(float x, int *signgamp)
gammaf_r(float x, int *signgamp)
{
return __ieee754_lgammaf_r(x,signgamp);
return lgammaf_r(x,signgamp);
}

4
libm/upstream-freebsd/lib/msun/src/e_hypot.c
View file

@ -14,7 +14,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_hypot(x,y)
/* hypot(x,y)
*
* Method :
* If (assume round-to-nearest) z=x*x+y*y
@ -52,7 +52,7 @@ __FBSDID("$FreeBSD$");
#include "math_private.h"
double
__ieee754_hypot(double x, double y)
hypot(double x, double y)
{
double a,b,t1,t2,y1,y2,w;
int32_t j,k,ha,hb;

6
libm/upstream-freebsd/lib/msun/src/e_hypotf.c
View file

@ -20,7 +20,7 @@ __FBSDID("$FreeBSD$");
#include "math_private.h"
float
__ieee754_hypotf(float x, float y)
hypotf(float x, float y)
{
float a,b,t1,t2,y1,y2,w;
int32_t j,k,ha,hb;
@ -67,14 +67,14 @@ __ieee754_hypotf(float x, float y)
if (w>b) {
SET_FLOAT_WORD(t1,ha&0xfffff000);
t2 = a-t1;
w = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
w = sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
SET_FLOAT_WORD(y1,hb&0xfffff000);
y2 = b - y1;
SET_FLOAT_WORD(t1,(ha+0x00800000)&0xfffff000);
t2 = a - t1;
w = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
w = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
SET_FLOAT_WORD(t1,(127+k)<<23);

10
libm/upstream-freebsd/lib/msun/src/e_j0.c
View file

@ -13,7 +13,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_j0(x), __ieee754_y0(x)
/* j0(x), y0(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j0(x):
* 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
@ -83,7 +83,7 @@ S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
static const double zero = 0, qrtr = 0.25;
double
__ieee754_j0(double x)
j0(double x)
{
double z, s,c,ss,cc,r,u,v;
int32_t hx,ix;
@ -143,7 +143,7 @@ v03 = 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
double
__ieee754_y0(double x)
y0(double x)
{
double z, s,c,ss,cc,u,v;
int32_t hx,ix,lx;
@ -192,12 +192,12 @@ __ieee754_y0(double x)
return z;
}
if(ix<=0x3e400000) { /* x < 2**-27 */
return(u00 + tpi*__ieee754_log(x));
return(u00 + tpi*log(x));
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = one+z*(v01+z*(v02+z*(v03+z*v04)));
return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
return(u/v + tpi*(j0(x)*log(x)));
}
/* The asymptotic expansions of pzero is

8
libm/upstream-freebsd/lib/msun/src/e_j0f.c
View file

@ -45,7 +45,7 @@ S04 = 1.1661400734e-09; /* 0x30a045e8 */
static const float zero = 0, qrtr = 0.25;
float
__ieee754_j0f(float x)
j0f(float x)
{
float z, s,c,ss,cc,r,u,v;
int32_t hx,ix;
@ -105,7 +105,7 @@ v03 = 2.5915085189e-07, /* 0x348b216c */
v04 = 4.4111031494e-10; /* 0x2ff280c2 */
float
__ieee754_y0f(float x)
y0f(float x)
{
float z, s,c,ss,cc,u,v;
int32_t hx,ix;
@ -147,12 +147,12 @@ __ieee754_y0f(float x)
return z;
}
if(ix<=0x39000000) { /* x < 2**-13 */
return(u00 + tpi*__ieee754_logf(x));
return(u00 + tpi*logf(x));
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = one+z*(v01+z*(v02+z*(v03+z*v04)));
return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x)));
return(u/v + tpi*(j0f(x)*logf(x)));
}
/* The asymptotic expansions of pzero is

8
libm/upstream-freebsd/lib/msun/src/e_j1.c
View file

@ -13,7 +13,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_j1(x), __ieee754_y1(x)
/* j1(x), y1(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j1(x):
* 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
@ -84,7 +84,7 @@ s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
static const double zero = 0.0;
double
__ieee754_j1(double x)
j1(double x)
{
double z, s,c,ss,cc,r,u,v,y;
int32_t hx,ix;
@ -140,7 +140,7 @@ static const double V0[5] = {
};
double
__ieee754_y1(double x)
y1(double x)
{
double z, s,c,ss,cc,u,v;
int32_t hx,ix,lx;
@ -190,7 +190,7 @@ __ieee754_y1(double x)
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
return(x*(u/v) + tpi*(j1(x)*log(x)-one/x));
}
/* For x >= 8, the asymptotic expansions of pone is

6
libm/upstream-freebsd/lib/msun/src/e_j1f.c
View file

@ -46,7 +46,7 @@ s05 = 1.2354227016e-11; /* 0x2d59567e */
static const float zero = 0.0;
float
__ieee754_j1f(float x)
j1f(float x)
{
float z, s,c,ss,cc,r,u,v,y;
int32_t hx,ix;
@ -102,7 +102,7 @@ static const float V0[5] = {
};
float
__ieee754_y1f(float x)
y1f(float x)
{
float z, s,c,ss,cc,u,v;
int32_t hx,ix;
@ -145,7 +145,7 @@ __ieee754_y1f(float x)
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return(x*(u/v) + tpi*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x));
return(x*(u/v) + tpi*(j1f(x)*logf(x)-one/x));
}
/* For x >= 8, the asymptotic expansions of pone is

28
libm/upstream-freebsd/lib/msun/src/e_jn.c
View file

@ -14,7 +14,7 @@
__FBSDID("$FreeBSD$");
/*
* __ieee754_jn(n, x), __ieee754_yn(n, x)
* jn(n, x), yn(n, x)
* floating point Bessel's function of the 1st and 2nd kind
* of order n
*
@ -51,7 +51,7 @@ one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
static const double zero = 0.00000000000000000000e+00;
double
__ieee754_jn(int n, double x)
jn(int n, double x)
{
int32_t i,hx,ix,lx, sgn;
double a, b, c, s, temp, di;
@ -69,8 +69,8 @@ __ieee754_jn(int n, double x)
x = -x;
hx ^= 0x80000000;
}
if(n==0) return(__ieee754_j0(x));
if(n==1) return(__ieee754_j1(x));
if(n==0) return(j0(x));
if(n==1) return(j1(x));
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabs(x);
if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
@ -100,8 +100,8 @@ __ieee754_jn(int n, double x)
}
b = invsqrtpi*temp/sqrt(x);
} else {
a = __ieee754_j0(x);
b = __ieee754_j1(x);
a = j0(x);
b = j1(x);
for(i=1;i<n;i++){
temp = b;
b = b*((double)(i+i)/x) - a; /* avoid underflow */
@ -177,7 +177,7 @@ __ieee754_jn(int n, double x)
*/
tmp = n;
v = two/x;
tmp = tmp*__ieee754_log(fabs(v*tmp));
tmp = tmp*log(fabs(v*tmp));
if(tmp<7.09782712893383973096e+02) {
for(i=n-1,di=(double)(i+i);i>0;i--){
temp = b;
@ -201,8 +201,8 @@ __ieee754_jn(int n, double x)
}
}
}
z = __ieee754_j0(x);
w = __ieee754_j1(x);
z = j0(x);
w = j1(x);
if (fabs(z) >= fabs(w))
b = (t*z/b);
else
@ -213,7 +213,7 @@ __ieee754_jn(int n, double x)
}
double
__ieee754_yn(int n, double x)
yn(int n, double x)
{
int32_t i,hx,ix,lx;
int32_t sign;
@ -232,8 +232,8 @@ __ieee754_yn(int n, double x)
n = -n;
sign = 1 - ((n&1)<<1);
}
if(n==0) return(__ieee754_y0(x));
if(n==1) return(sign*__ieee754_y1(x));
if(n==0) return(y0(x));
if(n==1) return(sign*y1(x));
if(ix==0x7ff00000) return zero;
if(ix>=0x52D00000) { /* x > 2**302 */
/* (x >> n**2)
@ -259,8 +259,8 @@ __ieee754_yn(int n, double x)
b = invsqrtpi*temp/sqrt(x);
} else {
u_int32_t high;
a = __ieee754_y0(x);
b = __ieee754_y1(x);
a = y0(x);
b = y1(x);
/* quit if b is -inf */
GET_HIGH_WORD(high,b);
for(i=1;i<n&&high!=0xfff00000;i++){

26
libm/upstream-freebsd/lib/msun/src/e_jnf.c
View file

@ -32,7 +32,7 @@ one = 1.0000000000e+00; /* 0x3F800000 */
static const float zero = 0.0000000000e+00;
float
__ieee754_jnf(int n, float x)
jnf(int n, float x)
{
int32_t i,hx,ix, sgn;
float a, b, temp, di;
@ -50,16 +50,16 @@ __ieee754_jnf(int n, float x)
x = -x;
hx ^= 0x80000000;
}
if(n==0) return(__ieee754_j0f(x));
if(n==1) return(__ieee754_j1f(x));
if(n==0) return(j0f(x));
if(n==1) return(j1f(x));
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x);
if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */
b = zero;
else if((float)n<=x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = __ieee754_j0f(x);
b = __ieee754_j1f(x);
a = j0f(x);
b = j1f(x);
for(i=1;i<n;i++){
temp = b;
b = b*((float)(i+i)/x) - a; /* avoid underflow */
@ -134,7 +134,7 @@ __ieee754_jnf(int n, float x)
*/
tmp = n;
v = two/x;
tmp = tmp*__ieee754_logf(fabsf(v*tmp));
tmp = tmp*logf(fabsf(v*tmp));
if(tmp<(float)8.8721679688e+01) {
for(i=n-1,di=(float)(i+i);i>0;i--){
temp = b;
@ -158,8 +158,8 @@ __ieee754_jnf(int n, float x)
}
}
}
z = __ieee754_j0f(x);
w = __ieee754_j1f(x);
z = j0f(x);
w = j1f(x);
if (fabsf(z) >= fabsf(w))
b = (t*z/b);
else
@ -170,7 +170,7 @@ __ieee754_jnf(int n, float x)
}
float
__ieee754_ynf(int n, float x)
ynf(int n, float x)
{
int32_t i,hx,ix,ib;
int32_t sign;
@ -186,12 +186,12 @@ __ieee754_ynf(int n, float x)
n = -n;
sign = 1 - ((n&1)<<1);
}
if(n==0) return(__ieee754_y0f(x));
if(n==1) return(sign*__ieee754_y1f(x));
if(n==0) return(y0f(x));
if(n==1) return(sign*y1f(x));
if(ix==0x7f800000) return zero;
a = __ieee754_y0f(x);
b = __ieee754_y1f(x);
a = y0f(x);
b = y1f(x);
/* quit if b is -inf */
GET_FLOAT_WORD(ib,b);
for(i=1;i<n&&ib!=0xff800000;i++){

8
libm/upstream-freebsd/lib/msun/src/e_lgamma.c
View file

@ -15,10 +15,10 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_lgamma(x)
/* lgamma(x)
* Return the logarithm of the Gamma function of x.
*
* Method: call __ieee754_lgamma_r
* Method: call lgamma_r
*/
#include <float.h>
@ -29,9 +29,9 @@ __FBSDID("$FreeBSD$");
extern int signgam;
double
__ieee754_lgamma(double x)
lgamma(double x)
{
return __ieee754_lgamma_r(x,&signgam);
return lgamma_r(x,&signgam);
}
#if (LDBL_MANT_DIG == 53)

16
libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c
View file

@ -13,7 +13,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_lgamma_r(x, signgamp)
/* lgamma_r(x, signgamp)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
@ -199,7 +199,7 @@ sin_pi(double x)
double
__ieee754_lgamma_r(double x, int *signgamp)
lgamma_r(double x, int *signgamp)
{
double nadj,p,p1,p2,p3,q,r,t,w,y,z;
int32_t hx;
@ -217,7 +217,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
if(ix<0x3c700000) { /* |x|<2**-56, return -log(|x|) */
if((ix|lx)==0)
return one/vzero;
return -__ieee754_log(fabs(x));
return -log(fabs(x));
}
/* purge negative integers and start evaluation for other x < 0 */
@ -227,7 +227,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
return one/vzero;
t = sin_pi(x);
if(t==zero) return one/vzero; /* -integer */
nadj = __ieee754_log(pi/fabs(t*x));
nadj = log(pi/fabs(t*x));
if(t<zero) *signgamp = -1;
x = -x;
}
@ -237,7 +237,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
/* for x < 2.0 */
else if(ix<0x40000000) {
if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -__ieee754_log(x);
r = -log(x);
if(ix>=0x3FE76944) {y = one-x; i= 0;}
else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
else {y = x; i=2;}
@ -282,18 +282,18 @@ __ieee754_lgamma_r(double x, int *signgamp)
case 5: z *= (y+4); /* FALLTHRU */
case 4: z *= (y+3); /* FALLTHRU */
case 3: z *= (y+2); /* FALLTHRU */
r += __ieee754_log(z); break;
r += log(z); break;
}
/* 8.0 <= x < 2**56 */
} else if (ix < 0x43700000) {
t = __ieee754_log(x);
t = log(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**56 <= x <= inf */
r = x*(__ieee754_log(x)-one);
r = x*(log(x)-one);
if(hx<0) r = nadj - r;
return r;
}

8
libm/upstream-freebsd/lib/msun/src/e_lgammaf.c
View file

@ -16,10 +16,10 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_lgammaf(x)
/* lgammaf(x)
* Return the logarithm of the Gamma function of x.
*
* Method: call __ieee754_lgammaf_r
* Method: call lgammaf_r
*/
#include "math.h"
@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
extern int signgam;
float
__ieee754_lgammaf(float x)
lgammaf(float x)
{
return __ieee754_lgammaf_r(x,&signgam);
return lgammaf_r(x,&signgam);
}

14
libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c
View file

@ -120,7 +120,7 @@ sin_pif(float x)
float
__ieee754_lgammaf_r(float x, int *signgamp)
lgammaf_r(float x, int *signgamp)
{
float nadj,p,p1,p2,q,r,t,w,y,z;
int32_t hx;
@ -138,7 +138,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
if(ix<0x32000000) { /* |x|<2**-27, return -log(|x|) */
if(ix==0)
return one/vzero;
return -__ieee754_logf(fabsf(x));
return -logf(fabsf(x));
}
/* purge negative integers and start evaluation for other x < 0 */
@ -148,7 +148,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
return one/vzero;
t = sin_pif(x);
if(t==zero) return one/vzero; /* -integer */
nadj = __ieee754_logf(pi/fabsf(t*x));
nadj = logf(pi/fabsf(t*x));
if(t<zero) *signgamp = -1;
x = -x;
}
@ -158,7 +158,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
/* for x < 2.0 */
else if(ix<0x40000000) {
if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -__ieee754_logf(x);
r = -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;}
@ -198,18 +198,18 @@ __ieee754_lgammaf_r(float x, int *signgamp)
case 5: z *= (y+4); /* FALLTHRU */
case 4: z *= (y+3); /* FALLTHRU */
case 3: z *= (y+2); /* FALLTHRU */
r += __ieee754_logf(z); break;
r += logf(z); break;
}
/* 8.0 <= x < 2**27 */
} else if (ix < 0x4d000000) {
t = __ieee754_logf(x);
t = logf(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*w2);
r = (x-half)*(t-one)+w;
} else
/* 2**27 <= x <= inf */
r = x*(__ieee754_logf(x)-one);
r = x*(logf(x)-one);
if(hx<0) r = nadj - r;
return r;
}

4
libm/upstream-freebsd/lib/msun/src/e_log.c
View file

@ -14,7 +14,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_log(x)
/* log(x)
* Return the logrithm of x
*
* Method :
@ -86,7 +86,7 @@ static const double zero = 0.0;
static volatile double vzero = 0.0;
double
__ieee754_log(double x)
log(double x)
{
double hfsq,f,s,z,R,w,t1,t2,dk;
int32_t k,hx,i,j;

2
libm/upstream-freebsd/lib/msun/src/e_log10.c
View file

@ -39,7 +39,7 @@ static const double zero = 0.0;
static volatile double vzero = 0.0;
double
__ieee754_log10(double x)
log10(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
int32_t i,k,hx;

2
libm/upstream-freebsd/lib/msun/src/e_log10f.c
View file

@ -31,7 +31,7 @@ static const float zero = 0.0;
static volatile float vzero = 0.0;
float
__ieee754_log10f(float x)
log10f(float x)
{
float f,hfsq,hi,lo,r,y;
int32_t i,k,hx;

2
libm/upstream-freebsd/lib/msun/src/e_log2.c
View file

@ -39,7 +39,7 @@ static const double zero = 0.0;
static volatile double vzero = 0.0;
double
__ieee754_log2(double x)
log2(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
int32_t i,k,hx;

2
libm/upstream-freebsd/lib/msun/src/e_log2f.c
View file

@ -29,7 +29,7 @@ static const float zero = 0.0;
static volatile float vzero = 0.0;
float
__ieee754_log2f(float x)
log2f(float x)
{
float f,hfsq,hi,lo,r,y;
int32_t i,k,hx;

2
libm/upstream-freebsd/lib/msun/src/e_logf.c
View file

@ -33,7 +33,7 @@ static const float zero = 0.0;
static volatile float vzero = 0.0;
float
__ieee754_logf(float x)
logf(float x)
{
float hfsq,f,s,z,R,w,t1,t2,dk;
int32_t k,ix,i,j;

4
libm/upstream-freebsd/lib/msun/src/e_pow.c
View file

@ -12,7 +12,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_pow(x,y) return x**y
/* pow(x,y) return x**y
*
* n
* Method: Let x = 2 * (1+f)
@ -98,7 +98,7 @@ ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
double
__ieee754_pow(double x, double y)
pow(double x, double y)
{
double z,ax,z_h,z_l,p_h,p_l;
double y1,t1,t2,r,s,t,u,v,w;

4
libm/upstream-freebsd/lib/msun/src/e_powf.c
View file

@ -56,7 +56,7 @@ ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
float
__ieee754_powf(float x, float y)
powf(float x, float y)
{
float z,ax,z_h,z_l,p_h,p_l;
float y1,t1,t2,r,s,sn,t,u,v,w;
@ -108,7 +108,7 @@ __ieee754_powf(float x, float y)
if(hy==0x40000000) return x*x; /* y is 2 */
if(hy==0x3f000000) { /* y is 0.5 */
if(hx>=0) /* x >= +0 */
return __ieee754_sqrtf(x);
return sqrtf(x);
}
ax = fabsf(x);

6
libm/upstream-freebsd/lib/msun/src/e_remainder.c
View file

@ -14,7 +14,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_remainder(x,p)
/* remainder(x,p)
* Return :
* returns x REM p = x - [x/p]*p as if in infinite
* precise arithmetic, where [x/p] is the (infinite bit)
@ -32,7 +32,7 @@ static const double zero = 0.0;
double
__ieee754_remainder(double x, double p)
remainder(double x, double p)
{
int32_t hx,hp;
u_int32_t sx,lx,lp;
@ -52,7 +52,7 @@ __ieee754_remainder(double x, double p)
return nan_mix_op(x, p, *)/nan_mix_op(x, p, *);
if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */
if (hp<=0x7fdfffff) x = fmod(x,p+p); /* now x < 2p */
if (((hx-hp)|(lx-lp))==0) return zero*x;
x = fabs(x);
p = fabs(p);

4
libm/upstream-freebsd/lib/msun/src/e_remainderf.c
View file

@ -23,7 +23,7 @@ static const float zero = 0.0;
float
__ieee754_remainderf(float x, float p)
remainderf(float x, float p)
{
int32_t hx,hp;
u_int32_t sx;
@ -42,7 +42,7 @@ __ieee754_remainderf(float x, float p)
return nan_mix_op(x, p, *)/nan_mix_op(x, p, *);
if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */
if (hp<=0x7effffff) x = fmodf(x,p+p); /* now x < 2p */
if ((hx-hp)==0) return zero*x;
x = fabsf(x);
p = fabsf(p);

6
libm/upstream-freebsd/lib/msun/src/e_scalb.c
View file

@ -15,7 +15,7 @@
__FBSDID("$FreeBSD$");
/*
* __ieee754_scalb(x, fn) is provide for
* scalb(x, fn) is provide for
* passing various standard test suite. One
* should use scalbn() instead.
*/
@ -25,10 +25,10 @@ __FBSDID("$FreeBSD$");
#ifdef _SCALB_INT
double
__ieee754_scalb(double x, int fn)
scalb(double x, int fn)
#else
double
__ieee754_scalb(double x, double fn)
scalb(double x, double fn)
#endif
{
#ifdef _SCALB_INT

4
libm/upstream-freebsd/lib/msun/src/e_scalbf.c
View file

@ -21,10 +21,10 @@ __FBSDID("$FreeBSD$");
#ifdef _SCALB_INT
float
__ieee754_scalbf(float x, int fn)
scalbf(float x, int fn)
#else
float
__ieee754_scalbf(float x, float fn)
scalbf(float x, float fn)
#endif
{
#ifdef _SCALB_INT

6
libm/upstream-freebsd/lib/msun/src/e_sinh.c
View file

@ -14,7 +14,7 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_sinh(x)
/* sinh(x)
* Method :
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
* 1. Replace x by |x| (sinh(-x) = -sinh(x)).
@ -40,7 +40,7 @@ __FBSDID("$FreeBSD$");
static const double one = 1.0, shuge = 1.0e307;
double
__ieee754_sinh(double x)
sinh(double x)
{
double t,h;
int32_t ix,jx;
@ -64,7 +64,7 @@ __ieee754_sinh(double x)
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x));
if (ix < 0x40862E42) return h*exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
if (ix<=0x408633CE)

4
libm/upstream-freebsd/lib/msun/src/e_sinhf.c
View file

@ -22,7 +22,7 @@ __FBSDID("$FreeBSD$");
static const float one = 1.0, shuge = 1.0e37;
float
__ieee754_sinhf(float x)
sinhf(float x)
{
float t,h;
int32_t ix,jx;
@ -45,7 +45,7 @@ __ieee754_sinhf(float x)
}
/* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */
if (ix < 0x42b17217) return h*__ieee754_expf(fabsf(x));
if (ix < 0x42b17217) return h*expf(fabsf(x));
/* |x| in [logf(maxfloat), overflowthresold] */
if (ix<=0x42b2d4fc)

231
libm/upstream-freebsd/lib/msun/src/math_private.h
View file

@ -624,7 +624,7 @@ rnintf(__float_t x)
* The complications for extra precision are smaller for rnintl() since it
* can safely assume that the rounding precision has been increased from
* its default to FP_PE on x86. We don't exploit that here to get small
* optimizations from limiting the rangle to double. We just need it for
* optimizations from limiting the range to double. We just need it for
* the magic number to work with long doubles. ld128 callers should use
* rnint() instead of this if possible. ld80 callers should prefer
* rnintl() since for amd64 this avoids swapping the register set, while
@ -688,6 +688,59 @@ irintl(long double x)
}
#endif
/*
* The following are fast floor macros for 0 <= |x| < 0x1p(N-1), where
* N is the precision of the type of x. These macros are used in the
* half-cycle trignometric functions (e.g., sinpi(x)).
*/
#define FFLOORF(x, j0, ix) do { \
(j0) = (((ix) >> 23) & 0xff) - 0x7f; \
(ix) &= ~(0x007fffff >> (j0)); \
SET_FLOAT_WORD((x), (ix)); \
} while (0)
#define FFLOOR(x, j0, ix, lx) do { \
(j0) = (((ix) >> 20) & 0x7ff) - 0x3ff; \
if ((j0) < 20) { \
(ix) &= ~(0x000fffff >> (j0)); \
(lx) = 0; \
} else { \
(lx) &= ~((uint32_t)0xffffffff >> ((j0) - 20)); \
} \
INSERT_WORDS((x), (ix), (lx)); \
} while (0)
#define FFLOORL80(x, j0, ix, lx) do { \
j0 = ix - 0x3fff + 1; \
if ((j0) < 32) { \
(lx) = ((lx) >> 32) << 32; \
(lx) &= ~((((lx) << 32)-1) >> (j0)); \
} else { \
uint64_t _m; \
_m = (uint64_t)-1 >> (j0); \
if ((lx) & _m) (lx) &= ~_m; \
} \
INSERT_LDBL80_WORDS((x), (ix), (lx)); \
} while (0)
#define FFLOORL128(x, ai, ar) do { \
union IEEEl2bits u; \
uint64_t m; \
int e; \
u.e = (x); \
e = u.bits.exp - 16383; \
if (e < 48) { \
m = ((1llu << 49) - 1) >> (e + 1); \
u.bits.manh &= ~m; \
u.bits.manl = 0; \
} else { \
m = (uint64_t)-1 >> (e - 48); \
u.bits.manl &= ~m; \
} \
(ai) = u.e; \
(ar) = (x) - (ai); \
} while (0)
#ifdef DEBUG
#if defined(__amd64__) || defined(__i386__)
#define breakpoint() asm("int $3")
@ -698,191 +751,25 @@ irintl(long double x)
#endif
#endif
/* Write a pari script to test things externally. */
#ifdef DOPRINT
#include <stdio.h>
#ifndef DOPRINT_SWIZZLE
#define DOPRINT_SWIZZLE 0
#endif
#ifdef DOPRINT_LD80
#define DOPRINT_START(xp) do { \
uint64_t __lx; \
uint16_t __hx; \
\
/* Hack to give more-problematic args. */ \
EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \
__lx ^= DOPRINT_SWIZZLE; \
INSERT_LDBL80_WORDS(*xp, __hx, __lx); \
printf("x = %.21Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#elif defined(DOPRINT_D64)
#define DOPRINT_START(xp) do { \
uint32_t __hx, __lx; \
\
EXTRACT_WORDS(__hx, __lx, *xp); \
__lx ^= DOPRINT_SWIZZLE; \
INSERT_WORDS(*xp, __hx, __lx); \
printf("x = %.21Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#elif defined(DOPRINT_F32)
#define DOPRINT_START(xp) do { \
uint32_t __hx; \
\
GET_FLOAT_WORD(__hx, *xp); \
__hx ^= DOPRINT_SWIZZLE; \
SET_FLOAT_WORD(*xp, __hx); \
printf("x = %.21Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
#ifndef DOPRINT_SWIZZLE_HIGH
#define DOPRINT_SWIZZLE_HIGH 0
#endif
#define DOPRINT_START(xp) do { \
uint64_t __lx, __llx; \
uint16_t __hx; \
\
EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \
__llx ^= DOPRINT_SWIZZLE; \
__lx ^= DOPRINT_SWIZZLE_HIGH; \
INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \
printf("x = %.36Lg; ", (long double)*xp); \
} while (0)
#define DOPRINT_END1(v) \
printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define DOPRINT_END2(hi, lo) \
printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \
(long double)(hi), (long double)(lo))
#endif /* DOPRINT_LD80 */
#else /* !DOPRINT */
#define DOPRINT_START(xp)
#define DOPRINT_END1(v)
#define DOPRINT_END2(hi, lo)
#endif /* DOPRINT */
#define RETURNP(x) do { \
DOPRINT_END1(x); \
RETURNF(x); \
} while (0)
#define RETURNPI(x) do { \
DOPRINT_END1(x); \
RETURNI(x); \
} while (0)
#define RETURN2P(x, y) do { \
DOPRINT_END2((x), (y)); \
RETURNF((x) + (y)); \
} while (0)
#define RETURN2PI(x, y) do { \
DOPRINT_END2((x), (y)); \
RETURNI((x) + (y)); \
} while (0)
#ifdef STRUCT_RETURN
#define RETURNSP(rp) do { \
if (!(rp)->lo_set) \
RETURNP((rp)->hi); \
RETURN2P((rp)->hi, (rp)->lo); \
RETURNF((rp)->hi); \
RETURNF((rp)->hi + (rp)->lo); \
} while (0)
#define RETURNSPI(rp) do { \
if (!(rp)->lo_set) \
RETURNPI((rp)->hi); \
RETURN2PI((rp)->hi, (rp)->lo); \
RETURNI((rp)->hi); \
RETURNI((rp)->hi + (rp)->lo); \
} while (0)
#endif
#define SUM2P(x, y) ({ \
const __typeof (x) __x = (x); \
const __typeof (y) __y = (y); \
\
DOPRINT_END2(__x, __y); \
__x + __y; \
})
/*
* ieee style elementary functions
*
* We rename functions here to improve other sources' diffability
* against fdlibm.
*/
#define __ieee754_sqrt sqrt
#define __ieee754_acos acos
#define __ieee754_acosh acosh
#define __ieee754_log log
#define __ieee754_log2 log2
#define __ieee754_atanh atanh
#define __ieee754_asin asin
#define __ieee754_atan2 atan2
#define __ieee754_exp exp
#define __ieee754_cosh cosh
#define __ieee754_fmod fmod
#define __ieee754_pow pow
#define __ieee754_lgamma lgamma
#define __ieee754_gamma gamma
#define __ieee754_lgamma_r lgamma_r
#define __ieee754_gamma_r gamma_r
#define __ieee754_log10 log10
#define __ieee754_sinh sinh
#define __ieee754_hypot hypot
#define __ieee754_j0 j0
#define __ieee754_j1 j1
#define __ieee754_y0 y0
#define __ieee754_y1 y1
#define __ieee754_jn jn
#define __ieee754_yn yn
#define __ieee754_remainder remainder
#define __ieee754_scalb scalb
#define __ieee754_sqrtf sqrtf
#define __ieee754_acosf acosf
#define __ieee754_acoshf acoshf
#define __ieee754_logf logf
#define __ieee754_atanhf atanhf
#define __ieee754_asinf asinf
#define __ieee754_atan2f atan2f
#define __ieee754_expf expf
#define __ieee754_coshf coshf
#define __ieee754_fmodf fmodf
#define __ieee754_powf powf
#define __ieee754_lgammaf lgammaf
#define __ieee754_gammaf gammaf
#define __ieee754_lgammaf_r lgammaf_r
#define __ieee754_gammaf_r gammaf_r
#define __ieee754_log10f log10f
#define __ieee754_log2f log2f
#define __ieee754_sinhf sinhf
#define __ieee754_hypotf hypotf
#define __ieee754_j0f j0f
#define __ieee754_j1f j1f
#define __ieee754_y0f y0f
#define __ieee754_y1f y1f
#define __ieee754_jnf jnf
#define __ieee754_ynf ynf
#define __ieee754_remainderf remainderf
#define __ieee754_scalbf scalbf
/* fdlibm kernel function */
int __kernel_rem_pio2(double*,double*,int,int,int);

6
libm/upstream-freebsd/lib/msun/src/s_asinh.c
View file

@ -46,13 +46,13 @@ asinh(double x)
if(huge+x>one) return x; /* return x inexact except 0 */
}
if(ix>0x41b00000) { /* |x| > 2**28 */
w = __ieee754_log(fabs(x))+ln2;
w = log(fabs(x))+ln2;
} else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabs(x);
w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
w = log(2.0*t+one/(sqrt(x*x+one)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
w =log1p(fabs(x)+t/(one+sqrt(one+t)));
}
if(hx>0) return w; else return -w;
}

6
libm/upstream-freebsd/lib/msun/src/s_asinhf.c
View file

@ -36,13 +36,13 @@ asinhf(float x)
if(huge+x>one) return x; /* return x inexact except 0 */
}
if(ix>0x4d800000) { /* |x| > 2**28 */
w = __ieee754_logf(fabsf(x))+ln2;
w = logf(fabsf(x))+ln2;
} else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabsf(x);
w = __ieee754_logf((float)2.0*t+one/(__ieee754_sqrtf(x*x+one)+t));
w = logf((float)2.0*t+one/(sqrtf(x*x+one)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1pf(fabsf(x)+t/(one+__ieee754_sqrtf(one+t)));
w =log1pf(fabsf(x)+t/(one+sqrtf(one+t)));
}
if(hx>0) return w; else return -w;
}

20
libm/upstream-freebsd/lib/msun/src/s_cospi.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -104,20 +104,10 @@ cospi(double x)
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
}
INSERT_WORDS(x, ix, lx);
FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
ax -= x;
EXTRACT_WORDS(ix, lx, ax);
if (ix < 0x3fe00000) { /* |x| < 0.5 */
if (ix < 0x3fd00000) /* |x| < 0.25 */
c = ix == 0 ? 1 : __kernel_cospi(ax);
@ -143,9 +133,11 @@ cospi(double x)
return (vzero / vzero);
/*
* |x| >= 0x1p52 is always an even integer, so return 1.
* For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even
* or odd integer to return +1 or -1.
* For |x| >= 0x1p53, it is always an even integer, so return 1.
*/
return (1);
return (ix < 0x43400000 ? ((lx & 1) ? -1 : 1) : 1);
}
#if LDBL_MANT_DIG == 53

4
libm/upstream-freebsd/lib/msun/src/s_erf.c
View file

@ -238,7 +238,7 @@ erf(double x)
}
z = x;
SET_LOW_WORD(z,0);
r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if(hx>=0) return one-r/x; else return r/x-one;
}
@ -297,7 +297,7 @@ erfc(double x)
}
z = x;
SET_LOW_WORD(z,0);
r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if(hx>0) return r/x; else return two-r/x;
} else {
if(hx>0) return tiny*tiny; else return two-tiny;

2
libm/upstream-freebsd/lib/msun/src/s_significand.c
View file

@ -25,5 +25,5 @@ __FBSDID("$FreeBSD$");
double
significand(double x)
{
return __ieee754_scalb(x,(double) -ilogb(x));
return scalb(x,(double) -ilogb(x));
}

2
libm/upstream-freebsd/lib/msun/src/s_significandf.c
View file

@ -22,5 +22,5 @@ __FBSDID("$FreeBSD$");
float
significandf(float x)
{
return __ieee754_scalbf(x,(float) -ilogbf(x));
return scalbf(x,(float) -ilogbf(x));
}

13
libm/upstream-freebsd/lib/msun/src/s_sinpi.c
View file

@ -1,5 +1,5 @@
/*-
* Copyright (c) 2017 Steven G. Kargl
* Copyright (c) 2017, 2023 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@ -118,16 +118,7 @@ sinpi(double x)
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
}
INSERT_WORDS(x, ix, lx);
FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
ax -= x;
EXTRACT_WORDS(ix, lx, ax);