Clean up libm use of the __ieee754_ prefix

This removes the __ieee754_ prefix from a number of the math functions.
msun/src/math_private.h contains the statement that

  /*
   * ieee style elementary functions
   *
   * We rename functions here to improve other sources' diffability
   * against fdlibm.
   */
   #define        __ieee754_sqrt  sqrt
   ...

Here, fdlibm refers to https://netlib.org/fdlibm. It is seen from
https://netlib.org/fdlibm/readme that this prefix was used to
differentiate between different standards:

   Wrapper functions will twist the result of the ieee754
   function to comply to the standard specified by the value
   of _LIB_VERSION
      if _LIB_VERSION = _IEEE_, return the ieee754 result;
      if _LIB_VERSION = _SVID_, return SVID result;
      if _LIB_VERSION = _XOPEN_, return XOPEN result;
      if _LIB_VERSION = _POSIX_, return POSIX/ANSI result.
   (These are macros, see fdlibm.h for their definition.)

AFAICT, FreeBSD has never supported these wrappers. In addition, as C99,
principally the long double, functions were added to libm, this
convention was not maintained. Given that only 148 of 324 files under
lib/msun contain a "Copyright (C) 1993 by Sun Microsystems" statement,
the removal of the __ieee754_ prefix provides consistency across all
source files.

The last time someone compared lib/msun to fdlibm appears to be

  commit 3f70824172
  Author: David Schultz <das@FreeBSD.org>
  Date:   Fri Feb 4 18:26:06 2005 +0000

  Reduce diffs against vendor source (Sun fdlibm 5.3).

The most recent fdlibm RCS string that appears in a Sun Microsystem
copyrighted file is date "95/01/18". With Oracle Corporation's
acquisition of Sun Microsystems in 2009, it is unlikely that fdlibm will
ever be updated. A search for fdlibm at https://opensource.oracle.com/
yields no hits.

Finally, OpenBSD removed the use of this prefix over 21 years ago. pSee
revision 1.6 of OpenBSD's math_private.h.

Note: this does not drop the __ieee754_ prefix from the trigonometric
argument reduction functions, e.g., __ieee754_rem_pio2. These functions
are internal to the libm and exported through Symbol.map; and thus,
reserved for the implementation.

PR:		272783
MFC after:	1 week
This commit is contained in:
Steve Kargl 2023-08-03 21:51:17 +02:00 committed by Dimitry Andric
parent 56b9145cd0
commit 99843eb899
54 changed files with 162 additions and 223 deletions

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;

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;

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));

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));
}
}

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;

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;

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;

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;

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;

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;

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)

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)

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;

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;

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;

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;

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);
}

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);
}

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);
}

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);
}

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;

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);

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

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

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

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

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++){

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++){

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)

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;
}

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);
}

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;
}

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;

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;

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;

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;

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;

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;

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;

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);

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);

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);

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

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

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)

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)

View file

@ -21,12 +21,12 @@ __FBSDID("$FreeBSD$");
#ifdef USE_BUILTIN_SQRT
double
__ieee754_sqrt(double x)
sqrt(double x)
{
return (__builtin_sqrt(x));
}
#else
/* __ieee754_sqrt(x)
/* sqrt(x)
* Return correctly rounded sqrt.
* ------------------------------------------
* | Use the hardware sqrt if you have one |
@ -99,7 +99,7 @@ __ieee754_sqrt(double x)
static const double one = 1.0, tiny=1.0e-300;
double
__ieee754_sqrt(double x)
sqrt(double x)
{
double z;
int32_t sign = (int)0x80000000;

View file

@ -22,7 +22,7 @@ static char rcsid[] = "$FreeBSD$";
#ifdef USE_BUILTIN_SQRTF
float
__ieee754_sqrtf(float x)
sqrtf(float x)
{
return (__builtin_sqrtf(x));
}
@ -30,7 +30,7 @@ __ieee754_sqrtf(float x)
static const float one = 1.0, tiny=1.0e-30;
float
__ieee754_sqrtf(float x)
sqrtf(float x)
{
float z;
int32_t sign = (int)0x80000000;

View file

@ -770,67 +770,6 @@ irintl(long double x)
__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);

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;
}

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;
}

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;

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));
}

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));
}