Fcalva/OpenLibm
1
0
Fork 0
mirror of https://git.planet-casio.com/Lephenixnoir/OpenLibm.git synced 2024-12-28 20:43:41 +01:00
Code Issues Projects Releases Packages Wiki Activity Actions

Import long double versions from OpenBSD.

Browse source
This commit is contained in:
Viral B. Shah 2014-12-04 23:11:16 +05:30
parent 691b989ba3
commit 52c901a68c
53 changed files with 9450 additions and 3 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

58
ld128/e_acoshl.c Normal file
View file

@ -0,0 +1,58 @@
/* @(#)e_acosh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* acoshl(x)
* Method :
* Based on
* acoshl(x) = logl [ x + sqrtl(x*x-1) ]
* we have
* acoshl(x) := logl(x)+ln2, if x is large; else
* acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
* acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
*
* Special cases:
* acoshl(x) is NaN with signal if x<1.
* acoshl(NaN) is NaN without signal.
*/
#include <math.h>
#include "math_private.h"
static const long double
one = 1.0,
ln2 = 0.6931471805599453094172321214581766L;
long double
acoshl(long double x)
{
long double t;
u_int64_t lx;
int64_t hx;
GET_LDOUBLE_WORDS64(hx,lx,x);
if(hx<0x3fff000000000000LL) { /* x < 1 */
return (x-x)/(x-x);
} else if(hx >=0x4035000000000000LL) { /* x > 2**54 */
if(hx >=0x7fff000000000000LL) { /* x is inf of NaN */
return x+x;
} else
return logl(x)+ln2; /* acoshl(huge)=logl(2x) */
} else if(((hx-0x3fff000000000000LL)|lx)==0) {
return 0.0L; /* acosh(1) = 0 */
} else if (hx > 0x4000000000000000LL) { /* 2**28 > x > 2 */
t=x*x;
return logl(2.0L*x-one/(x+sqrtl(t-one)));
} else { /* 1<x<2 */
t = x-one;
return log1pl(t+sqrtl(2.0L*t+t*t));
}
}

65
ld128/e_atanhl.c Normal file
View file

@ -0,0 +1,65 @@
/* @(#)e_atanh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* atanhl(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
* 1 2x x
* atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x
*
* For x<0.5
* atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
*
* Special cases:
* atanhl(x) is NaN if |x| > 1 with signal;
* atanhl(NaN) is that NaN with no signal;
* atanhl(+-1) is +-INF with signal.
*
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0L, huge = 1e4900L;
static const long double zero = 0.0L;
long double
atanhl(long double x)
{
long double t;
u_int32_t jx, ix;
ieee_quad_shape_type u;
u.value = x;
jx = u.parts32.mswhi;
ix = jx & 0x7fffffff;
u.parts32.mswhi = ix;
if (ix >= 0x3fff0000) /* |x| >= 1.0 or infinity or NaN */
{
if (u.value == one)
return x/zero;
else
return (x-x)/(x-x);
}
if(ix<0x3fc60000 && (huge+x)>zero) return x; /* x < 2^-57 */
if(ix<0x3ffe0000) { /* x < 0.5 */
t = u.value+u.value;
t = 0.5*log1pl(t+t*u.value/(one-u.value));
} else
t = 0.5*log1pl((u.value+u.value)/(one-u.value));
if(jx & 0x80000000) return -t; else return t;
}

105
ld128/e_coshl.c Normal file
View file

@ -0,0 +1,105 @@
/*
* ====================================================
* 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) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* coshl(x)
* Method :
* mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
* 1. Replace x by |x| (coshl(x) = coshl(-x)).
* 2.
* [ exp(x) - 1 ]^2
* 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
* 2*exp(x)
*
* exp(x) + 1/exp(x)
* ln2/2 <= x <= 22 : coshl(x) := -------------------
* 2
* 22 <= x <= lnovft : coshl(x) := expl(x)/2
* lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2)
* ln2ovft < x : coshl(x) := huge*huge (overflow)
*
* Special cases:
* coshl(x) is |x| if x is +INF, -INF, or NaN.
* only coshl(0)=1 is exact for finite x.
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0, half = 0.5, huge = 1.0e4900L,
ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
long double
coshl(long double x)
{
long double t, w;
int32_t ex;
ieee_quad_shape_type u;
u.value = x;
ex = u.parts32.mswhi & 0x7fffffff;
/* Absolute value of x. */
u.parts32.mswhi = ex;
/* x is INF or NaN */
if (ex >= 0x7fff0000)
return x * x;
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
{
t = expm1l (u.value);
w = one + t;
if (ex < 0x3fb80000) /* |x| < 2^-116 */
return w; /* cosh(tiny) = 1 */
return one + (t * t) / (w + w);
}
/* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
if (ex < 0x40044000)
{
t = expl (u.value);
return half * t + half / t;
}
/* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
if (ex <= 0x400c62e3) /* 11356.375 */
return half * expl (u.value);
/* |x| in [log(maxdouble), overflowthresold] */
if (u.value <= ovf_thresh)
{
w = expl (half * u.value);
t = half * w;
return t * w;
}
/* |x| > overflowthresold, cosh(x) overflow */
return huge * huge;
}

145
ld128/e_expl.c Normal file
View file

@ -0,0 +1,145 @@
/* $OpenBSD: e_expl.c,v 1.3 2013/11/12 20:35:18 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* expl.c
*
* Exponential function, 128-bit long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, expl();
*
* y = expl( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* A Pade' form of degree 2/3 is used to approximate exp(f) - 1
* in the basic range [-0.5 ln 2, 0.5 ln 2].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +-MAXLOG 100,000 2.6e-34 8.6e-35
*
*
* Error amplification in the exponential function can be
* a serious matter. The error propagation involves
* exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
* which shows that a 1 lsb error in representing X produces
* a relative error of X times 1 lsb in the function.
* While the routine gives an accurate result for arguments
* that are exactly represented by a long double precision
* computer number, the result contains amplified roundoff
* error for large arguments not exactly represented.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* exp underflow x < MINLOG 0.0
* exp overflow x > MAXLOG MAXNUM
*
*/
/* Exponential function */
#include <float.h>
#include <math.h>
#include "math_private.h"
/* Pade' coefficients for exp(x) - 1
Theoretical peak relative error = 2.2e-37,
relative peak error spread = 9.2e-38
*/
static long double P[5] = {
3.279723985560247033712687707263393506266E-10L,
6.141506007208645008909088812338454698548E-7L,
2.708775201978218837374512615596512792224E-4L,
3.508710990737834361215404761139478627390E-2L,
9.999999999999999999999999999999999998502E-1L
};
static long double Q[6] = {
2.980756652081995192255342779918052538681E-12L,
1.771372078166251484503904874657985291164E-8L,
1.504792651814944826817779302637284053660E-5L,
3.611828913847589925056132680618007270344E-3L,
2.368408864814233538909747618894558968880E-1L,
2.000000000000000000000000000000000000150E0L
};
/* C1 + C2 = ln 2 */
static const long double C1 = -6.93145751953125E-1L;
static const long double C2 = -1.428606820309417232121458176568075500134E-6L;
static const long double LOG2EL = 1.442695040888963407359924681001892137426646L;
static const long double MAXLOGL = 1.1356523406294143949491931077970764891253E4L;
static const long double MINLOGL = -1.143276959615573793352782661133116431383730e4L;
static const long double huge = 0x1p10000L;
#if 0 /* XXX Prevent gcc from erroneously constant folding this. */
static const long double twom10000 = 0x1p-10000L;
#else
static volatile long double twom10000 = 0x1p-10000L;
#endif
long double
expl(long double x)
{
long double px, xx;
int n;
if( x > MAXLOGL)
return (huge*huge); /* overflow */
if( x < MINLOGL )
return (twom10000*twom10000); /* underflow */
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */
n = px;
x += px * C1;
x += px * C2;
/* rational approximation for exponential
* of the fractional part:
* e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*/
xx = x * x;
px = x * __polevll( xx, P, 4 );
xx = __polevll( xx, Q, 5 );
x = px/( xx - px );
x = 1.0L + x + x;
x = ldexpl( x, n );
return(x);
}

129
ld128/e_fmodl.c Normal file
View file

@ -0,0 +1,129 @@
/* @(#)e_fmod.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* fmodl(x,y)
* Return x mod y in exact arithmetic
* Method: shift and subtract
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0, Zero[] = {0.0, -0.0,};
long double
fmodl(long double x, long double y)
{
int64_t n,hx,hy,hz,ix,iy,sx,i;
u_int64_t lx,ly,lz;
GET_LDOUBLE_WORDS64(hx,lx,x);
GET_LDOUBLE_WORDS64(hy,ly,y);
sx = hx&0x8000000000000000ULL; /* sign of x */
hx ^=sx; /* |x| */
hy &= 0x7fffffffffffffffLL; /* |y| */
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>63))>0x7fff000000000000LL)) /* or y is NaN */
return (x*y)/(x*y);
if(hx<=hy) {
if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
if(lx==ly)
return Zero[(u_int64_t)sx>>63]; /* |x|=|y| return x*0*/
}
/* determine ix = ilogb(x) */
if(hx<0x0001000000000000LL) { /* subnormal x */
if(hx==0) {
for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
} else {
for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
}
} else ix = (hx>>48)-0x3fff;
/* determine iy = ilogb(y) */
if(hy<0x0001000000000000LL) { /* subnormal y */
if(hy==0) {
for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
} else {
for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
}
} else iy = (hy>>48)-0x3fff;
/* set up {hx,lx}, {hy,ly} and align y to x */
if(ix >= -16382)
hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
else { /* subnormal x, shift x to normal */
n = -16382-ix;
if(n<=63) {
hx = (hx<<n)|(lx>>(64-n));
lx <<= n;
} else {
hx = lx<<(n-64);
lx = 0;
}
}
if(iy >= -16382)
hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
else { /* subnormal y, shift y to normal */
n = -16382-iy;
if(n<=63) {
hy = (hy<<n)|(ly>>(64-n));
ly <<= n;
} else {
hy = ly<<(n-64);
ly = 0;
}
}
/* fix point fmod */
n = ix - iy;
while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
else {
if((hz|lz)==0) /* return sign(x)*0 */
return Zero[(u_int64_t)sx>>63];
hx = hz+hz+(lz>>63); lx = lz+lz;
}
}
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;}
/* convert back to floating value and restore the sign */
if((hx|lx)==0) /* return sign(x)*0 */
return Zero[(u_int64_t)sx>>63];
while(hx<0x0001000000000000LL) { /* normalize x */
hx = hx+hx+(lx>>63); lx = lx+lx;
iy -= 1;
}
if(iy>= -16382) { /* normalize output */
hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
SET_LDOUBLE_WORDS64(x,hx|sx,lx);
} else { /* subnormal output */
n = -16382 - iy;
if(n<=48) {
lx = (lx>>n)|((u_int64_t)hx<<(64-n));
hx >>= n;
} else if (n<=63) {
lx = (hx<<(64-n))|(lx>>n); hx = sx;
} else {
lx = hx>>(n-64); hx = sx;
}
SET_LDOUBLE_WORDS64(x,hx|sx,lx);
x *= one; /* create necessary signal */
}
return x; /* exact output */
}

122
ld128/e_hypotl.c Normal file
View file

@ -0,0 +1,122 @@
/* @(#)e_hypot.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* hypotl(x,y)
*
* Method :
* If (assume round-to-nearest) z=x*x+y*y
* has error less than sqrtl(2)/2 ulp, than
* sqrtl(z) has error less than 1 ulp (exercise).
*
* So, compute sqrtl(x*x+y*y) with some care as
* follows to get the error below 1 ulp:
*
* Assume x>y>0;
* (if possible, set rounding to round-to-nearest)
* 1. if x > 2y use
* x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
* where x1 = x with lower 64 bits cleared, x2 = x-x1; else
* 2. if x <= 2y use
* t1*yy1+((x-y)*(x-y)+(t1*y2+t2*y))
* where t1 = 2x with lower 64 bits cleared, t2 = 2x-t1,
* yy1= y with lower 64 bits chopped, y2 = y-yy1.
*
* NOTE: scaling may be necessary if some argument is too
* large or too tiny
*
* Special cases:
* hypotl(x,y) is INF if x or y is +INF or -INF; else
* hypotl(x,y) is NAN if x or y is NAN.
*
* Accuracy:
* hypotl(x,y) returns sqrtl(x^2+y^2) with error less
* than 1 ulps (units in the last place)
*/
#include <math.h>
#include "math_private.h"
long double
hypotl(long double x, long double y)
{
long double a,b,t1,t2,yy1,y2,w;
int64_t j,k,ha,hb;
GET_LDOUBLE_MSW64(ha,x);
ha &= 0x7fffffffffffffffLL;
GET_LDOUBLE_MSW64(hb,y);
hb &= 0x7fffffffffffffffLL;
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
SET_LDOUBLE_MSW64(a,ha); /* a <- |a| */
SET_LDOUBLE_MSW64(b,hb); /* b <- |b| */
if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */
k=0;
if(ha > 0x5f3f000000000000LL) { /* a>2**8000 */
if(ha >= 0x7fff000000000000LL) { /* Inf or NaN */
u_int64_t low;
w = a+b; /* for sNaN */
GET_LDOUBLE_LSW64(low,a);
if(((ha&0xffffffffffffLL)|low)==0) w = a;
GET_LDOUBLE_LSW64(low,b);
if(((hb^0x7fff000000000000LL)|low)==0) w = b;
return w;
}
/* scale a and b by 2**-9600 */
ha -= 0x2580000000000000LL;
hb -= 0x2580000000000000LL; k += 9600;
SET_LDOUBLE_MSW64(a,ha);
SET_LDOUBLE_MSW64(b,hb);
}
if(hb < 0x20bf000000000000LL) { /* b < 2**-8000 */
if(hb <= 0x0000ffffffffffffLL) { /* subnormal b or 0 */
u_int64_t low;
GET_LDOUBLE_LSW64(low,b);
if((hb|low)==0) return a;
t1=0;
SET_LDOUBLE_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */
b *= t1;
a *= t1;
k -= 16382;
} else { /* scale a and b by 2^9600 */
ha += 0x2580000000000000LL; /* a *= 2^9600 */
hb += 0x2580000000000000LL; /* b *= 2^9600 */
k -= 9600;
SET_LDOUBLE_MSW64(a,ha);
SET_LDOUBLE_MSW64(b,hb);
}
}
/* medium size a and b */
w = a-b;
if (w>b) {
t1 = 0;
SET_LDOUBLE_MSW64(t1,ha);
t2 = a-t1;
w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
yy1 = 0;
SET_LDOUBLE_MSW64(yy1,hb);
y2 = b - yy1;
t1 = 0;
SET_LDOUBLE_MSW64(t1,ha+0x0001000000000000LL);
t2 = a - t1;
w = sqrtl(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
u_int64_t high;
t1 = 1.0L;
GET_LDOUBLE_MSW64(high,t1);
SET_LDOUBLE_MSW64(t1,high+(k<<48));
return t1*w;
} else return w;
}

1038
ld128/e_lgammal.c Normal file
View file

File diff suppressed because it is too large Load diff

255
ld128/e_log10l.c Normal file
View file

@ -0,0 +1,255 @@
/* $OpenBSD: e_log10l.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* log10l.c
*
* Common logarithm, 128-bit long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, log10l();
*
* y = log10l( x );
*
*
*
* DESCRIPTION:
*
* Returns the base 10 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z^3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 2.3e-34 4.9e-35
* IEEE exp(+-10000) 30000 1.0e-34 4.1e-35
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
*/
#include <math.h>
#include "math_private.h"
/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 5.3e-37,
* relative peak error spread = 2.3e-14
*/
static const long double P[13] =
{
1.313572404063446165910279910527789794488E4L,
7.771154681358524243729929227226708890930E4L,
2.014652742082537582487669938141683759923E5L,
3.007007295140399532324943111654767187848E5L,
2.854829159639697837788887080758954924001E5L,
1.797628303815655343403735250238293741397E5L,
7.594356839258970405033155585486712125861E4L,
2.128857716871515081352991964243375186031E4L,
3.824952356185897735160588078446136783779E3L,
4.114517881637811823002128927449878962058E2L,
2.321125933898420063925789532045674660756E1L,
4.998469661968096229986658302195402690910E-1L,
1.538612243596254322971797716843006400388E-6L
};
static const long double Q[12] =
{
3.940717212190338497730839731583397586124E4L,
2.626900195321832660448791748036714883242E5L,
7.777690340007566932935753241556479363645E5L,
1.347518538384329112529391120390701166528E6L,
1.514882452993549494932585972882995548426E6L,
1.158019977462989115839826904108208787040E6L,
6.132189329546557743179177159925690841200E5L,
2.248234257620569139969141618556349415120E5L,
5.605842085972455027590989944010492125825E4L,
9.147150349299596453976674231612674085381E3L,
9.104928120962988414618126155557301584078E2L,
4.839208193348159620282142911143429644326E1L
/* 1.000000000000000000000000000000000000000E0L, */
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 1.1e-35,
* relative peak error spread 1.1e-9
*/
static const long double R[6] =
{
1.418134209872192732479751274970992665513E5L,
-8.977257995689735303686582344659576526998E4L,
2.048819892795278657810231591630928516206E4L,
-2.024301798136027039250415126250455056397E3L,
8.057002716646055371965756206836056074715E1L,
-8.828896441624934385266096344596648080902E-1L
};
static const long double S[6] =
{
1.701761051846631278975701529965589676574E6L,
-1.332535117259762928288745111081235577029E6L,
4.001557694070773974936904547424676279307E5L,
-5.748542087379434595104154610899551484314E4L,
3.998526750980007367835804959888064681098E3L,
-1.186359407982897997337150403816839480438E2L
/* 1.000000000000000000000000000000000000000E0L, */
};
static const long double
/* log10(2) */
L102A = 0.3125L,
L102B = -1.14700043360188047862611052755069732318101185E-2L,
/* log10(e) */
L10EA = 0.5L,
L10EB = -6.570551809674817234887108108339491770560299E-2L,
/* sqrt(2)/2 */
SQRTH = 7.071067811865475244008443621048490392848359E-1L;
/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
static long double
neval (long double x, const long double *p, int n)
{
long double y;
p += n;
y = *p--;
do
{
y = y * x + *p--;
}
while (--n > 0);
return y;
}
/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
static long double
deval (long double x, const long double *p, int n)
{
long double y;
p += n;
y = x + *p--;
do
{
y = y * x + *p--;
}
while (--n > 0);
return y;
}
long double
log10l(long double x)
{
long double z;
long double y;
int e;
int64_t hx, lx;
/* Test for domain */
GET_LDOUBLE_WORDS64 (hx, lx, x);
if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
return (-1.0L / (x - x));
if (hx < 0)
return (x - x) / (x - x);
if (hx >= 0x7fff000000000000LL)
return (x + x);
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl (x, &e);
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if ((e > 2) || (e < -2))
{
if (x < SQRTH)
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x * x;
y = x * (z * neval (z, R, 5) / deval (z, S, 5));
goto done;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH)
{
e -= 1;
x = 2.0 * x - 1.0L; /* 2x - 1 */
}
else
{
x = x - 1.0L;
}
z = x * x;
y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
y = y - 0.5 * z;
done:
/* Multiply log of fraction by log10(e)
* and base 2 exponent by log10(2).
*/
z = y * L10EB;
z += x * L10EB;
z += e * L102B;
z += y * L10EA;
z += x * L10EA;
z += e * L102A;
return (z);
}

248
ld128/e_log2l.c Normal file
View file

@ -0,0 +1,248 @@
/* $OpenBSD: e_log2l.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* log2l.c
* Base 2 logarithm, 128-bit long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, log2l();
*
* y = log2l( x );
*
*
*
* DESCRIPTION:
*
* Returns the base 2 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the (natural)
* logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z^3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 100,000 2.6e-34 4.9e-35
* IEEE exp(+-10000) 100,000 9.6e-35 4.0e-35
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
*/
#include <math.h>
#include "math_private.h"
/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 5.3e-37,
* relative peak error spread = 2.3e-14
*/
static const long double P[13] =
{
1.313572404063446165910279910527789794488E4L,
7.771154681358524243729929227226708890930E4L,
2.014652742082537582487669938141683759923E5L,
3.007007295140399532324943111654767187848E5L,
2.854829159639697837788887080758954924001E5L,
1.797628303815655343403735250238293741397E5L,
7.594356839258970405033155585486712125861E4L,
2.128857716871515081352991964243375186031E4L,
3.824952356185897735160588078446136783779E3L,
4.114517881637811823002128927449878962058E2L,
2.321125933898420063925789532045674660756E1L,
4.998469661968096229986658302195402690910E-1L,
1.538612243596254322971797716843006400388E-6L
};
static const long double Q[12] =
{
3.940717212190338497730839731583397586124E4L,
2.626900195321832660448791748036714883242E5L,
7.777690340007566932935753241556479363645E5L,
1.347518538384329112529391120390701166528E6L,
1.514882452993549494932585972882995548426E6L,
1.158019977462989115839826904108208787040E6L,
6.132189329546557743179177159925690841200E5L,
2.248234257620569139969141618556349415120E5L,
5.605842085972455027590989944010492125825E4L,
9.147150349299596453976674231612674085381E3L,
9.104928120962988414618126155557301584078E2L,
4.839208193348159620282142911143429644326E1L
/* 1.000000000000000000000000000000000000000E0L, */
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 1.1e-35,
* relative peak error spread 1.1e-9
*/
static const long double R[6] =
{
1.418134209872192732479751274970992665513E5L,
-8.977257995689735303686582344659576526998E4L,
2.048819892795278657810231591630928516206E4L,
-2.024301798136027039250415126250455056397E3L,
8.057002716646055371965756206836056074715E1L,
-8.828896441624934385266096344596648080902E-1L
};
static const long double S[6] =
{
1.701761051846631278975701529965589676574E6L,
-1.332535117259762928288745111081235577029E6L,
4.001557694070773974936904547424676279307E5L,
-5.748542087379434595104154610899551484314E4L,
3.998526750980007367835804959888064681098E3L,
-1.186359407982897997337150403816839480438E2L
/* 1.000000000000000000000000000000000000000E0L, */
};
static const long double
/* log2(e) - 1 */
LOG2EA = 4.4269504088896340735992468100189213742664595E-1L,
/* sqrt(2)/2 */
SQRTH = 7.071067811865475244008443621048490392848359E-1L;
/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
static long double
neval (long double x, const long double *p, int n)
{
long double y;
p += n;
y = *p--;
do
{
y = y * x + *p--;
}
while (--n > 0);
return y;
}
/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
static long double
deval (long double x, const long double *p, int n)
{
long double y;
p += n;
y = x + *p--;
do
{
y = y * x + *p--;
}
while (--n > 0);
return y;
}
long double
log2l(long double x)
{
long double z;
long double y;
int e;
int64_t hx, lx;
/* Test for domain */
GET_LDOUBLE_WORDS64 (hx, lx, x);
if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
return (-1.0L / (x - x));
if (hx < 0)
return (x - x) / (x - x);
if (hx >= 0x7fff000000000000LL)
return (x + x);
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl (x, &e);
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if ((e > 2) || (e < -2))
{
if (x < SQRTH)
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x * x;
y = x * (z * neval (z, R, 5) / deval (z, S, 5));
goto done;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH)
{
e -= 1;
x = 2.0 * x - 1.0L; /* 2x - 1 */
}
else
{
x = x - 1.0L;
}
z = x * x;
y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
y = y - 0.5 * z;
done:
/* Multiply log of fraction by log2(e)
* and base 2 exponent by 1
*/
z = y * LOG2EA;
z += x * LOG2EA;
z += y;
z += x;
z += e;
return (z);
}

283
ld128/e_logl.c Normal file
View file

@ -0,0 +1,283 @@
/* $OpenBSD: e_logl.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* logl.c
*
* Natural logarithm for 128-bit long double precision.
*
*
*
* SYNOPSIS:
*
* long double x, y, logl();
*
* y = logl( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. Use of a lookup table increases the speed of the routine.
* The program uses logarithms tabulated at intervals of 1/128 to
* cover the domain from approximately 0.7 to 1.4.
*
* On the interval [-1/128, +1/128] the logarithm of 1+x is approximated by
* log(1+x) = x - 0.5 x^2 + x^3 P(x) .
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.875, 1.125 100000 1.2e-34 4.1e-35
* IEEE 0.125, 8 100000 1.2e-34 4.1e-35
*
*
* WARNING:
*
* This program uses integer operations on bit fields of floating-point
* numbers. It does not work with data structures other than the
* structure assumed.
*
*/
#include <math.h>
#include "math_private.h"
/* log(1+x) = x - .5 x^2 + x^3 l(x)
-.0078125 <= x <= +.0078125
peak relative error 1.2e-37 */
static const long double
l3 = 3.333333333333333333333333333333336096926E-1L,
l4 = -2.499999999999999999999999999486853077002E-1L,
l5 = 1.999999999999999999999999998515277861905E-1L,
l6 = -1.666666666666666666666798448356171665678E-1L,
l7 = 1.428571428571428571428808945895490721564E-1L,
l8 = -1.249999999999999987884655626377588149000E-1L,
l9 = 1.111111111111111093947834982832456459186E-1L,
l10 = -1.000000000000532974938900317952530453248E-1L,
l11 = 9.090909090915566247008015301349979892689E-2L,
l12 = -8.333333211818065121250921925397567745734E-2L,
l13 = 7.692307559897661630807048686258659316091E-2L,
l14 = -7.144242754190814657241902218399056829264E-2L,
l15 = 6.668057591071739754844678883223432347481E-2L;
/* Lookup table of ln(t) - (t-1)
t = 0.5 + (k+26)/128)
k = 0, ..., 91 */
static const long double logtbl[92] = {
-5.5345593589352099112142921677820359632418E-2L,
-5.2108257402767124761784665198737642086148E-2L,
-4.8991686870576856279407775480686721935120E-2L,
-4.5993270766361228596215288742353061431071E-2L,
-4.3110481649613269682442058976885699556950E-2L,
-4.0340872319076331310838085093194799765520E-2L,
-3.7682072451780927439219005993827431503510E-2L,
-3.5131785416234343803903228503274262719586E-2L,
-3.2687785249045246292687241862699949178831E-2L,
-3.0347913785027239068190798397055267411813E-2L,
-2.8110077931525797884641940838507561326298E-2L,
-2.5972247078357715036426583294246819637618E-2L,
-2.3932450635346084858612873953407168217307E-2L,
-2.1988775689981395152022535153795155900240E-2L,
-2.0139364778244501615441044267387667496733E-2L,
-1.8382413762093794819267536615342902718324E-2L,
-1.6716169807550022358923589720001638093023E-2L,
-1.5138929457710992616226033183958974965355E-2L,
-1.3649036795397472900424896523305726435029E-2L,
-1.2244881690473465543308397998034325468152E-2L,
-1.0924898127200937840689817557742469105693E-2L,
-9.6875626072830301572839422532631079809328E-3L,
-8.5313926245226231463436209313499745894157E-3L,
-7.4549452072765973384933565912143044991706E-3L,
-6.4568155251217050991200599386801665681310E-3L,
-5.5356355563671005131126851708522185605193E-3L,
-4.6900728132525199028885749289712348829878E-3L,
-3.9188291218610470766469347968659624282519E-3L,
-3.2206394539524058873423550293617843896540E-3L,
-2.5942708080877805657374888909297113032132E-3L,
-2.0385211375711716729239156839929281289086E-3L,
-1.5522183228760777967376942769773768850872E-3L,
-1.1342191863606077520036253234446621373191E-3L,
-7.8340854719967065861624024730268350459991E-4L,
-4.9869831458030115699628274852562992756174E-4L,
-2.7902661731604211834685052867305795169688E-4L,
-1.2335696813916860754951146082826952093496E-4L,
-3.0677461025892873184042490943581654591817E-5L,
#define ZERO logtbl[38]
0.0000000000000000000000000000000000000000E0L,
-3.0359557945051052537099938863236321874198E-5L,
-1.2081346403474584914595395755316412213151E-4L,
-2.7044071846562177120083903771008342059094E-4L,
-4.7834133324631162897179240322783590830326E-4L,
-7.4363569786340080624467487620270965403695E-4L,
-1.0654639687057968333207323853366578860679E-3L,
-1.4429854811877171341298062134712230604279E-3L,
-1.8753781835651574193938679595797367137975E-3L,
-2.3618380914922506054347222273705859653658E-3L,
-2.9015787624124743013946600163375853631299E-3L,
-3.4938307889254087318399313316921940859043E-3L,
-4.1378413103128673800485306215154712148146E-3L,
-4.8328735414488877044289435125365629849599E-3L,
-5.5782063183564351739381962360253116934243E-3L,
-6.3731336597098858051938306767880719015261E-3L,
-7.2169643436165454612058905294782949315193E-3L,
-8.1090214990427641365934846191367315083867E-3L,
-9.0486422112807274112838713105168375482480E-3L,
-1.0035177140880864314674126398350812606841E-2L,
-1.1067990155502102718064936259435676477423E-2L,
-1.2146457974158024928196575103115488672416E-2L,
-1.3269969823361415906628825374158424754308E-2L,
-1.4437927104692837124388550722759686270765E-2L,
-1.5649743073340777659901053944852735064621E-2L,
-1.6904842527181702880599758489058031645317E-2L,
-1.8202661505988007336096407340750378994209E-2L,
-1.9542647000370545390701192438691126552961E-2L,
-2.0924256670080119637427928803038530924742E-2L,
-2.2346958571309108496179613803760727786257E-2L,
-2.3810230892650362330447187267648486279460E-2L,
-2.5313561699385640380910474255652501521033E-2L,
-2.6856448685790244233704909690165496625399E-2L,
-2.8438398935154170008519274953860128449036E-2L,
-3.0058928687233090922411781058956589863039E-2L,
-3.1717563112854831855692484086486099896614E-2L,
-3.3413836095418743219397234253475252001090E-2L,
-3.5147290019036555862676702093393332533702E-2L,
-3.6917475563073933027920505457688955423688E-2L,
-3.8723951502862058660874073462456610731178E-2L,
-4.0566284516358241168330505467000838017425E-2L,
-4.2444048996543693813649967076598766917965E-2L,
-4.4356826869355401653098777649745233339196E-2L,
-4.6304207416957323121106944474331029996141E-2L,
-4.8285787106164123613318093945035804818364E-2L,
-5.0301169421838218987124461766244507342648E-2L,
-5.2349964705088137924875459464622098310997E-2L,
-5.4431789996103111613753440311680967840214E-2L,
-5.6546268881465384189752786409400404404794E-2L,
-5.8693031345788023909329239565012647817664E-2L,
-6.0871713627532018185577188079210189048340E-2L,
-6.3081958078862169742820420185833800925568E-2L,
-6.5323413029406789694910800219643791556918E-2L,
-6.7595732653791419081537811574227049288168E-2L
};
/* ln(2) = ln2a + ln2b with extended precision. */
static const long double
ln2a = 6.93145751953125e-1L,
ln2b = 1.4286068203094172321214581765680755001344E-6L;
long double
logl(long double x)
{
long double z, y, w;
ieee_quad_shape_type u, t;
unsigned int m;
int k, e;
u.value = x;
m = u.parts32.mswhi;
/* Check for IEEE special cases. */
k = m & 0x7fffffff;
/* log(0) = -infinity. */
if ((k | u.parts32.mswlo | u.parts32.lswhi | u.parts32.lswlo) == 0)
{
return -0.5L / ZERO;
}
/* log ( x < 0 ) = NaN */
if (m & 0x80000000)
{
return (x - x) / ZERO;
}
/* log (infinity or NaN) */
if (k >= 0x7fff0000)
{
return x + x;
}
/* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */
e = (int) (m >> 16) - (int) 0x3ffe;
m &= 0xffff;
u.parts32.mswhi = m | 0x3ffe0000;
m |= 0x10000;
/* Find lookup table index k from high order bits of the significand. */
if (m < 0x16800)
{
k = (m - 0xff00) >> 9;
/* t is the argument 0.5 + (k+26)/128
of the nearest item to u in the lookup table. */
t.parts32.mswhi = 0x3fff0000 + (k << 9);
t.parts32.mswlo = 0;
t.parts32.lswhi = 0;
t.parts32.lswlo = 0;
u.parts32.mswhi += 0x10000;
e -= 1;
k += 64;
}
else
{
k = (m - 0xfe00) >> 10;
t.parts32.mswhi = 0x3ffe0000 + (k << 10);
t.parts32.mswlo = 0;
t.parts32.lswhi = 0;
t.parts32.lswlo = 0;
}
/* On this interval the table is not used due to cancellation error. */
if ((x <= 1.0078125L) && (x >= 0.9921875L))
{
z = x - 1.0L;
k = 64;
t.value = 1.0L;
e = 0;
}
else
{
/* log(u) = log( t u/t ) = log(t) + log(u/t)
log(t) is tabulated in the lookup table.
Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t.
cf. Cody & Waite. */
z = (u.value - t.value) / t.value;
}
/* Series expansion of log(1+z). */
w = z * z;
y = ((((((((((((l15 * z
+ l14) * z
+ l13) * z
+ l12) * z
+ l11) * z
+ l10) * z
+ l9) * z
+ l8) * z
+ l7) * z
+ l6) * z
+ l5) * z
+ l4) * z
+ l3) * z * w;
y -= 0.5 * w;
y += e * ln2b; /* Base 2 exponent offset times ln(2). */
y += z;
y += logtbl[k-26]; /* log(t) - (t-1) */
y += (t.value - 1.0L);
y += e * ln2a;
return y;
}

439
ld128/e_powl.c Normal file
View file

@ -0,0 +1,439 @@
/*
* ====================================================
* 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) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* powl(x,y) return x**y
*
* n
* Method: Let x = 2 * (1+f)
* 1. Compute and return log2(x) in two pieces:
* log2(x) = w1 + w2,
* where w1 has 113-53 = 60 bit trailing zeros.
* 2. Perform y*log2(x) = n+y' by simulating muti-precision
* arithmetic, where |y'|<=0.5.
* 3. Return x**y = 2**n*exp(y'*log2)
*
* Special cases:
* 1. (anything) ** 0 is 1
* 2. (anything) ** 1 is itself
* 3. (anything) ** NAN is NAN
* 4. NAN ** (anything except 0) is NAN
* 5. +-(|x| > 1) ** +INF is +INF
* 6. +-(|x| > 1) ** -INF is +0
* 7. +-(|x| < 1) ** +INF is +0
* 8. +-(|x| < 1) ** -INF is +INF
* 9. +-1 ** +-INF is NAN
* 10. +0 ** (+anything except 0, NAN) is +0
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
* 12. +0 ** (-anything except 0, NAN) is +INF
* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
* 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
* 15. +INF ** (+anything except 0,NAN) is +INF
* 16. +INF ** (-anything except 0,NAN) is +0
* 17. -INF ** (anything) = -0 ** (-anything)
* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
* 19. (-anything except 0 and inf) ** (non-integer) is NAN
*
*/
#include <math.h>
#include "math_private.h"
static const long double bp[] = {
1.0L,
1.5L,
};
/* log_2(1.5) */
static const long double dp_h[] = {
0.0,
5.8496250072115607565592654282227158546448E-1L
};
/* Low part of log_2(1.5) */
static const long double dp_l[] = {
0.0,
1.0579781240112554492329533686862998106046E-16L
};
static const long double zero = 0.0L,
one = 1.0L,
two = 2.0L,
two113 = 1.0384593717069655257060992658440192E34L,
huge = 1.0e3000L,
tiny = 1.0e-3000L;
/* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2))
z = (x-1)/(x+1)
1 <= x <= 1.25
Peak relative error 2.3e-37 */
static const long double LN[] =
{
-3.0779177200290054398792536829702930623200E1L,
6.5135778082209159921251824580292116201640E1L,
-4.6312921812152436921591152809994014413540E1L,
1.2510208195629420304615674658258363295208E1L,
-9.9266909031921425609179910128531667336670E-1L
};
static const long double LD[] =
{
-5.129862866715009066465422805058933131960E1L,
1.452015077564081884387441590064272782044E2L,
-1.524043275549860505277434040464085593165E2L,
7.236063513651544224319663428634139768808E1L,
-1.494198912340228235853027849917095580053E1L
/* 1.0E0 */
};
/* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2)))
0 <= x <= 0.5
Peak relative error 5.7e-38 */
static const long double PN[] =
{
5.081801691915377692446852383385968225675E8L,
9.360895299872484512023336636427675327355E6L,
4.213701282274196030811629773097579432957E4L,
5.201006511142748908655720086041570288182E1L,
9.088368420359444263703202925095675982530E-3L,
};
static const long double PD[] =
{
3.049081015149226615468111430031590411682E9L,
1.069833887183886839966085436512368982758E8L,
8.259257717868875207333991924545445705394E5L,
1.872583833284143212651746812884298360922E3L,
/* 1.0E0 */
};
static const long double
/* ln 2 */
lg2 = 6.9314718055994530941723212145817656807550E-1L,
lg2_h = 6.9314718055994528622676398299518041312695E-1L,
lg2_l = 2.3190468138462996154948554638754786504121E-17L,
ovt = 8.0085662595372944372e-0017L,
/* 2/(3*log(2)) */
cp = 9.6179669392597560490661645400126142495110E-1L,
cp_h = 9.6179669392597555432899980587535537779331E-1L,
cp_l = 5.0577616648125906047157785230014751039424E-17L;
long double
powl(long double x, long double y)
{
long double z, ax, z_h, z_l, p_h, p_l;
long double yy1, t1, t2, r, s, t, u, v, w;
long double s2, s_h, s_l, t_h, t_l;
int32_t i, j, k, yisint, n;
u_int32_t ix, iy;
int32_t hx, hy;
ieee_quad_shape_type o, p, q;
p.value = x;
hx = p.parts32.mswhi;
ix = hx & 0x7fffffff;
q.value = y;
hy = q.parts32.mswhi;
iy = hy & 0x7fffffff;
/* y==zero: x**0 = 1 */
if ((iy | q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0)
return one;
/* 1.0**y = 1; -1.0**+-Inf = 1 */
if (x == one)
return one;
if (x == -1.0L && iy == 0x7fff0000
&& (q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0)
return one;
/* +-NaN return x+y */
if ((ix > 0x7fff0000)
|| ((ix == 0x7fff0000)
&& ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) != 0))
|| (iy > 0x7fff0000)
|| ((iy == 0x7fff0000)
&& ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) != 0)))
return x + y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if (hx < 0)
{
if (iy >= 0x40700000) /* 2^113 */
yisint = 2; /* even integer y */
else if (iy >= 0x3fff0000) /* 1.0 */
{
if (floorl (y) == y)
{
z = 0.5 * y;
if (floorl (z) == z)
yisint = 2;
else
yisint = 1;
}
}
}
/* special value of y */
if ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0)
{
if (iy == 0x7fff0000) /* y is +-inf */
{
if (((ix - 0x3fff0000) | p.parts32.mswlo | p.parts32.lswhi |
p.parts32.lswlo) == 0)
return y - y; /* +-1**inf is NaN */
else if (ix >= 0x3fff0000) /* (|x|>1)**+-inf = inf,0 */
return (hy >= 0) ? y : zero;
else /* (|x|<1)**-,+inf = inf,0 */
return (hy < 0) ? -y : zero;
}
if (iy == 0x3fff0000)
{ /* y is +-1 */
if (hy < 0)
return one / x;
else
return x;
}
if (hy == 0x40000000)
return x * x; /* y is 2 */
if (hy == 0x3ffe0000)
{ /* y is 0.5 */
if (hx >= 0) /* x >= +0 */
return sqrtl (x);
}
}
ax = fabsl (x);
/* special value of x */
if ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) == 0)
{
if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000)
{
z = ax; /*x is +-0,+-inf,+-1 */
if (hy < 0)
z = one / z; /* z = (1/|x|) */
if (hx < 0)
{
if (((ix - 0x3fff0000) | yisint) == 0)
{
z = (z - z) / (z - z); /* (-1)**non-int is NaN */
}
else if (yisint == 1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
return z;
}
}
/* (x<0)**(non-int) is NaN */
if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
return (x - x) / (x - x);
/* |y| is huge.
2^-16495 = 1/2 of smallest representable value.
If (1 - 1/131072)^y underflows, y > 1.4986e9 */
if (iy > 0x401d654b)
{
/* if (1 - 2^-113)^y underflows, y > 1.1873e38 */
if (iy > 0x407d654b)
{
if (ix <= 0x3ffeffff)
return (hy < 0) ? huge * huge : tiny * tiny;
if (ix >= 0x3fff0000)
return (hy > 0) ? huge * huge : tiny * tiny;
}
/* over/underflow if x is not close to one */
if (ix < 0x3ffeffff)
return (hy < 0) ? huge * huge : tiny * tiny;
if (ix > 0x3fff0000)
return (hy > 0) ? huge * huge : tiny * tiny;
}
n = 0;
/* take care subnormal number */
if (ix < 0x00010000)
{
ax *= two113;
n -= 113;
o.value = ax;
ix = o.parts32.mswhi;
}
n += ((ix) >> 16) - 0x3fff;
j = ix & 0x0000ffff;
/* determine interval */
ix = j | 0x3fff0000; /* normalize ix */
if (j <= 0x3988)
k = 0; /* |x|<sqrt(3/2) */
else if (j < 0xbb67)
k = 1; /* |x|<sqrt(3) */
else
{
k = 0;
n += 1;
ix -= 0x00010000;
}
o.value = ax;
o.parts32.mswhi = ix;
ax = o.value;
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one / (ax + bp[k]);
s = u * v;
s_h = s;
o.value = s_h;
o.parts32.lswlo = 0;
o.parts32.lswhi &= 0xf8000000;
s_h = o.value;
/* t_h=ax+bp[k] High */
t_h = ax + bp[k];
o.value = t_h;
o.parts32.lswlo = 0;
o.parts32.lswhi &= 0xf8000000;
t_h = o.value;
t_l = ax - (t_h - bp[k]);
s_l = v * ((u - s_h * t_h) - s_h * t_l);
/* compute log(ax) */
s2 = s * s;
u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4])));
v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2))));
r = s2 * s2 * u / v;
r += s_l * (s_h + s);
s2 = s_h * s_h;
t_h = 3.0 + s2 + r;
o.value = t_h;
o.parts32.lswlo = 0;
o.parts32.lswhi &= 0xf8000000;
t_h = o.value;
t_l = r - ((t_h - 3.0) - s2);
/* u+v = s*(1+...) */
u = s_h * t_h;
v = s_l * t_h + t_l * s;
/* 2/(3log2)*(s+...) */
p_h = u + v;
o.value = p_h;
o.parts32.lswlo = 0;
o.parts32.lswhi &= 0xf8000000;
p_h = o.value;
p_l = v - (p_h - u);
z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l * p_h + p_l * cp + dp_l[k];
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (long double) n;
t1 = (((z_h + z_l) + dp_h[k]) + t);
o.value = t1;
o.parts32.lswlo = 0;
o.parts32.lswhi &= 0xf8000000;
t1 = o.value;
t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
/* s (sign of result -ve**odd) = -1 else = 1 */
s = one;
if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
s = -one; /* (-ve)**(odd int) */
/* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */
yy1 = y;
o.value = yy1;
o.parts32.lswlo = 0;
o.parts32.lswhi &= 0xf8000000;
yy1 = o.value;
p_l = (y - yy1) * t1 + y * t2;
p_h = yy1 * t1;
z = p_l + p_h;
o.value = z;
j = o.parts32.mswhi;
if (j >= 0x400d0000) /* z >= 16384 */
{
/* if z > 16384 */
if (((j - 0x400d0000) | o.parts32.mswlo | o.parts32.lswhi |
o.parts32.lswlo) != 0)
return s * huge * huge; /* overflow */
else
{
if (p_l + ovt > z - p_h)
return s * huge * huge; /* overflow */
}
}
else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */
{
/* z < -16495 */
if (((j - 0xc00d01bc) | o.parts32.mswlo | o.parts32.lswhi |
o.parts32.lswlo)
!= 0)
return s * tiny * tiny; /* underflow */
else
{
if (p_l <= z - p_h)
return s * tiny * tiny; /* underflow */
}
}
/* compute 2**(p_h+p_l) */
i = j & 0x7fffffff;
k = (i >> 16) - 0x3fff;
n = 0;
if (i > 0x3ffe0000)
{ /* if |z| > 0.5, set n = [z+0.5] */
n = floorl (z + 0.5L);
t = n;
p_h -= t;
}
t = p_l + p_h;
o.value = t;
o.parts32.lswlo = 0;
o.parts32.lswhi &= 0xf8000000;
t = o.value;
u = t * lg2_h;
v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
z = u + v;
w = v - (z - u);
/* exp(z) */
t = z * z;
u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4])));
v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t)));
t1 = z - t * u / v;
r = (z * t1) / (t1 - two) - (w + z * w);
z = one - (r - z);
o.value = z;
j = o.parts32.mswhi;
j += (n << 16);
if ((j >> 16) <= 0)
z = scalbnl (z, n); /* subnormal output */
else
{
o.parts32.mswhi = j;
z = o.value;
}
return s * z;
}

104
ld128/e_sinhl.c Normal file
View file

@ -0,0 +1,104 @@
/* @(#)e_sinh.c 5.1 93/09/24 */
/*
* ====================================================
* 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) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* sinhl(x)
* Method :
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
* 1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
* 2.
* E + E/(E+1)
* 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x)
* 2
*
* 25 <= x <= lnovft : sinhl(x) := expl(x)/2
* lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2)
* ln2ovft < x : sinhl(x) := x*shuge (overflow)
*
* Special cases:
* sinhl(x) is |x| if x is +INF, -INF, or NaN.
* only sinhl(0)=0 is exact for finite x.
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0, shuge = 1.0e4931L,
ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
long double
sinhl(long double x)
{
long double t, w, h;
u_int32_t jx, ix;
ieee_quad_shape_type u;
/* Words of |x|. */
u.value = x;
jx = u.parts32.mswhi;
ix = jx & 0x7fffffff;
/* x is INF or NaN */
if (ix >= 0x7fff0000)
return x + x;
h = 0.5;
if (jx & 0x80000000)
h = -h;
/* Absolute value of x. */
u.parts32.mswhi = ix;
/* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */
if (ix <= 0x40044000)
{
if (ix < 0x3fc60000) /* |x| < 2^-57 */
if (shuge + x > one)
return x; /* sinh(tiny) = tiny with inexact */
t = expm1l (u.value);
if (ix < 0x3fff0000)
return h * (2.0 * t - t * t / (t + one));
return h * (t + t / (t + one));
}
/* |x| in [40, log(maxdouble)] return 0.5*exp(|x|) */
if (ix <= 0x400c62e3) /* 11356.375 */
return h * expl (u.value);
/* |x| in [log(maxdouble), overflowthreshold]
Overflow threshold is log(2 * maxdouble). */
if (u.value <= ovf_thresh)
{
w = expl (0.5 * u.value);
t = h * w;
return t * w;
}
/* |x| > overflowthreshold, sinhl(x) overflow */
return x * shuge;
}

45
ld128/e_tgammal.c Normal file
View file

@ -0,0 +1,45 @@
/* $OpenBSD: e_tgammal.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */
/*
* Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
#include <math.h>
#include "math_private.h"
long double
tgammal(long double x)
{
int64_t i0,i1;
GET_LDOUBLE_WORDS64(i0,i1,x);
if (((i0&0x7fffffffffffffffLL)|i1) == 0) {
signgam = 0;
return (1.0/x);
}
if (i0<0 && (u_int64_t)i0<0xffff000000000000ULL && rintl(x)==x) {
signgam = 0;
return (x-x)/(x-x);
}
if (i0==0xffff000000000000ULL && i1==0) {
signgam = 0;
return (x-x);
}
return expl(lgammal(x));
}

69
ld128/s_asinhl.c Normal file
View file

@ -0,0 +1,69 @@
/* @(#)s_asinh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* asinhl(x)
* Method :
* Based on
* asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ]
* we have
* asinhl(x) := x if 1+x*x=1,
* := signl(x)*(logl(x)+ln2)) for large |x|, else
* := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else
* := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2)))
*/
#include <math.h>
#include "math_private.h"
static const long double
one = 1.0L,
ln2 = 6.931471805599453094172321214581765681e-1L,
huge = 1.0e+4900L;
long double
asinhl(long double x)
{
long double t, w;
int32_t ix, sign;
ieee_quad_shape_type u;
u.value = x;
sign = u.parts32.mswhi;
ix = sign & 0x7fffffff;
if (ix == 0x7fff0000)
return x + x; /* x is inf or NaN */
if (ix < 0x3fc70000)
{ /* |x| < 2^ -56 */
if (huge + x > one)
return x; /* return x inexact except 0 */
}
u.parts32.mswhi = ix;
if (ix > 0x40350000)
{ /* |x| > 2 ^ 54 */
w = logl (u.value) + ln2;
}
else if (ix >0x40000000)
{ /* 2^ 54 > |x| > 2.0 */
t = u.value;
w = logl (2.0 * t + one / (sqrtl (x * x + one) + t));
}
else
{ /* 2.0 > |x| > 2 ^ -56 */
t = x * x;
w = log1pl (u.value + t / (one + sqrtl (one + t)));
}
if (sign & 0x80000000)
return -w;
else
return w;
}

69
ld128/s_ceill.c Normal file
View file

@ -0,0 +1,69 @@
/* @(#)s_ceil.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* ceill(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to ceil(x).
*/
#include <math.h>
#include "math_private.h"
static const long double huge = 1.0e4930L;
long double
ceill(long double x)
{
int64_t i0,i1,jj0;
u_int64_t i,j;
GET_LDOUBLE_WORDS64(i0,i1,x);
jj0 = ((i0>>48)&0x7fff)-0x3fff;
if(jj0<48) {
if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
if(i0<0) {i0=0x8000000000000000ULL;i1=0;}
else if((i0|i1)!=0) { i0=0x3fff000000000000ULL;i1=0;}
}
} else {
i = (0x0000ffffffffffffULL)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(i0>0) i0 += (0x0001000000000000LL)>>jj0;
i0 &= (~i); i1=0;
}
}
} else if (jj0>111) {
if(jj0==0x4000) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = -1ULL>>(jj0-48);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(i0>0) {
if(jj0==48) i0+=1;
else {
j = i1+(1LL<<(112-jj0));
if(j<i1) i0 +=1 ; /* got a carry */
i1=j;
}
}
i1 &= (~i);
}
}
SET_LDOUBLE_WORDS64(x,i0,i1);
return x;
}

926
ld128/s_erfl.c Normal file
View file

@ -0,0 +1,926 @@
/*
* ====================================================
* 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) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* double erf(double x)
* double erfc(double x)
* x
* 2 |\
* erf(x) = --------- | exp(-t*t)dt
* sqrt(pi) \|
* 0
*
* erfc(x) = 1-erf(x)
* Note that
* erf(-x) = -erf(x)
* erfc(-x) = 2 - erfc(x)
*
* Method:
* 1. erf(x) = x + x*R(x^2) for |x| in [0, 7/8]
* Remark. The formula is derived by noting
* erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
* and that
* 2/sqrt(pi) = 1.128379167095512573896158903121545171688
* is close to one.
*
* 1a. erf(x) = 1 - erfc(x), for |x| > 1.0
* erfc(x) = 1 - erf(x) if |x| < 1/4
*
* 2. For |x| in [7/8, 1], let s = |x| - 1, and
* c = 0.84506291151 rounded to single (24 bits)
* erf(s + c) = sign(x) * (c + P1(s)/Q1(s))
* Remark: here we use the taylor series expansion at x=1.
* erf(1+s) = erf(1) + s*Poly(s)
* = 0.845.. + P1(s)/Q1(s)
* Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
*
* 3. For x in [1/4, 5/4],
* erfc(s + const) = erfc(const) + s P1(s)/Q1(s)
* for const = 1/4, 3/8, ..., 9/8
* and 0 <= s <= 1/8 .
*
* 4. For x in [5/4, 107],
* erfc(x) = (1/x)*exp(-x*x-0.5625 + R(z))
* z=1/x^2
* The interval is partitioned into several segments
* of width 1/8 in 1/x.
*
* Note1:
* To compute exp(-x*x-0.5625+R/S), let s be a single
* precision number and s := x; then
* -x*x = -s*s + (s-x)*(s+x)
* exp(-x*x-0.5626+R/S) =
* exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
* Note2:
* Here 4 and 5 make use of the asymptotic series
* exp(-x*x)
* erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
* x*sqrt(pi)
*
* 5. For inf > x >= 107
* erf(x) = sign(x) *(1 - tiny) (raise inexact)
* erfc(x) = tiny*tiny (raise underflow) if x > 0
* = 2 - tiny if x<0
*
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*/
#include <math.h>
#include "math_private.h"
/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
static long double
neval (long double x, const long double *p, int n)
{
long double y;
p += n;
y = *p--;
do
{
y = y * x + *p--;
}
while (--n > 0);
return y;
}
/* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
static long double
deval (long double x, const long double *p, int n)
{
long double y;
p += n;
y = x + *p--;
do
{
y = y * x + *p--;
}
while (--n > 0);
return y;
}
static const long double
tiny = 1e-4931L,
one = 1.0L,
two = 2.0L,
/* 2/sqrt(pi) - 1 */
efx = 1.2837916709551257389615890312154517168810E-1L,
/* 8 * (2/sqrt(pi) - 1) */
efx8 = 1.0270333367641005911692712249723613735048E0L;
/* erf(x) = x + x R(x^2)
0 <= x <= 7/8
Peak relative error 1.8e-35 */
#define NTN1 8
static const long double TN1[NTN1 + 1] =
{
-3.858252324254637124543172907442106422373E10L,
9.580319248590464682316366876952214879858E10L,
1.302170519734879977595901236693040544854E10L,
2.922956950426397417800321486727032845006E9L,
1.764317520783319397868923218385468729799E8L,
1.573436014601118630105796794840834145120E7L,
4.028077380105721388745632295157816229289E5L,
1.644056806467289066852135096352853491530E4L,
3.390868480059991640235675479463287886081E1L
};
#define NTD1 8
static const long double TD1[NTD1 + 1] =
{
-3.005357030696532927149885530689529032152E11L,
-1.342602283126282827411658673839982164042E11L,
-2.777153893355340961288511024443668743399E10L,
-3.483826391033531996955620074072768276974E9L,
-2.906321047071299585682722511260895227921E8L,
-1.653347985722154162439387878512427542691E7L,
-6.245520581562848778466500301865173123136E5L,
-1.402124304177498828590239373389110545142E4L,
-1.209368072473510674493129989468348633579E2L
/* 1.0E0 */
};
/* erf(z+1) = erf_const + P(z)/Q(z)
-.125 <= z <= 0
Peak relative error 7.3e-36 */
static const long double erf_const = 0.845062911510467529296875L;
#define NTN2 8
static const long double TN2[NTN2 + 1] =
{
-4.088889697077485301010486931817357000235E1L,
7.157046430681808553842307502826960051036E3L,
-2.191561912574409865550015485451373731780E3L,
2.180174916555316874988981177654057337219E3L,
2.848578658049670668231333682379720943455E2L,
1.630362490952512836762810462174798925274E2L,
6.317712353961866974143739396865293596895E0L,
2.450441034183492434655586496522857578066E1L,
5.127662277706787664956025545897050896203E-1L
};
#define NTD2 8
static const long double TD2[NTD2 + 1] =
{
1.731026445926834008273768924015161048885E4L,
1.209682239007990370796112604286048173750E4L,
1.160950290217993641320602282462976163857E4L,
5.394294645127126577825507169061355698157E3L,
2.791239340533632669442158497532521776093E3L,
8.989365571337319032943005387378993827684E2L,
2.974016493766349409725385710897298069677E2L,
6.148192754590376378740261072533527271947E1L,
1.178502892490738445655468927408440847480E1L
/* 1.0E0 */
};
/* erfc(x + 0.25) = erfc(0.25) + x R(x)
0 <= x < 0.125
Peak relative error 1.4e-35 */
#define NRNr13 8
static const long double RNr13[NRNr13 + 1] =
{
-2.353707097641280550282633036456457014829E3L,
3.871159656228743599994116143079870279866E2L,
-3.888105134258266192210485617504098426679E2L,
-2.129998539120061668038806696199343094971E1L,
-8.125462263594034672468446317145384108734E1L,
8.151549093983505810118308635926270319660E0L,
-5.033362032729207310462422357772568553670E0L,
-4.253956621135136090295893547735851168471E-2L,
-8.098602878463854789780108161581050357814E-2L
};
#define NRDr13 7
static const long double RDr13[NRDr13 + 1] =
{
2.220448796306693503549505450626652881752E3L,
1.899133258779578688791041599040951431383E2L,
1.061906712284961110196427571557149268454E3L,
7.497086072306967965180978101974566760042E1L,
2.146796115662672795876463568170441327274E2L,
1.120156008362573736664338015952284925592E1L,
2.211014952075052616409845051695042741074E1L,
6.469655675326150785692908453094054988938E-1L
/* 1.0E0 */
};
/* erfc(0.25) = C13a + C13b to extra precision. */
static const long double C13a = 0.723663330078125L;
static const long double C13b = 1.0279753638067014931732235184287934646022E-5L;
/* erfc(x + 0.375) = erfc(0.375) + x R(x)
0 <= x < 0.125
Peak relative error 1.2e-35 */
#define NRNr14 8
static const long double RNr14[NRNr14 + 1] =
{
-2.446164016404426277577283038988918202456E3L,
6.718753324496563913392217011618096698140E2L,
-4.581631138049836157425391886957389240794E2L,
-2.382844088987092233033215402335026078208E1L,
-7.119237852400600507927038680970936336458E1L,
1.313609646108420136332418282286454287146E1L,
-6.188608702082264389155862490056401365834E0L,
-2.787116601106678287277373011101132659279E-2L,
-2.230395570574153963203348263549700967918E-2L
};
#define NRDr14 7
static const long double RDr14[NRDr14 + 1] =
{
2.495187439241869732696223349840963702875E3L,
2.503549449872925580011284635695738412162E2L,
1.159033560988895481698051531263861842461E3L,
9.493751466542304491261487998684383688622E1L,
2.276214929562354328261422263078480321204E2L,
1.367697521219069280358984081407807931847E1L,
2.276988395995528495055594829206582732682E1L,
7.647745753648996559837591812375456641163E-1L
/* 1.0E0 */
};
/* erfc(0.375) = C14a + C14b to extra precision. */
static const long double C14a = 0.5958709716796875L;
static const long double C14b = 1.2118885490201676174914080878232469565953E-5L;
/* erfc(x + 0.5) = erfc(0.5) + x R(x)
0 <= x < 0.125
Peak relative error 4.7e-36 */
#define NRNr15 8
static const long double RNr15[NRNr15 + 1] =
{
-2.624212418011181487924855581955853461925E3L,
8.473828904647825181073831556439301342756E2L,
-5.286207458628380765099405359607331669027E2L,
-3.895781234155315729088407259045269652318E1L,
-6.200857908065163618041240848728398496256E1L,
1.469324610346924001393137895116129204737E1L,
-6.961356525370658572800674953305625578903E0L,
5.145724386641163809595512876629030548495E-3L,
1.990253655948179713415957791776180406812E-2L
};
#define NRDr15 7
static const long double RDr15[NRDr15 + 1] =
{
2.986190760847974943034021764693341524962E3L,
5.288262758961073066335410218650047725985E2L,
1.363649178071006978355113026427856008978E3L,
1.921707975649915894241864988942255320833E2L,
2.588651100651029023069013885900085533226E2L,
2.628752920321455606558942309396855629459E1L,
2.455649035885114308978333741080991380610E1L,
1.378826653595128464383127836412100939126E0L
/* 1.0E0 */
};
/* erfc(0.5) = C15a + C15b to extra precision. */
static const long double C15a = 0.4794921875L;
static const long double C15b = 7.9346869534623172533461080354712635484242E-6L;
/* erfc(x + 0.625) = erfc(0.625) + x R(x)
0 <= x < 0.125
Peak relative error 5.1e-36 */
#define NRNr16 8
static const long double RNr16[NRNr16 + 1] =
{
-2.347887943200680563784690094002722906820E3L,
8.008590660692105004780722726421020136482E2L,
-5.257363310384119728760181252132311447963E2L,
-4.471737717857801230450290232600243795637E1L,
-4.849540386452573306708795324759300320304E1L,
1.140885264677134679275986782978655952843E1L,
-6.731591085460269447926746876983786152300E0L,
1.370831653033047440345050025876085121231E-1L,
2.022958279982138755020825717073966576670E-2L,
};
#define NRDr16 7
static const long double RDr16[NRDr16 + 1] =
{
3.075166170024837215399323264868308087281E3L,
8.730468942160798031608053127270430036627E2L,
1.458472799166340479742581949088453244767E3L,
3.230423687568019709453130785873540386217E2L,
2.804009872719893612081109617983169474655E2L,
4.465334221323222943418085830026979293091E1L,
2.612723259683205928103787842214809134746E1L,
2.341526751185244109722204018543276124997E0L,
/* 1.0E0 */
};
/* erfc(0.625) = C16a + C16b to extra precision. */
static const long double C16a = 0.3767547607421875L;
static const long double C16b = 4.3570693945275513594941232097252997287766E-6L;
/* erfc(x + 0.75) = erfc(0.75) + x R(x)
0 <= x < 0.125
Peak relative error 1.7e-35 */
#define NRNr17 8
static const long double RNr17[NRNr17 + 1] =
{
-1.767068734220277728233364375724380366826E3L,
6.693746645665242832426891888805363898707E2L,
-4.746224241837275958126060307406616817753E2L,
-2.274160637728782675145666064841883803196E1L,
-3.541232266140939050094370552538987982637E1L,
6.988950514747052676394491563585179503865E0L,
-5.807687216836540830881352383529281215100E0L,
3.631915988567346438830283503729569443642E-1L,
-1.488945487149634820537348176770282391202E-2L
};
#define NRDr17 7
static const long double RDr17[NRDr17 + 1] =
{
2.748457523498150741964464942246913394647E3L,
1.020213390713477686776037331757871252652E3L,
1.388857635935432621972601695296561952738E3L,
3.903363681143817750895999579637315491087E2L,
2.784568344378139499217928969529219886578E2L,
5.555800830216764702779238020065345401144E1L,
2.646215470959050279430447295801291168941E1L,
2.984905282103517497081766758550112011265E0L,
/* 1.0E0 */
};
/* erfc(0.75) = C17a + C17b to extra precision. */
static const long double C17a = 0.2888336181640625L;
static const long double C17b = 1.0748182422368401062165408589222625794046E-5L;
/* erfc(x + 0.875) = erfc(0.875) + x R(x)
0 <= x < 0.125
Peak relative error 2.2e-35 */
#define NRNr18 8
static const long double RNr18[NRNr18 + 1] =
{
-1.342044899087593397419622771847219619588E3L,
6.127221294229172997509252330961641850598E2L,
-4.519821356522291185621206350470820610727E2L,
1.223275177825128732497510264197915160235E1L,
-2.730789571382971355625020710543532867692E1L,
4.045181204921538886880171727755445395862E0L,
-4.925146477876592723401384464691452700539E0L,
5.933878036611279244654299924101068088582E-1L,
-5.557645435858916025452563379795159124753E-2L
};
#define NRDr18 7
static const long double RDr18[NRDr18 + 1] =
{
2.557518000661700588758505116291983092951E3L,
1.070171433382888994954602511991940418588E3L,
1.344842834423493081054489613250688918709E3L,
4.161144478449381901208660598266288188426E2L,
2.763670252219855198052378138756906980422E2L,
5.998153487868943708236273854747564557632E1L,
2.657695108438628847733050476209037025318E1L,
3.252140524394421868923289114410336976512E0L,
/* 1.0E0 */
};
/* erfc(0.875) = C18a + C18b to extra precision. */
static const long double C18a = 0.215911865234375L;
static const long double C18b = 1.3073705765341685464282101150637224028267E-5L;
/* erfc(x + 1.0) = erfc(1.0) + x R(x)
0 <= x < 0.125
Peak relative error 1.6e-35 */
#define NRNr19 8
static const long double RNr19[NRNr19 + 1] =
{
-1.139180936454157193495882956565663294826E3L,
6.134903129086899737514712477207945973616E2L,
-4.628909024715329562325555164720732868263E2L,
4.165702387210732352564932347500364010833E1L,
-2.286979913515229747204101330405771801610E1L,
1.870695256449872743066783202326943667722E0L,
-4.177486601273105752879868187237000032364E0L,
7.533980372789646140112424811291782526263E-1L,
-8.629945436917752003058064731308767664446E-2L
};
#define NRDr19 7
static const long double RDr19[NRDr19 + 1] =
{
2.744303447981132701432716278363418643778E3L,
1.266396359526187065222528050591302171471E3L,
1.466739461422073351497972255511919814273E3L,
4.868710570759693955597496520298058147162E2L,
2.993694301559756046478189634131722579643E2L,
6.868976819510254139741559102693828237440E1L,
2.801505816247677193480190483913753613630E1L,
3.604439909194350263552750347742663954481E0L,
/* 1.0E0 */
};
/* erfc(1.0) = C19a + C19b to extra precision. */
static const long double C19a = 0.15728759765625L;
static const long double C19b = 1.1609394035130658779364917390740703933002E-5L;
/* erfc(x + 1.125) = erfc(1.125) + x R(x)
0 <= x < 0.125
Peak relative error 3.6e-36 */
#define NRNr20 8
static const long double RNr20[NRNr20 + 1] =
{
-9.652706916457973956366721379612508047640E2L,
5.577066396050932776683469951773643880634E2L,
-4.406335508848496713572223098693575485978E2L,
5.202893466490242733570232680736966655434E1L,
-1.931311847665757913322495948705563937159E1L,
-9.364318268748287664267341457164918090611E-2L,
-3.306390351286352764891355375882586201069E0L,
7.573806045289044647727613003096916516475E-1L,
-9.611744011489092894027478899545635991213E-2L
};
#define NRDr20 7
static const long double RDr20[NRDr20 + 1] =
{
3.032829629520142564106649167182428189014E3L,
1.659648470721967719961167083684972196891E3L,
1.703545128657284619402511356932569292535E3L,
6.393465677731598872500200253155257708763E2L,
3.489131397281030947405287112726059221934E2L,
8.848641738570783406484348434387611713070E1L,
3.132269062552392974833215844236160958502E1L,
4.430131663290563523933419966185230513168E0L
/* 1.0E0 */
};
/* erfc(1.125) = C20a + C20b to extra precision. */
static const long double C20a = 0.111602783203125L;
static const long double C20b = 8.9850951672359304215530728365232161564636E-6L;
/* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
7/8 <= 1/x < 1
Peak relative error 1.4e-35 */
#define NRNr8 9
static const long double RNr8[NRNr8 + 1] =
{
3.587451489255356250759834295199296936784E1L,
5.406249749087340431871378009874875889602E2L,
2.931301290625250886238822286506381194157E3L,
7.359254185241795584113047248898753470923E3L,
9.201031849810636104112101947312492532314E3L,
5.749697096193191467751650366613289284777E3L,
1.710415234419860825710780802678697889231E3L,
2.150753982543378580859546706243022719599E2L,
8.740953582272147335100537849981160931197E0L,
4.876422978828717219629814794707963640913E-2L
};
#define NRDr8 8
static const long double RDr8[NRDr8 + 1] =
{
6.358593134096908350929496535931630140282E1L,
9.900253816552450073757174323424051765523E2L,
5.642928777856801020545245437089490805186E3L,
1.524195375199570868195152698617273739609E4L,
2.113829644500006749947332935305800887345E4L,
1.526438562626465706267943737310282977138E4L,
5.561370922149241457131421914140039411782E3L,
9.394035530179705051609070428036834496942E2L,
6.147019596150394577984175188032707343615E1L
/* 1.0E0 */
};
/* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
0.75 <= 1/x <= 0.875
Peak relative error 2.0e-36 */
#define NRNr7 9
static const long double RNr7[NRNr7 + 1] =
{
1.686222193385987690785945787708644476545E1L,
1.178224543567604215602418571310612066594E3L,
1.764550584290149466653899886088166091093E4L,
1.073758321890334822002849369898232811561E5L,
3.132840749205943137619839114451290324371E5L,
4.607864939974100224615527007793867585915E5L,
3.389781820105852303125270837910972384510E5L,
1.174042187110565202875011358512564753399E5L,
1.660013606011167144046604892622504338313E4L,
6.700393957480661937695573729183733234400E2L
};
#define NRDr7 9
static const long double RDr7[NRDr7 + 1] =
{
-1.709305024718358874701575813642933561169E3L,
-3.280033887481333199580464617020514788369E4L,
-2.345284228022521885093072363418750835214E5L,
-8.086758123097763971926711729242327554917E5L,
-1.456900414510108718402423999575992450138E6L,
-1.391654264881255068392389037292702041855E6L,
-6.842360801869939983674527468509852583855E5L,
-1.597430214446573566179675395199807533371E5L,
-1.488876130609876681421645314851760773480E4L,
-3.511762950935060301403599443436465645703E2L
/* 1.0E0 */
};
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
5/8 <= 1/x < 3/4
Peak relative error 1.9e-35 */
#define NRNr6 9
static const long double RNr6[NRNr6 + 1] =
{
1.642076876176834390623842732352935761108E0L,
1.207150003611117689000664385596211076662E2L,
2.119260779316389904742873816462800103939E3L,
1.562942227734663441801452930916044224174E4L,
5.656779189549710079988084081145693580479E4L,
1.052166241021481691922831746350942786299E5L,
9.949798524786000595621602790068349165758E4L,
4.491790734080265043407035220188849562856E4L,
8.377074098301530326270432059434791287601E3L,
4.506934806567986810091824791963991057083E2L
};
#define NRDr6 9
static const long double RDr6[NRDr6 + 1] =
{
-1.664557643928263091879301304019826629067E2L,
-3.800035902507656624590531122291160668452E3L,
-3.277028191591734928360050685359277076056E4L,
-1.381359471502885446400589109566587443987E5L,
-3.082204287382581873532528989283748656546E5L,
-3.691071488256738343008271448234631037095E5L,
-2.300482443038349815750714219117566715043E5L,
-6.873955300927636236692803579555752171530E4L,
-8.262158817978334142081581542749986845399E3L,
-2.517122254384430859629423488157361983661E2L
/* 1.00 */
};
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
1/2 <= 1/x < 5/8
Peak relative error 4.6e-36 */
#define NRNr5 10
static const long double RNr5[NRNr5 + 1] =
{
-3.332258927455285458355550878136506961608E-3L,
-2.697100758900280402659586595884478660721E-1L,
-6.083328551139621521416618424949137195536E0L,
-6.119863528983308012970821226810162441263E1L,
-3.176535282475593173248810678636522589861E2L,
-8.933395175080560925809992467187963260693E2L,
-1.360019508488475978060917477620199499560E3L,
-1.075075579828188621541398761300910213280E3L,
-4.017346561586014822824459436695197089916E2L,
-5.857581368145266249509589726077645791341E1L,
-2.077715925587834606379119585995758954399E0L
};
#define NRDr5 9
static const long double RDr5[NRDr5 + 1] =
{
3.377879570417399341550710467744693125385E-1L,
1.021963322742390735430008860602594456187E1L,
1.200847646592942095192766255154827011939E2L,
7.118915528142927104078182863387116942836E2L,
2.318159380062066469386544552429625026238E3L,
4.238729853534009221025582008928765281620E3L,
4.279114907284825886266493994833515580782E3L,
2.257277186663261531053293222591851737504E3L,
5.570475501285054293371908382916063822957E2L,
5.142189243856288981145786492585432443560E1L
/* 1.0E0 */
};
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
3/8 <= 1/x < 1/2
Peak relative error 2.0e-36 */
#define NRNr4 10
static const long double RNr4[NRNr4 + 1] =
{
3.258530712024527835089319075288494524465E-3L,
2.987056016877277929720231688689431056567E-1L,
8.738729089340199750734409156830371528862E0L,
1.207211160148647782396337792426311125923E2L,
8.997558632489032902250523945248208224445E2L,
3.798025197699757225978410230530640879762E3L,
9.113203668683080975637043118209210146846E3L,
1.203285891339933238608683715194034900149E4L,
8.100647057919140328536743641735339740855E3L,
2.383888249907144945837976899822927411769E3L,
2.127493573166454249221983582495245662319E2L
};
#define NRDr4 10
static const long double RDr4[NRDr4 + 1] =
{
-3.303141981514540274165450687270180479586E-1L,
-1.353768629363605300707949368917687066724E1L,
-2.206127630303621521950193783894598987033E2L,
-1.861800338758066696514480386180875607204E3L,
-8.889048775872605708249140016201753255599E3L,
-2.465888106627948210478692168261494857089E4L,
-3.934642211710774494879042116768390014289E4L,
-3.455077258242252974937480623730228841003E4L,
-1.524083977439690284820586063729912653196E4L,
-2.810541887397984804237552337349093953857E3L,
-1.343929553541159933824901621702567066156E2L
/* 1.0E0 */
};
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
1/4 <= 1/x < 3/8
Peak relative error 8.4e-37 */
#define NRNr3 11
static const long double RNr3[NRNr3 + 1] =
{
-1.952401126551202208698629992497306292987E-6L,
-2.130881743066372952515162564941682716125E-4L,
-8.376493958090190943737529486107282224387E-3L,
-1.650592646560987700661598877522831234791E-1L,
-1.839290818933317338111364667708678163199E0L,
-1.216278715570882422410442318517814388470E1L,
-4.818759344462360427612133632533779091386E1L,
-1.120994661297476876804405329172164436784E2L,
-1.452850765662319264191141091859300126931E2L,
-9.485207851128957108648038238656777241333E1L,
-2.563663855025796641216191848818620020073E1L,
-1.787995944187565676837847610706317833247E0L
};
#define NRDr3 10
static const long double RDr3[NRDr3 + 1] =
{
1.979130686770349481460559711878399476903E-4L,
1.156941716128488266238105813374635099057E-2L,
2.752657634309886336431266395637285974292E-1L,
3.482245457248318787349778336603569327521E0L,
2.569347069372696358578399521203959253162E1L,
1.142279000180457419740314694631879921561E2L,
3.056503977190564294341422623108332700840E2L,
4.780844020923794821656358157128719184422E2L,
4.105972727212554277496256802312730410518E2L,
1.724072188063746970865027817017067646246E2L,
2.815939183464818198705278118326590370435E1L
/* 1.0E0 */
};
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
1/8 <= 1/x < 1/4
Peak relative error 1.5e-36 */
#define NRNr2 11
static const long double RNr2[NRNr2 + 1] =
{
-2.638914383420287212401687401284326363787E-8L,
-3.479198370260633977258201271399116766619E-6L,
-1.783985295335697686382487087502222519983E-4L,
-4.777876933122576014266349277217559356276E-3L,
-7.450634738987325004070761301045014986520E-2L,
-7.068318854874733315971973707247467326619E-1L,
-4.113919921935944795764071670806867038732E0L,
-1.440447573226906222417767283691888875082E1L,
-2.883484031530718428417168042141288943905E1L,
-2.990886974328476387277797361464279931446E1L,
-1.325283914915104866248279787536128997331E1L,
-1.572436106228070195510230310658206154374E0L
};
#define NRDr2 10
static const long double RDr2[NRDr2 + 1] =
{
2.675042728136731923554119302571867799673E-6L,
2.170997868451812708585443282998329996268E-4L,
7.249969752687540289422684951196241427445E-3L,
1.302040375859768674620410563307838448508E-1L,
1.380202483082910888897654537144485285549E0L,
8.926594113174165352623847870299170069350E0L,
3.521089584782616472372909095331572607185E1L,
8.233547427533181375185259050330809105570E1L,
1.072971579885803033079469639073292840135E2L,
6.943803113337964469736022094105143158033E1L,
1.775695341031607738233608307835017282662E1L
/* 1.0E0 */
};
/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
1/128 <= 1/x < 1/8
Peak relative error 2.2e-36 */
#define NRNr1 9
static const long double RNr1[NRNr1 + 1] =
{
-4.250780883202361946697751475473042685782E-8L,
-5.375777053288612282487696975623206383019E-6L,
-2.573645949220896816208565944117382460452E-4L,
-6.199032928113542080263152610799113086319E-3L,
-8.262721198693404060380104048479916247786E-2L,
-6.242615227257324746371284637695778043982E-1L,
-2.609874739199595400225113299437099626386E0L,
-5.581967563336676737146358534602770006970E0L,
-5.124398923356022609707490956634280573882E0L,
-1.290865243944292370661544030414667556649E0L
};
#define NRDr1 8
static const long double RDr1[NRDr1 + 1] =
{
4.308976661749509034845251315983612976224E-6L,
3.265390126432780184125233455960049294580E-4L,
9.811328839187040701901866531796570418691E-3L,
1.511222515036021033410078631914783519649E-1L,
1.289264341917429958858379585970225092274E0L,
6.147640356182230769548007536914983522270E0L,
1.573966871337739784518246317003956180750E1L,
1.955534123435095067199574045529218238263E1L,
9.472613121363135472247929109615785855865E0L
/* 1.0E0 */
};
long double
erfl(long double x)
{
long double a, y, z;
int32_t i, ix, sign;
ieee_quad_shape_type u;
u.value = x;
sign = u.parts32.mswhi;
ix = sign & 0x7fffffff;
if (ix >= 0x7fff0000)
{ /* erf(nan)=nan */
i = ((sign & 0xffff0000) >> 31) << 1;
return (long double) (1 - i) + one / x; /* erf(+-inf)=+-1 */
}
if (ix >= 0x3fff0000) /* |x| >= 1.0 */
{
y = erfcl (x);
return (one - y);
/* return (one - erfcl (x)); */
}
u.parts32.mswhi = ix;
a = u.value;
z = x * x;
if (ix < 0x3ffec000) /* a < 0.875 */
{
if (ix < 0x3fc60000) /* |x|<2**-57 */
{
if (ix < 0x00080000)
return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */
return x + efx * x;
}
y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1);
}
else
{
a = a - one;
y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2);
}
if (sign & 0x80000000) /* x < 0 */
y = -y;
return( y );
}
long double
erfcl(long double x)
{
long double y, z, p, r;
int32_t i, ix, sign;
ieee_quad_shape_type u;
u.value = x;
sign = u.parts32.mswhi;
ix = sign & 0x7fffffff;
u.parts32.mswhi = ix;
if (ix >= 0x7fff0000)
{ /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (long double) (((u_int32_t) sign >> 31) << 1) + one / x;
}
if (ix < 0x3ffd0000) /* |x| <1/4 */
{
if (ix < 0x3f8d0000) /* |x|<2**-114 */
return one - x;
return one - erfl (x);
}
if (ix < 0x3fff4000) /* 1.25 */
{
x = u.value;
i = 8.0 * x;
switch (i)
{
case 2:
z = x - 0.25L;
y = C13b + z * neval (z, RNr13, NRNr13) / deval (z, RDr13, NRDr13);
y += C13a;
break;
case 3:
z = x - 0.375L;
y = C14b + z * neval (z, RNr14, NRNr14) / deval (z, RDr14, NRDr14);
y += C14a;
break;
case 4:
z = x - 0.5L;
y = C15b + z * neval (z, RNr15, NRNr15) / deval (z, RDr15, NRDr15);
y += C15a;
break;
case 5:
z = x - 0.625L;
y = C16b + z * neval (z, RNr16, NRNr16) / deval (z, RDr16, NRDr16);
y += C16a;
break;
case 6:
z = x - 0.75L;
y = C17b + z * neval (z, RNr17, NRNr17) / deval (z, RDr17, NRDr17);
y += C17a;
break;
case 7:
z = x - 0.875L;
y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18);
y += C18a;
break;
case 8:
z = x - 1.0L;
y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
y += C19a;
break;
case 9:
z = x - 1.125L;
y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20);
y += C20a;
break;
}
if (sign & 0x80000000)
y = 2.0L - y;
return y;
}
/* 1.25 < |x| < 107 */
if (ix < 0x4005ac00)
{
/* x < -9 */
if ((ix >= 0x40022000) && (sign & 0x80000000))
return two - tiny;
x = fabsl (x);
z = one / (x * x);
i = 8.0 / x;
switch (i)
{
default:
case 0:
p = neval (z, RNr1, NRNr1) / deval (z, RDr1, NRDr1);
break;
case 1:
p = neval (z, RNr2, NRNr2) / deval (z, RDr2, NRDr2);
break;
case 2:
p = neval (z, RNr3, NRNr3) / deval (z, RDr3, NRDr3);
break;
case 3:
p = neval (z, RNr4, NRNr4) / deval (z, RDr4, NRDr4);
break;
case 4:
p = neval (z, RNr5, NRNr5) / deval (z, RDr5, NRDr5);
break;
case 5:
p = neval (z, RNr6, NRNr6) / deval (z, RDr6, NRDr6);
break;
case 6:
p = neval (z, RNr7, NRNr7) / deval (z, RDr7, NRDr7);
break;
case 7:
p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8);
break;
}
u.value = x;
u.parts32.lswlo = 0;
u.parts32.lswhi &= 0xfe000000;
z = u.value;
r = expl (-z * z - 0.5625) *
expl ((z - x) * (z + x) + p);
if ((sign & 0x80000000) == 0)
return r / x;
else
return two - r / x;
}
else
{
if ((sign & 0x80000000) == 0)
return tiny * tiny;
else
return two - tiny;
}
}

162
ld128/s_expm1l.c Normal file
View file

@ -0,0 +1,162 @@
/* $OpenBSD: s_expm1l.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* expm1l.c
*
* Exponential function, minus 1
* 128-bit long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, expm1l();
*
* y = expm1l( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power, minus one.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
* in the basic range [-0.5 ln 2, 0.5 ln 2].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -79,+MAXLOG 100,000 1.7e-34 4.5e-35
*
*/
#include <errno.h>
#include <math.h>
#include "math_private.h"
/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
-.5 ln 2 < x < .5 ln 2
Theoretical peak relative error = 8.1e-36 */
static const long double
P0 = 2.943520915569954073888921213330863757240E8L,
P1 = -5.722847283900608941516165725053359168840E7L,
P2 = 8.944630806357575461578107295909719817253E6L,
P3 = -7.212432713558031519943281748462837065308E5L,
P4 = 4.578962475841642634225390068461943438441E4L,
P5 = -1.716772506388927649032068540558788106762E3L,
P6 = 4.401308817383362136048032038528753151144E1L,
P7 = -4.888737542888633647784737721812546636240E-1L,
Q0 = 1.766112549341972444333352727998584753865E9L,
Q1 = -7.848989743695296475743081255027098295771E8L,
Q2 = 1.615869009634292424463780387327037251069E8L,
Q3 = -2.019684072836541751428967854947019415698E7L,
Q4 = 1.682912729190313538934190635536631941751E6L,
Q5 = -9.615511549171441430850103489315371768998E4L,
Q6 = 3.697714952261803935521187272204485251835E3L,
Q7 = -8.802340681794263968892934703309274564037E1L,
/* Q8 = 1.000000000000000000000000000000000000000E0 */
/* C1 + C2 = ln 2 */
C1 = 6.93145751953125E-1L,
C2 = 1.428606820309417232121458176568075500134E-6L,
/* ln (2^16384 * (1 - 2^-113)) */
maxlog = 1.1356523406294143949491931077970764891253E4L,
/* ln 2^-114 */
minarg = -7.9018778583833765273564461846232128760607E1L, big = 1e4932L;
long double
expm1l(long double x)
{
long double px, qx, xx;
int32_t ix, sign;
ieee_quad_shape_type u;
int k;
/* Detect infinity and NaN. */
u.value = x;
ix = u.parts32.mswhi;
sign = ix & 0x80000000;
ix &= 0x7fffffff;
if (ix >= 0x7fff0000)
{
/* Infinity. */
if (((ix & 0xffff) | u.parts32.mswlo | u.parts32.lswhi |
u.parts32.lswlo) == 0)
{
if (sign)
return -1.0L;
else
return x;
}
/* NaN. No invalid exception. */
return x;
}
/* expm1(+- 0) = +- 0. */
if ((ix == 0) && (u.parts32.mswlo | u.parts32.lswhi | u.parts32.lswlo) == 0)
return x;
/* Overflow. */
if (x > maxlog)
return (big * big);
/* Minimum value. */
if (x < minarg)
return (4.0/big - 1.0L);
/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
xx = C1 + C2; /* ln 2. */
px = floorl (0.5 + x / xx);
k = px;
/* remainder times ln 2 */
x -= px * C1;
x -= px * C2;
/* Approximate exp(remainder ln 2). */
px = (((((((P7 * x
+ P6) * x
+ P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
qx = (((((((x
+ Q7) * x
+ Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
xx = x * x;
qx = x + (0.5 * xx + xx * px / qx);
/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
We have qx = exp(remainder ln 2) - 1, so
exp(x) - 1 = 2^k (qx + 1) - 1
= 2^k qx + 2^k - 1. */
px = ldexpl (1.0L, k);
x = px * qx + (px - 1.0);
return x;
}

71
ld128/s_floorl.c Normal file
View file

@ -0,0 +1,71 @@
/* @(#)s_floor.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* floorl(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to floor(x).
*/
#include <math.h>
#include "math_private.h"
static const long double huge = 1.0e4930L;
long double
floorl(long double x)
{
int64_t i0,i1,jj0;
u_int64_t i,j;
GET_LDOUBLE_WORDS64(i0,i1,x);
jj0 = ((i0>>48)&0x7fff)-0x3fff;
if(jj0<48) {
if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {
if(i0>=0)
return 0.0L;
else if(((i0&0x7fffffffffffffffLL)|i1)!=0)
return -1.0L;
}
} else {
i = (0x0000ffffffffffffULL)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(i0<0) i0 += (0x0001000000000000LL)>>jj0;
i0 &= (~i); i1=0;
}
}
} else if (jj0>111) {
if(jj0==0x4000) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = -1ULL>>(jj0-48);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(i0<0) {
if(jj0==48) i0+=1;
else {
j = i1+(1LL<<(112-jj0));
if(j<i1) i0 +=1 ; /* got a carry */
i1=j;
}
}
i1 &= (~i);
}
}
SET_LDOUBLE_WORDS64(x,i0,i1);
return x;
}

247
ld128/s_log1pl.c Normal file
View file

@ -0,0 +1,247 @@
/* $OpenBSD: s_log1pl.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* log1pl.c
*
* Relative error logarithm
* Natural logarithm of 1+x, 128-bit long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, log1pl();
*
* y = log1pl( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of 1+x.
*
* The argument 1+x is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
*
* Otherwise, setting z = 2(w-1)/(w+1),
*
* log(w) = z + z^3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1, 8 100000 1.9e-34 4.3e-35
*/
#include <math.h>
#include "math_private.h"
/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
* 1/sqrt(2) <= 1+x < sqrt(2)
* Theoretical peak relative error = 5.3e-37,
* relative peak error spread = 2.3e-14
*/
static const long double
P12 = 1.538612243596254322971797716843006400388E-6L,
P11 = 4.998469661968096229986658302195402690910E-1L,
P10 = 2.321125933898420063925789532045674660756E1L,
P9 = 4.114517881637811823002128927449878962058E2L,
P8 = 3.824952356185897735160588078446136783779E3L,
P7 = 2.128857716871515081352991964243375186031E4L,
P6 = 7.594356839258970405033155585486712125861E4L,
P5 = 1.797628303815655343403735250238293741397E5L,
P4 = 2.854829159639697837788887080758954924001E5L,
P3 = 3.007007295140399532324943111654767187848E5L,
P2 = 2.014652742082537582487669938141683759923E5L,
P1 = 7.771154681358524243729929227226708890930E4L,
P0 = 1.313572404063446165910279910527789794488E4L,
/* Q12 = 1.000000000000000000000000000000000000000E0L, */
Q11 = 4.839208193348159620282142911143429644326E1L,
Q10 = 9.104928120962988414618126155557301584078E2L,
Q9 = 9.147150349299596453976674231612674085381E3L,
Q8 = 5.605842085972455027590989944010492125825E4L,
Q7 = 2.248234257620569139969141618556349415120E5L,
Q6 = 6.132189329546557743179177159925690841200E5L,
Q5 = 1.158019977462989115839826904108208787040E6L,
Q4 = 1.514882452993549494932585972882995548426E6L,
Q3 = 1.347518538384329112529391120390701166528E6L,
Q2 = 7.777690340007566932935753241556479363645E5L,
Q1 = 2.626900195321832660448791748036714883242E5L,
Q0 = 3.940717212190338497730839731583397586124E4L;
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 1.1e-35,
* relative peak error spread 1.1e-9
*/
static const long double
R5 = -8.828896441624934385266096344596648080902E-1L,
R4 = 8.057002716646055371965756206836056074715E1L,
R3 = -2.024301798136027039250415126250455056397E3L,
R2 = 2.048819892795278657810231591630928516206E4L,
R1 = -8.977257995689735303686582344659576526998E4L,
R0 = 1.418134209872192732479751274970992665513E5L,
/* S6 = 1.000000000000000000000000000000000000000E0L, */
S5 = -1.186359407982897997337150403816839480438E2L,
S4 = 3.998526750980007367835804959888064681098E3L,
S3 = -5.748542087379434595104154610899551484314E4L,
S2 = 4.001557694070773974936904547424676279307E5L,
S1 = -1.332535117259762928288745111081235577029E6L,
S0 = 1.701761051846631278975701529965589676574E6L;
/* C1 + C2 = ln 2 */
static const long double C1 = 6.93145751953125E-1L;
static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
static const long double sqrth = 0.7071067811865475244008443621048490392848L;
/* ln (2^16384 * (1 - 2^-113)) */
static const long double zero = 0.0L;
long double
log1pl(long double xm1)
{
long double x, y, z, r, s;
ieee_quad_shape_type u;
int32_t hx;
int e;
/* Test for NaN or infinity input. */
u.value = xm1;
hx = u.parts32.mswhi;
if (hx >= 0x7fff0000)
return xm1;
/* log1p(+- 0) = +- 0. */
if (((hx & 0x7fffffff) == 0)
&& (u.parts32.mswlo | u.parts32.lswhi | u.parts32.lswlo) == 0)
return xm1;
x = xm1 + 1.0L;
/* log1p(-1) = -inf */
if (x <= 0.0L)
{
if (x == 0.0L)
return (-1.0L / (x - x));
else
return (zero / (x - x));
}
/* Separate mantissa from exponent. */
/* Use frexp used so that denormal numbers will be handled properly. */
x = frexpl (x, &e);
/* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
where z = 2(x-1)/x+1). */
if ((e > 2) || (e < -2))
{
if (x < sqrth)
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x * x;
r = ((((R5 * z
+ R4) * z
+ R3) * z
+ R2) * z
+ R1) * z
+ R0;
s = (((((z
+ S5) * z
+ S4) * z
+ S3) * z
+ S2) * z
+ S1) * z
+ S0;
z = x * (z * r / s);
z = z + e * C2;
z = z + x;
z = z + e * C1;
return (z);
}
/* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
if (x < sqrth)
{
e -= 1;
if (e != 0)
x = 2.0L * x - 1.0L; /* 2x - 1 */
else
x = xm1;
}
else
{
if (e != 0)
x = x - 1.0L;
else
x = xm1;
}
z = x * x;
r = (((((((((((P12 * x
+ P11) * x
+ P10) * x
+ P9) * x
+ P8) * x
+ P7) * x
+ P6) * x
+ P5) * x
+ P4) * x
+ P3) * x
+ P2) * x
+ P1) * x
+ P0;
s = (((((((((((x
+ Q11) * x
+ Q10) * x
+ Q9) * x
+ Q8) * x
+ Q7) * x
+ Q6) * x
+ Q5) * x
+ Q4) * x
+ Q3) * x
+ Q2) * x
+ Q1) * x
+ Q0;
y = x * (z * r / s);
y = y + e * C2;
z = y - 0.5L * z;
z = z + x;
z = z + e * C1;
return (z);
}

73
ld128/s_modfl.c Normal file
View file

@ -0,0 +1,73 @@
/* @(#)s_modf.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* modfl(long double x, long double *iptr)
* return fraction part of x, and return x's integral part in *iptr.
* Method:
* Bit twiddling.
*
* Exception:
* No exception.
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0;
long double
modfl(long double x, long double *iptr)
{
int64_t i0,i1,jj0;
u_int64_t i;
GET_LDOUBLE_WORDS64(i0,i1,x);
jj0 = ((i0>>48)&0x7fff)-0x3fff; /* exponent of x */
if(jj0<48) { /* integer part in high x */
if(jj0<0) { /* |x|<1 */
/* *iptr = +-0 */
SET_LDOUBLE_WORDS64(*iptr,i0&0x8000000000000000ULL,0);
return x;
} else {
i = (0x0000ffffffffffffLL)>>jj0;
if(((i0&i)|i1)==0) { /* x is integral */
*iptr = x;
/* return +-0 */
SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
return x;
} else {
SET_LDOUBLE_WORDS64(*iptr,i0&(~i),0);
return x - *iptr;
}
}
} else if (jj0>111) { /* no fraction part */
*iptr = x*one;
/* We must handle NaNs separately. */
if (jj0 == 0x4000 && ((i0 & 0x0000ffffffffffffLL) | i1))
return x*one;
/* return +-0 */
SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
return x;
} else { /* fraction part in low x */
i = -1ULL>>(jj0-48);
if((i1&i)==0) { /* x is integral */
*iptr = x;
/* return +-0 */
SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
return x;
} else {
SET_LDOUBLE_WORDS64(*iptr,i0,i1&(~i));
return x - *iptr;
}
}
}

72
ld128/s_nextafterl.c Normal file
View file

@ -0,0 +1,72 @@
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* IEEE functions
* nextafterl(x,y)
* return the next machine floating-point number of x in the
* direction toward y.
* Special cases:
*/
#include <math.h>
#include "math_private.h"
long double
nextafterl(long double x, long double y)
{
int64_t hx,hy,ix,iy;
u_int64_t lx,ly;
GET_LDOUBLE_WORDS64(hx,lx,x);
GET_LDOUBLE_WORDS64(hy,ly,y);
ix = hx&0x7fffffffffffffffLL; /* |x| */
iy = hy&0x7fffffffffffffffLL; /* |y| */
if(((ix>=0x7fff000000000000LL)&&((ix-0x7fff000000000000LL)|lx)!=0) || /* x is nan */
((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0)) /* y is nan */
return x+y;
if(x==y) return y; /* x=y, return y */
if((ix|lx)==0) { /* x == 0 */
volatile long double u;
SET_LDOUBLE_WORDS64(x,hy&0x8000000000000000ULL,1);/* return +-minsubnormal */
u = x;
u = u * u; /* raise underflow flag */
return x;
}
if(hx>=0) { /* x > 0 */
if(hx>hy||((hx==hy)&&(lx>ly))) { /* x > y, x -= ulp */
if(lx==0) hx--;
lx--;
} else { /* x < y, x += ulp */
lx++;
if(lx==0) hx++;
}
} else { /* x < 0 */
if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */
if(lx==0) hx--;
lx--;
} else { /* x > y, x += ulp */
lx++;
if(lx==0) hx++;
}
}
hy = hx&0x7fff000000000000LL;
if(hy==0x7fff000000000000LL) return x+x;/* overflow */
if(hy==0) {
volatile long double u = x*x; /* underflow */
}
SET_LDOUBLE_WORDS64(x,hx,lx);
return x;
}
__strong_alias(nexttowardl, nextafterl);

85
ld128/s_nexttoward.c Normal file
View file

@ -0,0 +1,85 @@
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* IEEE functions
* nexttoward(x,y)
* return the next machine floating-point number of x in the
* direction toward y.
* Special cases:
*/
#include <math.h>
#include <float.h>
#include "math_private.h"
double
nexttoward(double x, long double y)
{
int32_t hx,ix;
int64_t hy,iy;
u_int32_t lx;
u_int64_t ly;
EXTRACT_WORDS(hx,lx,x);
GET_LDOUBLE_WORDS64(hy,ly,y);
ix = hx&0x7fffffff; /* |x| */
iy = hy&0x7fffffffffffffffLL; /* |y| */
if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */
((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0))
/* y is nan */
return x+y;
if((long double) x==y) return y; /* x=y, return y */
if((ix|lx)==0) { /* x == 0 */
volatile double u;
INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */
u = x;
u = u * u; /* raise underflow flag */
return x;
}
if(hx>=0) { /* x > 0 */
if (hy<0||(ix>>20)>(iy>>48)-0x3c00
|| ((ix>>20)==(iy>>48)-0x3c00
&& (((((int64_t)hx)<<28)|(lx>>4))>(hy&0x0000ffffffffffffLL)
|| (((((int64_t)hx)<<28)|(lx>>4))==(hy&0x0000ffffffffffffLL)
&& (lx&0xf)>(ly>>60))))) { /* x > y, x -= ulp */
if(lx==0) hx -= 1;
lx -= 1;
} else { /* x < y, x += ulp */
lx += 1;
if(lx==0) hx += 1;
}
} else { /* x < 0 */
if (hy>=0||(ix>>20)>(iy>>48)-0x3c00
|| ((ix>>20)==(iy>>48)-0x3c00
&& (((((int64_t)hx)<<28)|(lx>>4))>(hy&0x0000ffffffffffffLL)
|| (((((int64_t)hx)<<28)|(lx>>4))==(hy&0x0000ffffffffffffLL)
&& (lx&0xf)>(ly>>60))))) { /* x < y, x -= ulp */
if(lx==0) hx -= 1;
lx -= 1;
} else { /* x > y, x += ulp */
lx += 1;
if(lx==0) hx += 1;
}
}
hy = hx&0x7ff00000;
if(hy>=0x7ff00000) {
x = x+x; /* overflow */
return x;
}
if(hy<0x00100000) {
volatile double u = x*x; /* underflow */
}
INSERT_WORDS(x,hx,lx);
return x;
}

65
ld128/s_nexttowardf.c Normal file
View file

@ -0,0 +1,65 @@
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
#include <math.h>
#include "math_private.h"
float
nexttowardf(float x, long double y)
{
int32_t hx,ix;
int64_t hy,iy;
u_int64_t ly;
GET_FLOAT_WORD(hx,x);
GET_LDOUBLE_WORDS64(hy,ly,y);
ix = hx&0x7fffffff; /* |x| */
iy = hy&0x7fffffffffffffffLL; /* |y| */
if((ix>0x7f800000) || /* x is nan */
((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0))
/* y is nan */
return x+y;
if((long double) x==y) return y; /* x=y, return y */
if(ix==0) { /* x == 0 */
volatile float u;
SET_FLOAT_WORD(x,(u_int32_t)((hy>>32)&0x80000000)|1);/* return +-minsub*/
u = x;
u = u * u; /* raise underflow flag */
return x;
}
if(hx>=0) { /* x > 0 */
if(hy<0||(ix>>23)>(iy>>48)-0x3f80
|| ((ix>>23)==(iy>>48)-0x3f80
&& (ix&0x7fffff)>((hy>>25)&0x7fffff))) {/* x > y, x -= ulp */
hx -= 1;
} else { /* x < y, x += ulp */
hx += 1;
}
} else { /* x < 0 */
if(hy>=0||(ix>>23)>(iy>>48)-0x3f80
|| ((ix>>23)==(iy>>48)-0x3f80
&& (ix&0x7fffff)>((hy>>25)&0x7fffff))) {/* x < y, x -= ulp */
hx -= 1;
} else { /* x > y, x += ulp */
hx += 1;
}
}
hy = hx&0x7f800000;
if(hy>=0x7f800000) return x+x; /* overflow */
if(hy<0x00800000) {
volatile float u = x*x; /* underflow */
}
SET_FLOAT_WORD(x,hx);
return x;
}

168
ld128/s_remquol.c Normal file
View file

@ -0,0 +1,168 @@
/* @(#)e_fmod.c 1.3 95/01/18 */
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/types.h>
#include <machine/ieee.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include "math_private.h"
#define BIAS (LDBL_MAX_EXP - 1)
/*
* These macros add and remove an explicit integer bit in front of the
* fractional mantissa, if the architecture doesn't have such a bit by
* default already.
*/
#ifdef LDBL_IMPLICIT_NBIT
#define LDBL_NBIT 0
#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
#define HFRAC_BITS (EXT_FRACHBITS + EXT_FRACHMBITS)
#else
#define LDBL_NBIT 0x80000000
#define SET_NBIT(hx) (hx)
#define HFRAC_BITS (EXT_FRACHBITS + EXT_FRACHMBITS - 1)
#endif
#define MANL_SHIFT (EXT_FRACLMBITS + EXT_FRACLBITS - 1)
static const long double Zero[] = {0.0L, -0.0L};
/*
* Return the IEEE remainder and set *quo to the last n bits of the
* quotient, rounded to the nearest integer. We choose n=31 because
* we wind up computing all the integer bits of the quotient anyway as
* a side-effect of computing the remainder by the shift and subtract
* method. In practice, this is far more bits than are needed to use
* remquo in reduction algorithms.
*
* Assumptions:
* - The low part of the mantissa fits in a manl_t exactly.
* - The high part of the mantissa fits in an int64_t with enough room
* for an explicit integer bit in front of the fractional bits.
*/
long double
remquol(long double x, long double y, int *quo)
{
int64_t hx,hz,hy,_hx;
uint64_t lx,ly,lz;
uint64_t sx,sxy;
int ix,iy,n,q;
GET_LDOUBLE_WORDS64(hx,lx,x);
GET_LDOUBLE_WORDS64(hy,ly,y);
sx = (hx>>48)&0x8000;
sxy = sx ^ ((hy>>48)&0x8000);
hx &= 0x7fffffffffffffffLL; /* |x| */
hy &= 0x7fffffffffffffffLL; /* |y| */
SET_LDOUBLE_WORDS64(x,hx,lx);
SET_LDOUBLE_WORDS64(y,hy,ly);
/* purge off exception values */
if((hy|ly)==0 || /* y=0 */
((hx>>48) == BIAS + LDBL_MAX_EXP) || /* or x not finite */
((hy>>48) == BIAS + LDBL_MAX_EXP &&
(((hy&0x0000ffffffffffffLL)&~LDBL_NBIT)|ly)!=0)) /* or y is NaN */
return (x*y)/(x*y);
if((hx>>48)<=(hy>>48)) {
if(((hx>>48)<(hy>>48)) ||
((hx&0x0000ffffffffffffLL)<=(hy&0x0000ffffffffffffLL) &&
((hx&0x0000ffffffffffffLL)<(hy&0x0000ffffffffffffLL) ||
lx<ly))) {
q = 0;
goto fixup; /* |x|<|y| return x or x-y */
}
if((hx&0x0000ffffffffffffLL)==(hy&0x0000ffffffffffffLL) &&
lx==ly) {
*quo = 1;
return Zero[sx!=0]; /* |x|=|y| return x*0*/
}
}
/* determine ix = ilogb(x) */
if((hx>>48) == 0) { /* subnormal x */
x *= 0x1.0p512;
GET_LDOUBLE_WORDS64(hx,lx,x);
ix = (hx>>48) - (BIAS + 512);
} else {
ix = (hx>>48) - BIAS;
}
/* determine iy = ilogb(y) */
if((hy>>48) == 0) { /* subnormal y */
y *= 0x1.0p512;
GET_LDOUBLE_WORDS64(hy,ly,y);
iy = (hy>>48) - (BIAS + 512);
} else {
iy = (hy>>48) - BIAS;
}
/* set up {hx,lx}, {hy,ly} and align y to x */
_hx = SET_NBIT(hx) & 0x0000ffffffffffffLL;
hy = SET_NBIT(hy);
/* fix point fmod */
n = ix - iy;
q = 0;
while(n--) {
hz=_hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){_hx = _hx+_hx+(lx>>MANL_SHIFT); lx = lx+lx;}
else {_hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; q++;}
q <<= 1;
}
hz=_hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {_hx=hz;lx=lz;q++;}
/* convert back to floating value and restore the sign */
if((_hx|lx)==0) { /* return sign(x)*0 */
*quo = (sxy ? -q : q);
return Zero[sx!=0];
}
while(_hx<(1ULL<<HFRAC_BITS)) { /* normalize x */
_hx = _hx+_hx+(lx>>MANL_SHIFT); lx = lx+lx;
iy -= 1;
}
hx = (hx&0xffff000000000000LL) | (_hx&0x0000ffffffffffffLL);
if (iy < LDBL_MIN_EXP) {
hx = (hx&0x0000ffffffffffffLL) | (uint64_t)(iy + BIAS + 512)<<48;
SET_LDOUBLE_WORDS64(x,hx,lx);
x *= 0x1p-512;
GET_LDOUBLE_WORDS64(hx,lx,x);
} else {
hx = (hx&0x0000ffffffffffffLL) | (uint64_t)(iy + BIAS)<<48;
}
hx &= 0x7fffffffffffffffLL;
SET_LDOUBLE_WORDS64(x,hx,lx);
fixup:
y = fabsl(y);
if (y < LDBL_MIN * 2) {
if (x+x>y || (x+x==y && (q & 1))) {
q++;
x-=y;
}
} else if (x>0.5*y || (x==0.5*y && (q & 1))) {
q++;
x-=y;
}
GET_LDOUBLE_MSW64(hx,x);
hx ^= sx;
SET_LDOUBLE_MSW64(x,hx);
q &= 0x7fffffff;
*quo = (sxy ? -q : q);
return x;
}

104
ld128/s_tanhl.c Normal file
View file

@ -0,0 +1,104 @@
/* @(#)s_tanh.c 5.1 93/09/24 */
/*
* ====================================================
* 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) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* tanhl(x)
* Return the Hyperbolic Tangent of x
*
* Method :
* x -x
* e - e
* 0. tanhl(x) is defined to be -----------
* x -x
* e + e
* 1. reduce x to non-negative by tanhl(-x) = -tanhl(x).
* 2. 0 <= x <= 2**-57 : tanhl(x) := x*(one+x)
* -t
* 2**-57 < x <= 1 : tanhl(x) := -----; t = expm1l(-2x)
* t + 2
* 2
* 1 <= x <= 40.0 : tanhl(x) := 1- ----- ; t=expm1l(2x)
* t + 2
* 40.0 < x <= INF : tanhl(x) := 1.
*
* Special cases:
* tanhl(NaN) is NaN;
* only tanhl(0)=0 is exact for finite argument.
*/
#include "math.h"
#include "math_private.h"
static const long double one = 1.0, two = 2.0, tiny = 1.0e-4900L;
long double
tanhl(long double x)
{
long double t, z;
u_int32_t jx, ix;
ieee_quad_shape_type u;
/* Words of |x|. */
u.value = x;
jx = u.parts32.mswhi;
ix = jx & 0x7fffffff;
/* x is INF or NaN */
if (ix >= 0x7fff0000)
{
/* for NaN it's not important which branch: tanhl(NaN) = NaN */
if (jx & 0x80000000)
return one / x - one; /* tanhl(-inf)= -1; */
else
return one / x + one; /* tanhl(+inf)=+1 */
}
/* |x| < 40 */
if (ix < 0x40044000)
{
if (u.value == 0)
return x; /* x == +- 0 */
if (ix < 0x3fc60000) /* |x| < 2^-57 */
return x * (one + tiny); /* tanh(small) = small */
u.parts32.mswhi = ix; /* Absolute value of x. */
if (ix >= 0x3fff0000)
{ /* |x| >= 1 */
t = expm1l (two * u.value);
z = one - two / (t + two);
}
else
{
t = expm1l (-two * u.value);
z = -t / (t + two);
}
/* |x| > 40, return +-1 */
}
else
{
z = one - tiny; /* raised inexact flag */
}
return (jx & 0x80000000) ? -z : z;
}

72
ld128/s_truncl.c Normal file
View file

@ -0,0 +1,72 @@
/*
* ====================================================
* 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.
* ====================================================
*
* From: @(#)s_floor.c 5.1 93/09/24
*/
/*
* truncl(x)
* Return x rounded toward 0 to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to truncl(x).
*/
#include <sys/types.h>
#include <machine/ieee.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include "math_private.h"
#ifdef LDBL_IMPLICIT_NBIT
#define MANH_SIZE (EXT_FRACHBITS + EXT_FRACHMBITS + 1)
#else
#define MANH_SIZE (EXT_FRACHBITS + EXT_FRACHMBITS)
#endif
static const long double huge = 1.0e300;
static const float zero[] = { 0.0, -0.0 };
long double
truncl(long double x)
{
int e;
int64_t ix0, ix1;
GET_LDOUBLE_WORDS64(ix0,ix1,x);
e = ((ix0>>48)&0x7fff) - LDBL_MAX_EXP + 1;
if (e < MANH_SIZE - 1) {
if (e < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0)
return (zero[((ix0>>48)&0x8000)!=0]);
} else {
uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
if (((ix0 & m) | ix1) == 0)
return (x); /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
ix0 &= ~m;
ix1 = 0;
}
}
} else if (e < LDBL_MANT_DIG - 1) {
uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
if ((ix1 & m) == 0)
return (x); /* x is integral */
if (huge + x > 0.0) /* raise inexact flag */
ix1 &= ~m;
}
SET_LDOUBLE_WORDS64(x,ix0,ix1);
return (x);
}

11
ld80/Make.files
View file

@ -1,6 +1,11 @@
$(CUR_SRCS) += invtrig.c k_cosl.c k_sinl.c\
k_tanl.c s_exp2l.c
$(CUR_SRCS) += invtrig.c \
e_acoshl.c e_hypotl.c e_powl.c k_tanl.c s_exp2l.c s_nanl.c \
e_atanhl.c e_lgammal.c e_sinhl.c s_asinhl.c s_expm1l.c \
e_coshl.c e_log10l.c e_tgammal.c s_floorl.c s_nextafterl.c \
e_expl.c e_log2l.c k_cosl.c s_ceill.c s_log1pl.c s_tanhl.c \
e_fmodl.c e_logl.c k_sinl.c s_erfl.c s_modfl.c s_truncl.c
# s_remquol.c
ifneq ($(OS), WINNT)
$(CUR_SRCS) += s_nanl.c
endif
endif

57
ld80/e_acoshl.c Normal file
View file

@ -0,0 +1,57 @@
/* @(#)e_acosh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* acoshl(x)
* Method :
* Based on
* acoshl(x) = logl [ x + sqrtl(x*x-1) ]
* we have
* acoshl(x) := logl(x)+ln2, if x is large; else
* acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
* acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
*
* Special cases:
* acoshl(x) is NaN with signal if x<1.
* acoshl(NaN) is NaN without signal.
*/
#include <math.h>
#include "math_private.h"
static const long double
one = 1.0,
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
long double
acoshl(long double x)
{
long double t;
u_int32_t se,i0,i1;
GET_LDOUBLE_WORDS(se,i0,i1,x);
if(se<0x3fff || se & 0x8000) { /* x < 1 */
return (x-x)/(x-x);
} else if(se >=0x401d) { /* x > 2**30 */
if(se >=0x7fff) { /* x is inf of NaN */
return x+x;
} else
return logl(x)+ln2; /* acoshl(huge)=logl(2x) */
} else if(((se-0x3fff)|i0|i1)==0) {
return 0.0; /* acosh(1) = 0 */
} else if (se > 0x4000) { /* 2**28 > x > 2 */
t=x*x;
return logl(2.0*x-one/(x+sqrtl(t-one)));
} else { /* 1<x<2 */
t = x-one;
return log1pl(t+sqrtl(2.0*t+t*t));
}
}

60
ld80/e_atanhl.c Normal file
View file

@ -0,0 +1,60 @@
/* @(#)e_atanh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* atanhl(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
* 1 2x x
* atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x
*
* For x<0.5
* atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
*
* Special cases:
* atanhl(x) is NaN if |x| > 1 with signal;
* atanhl(NaN) is that NaN with no signal;
* atanhl(+-1) is +-INF with signal.
*
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0, huge = 1e4900L;
static const long double zero = 0.0;
long double
atanhl(long double x)
{
long double t;
int32_t ix;
u_int32_t se,i0,i1;
GET_LDOUBLE_WORDS(se,i0,i1,x);
ix = se&0x7fff;
if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff)
/* |x|>1 */
return (x-x)/(x-x);
if(ix==0x3fff)
return x/zero;
if(ix<0x3fe3&&(huge+x)>zero) return x; /* x<2**-28 */
SET_LDOUBLE_EXP(x,ix);
if(ix<0x3ffe) { /* x < 0.5 */
t = x+x;
t = 0.5*log1pl(t+t*x/(one-x));
} else
t = 0.5*log1pl((x+x)/(one-x));
if(se<=0x7fff) return t; else return -t;
}

82
ld80/e_coshl.c Normal file
View file

@ -0,0 +1,82 @@
/* @(#)e_cosh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* coshl(x)
* Method :
* mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
* 1. Replace x by |x| (coshl(x) = coshl(-x)).
* 2.
* [ exp(x) - 1 ]^2
* 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
* 2*exp(x)
*
* exp(x) + 1/exp(x)
* ln2/2 <= x <= 22 : coshl(x) := -------------------
* 2
* 22 <= x <= lnovft : coshl(x) := expl(x)/2
* lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2)
* ln2ovft < x : coshl(x) := huge*huge (overflow)
*
* Special cases:
* coshl(x) is |x| if x is +INF, -INF, or NaN.
* only coshl(0)=1 is exact for finite x.
*/
#include "math.h"
#include "math_private.h"
static const long double one = 1.0, half=0.5, huge = 1.0e4900L;
long double
coshl(long double x)
{
long double t,w;
int32_t ex;
u_int32_t mx,lx;
/* High word of |x|. */
GET_LDOUBLE_WORDS(ex,mx,lx,x);
ex &= 0x7fff;
/* x is INF or NaN */
if(ex==0x7fff) return x*x;
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
t = expm1l(fabsl(x));
w = one+t;
if (ex<0x3fbc) return w; /* cosh(tiny) = 1 */
return one+(t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
t = expl(fabsl(x));
return half*t+half/t;
}
/* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
return half*expl(fabsl(x));
/* |x| in [log(maxdouble), log(2*maxdouble)) */
if (ex == 0x400c && (mx < 0xb174ddc0u
|| (mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
{
w = expl(half*fabsl(x));
t = half*w;
return t*w;
}
/* |x| >= log(2*maxdouble), cosh(x) overflow */
return huge*huge;
}

131
ld80/e_expl.c Normal file
View file

@ -0,0 +1,131 @@
/* $OpenBSD: e_expl.c,v 1.3 2013/11/12 20:35:19 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* expl.c
*
* Exponential function, long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, expl();
*
* y = expl( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* A Pade' form of degree 2/3 is used to approximate exp(f) - 1
* in the basic range [-0.5 ln 2, 0.5 ln 2].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +-10000 50000 1.12e-19 2.81e-20
*
*
* Error amplification in the exponential function can be
* a serious matter. The error propagation involves
* exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
* which shows that a 1 lsb error in representing X produces
* a relative error of X times 1 lsb in the function.
* While the routine gives an accurate result for arguments
* that are exactly represented by a long double precision
* computer number, the result contains amplified roundoff
* error for large arguments not exactly represented.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* exp underflow x < MINLOG 0.0
* exp overflow x > MAXLOG MAXNUM
*
*/
/* Exponential function */
#include <math.h>
#include "math_private.h"
static long double P[3] = {
1.2617719307481059087798E-4L,
3.0299440770744196129956E-2L,
9.9999999999999999991025E-1L,
};
static long double Q[4] = {
3.0019850513866445504159E-6L,
2.5244834034968410419224E-3L,
2.2726554820815502876593E-1L,
2.0000000000000000000897E0L,
};
static const long double C1 = 6.9314575195312500000000E-1L;
static const long double C2 = 1.4286068203094172321215E-6L;
static const long double MAXLOGL = 1.1356523406294143949492E4L;
static const long double MINLOGL = -1.13994985314888605586758E4L;
static const long double LOG2EL = 1.4426950408889634073599E0L;
long double
expl(long double x)
{
long double px, xx;
int n;
if( isnan(x) )
return(x);
if( x > MAXLOGL)
return( INFINITY );
if( x < MINLOGL )
return(0.0L);
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */
n = px;
x -= px * C1;
x -= px * C2;
/* rational approximation for exponential
* of the fractional part:
* e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*/
xx = x * x;
px = x * __polevll( xx, P, 2 );
x = px/( __polevll( xx, Q, 3 ) - px );
x = 1.0L + ldexpl( x, 1 );
x = ldexpl( x, n );
return(x);
}

142
ld80/e_fmodl.c Normal file
View file

@ -0,0 +1,142 @@
/* @(#)e_fmod.c 1.3 95/01/18 */
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/types.h>
//#include <machine/ieee.h>
#include <float.h>
#include <openlibm.h>
#include <stdint.h>
#include "math_private.h"
#define BIAS (LDBL_MAX_EXP - 1)
/*
* These macros add and remove an explicit integer bit in front of the
* fractional mantissa, if the architecture doesn't have such a bit by
* default already.
*/
#ifdef LDBL_IMPLICIT_NBIT
#define LDBL_NBIT 0
#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
#define HFRAC_BITS EXT_FRACHBITS
#else
#define LDBL_NBIT 0x80000000
#define SET_NBIT(hx) (hx)
#define HFRAC_BITS (EXT_FRACHBITS - 1)
#endif
#define MANL_SHIFT (EXT_FRACLBITS - 1)
static const long double one = 1.0, Zero[] = {0.0, -0.0,};
/*
* fmodl(x,y)
* Return x mod y in exact arithmetic
* Method: shift and subtract
*
* Assumptions:
* - The low part of the mantissa fits in a manl_t exactly.
* - The high part of the mantissa fits in an int64_t with enough room
* for an explicit integer bit in front of the fractional bits.
*/
long double
fmodl(long double x, long double y)
{
union {
long double e;
struct ieee_ext bits;
} ux, uy;
int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */
uint32_t hy;
uint32_t lx,ly,lz;
int ix,iy,n,sx;
ux.e = x;
uy.e = y;
sx = ux.bits.ext_sign;
/* purge off exception values */
if((uy.bits.ext_exp|uy.bits.ext_frach|uy.bits.ext_fracl)==0 || /* y=0 */
(ux.bits.ext_exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
(uy.bits.ext_exp == BIAS + LDBL_MAX_EXP &&
((uy.bits.ext_frach&~LDBL_NBIT)|uy.bits.ext_fracl)!=0)) /* or y is NaN */
return (x*y)/(x*y);
if(ux.bits.ext_exp<=uy.bits.ext_exp) {
if((ux.bits.ext_exp<uy.bits.ext_exp) ||
(ux.bits.ext_frach<=uy.bits.ext_frach &&
(ux.bits.ext_frach<uy.bits.ext_frach ||
ux.bits.ext_fracl<uy.bits.ext_fracl))) {
return x; /* |x|<|y| return x or x-y */
}
if(ux.bits.ext_frach==uy.bits.ext_frach &&
ux.bits.ext_fracl==uy.bits.ext_fracl) {
return Zero[sx]; /* |x|=|y| return x*0*/
}
}
/* determine ix = ilogb(x) */
if(ux.bits.ext_exp == 0) { /* subnormal x */
ux.e *= 0x1.0p512;
ix = ux.bits.ext_exp - (BIAS + 512);
} else {
ix = ux.bits.ext_exp - BIAS;
}
/* determine iy = ilogb(y) */
if(uy.bits.ext_exp == 0) { /* subnormal y */
uy.e *= 0x1.0p512;
iy = uy.bits.ext_exp - (BIAS + 512);
} else {
iy = uy.bits.ext_exp - BIAS;
}
/* set up {hx,lx}, {hy,ly} and align y to x */
hx = SET_NBIT(ux.bits.ext_frach);
hy = SET_NBIT(uy.bits.ext_frach);
lx = ux.bits.ext_fracl;
ly = uy.bits.ext_fracl;
/* fix point fmod */
n = ix - iy;
while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
else {
if ((hz|lz)==0) /* return sign(x)*0 */
return Zero[sx];
hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz;
}
}
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;}
/* convert back to floating value and restore the sign */
if((hx|lx)==0) /* return sign(x)*0 */
return Zero[sx];
while(hx<(1ULL<<HFRAC_BITS)) { /* normalize x */
hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
iy -= 1;
}
ux.bits.ext_frach = hx; /* The mantissa is truncated here if needed. */
ux.bits.ext_fracl = lx;
if (iy < LDBL_MIN_EXP) {
ux.bits.ext_exp = iy + (BIAS + 512);
ux.e *= 0x1p-512;
} else {
ux.bits.ext_exp = iy + BIAS;
}
x = ux.e * one; /* create necessary signal */
return x; /* exact output */
}

122
ld80/e_hypotl.c Normal file
View file

@ -0,0 +1,122 @@
/* @(#)e_hypot.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* hypotl(x,y)
*
* Method :
* If (assume round-to-nearest) z=x*x+y*y
* has error less than sqrt(2)/2 ulp, than
* sqrt(z) has error less than 1 ulp (exercise).
*
* So, compute sqrt(x*x+y*y) with some care as
* follows to get the error below 1 ulp:
*
* Assume x>y>0;
* (if possible, set rounding to round-to-nearest)
* 1. if x > 2y use
* x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
* where x1 = x with lower 32 bits cleared, x2 = x-x1; else
* 2. if x <= 2y use
* t1*yy1+((x-y)*(x-y)+(t1*y2+t2*y))
* where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
* yy1= y with lower 32 bits chopped, y2 = y-yy1.
*
* NOTE: scaling may be necessary if some argument is too
* large or too tiny
*
* Special cases:
* hypot(x,y) is INF if x or y is +INF or -INF; else
* hypot(x,y) is NAN if x or y is NAN.
*
* Accuracy:
* hypot(x,y) returns sqrt(x^2+y^2) with error less
* than 1 ulps (units in the last place)
*/
#include <math.h>
#include "math_private.h"
long double
hypotl(long double x, long double y)
{
long double a,b,t1,t2,yy1,y2,w;
u_int32_t j,k,ea,eb;
GET_LDOUBLE_EXP(ea,x);
ea &= 0x7fff;
GET_LDOUBLE_EXP(eb,y);
eb &= 0x7fff;
if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;}
SET_LDOUBLE_EXP(a,ea); /* a <- |a| */
SET_LDOUBLE_EXP(b,eb); /* b <- |b| */
if((ea-eb)>0x46) {return a+b;} /* x/y > 2**70 */
k=0;
if(ea > 0x5f3f) { /* a>2**8000 */
if(ea == 0x7fff) { /* Inf or NaN */
u_int32_t es,high,low;
w = a+b; /* for sNaN */
GET_LDOUBLE_WORDS(es,high,low,a);
if(((high&0x7fffffff)|low)==0) w = a;
GET_LDOUBLE_WORDS(es,high,low,b);
if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b;
return w;
}
/* scale a and b by 2**-9600 */
ea -= 0x2580; eb -= 0x2580; k += 9600;
SET_LDOUBLE_EXP(a,ea);
SET_LDOUBLE_EXP(b,eb);
}
if(eb < 0x20bf) { /* b < 2**-8000 */
if(eb == 0) { /* subnormal b or 0 */
u_int32_t es,high,low;
GET_LDOUBLE_WORDS(es,high,low,b);
if((high|low)==0) return a;
SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0); /* t1=2^16382 */
b *= t1;
a *= t1;
k -= 16382;
} else { /* scale a and b by 2^9600 */
ea += 0x2580; /* a *= 2^9600 */
eb += 0x2580; /* b *= 2^9600 */
k -= 9600;
SET_LDOUBLE_EXP(a,ea);
SET_LDOUBLE_EXP(b,eb);
}
}
/* medium size a and b */
w = a-b;
if (w>b) {
u_int32_t high;
GET_LDOUBLE_MSW(high,a);
SET_LDOUBLE_WORDS(t1,ea,high,0);
t2 = a-t1;
w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
u_int32_t high;
GET_LDOUBLE_MSW(high,b);
a = a+a;
SET_LDOUBLE_WORDS(yy1,eb,high,0);
y2 = b - yy1;
GET_LDOUBLE_MSW(high,a);
SET_LDOUBLE_WORDS(t1,ea+1,high,0);
t2 = a - t1;
w = sqrtl(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
u_int32_t es;
t1 = 1.0;
GET_LDOUBLE_EXP(es,t1);
SET_LDOUBLE_EXP(t1,es+k);
return t1*w;
} else return w;
}

425
ld80/e_lgammal.c Normal file
View file

@ -0,0 +1,425 @@
/*
* ====================================================
* 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) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* lgammal(x)
* Reentrant version of the logarithm of the Gamma function
* with user provide pointer for the sign of Gamma(x).
*
* Method:
* 1. Argument Reduction for 0 < x <= 8
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* for example,
* lgamma(7.3) = log(6.3) + lgamma(6.3)
* = log(6.3*5.3) + lgamma(5.3)
* = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
* 2. Polynomial approximation of lgamma around its
* minimun ymin=1.461632144968362245 to maintain monotonicity.
* On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
* Let z = x-ymin;
* lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
* 2. Rational approximation in the primary interval [2,3]
* We use the following approximation:
* s = x-2.0;
* lgamma(x) = 0.5*s + s*P(s)/Q(s)
* Our algorithms are based on the following observation
*
* zeta(2)-1 2 zeta(3)-1 3
* lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
* 2 3
*
* where Euler = 0.5771... is the Euler constant, which is very
* close to 0.5.
*
* 3. For x>=8, we have
* lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
* (better formula:
* lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
* Let z = 1/x, then we approximation
* f(z) = lgamma(x) - (x-0.5)(log(x)-1)
* by
* 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
*
* 4. For negative x, since (G is gamma function)
* -x*G(-x)*G(x) = pi/sin(pi*x),
* we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
* Hence, for x<0, signgam = sign(sin(pi*x)) and
* lgamma(x) = log(|Gamma(x)|)
* = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
* Note: one should avoid compute pi*(-x) directly in the
* computation of sin(pi*(-x)).
*
* 5. Special Cases
* lgamma(2+s) ~ s*(1-Euler) for tiny s
* lgamma(1)=lgamma(2)=0
* lgamma(x) ~ -log(x) for tiny x
* lgamma(0) = lgamma(inf) = inf
* lgamma(-integer) = +-inf
*
*/
#include <math.h>
#include "math_private.h"
static const long double
half = 0.5L,
one = 1.0L,
pi = 3.14159265358979323846264L,
two63 = 9.223372036854775808e18L,
/* lgam(1+x) = 0.5 x + x a(x)/b(x)
-0.268402099609375 <= x <= 0
peak relative error 6.6e-22 */
a0 = -6.343246574721079391729402781192128239938E2L,
a1 = 1.856560238672465796768677717168371401378E3L,
a2 = 2.404733102163746263689288466865843408429E3L,
a3 = 8.804188795790383497379532868917517596322E2L,
a4 = 1.135361354097447729740103745999661157426E2L,
a5 = 3.766956539107615557608581581190400021285E0L,
b0 = 8.214973713960928795704317259806842490498E3L,
b1 = 1.026343508841367384879065363925870888012E4L,
b2 = 4.553337477045763320522762343132210919277E3L,
b3 = 8.506975785032585797446253359230031874803E2L,
b4 = 6.042447899703295436820744186992189445813E1L,
/* b5 = 1.000000000000000000000000000000000000000E0 */
tc = 1.4616321449683623412626595423257213284682E0L,
tf = -1.2148629053584961146050602565082954242826E-1,/* double precision */
/* tt = (tail of tf), i.e. tf + tt has extended precision. */
tt = 3.3649914684731379602768989080467587736363E-18L,
/* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
-1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
/* lgam (x + tc) = tf + tt + x g(x)/h(x)
- 0.230003726999612341262659542325721328468 <= x
<= 0.2699962730003876587373404576742786715318
peak relative error 2.1e-21 */
g0 = 3.645529916721223331888305293534095553827E-18L,
g1 = 5.126654642791082497002594216163574795690E3L,
g2 = 8.828603575854624811911631336122070070327E3L,
g3 = 5.464186426932117031234820886525701595203E3L,
g4 = 1.455427403530884193180776558102868592293E3L,
g5 = 1.541735456969245924860307497029155838446E2L,
g6 = 4.335498275274822298341872707453445815118E0L,
h0 = 1.059584930106085509696730443974495979641E4L,
h1 = 2.147921653490043010629481226937850618860E4L,
h2 = 1.643014770044524804175197151958100656728E4L,
h3 = 5.869021995186925517228323497501767586078E3L,
h4 = 9.764244777714344488787381271643502742293E2L,
h5 = 6.442485441570592541741092969581997002349E1L,
/* h6 = 1.000000000000000000000000000000000000000E0 */
/* lgam (x+1) = -0.5 x + x u(x)/v(x)
-0.100006103515625 <= x <= 0.231639862060546875
peak relative error 1.3e-21 */
u0 = -8.886217500092090678492242071879342025627E1L,
u1 = 6.840109978129177639438792958320783599310E2L,
u2 = 2.042626104514127267855588786511809932433E3L,
u3 = 1.911723903442667422201651063009856064275E3L,
u4 = 7.447065275665887457628865263491667767695E2L,
u5 = 1.132256494121790736268471016493103952637E2L,
u6 = 4.484398885516614191003094714505960972894E0L,
v0 = 1.150830924194461522996462401210374632929E3L,
v1 = 3.399692260848747447377972081399737098610E3L,
v2 = 3.786631705644460255229513563657226008015E3L,
v3 = 1.966450123004478374557778781564114347876E3L,
v4 = 4.741359068914069299837355438370682773122E2L,
v5 = 4.508989649747184050907206782117647852364E1L,
/* v6 = 1.000000000000000000000000000000000000000E0 */
/* lgam (x+2) = .5 x + x s(x)/r(x)
0 <= x <= 1
peak relative error 7.2e-22 */
s0 = 1.454726263410661942989109455292824853344E6L,
s1 = -3.901428390086348447890408306153378922752E6L,
s2 = -6.573568698209374121847873064292963089438E6L,
s3 = -3.319055881485044417245964508099095984643E6L,
s4 = -7.094891568758439227560184618114707107977E5L,
s5 = -6.263426646464505837422314539808112478303E4L,
s6 = -1.684926520999477529949915657519454051529E3L,
r0 = -1.883978160734303518163008696712983134698E7L,
r1 = -2.815206082812062064902202753264922306830E7L,
r2 = -1.600245495251915899081846093343626358398E7L,
r3 = -4.310526301881305003489257052083370058799E6L,
r4 = -5.563807682263923279438235987186184968542E5L,
r5 = -3.027734654434169996032905158145259713083E4L,
r6 = -4.501995652861105629217250715790764371267E2L,
/* r6 = 1.000000000000000000000000000000000000000E0 */
/* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
x >= 8
Peak relative error 1.51e-21
w0 = LS2PI - 0.5 */
w0 = 4.189385332046727417803e-1L,
w1 = 8.333333333333331447505E-2L,
w2 = -2.777777777750349603440E-3L,
w3 = 7.936507795855070755671E-4L,
w4 = -5.952345851765688514613E-4L,
w5 = 8.412723297322498080632E-4L,
w6 = -1.880801938119376907179E-3L,
w7 = 4.885026142432270781165E-3L;
static const long double zero = 0.0L;
static long double
sin_pi(long double x)
{
long double y, z;
int n, ix;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffd8000) /* 0.25 */
return sinl (pi * x);
y = -x; /* x is assume negative */
/*
* argument reduction, make sure inexact flag not raised if input
* is an integer
*/
z = floorl (y);
if (z != y)
{ /* inexact anyway */
y *= 0.5;
y = 2.0*(y - floorl(y)); /* y = |x| mod 2.0 */
n = (int) (y*4.0);
}
else
{
if (ix >= 0x403f8000) /* 2^64 */
{
y = zero; n = 0; /* y must be even */
}
else
{
if (ix < 0x403e8000) /* 2^63 */
z = y + two63; /* exact */
GET_LDOUBLE_WORDS (se, i0, i1, z);
n = i1 & 1;
y = n;
n <<= 2;
}
}
switch (n)
{
case 0:
y = sinl (pi * y);
break;
case 1:
case 2:
y = cosl (pi * (half - y));
break;
case 3:
case 4:
y = sinl (pi * (one - y));
break;
case 5:
case 6:
y = -cosl (pi * (y - 1.5));
break;
default:
y = sinl (pi * (y - 2.0));
break;
}
return -y;
}
long double
lgammal(long double x)
{
long double t, y, z, nadj, p, p1, p2, q, r, w;
int i, ix;
u_int32_t se, i0, i1;
signgam = 1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if ((ix | i0 | i1) == 0)
{
if (se & 0x8000)
signgam = -1;
return one / fabsl (x);
}
ix = (ix << 16) | (i0 >> 16);
/* purge off +-inf, NaN, +-0, and negative arguments */
if (ix >= 0x7fff0000)
return x * x;
if (ix < 0x3fc08000) /* 2^-63 */
{ /* |x|<2**-63, return -log(|x|) */
if (se & 0x8000)
{
signgam = -1;
return -logl (-x);
}
else
return -logl (x);
}
if (se & 0x8000)
{
t = sin_pi (x);
if (t == zero)
return one / fabsl (t); /* -integer */
nadj = logl (pi / fabsl (t * x));
if (t < zero)
signgam = -1;
x = -x;
}
/* purge off 1 and 2 */
if ((((ix - 0x3fff8000) | i0 | i1) == 0)
|| (((ix - 0x40008000) | i0 | i1) == 0))
r = 0;
else if (ix < 0x40008000) /* 2.0 */
{
/* x < 2.0 */
if (ix <= 0x3ffee666) /* 8.99993896484375e-1 */
{
/* lgamma(x) = lgamma(x+1) - log(x) */
r = -logl (x);
if (ix >= 0x3ffebb4a) /* 7.31597900390625e-1 */
{
y = x - one;
i = 0;
}
else if (ix >= 0x3ffced33)/* 2.31639862060546875e-1 */
{
y = x - (tc - one);
i = 1;
}
else
{
/* x < 0.23 */
y = x;
i = 2;
}
}
else
{
r = zero;
if (ix >= 0x3fffdda6) /* 1.73162841796875 */
{
/* [1.7316,2] */
y = x - 2.0;
i = 0;
}
else if (ix >= 0x3fff9da6)/* 1.23162841796875 */
{
/* [1.23,1.73] */
y = x - tc;
i = 1;
}
else
{
/* [0.9, 1.23] */
y = x - one;
i = 2;
}
}
switch (i)
{
case 0:
p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
r += half * y + y * p1/p2;
break;
case 1:
p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
p2 = h0 + y * (h1 + y * (h2 + y * (h3 + y * (h4 + y * (h5 + y)))));
p = tt + y * p1/p2;
r += (tf + p);
break;
case 2:
p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
r += (-half * y + p1 / p2);
}
}
else if (ix < 0x40028000) /* 8.0 */
{
/* x < 8.0 */
i = (int) x;
t = zero;
y = x - (double) i;
p = y *
(s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
r = half * y + p / q;
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i)
{
case 7:
z *= (y + 6.0); /* FALLTHRU */
case 6:
z *= (y + 5.0); /* FALLTHRU */
case 5:
z *= (y + 4.0); /* FALLTHRU */
case 4:
z *= (y + 3.0); /* FALLTHRU */
case 3:
z *= (y + 2.0); /* FALLTHRU */
r += logl (z);
break;
}
}
else if (ix < 0x40418000) /* 2^66 */
{
/* 8.0 <= x < 2**66 */
t = logl (x);
z = one / x;
y = z * z;
w = w0 + z * (w1
+ y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
r = (x - half) * (t - one) + w;
}
else
/* 2**66 <= x <= inf */
r = x * (logl (x) - one);
if (se & 0x8000)
r = nadj - r;
return r;
}

205
ld80/e_log10l.c Normal file
View file

@ -0,0 +1,205 @@
/* $OpenBSD: e_log10l.c,v 1.2 2013/11/12 20:35:19 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* log10l.c
*
* Common logarithm, long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, log10l();
*
* y = log10l( x );
*
*
*
* DESCRIPTION:
*
* Returns the base 10 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 9.0e-20 2.6e-20
* IEEE exp(+-10000) 30000 6.0e-20 2.3e-20
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns MINLOG
* log domain: x < 0; returns MINLOG
*/
#include <math.h>
#include "math_private.h"
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.2e-22
*/
static long double P[] = {
4.9962495940332550844739E-1L,
1.0767376367209449010438E1L,
7.7671073698359539859595E1L,
2.5620629828144409632571E2L,
4.2401812743503691187826E2L,
3.4258224542413922935104E2L,
1.0747524399916215149070E2L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
2.3479774160285863271658E1L,
1.9444210022760132894510E2L,
7.7952888181207260646090E2L,
1.6911722418503949084863E3L,
2.0307734695595183428202E3L,
1.2695660352705325274404E3L,
3.2242573199748645407652E2L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
/* log10(2) */
#define L102A 0.3125L
#define L102B -1.1470004336018804786261e-2L
/* log10(e) */
#define L10EA 0.5L
#define L10EB -6.5705518096748172348871e-2L
#define SQRTH 0.70710678118654752440L
long double
log10l(long double x)
{
long double y;
volatile long double z;
int e;
if( isnan(x) )
return(x);
/* Test for domain */
if( x <= 0.0L )
{
if( x == 0.0L )
return (-1.0L / (x - x));
else
return (x - x) / (x - x);
}
if( x == INFINITY )
return(INFINITY);
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl( x, &e );
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
y = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
goto done;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexpl( x, 1 ) - 1.0L; /* 2x - 1 */
}
else
{
x = x - 1.0L;
}
z = x*x;
y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 7 ) );
y = y - ldexpl( z, -1 ); /* -0.5x^2 + ... */
done:
/* Multiply log of fraction by log10(e)
* and base 2 exponent by log10(2).
*
* ***CAUTION***
*
* This sequence of operations is critical and it may
* be horribly defeated by some compiler optimizers.
*/
z = y * (L10EB);
z += x * (L10EB);
z += e * (L102B);
z += y * (L10EA);
z += x * (L10EA);
z += e * (L102A);
return( z );
}

199
ld80/e_log2l.c Normal file
View file

@ -0,0 +1,199 @@
/* $OpenBSD: e_log2l.c,v 1.2 2013/11/12 20:35:19 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* log2l.c
*
* Base 2 logarithm, long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, log2l();
*
* y = log2l( x );
*
*
*
* DESCRIPTION:
*
* Returns the base 2 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the (natural)
* logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 9.8e-20 2.7e-20
* IEEE exp(+-10000) 70000 5.4e-20 2.3e-20
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns -INFINITY
* log domain: x < 0; returns NAN
*/
#include <math.h>
#include "math_private.h"
/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.2e-22
*/
static long double P[] = {
4.9962495940332550844739E-1L,
1.0767376367209449010438E1L,
7.7671073698359539859595E1L,
2.5620629828144409632571E2L,
4.2401812743503691187826E2L,
3.4258224542413922935104E2L,
1.0747524399916215149070E2L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
2.3479774160285863271658E1L,
1.9444210022760132894510E2L,
7.7952888181207260646090E2L,
1.6911722418503949084863E3L,
2.0307734695595183428202E3L,
1.2695660352705325274404E3L,
3.2242573199748645407652E2L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
/* log2(e) - 1 */
#define LOG2EA 4.4269504088896340735992e-1L
#define SQRTH 0.70710678118654752440L
long double
log2l(long double x)
{
volatile long double z;
long double y;
int e;
if( isnan(x) )
return(x);
if( x == INFINITY )
return(x);
/* Test for domain */
if( x <= 0.0L )
{
if( x == 0.0L )
return( -INFINITY );
else
return( NAN );
}
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl( x, &e );
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
y = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
goto done;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexpl( x, 1 ) - 1.0L; /* 2x - 1 */
}
else
{
x = x - 1.0L;
}
z = x*x;
y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 7 ) );
y = y - ldexpl( z, -1 ); /* -0.5x^2 + ... */
done:
/* Multiply log of fraction by log2(e)
* and base 2 exponent by 1
*
* ***CAUTION***
*
* This sequence of operations is critical and it may
* be horribly defeated by some compiler optimizers.
*/
z = y * LOG2EA;
z += x * LOG2EA;
z += y;
z += x;
z += e;
return( z );
}

190
ld80/e_logl.c Normal file
View file

@ -0,0 +1,190 @@
/* $OpenBSD: e_logl.c,v 1.3 2013/11/12 20:35:19 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* logl.c
*
* Natural logarithm, long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, logl();
*
* y = logl( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 150000 8.71e-20 2.75e-20
* IEEE exp(+-10000) 100000 5.39e-20 2.34e-20
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns -INFINITY
* log domain: x < 0; returns NAN
*/
#include <math.h>
#include "math_private.h"
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 2.32e-20
*/
static long double P[] = {
4.5270000862445199635215E-5L,
4.9854102823193375972212E-1L,
6.5787325942061044846969E0L,
2.9911919328553073277375E1L,
6.0949667980987787057556E1L,
5.7112963590585538103336E1L,
2.0039553499201281259648E1L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
1.5062909083469192043167E1L,
8.3047565967967209469434E1L,
2.2176239823732856465394E2L,
3.0909872225312059774938E2L,
2.1642788614495947685003E2L,
6.0118660497603843919306E1L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
static const long double C1 = 6.9314575195312500000000E-1L;
static const long double C2 = 1.4286068203094172321215E-6L;
#define SQRTH 0.70710678118654752440L
long double
logl(long double x)
{
long double y, z;
int e;
if( isnan(x) )
return(x);
if( x == INFINITY )
return(x);
/* Test for domain */
if( x <= 0.0L )
{
if( x == 0.0L )
return( -INFINITY );
else
return( NAN );
}
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl( x, &e );
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
z = z + e * C2;
z = z + x;
z = z + e * C1;
return( z );
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexpl( x, 1 ) - 1.0L; /* 2x - 1 */
}
else
{
x = x - 1.0L;
}
z = x*x;
y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) );
y = y + e * C2;
z = y - ldexpl( z, -1 ); /* y - 0.5 * z */
/* Note, the sum of above terms does not exceed x/4,
* so it contributes at most about 1/4 lsb to the error.
*/
z = z + x;
z = z + e * C1; /* This sum has an error of 1/2 lsb. */
return( z );
}

615
ld80/e_powl.c Normal file
View file

@ -0,0 +1,615 @@
/* $OpenBSD: e_powl.c,v 1.5 2013/11/12 20:35:19 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* powl.c
*
* Power function, long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, z, powl();
*
* z = powl( x, y );
*
*
*
* DESCRIPTION:
*
* Computes x raised to the yth power. Analytically,
*
* x**y = exp( y log(x) ).
*
* Following Cody and Waite, this program uses a lookup table
* of 2**-i/32 and pseudo extended precision arithmetic to
* obtain several extra bits of accuracy in both the logarithm
* and the exponential.
*
*
*
* ACCURACY:
*
* The relative error of pow(x,y) can be estimated
* by y dl ln(2), where dl is the absolute error of
* the internally computed base 2 logarithm. At the ends
* of the approximation interval the logarithm equal 1/32
* and its relative error is about 1 lsb = 1.1e-19. Hence
* the predicted relative error in the result is 2.3e-21 y .
*
* Relative error:
* arithmetic domain # trials peak rms
*
* IEEE +-1000 40000 2.8e-18 3.7e-19
* .001 < x < 1000, with log(x) uniformly distributed.
* -1000 < y < 1000, y uniformly distributed.
*
* IEEE 0,8700 60000 6.5e-18 1.0e-18
* 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* pow overflow x**y > MAXNUM INFINITY
* pow underflow x**y < 1/MAXNUM 0.0
* pow domain x<0 and y noninteger 0.0
*
*/
#include <float.h>
#include <math.h>
#include "math_private.h"
/* Table size */
#define NXT 32
/* log2(Table size) */
#define LNXT 5
/* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z)
* on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1
*/
static long double P[] = {
8.3319510773868690346226E-4L,
4.9000050881978028599627E-1L,
1.7500123722550302671919E0L,
1.4000100839971580279335E0L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0L,*/
5.2500282295834889175431E0L,
8.4000598057587009834666E0L,
4.2000302519914740834728E0L,
};
/* A[i] = 2^(-i/32), rounded to IEEE long double precision.
* If i is even, A[i] + B[i/2] gives additional accuracy.
*/
static long double A[33] = {
1.0000000000000000000000E0L,
9.7857206208770013448287E-1L,
9.5760328069857364691013E-1L,
9.3708381705514995065011E-1L,
9.1700404320467123175367E-1L,
8.9735453750155359320742E-1L,
8.7812608018664974155474E-1L,
8.5930964906123895780165E-1L,
8.4089641525371454301892E-1L,
8.2287773907698242225554E-1L,
8.0524516597462715409607E-1L,
7.8799042255394324325455E-1L,
7.7110541270397041179298E-1L,
7.5458221379671136985669E-1L,
7.3841307296974965571198E-1L,
7.2259040348852331001267E-1L,
7.0710678118654752438189E-1L,
6.9195494098191597746178E-1L,
6.7712777346844636413344E-1L,
6.6261832157987064729696E-1L,
6.4841977732550483296079E-1L,
6.3452547859586661129850E-1L,
6.2092890603674202431705E-1L,
6.0762367999023443907803E-1L,
5.9460355750136053334378E-1L,
5.8186242938878875689693E-1L,
5.6939431737834582684856E-1L,
5.5719337129794626814472E-1L,
5.4525386633262882960438E-1L,
5.3357020033841180906486E-1L,
5.2213689121370692017331E-1L,
5.1094857432705833910408E-1L,
5.0000000000000000000000E-1L,
};
static long double B[17] = {
0.0000000000000000000000E0L,
2.6176170809902549338711E-20L,
-1.0126791927256478897086E-20L,
1.3438228172316276937655E-21L,
1.2207982955417546912101E-20L,
-6.3084814358060867200133E-21L,
1.3164426894366316434230E-20L,
-1.8527916071632873716786E-20L,
1.8950325588932570796551E-20L,
1.5564775779538780478155E-20L,
6.0859793637556860974380E-21L,
-2.0208749253662532228949E-20L,
1.4966292219224761844552E-20L,
3.3540909728056476875639E-21L,
-8.6987564101742849540743E-22L,
-1.2327176863327626135542E-20L,
0.0000000000000000000000E0L,
};
/* 2^x = 1 + x P(x),
* on the interval -1/32 <= x <= 0
*/
static long double R[] = {
1.5089970579127659901157E-5L,
1.5402715328927013076125E-4L,
1.3333556028915671091390E-3L,
9.6181291046036762031786E-3L,
5.5504108664798463044015E-2L,
2.4022650695910062854352E-1L,
6.9314718055994530931447E-1L,
};
#define douba(k) A[k]
#define doubb(k) B[k]
#define MEXP (NXT*16384.0L)
/* The following if denormal numbers are supported, else -MEXP: */
#define MNEXP (-NXT*(16384.0L+64.0L))
/* log2(e) - 1 */
#define LOG2EA 0.44269504088896340735992L
#define F W
#define Fa Wa
#define Fb Wb
#define G W
#define Ga Wa
#define Gb u
#define H W
#define Ha Wb
#define Hb Wb
static const long double MAXLOGL = 1.1356523406294143949492E4L;
static const long double MINLOGL = -1.13994985314888605586758E4L;
static const long double LOGE2L = 6.9314718055994530941723E-1L;
static volatile long double z;
static long double w, W, Wa, Wb, ya, yb, u;
static const long double huge = 0x1p10000L;
#if 0 /* XXX Prevent gcc from erroneously constant folding this. */
static const long double twom10000 = 0x1p-10000L;
#else
static volatile long double twom10000 = 0x1p-10000L;
#endif
static long double reducl( long double );
static long double powil ( long double, int );
long double
powl(long double x, long double y)
{
/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
int i, nflg, iyflg, yoddint;
long e;
if( y == 0.0L )
return( 1.0L );
if( x == 1.0L )
return( 1.0L );
if( isnan(x) )
return( x );
if( isnan(y) )
return( y );
if( y == 1.0L )
return( x );
if( !isfinite(y) && x == -1.0L )
return( 1.0L );
if( y >= LDBL_MAX )
{
if( x > 1.0L )
return( INFINITY );
if( x > 0.0L && x < 1.0L )
return( 0.0L );
if( x < -1.0L )
return( INFINITY );
if( x > -1.0L && x < 0.0L )
return( 0.0L );
}
if( y <= -LDBL_MAX )
{
if( x > 1.0L )
return( 0.0L );
if( x > 0.0L && x < 1.0L )
return( INFINITY );
if( x < -1.0L )
return( 0.0L );
if( x > -1.0L && x < 0.0L )
return( INFINITY );
}
if( x >= LDBL_MAX )
{
if( y > 0.0L )
return( INFINITY );
return( 0.0L );
}
w = floorl(y);
/* Set iyflg to 1 if y is an integer. */
iyflg = 0;
if( w == y )
iyflg = 1;
/* Test for odd integer y. */
yoddint = 0;
if( iyflg )
{
ya = fabsl(y);
ya = floorl(0.5L * ya);
yb = 0.5L * fabsl(w);
if( ya != yb )
yoddint = 1;
}
if( x <= -LDBL_MAX )
{
if( y > 0.0L )
{
if( yoddint )
return( -INFINITY );
return( INFINITY );
}
if( y < 0.0L )
{
if( yoddint )
return( -0.0L );
return( 0.0 );
}
}
nflg = 0; /* flag = 1 if x<0 raised to integer power */
if( x <= 0.0L )
{
if( x == 0.0L )
{
if( y < 0.0 )
{
if( signbit(x) && yoddint )
return( -INFINITY );
return( INFINITY );
}
if( y > 0.0 )
{
if( signbit(x) && yoddint )
return( -0.0L );
return( 0.0 );
}
if( y == 0.0L )
return( 1.0L ); /* 0**0 */
else
return( 0.0L ); /* 0**y */
}
else
{
if( iyflg == 0 )
return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
nflg = 1;
}
}
/* Integer power of an integer. */
if( iyflg )
{
i = w;
w = floorl(x);
if( (w == x) && (fabsl(y) < 32768.0) )
{
w = powil( x, (int) y );
return( w );
}
}
if( nflg )
x = fabsl(x);
/* separate significand from exponent */
x = frexpl( x, &i );
e = i;
/* find significand in antilog table A[] */
i = 1;
if( x <= douba(17) )
i = 17;
if( x <= douba(i+8) )
i += 8;
if( x <= douba(i+4) )
i += 4;
if( x <= douba(i+2) )
i += 2;
if( x >= douba(1) )
i = -1;
i += 1;
/* Find (x - A[i])/A[i]
* in order to compute log(x/A[i]):
*
* log(x) = log( a x/a ) = log(a) + log(x/a)
*
* log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a
*/
x -= douba(i);
x -= doubb(i/2);
x /= douba(i);
/* rational approximation for log(1+v):
*
* log(1+v) = v - v**2/2 + v**3 P(v) / Q(v)
*/
z = x*x;
w = x * ( z * __polevll( x, P, 3 ) / __p1evll( x, Q, 3 ) );
w = w - ldexpl( z, -1 ); /* w - 0.5 * z */
/* Convert to base 2 logarithm:
* multiply by log2(e) = 1 + LOG2EA
*/
z = LOG2EA * w;
z += w;
z += LOG2EA * x;
z += x;
/* Compute exponent term of the base 2 logarithm. */
w = -i;
w = ldexpl( w, -LNXT ); /* divide by NXT */
w += e;
/* Now base 2 log of x is w + z. */
/* Multiply base 2 log by y, in extended precision. */
/* separate y into large part ya
* and small part yb less than 1/NXT
*/
ya = reducl(y);
yb = y - ya;
/* (w+z)(ya+yb)
* = w*ya + w*yb + z*y
*/
F = z * y + w * yb;
Fa = reducl(F);
Fb = F - Fa;
G = Fa + w * ya;
Ga = reducl(G);
Gb = G - Ga;
H = Fb + Gb;
Ha = reducl(H);
w = ldexpl( Ga+Ha, LNXT );
/* Test the power of 2 for overflow */
if( w > MEXP )
return (huge * huge); /* overflow */
if( w < MNEXP )
return (twom10000 * twom10000); /* underflow */
e = w;
Hb = H - Ha;
if( Hb > 0.0L )
{
e += 1;
Hb -= (1.0L/NXT); /*0.0625L;*/
}
/* Now the product y * log2(x) = Hb + e/NXT.
*
* Compute base 2 exponential of Hb,
* where -0.0625 <= Hb <= 0.
*/
z = Hb * __polevll( Hb, R, 6 ); /* z = 2**Hb - 1 */
/* Express e/NXT as an integer plus a negative number of (1/NXT)ths.
* Find lookup table entry for the fractional power of 2.
*/
if( e < 0 )
i = 0;
else
i = 1;
i = e/NXT + i;
e = NXT*i - e;
w = douba( e );
z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = z + w;
z = ldexpl( z, i ); /* multiply by integer power of 2 */
if( nflg )
{
/* For negative x,
* find out if the integer exponent
* is odd or even.
*/
w = ldexpl( y, -1 );
w = floorl(w);
w = ldexpl( w, 1 );
if( w != y )
z = -z; /* odd exponent */
}
return( z );
}
/* Find a multiple of 1/NXT that is within 1/NXT of x. */
static long double
reducl(long double x)
{
long double t;
t = ldexpl( x, LNXT );
t = floorl( t );
t = ldexpl( t, -LNXT );
return(t);
}
/* powil.c
*
* Real raised to integer power, long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, powil();
* int n;
*
* y = powil( x, n );
*
*
*
* DESCRIPTION:
*
* Returns argument x raised to the nth power.
* The routine efficiently decomposes n as a sum of powers of
* two. The desired power is a product of two-to-the-kth
* powers of x. Thus to compute the 32767 power of x requires
* 28 multiplications instead of 32767 multiplications.
*
*
*
* ACCURACY:
*
*
* Relative error:
* arithmetic x domain n domain # trials peak rms
* IEEE .001,1000 -1022,1023 50000 4.3e-17 7.8e-18
* IEEE 1,2 -1022,1023 20000 3.9e-17 7.6e-18
* IEEE .99,1.01 0,8700 10000 3.6e-16 7.2e-17
*
* Returns MAXNUM on overflow, zero on underflow.
*
*/
static long double
powil(long double x, int nn)
{
long double ww, y;
long double s;
int n, e, sign, asign, lx;
if( x == 0.0L )
{
if( nn == 0 )
return( 1.0L );
else if( nn < 0 )
return( LDBL_MAX );
else
return( 0.0L );
}
if( nn == 0 )
return( 1.0L );
if( x < 0.0L )
{
asign = -1;
x = -x;
}
else
asign = 0;
if( nn < 0 )
{
sign = -1;
n = -nn;
}
else
{
sign = 1;
n = nn;
}
/* Overflow detection */
/* Calculate approximate logarithm of answer */
s = x;
s = frexpl( s, &lx );
e = (lx - 1)*n;
if( (e == 0) || (e > 64) || (e < -64) )
{
s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L);
s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
}
else
{
s = LOGE2L * e;
}
if( s > MAXLOGL )
return (huge * huge); /* overflow */
if( s < MINLOGL )
return (twom10000 * twom10000); /* underflow */
/* Handle tiny denormal answer, but with less accuracy
* since roundoff error in 1.0/x will be amplified.
* The precise demarcation should be the gradual underflow threshold.
*/
if( s < (-MAXLOGL+2.0L) )
{
x = 1.0L/x;
sign = -sign;
}
/* First bit of the power */
if( n & 1 )
y = x;
else
{
y = 1.0L;
asign = 0;
}
ww = x;
n >>= 1;
while( n )
{
ww = ww * ww; /* arg to the 2-to-the-kth power */
if( n & 1 ) /* if that bit is set, then include in product */
y *= ww;
n >>= 1;
}
if( asign )
y = -y; /* odd power of negative number */
if( sign < 0 )
y = 1.0L/y;
return(y);
}

76
ld80/e_sinhl.c Normal file
View file

@ -0,0 +1,76 @@
/* @(#)e_sinh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* sinhl(x)
* Method :
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
* 1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
* 2.
* E + E/(E+1)
* 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x)
* 2
*
* 25 <= x <= lnovft : sinhl(x) := expl(x)/2
* lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2)
* ln2ovft < x : sinhl(x) := x*shuge (overflow)
*
* Special cases:
* sinhl(x) is |x| if x is +INF, -INF, or NaN.
* only sinhl(0)=0 is exact for finite x.
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0, shuge = 1.0e4931L;
long double
sinhl(long double x)
{
long double t,w,h;
u_int32_t jx,ix,i0,i1;
/* Words of |x|. */
GET_LDOUBLE_WORDS(jx,i0,i1,x);
ix = jx&0x7fff;
/* x is INF or NaN */
if(ix==0x7fff) return x+x;
h = 0.5;
if (jx & 0x8000) h = -h;
/* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */
if (ix<0x3fdf) /* |x|<2**-32 */
if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
t = expm1l(fabsl(x));
if(ix<0x3fff) return h*(2.0*t-t*t/(t+one));
return h*(t+t/(t+one));
}
/* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7))
return h*expl(fabsl(x));
/* |x| in [log(maxdouble), overflowthreshold] */
if (ix<0x400c || (ix == 0x400c && (i0 < 0xb174ddc0
|| (i0 == 0xb174ddc0
&& i1 <= 0x31aec0ea)))) {
w = expl(0.5*fabsl(x));
t = h*w;
return t*w;
}
/* |x| > overflowthreshold, sinhl(x) overflow */
return x*shuge;
}

319
ld80/e_tgammal.c Normal file
View file

@ -0,0 +1,319 @@
/* $OpenBSD: e_tgammal.c,v 1.4 2013/11/12 20:35:19 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* tgammal.c
*
* Gamma function
*
*
*
* SYNOPSIS:
*
* long double x, y, tgammal();
* extern int signgam;
*
* y = tgammal( x );
*
*
*
* DESCRIPTION:
*
* Returns gamma function of the argument. The result is
* correctly signed, and the sign (+1 or -1) is also
* returned in a global (extern) variable named signgam.
* This variable is also filled in by the logarithmic gamma
* function lgamma().
*
* Arguments |x| <= 13 are reduced by recurrence and the function
* approximated by a rational function of degree 7/8 in the
* interval (2,3). Large arguments are handled by Stirling's
* formula. Large negative arguments are made positive using
* a reflection formula.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -40,+40 10000 3.6e-19 7.9e-20
* IEEE -1755,+1755 10000 4.8e-18 6.5e-19
*
* Accuracy for large arguments is dominated by error in powl().
*
*/
#include <float.h>
#include <math.h>
#include "math_private.h"
/*
tgamma(x+2) = tgamma(x+2) P(x)/Q(x)
0 <= x <= 1
Relative error
n=7, d=8
Peak error = 1.83e-20
Relative error spread = 8.4e-23
*/
static long double P[8] = {
4.212760487471622013093E-5L,
4.542931960608009155600E-4L,
4.092666828394035500949E-3L,
2.385363243461108252554E-2L,
1.113062816019361559013E-1L,
3.629515436640239168939E-1L,
8.378004301573126728826E-1L,
1.000000000000000000009E0L,
};
static long double Q[9] = {
-1.397148517476170440917E-5L,
2.346584059160635244282E-4L,
-1.237799246653152231188E-3L,
-7.955933682494738320586E-4L,
2.773706565840072979165E-2L,
-4.633887671244534213831E-2L,
-2.243510905670329164562E-1L,
4.150160950588455434583E-1L,
9.999999999999999999908E-1L,
};
/*
static long double P[] = {
-3.01525602666895735709e0L,
-3.25157411956062339893e1L,
-2.92929976820724030353e2L,
-1.70730828800510297666e3L,
-7.96667499622741999770e3L,
-2.59780216007146401957e4L,
-5.99650230220855581642e4L,
-7.15743521530849602425e4L
};
static long double Q[] = {
1.00000000000000000000e0L,
-1.67955233807178858919e1L,
8.85946791747759881659e1L,
5.69440799097468430177e1L,
-1.98526250512761318471e3L,
3.31667508019495079814e3L,
1.60577839621734713377e4L,
-2.97045081369399940529e4L,
-7.15743521530849602412e4L
};
*/
#define MAXGAML 1755.455L
/*static const long double LOGPI = 1.14472988584940017414L;*/
/* Stirling's formula for the gamma function
tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
z(x) = x
13 <= x <= 1024
Relative error
n=8, d=0
Peak error = 9.44e-21
Relative error spread = 8.8e-4
*/
static long double STIR[9] = {
7.147391378143610789273E-4L,
-2.363848809501759061727E-5L,
-5.950237554056330156018E-4L,
6.989332260623193171870E-5L,
7.840334842744753003862E-4L,
-2.294719747873185405699E-4L,
-2.681327161876304418288E-3L,
3.472222222230075327854E-3L,
8.333333333333331800504E-2L,
};
#define MAXSTIR 1024.0L
static const long double SQTPI = 2.50662827463100050242E0L;
/* 1/tgamma(x) = z P(z)
* z(x) = 1/x
* 0 < x < 0.03125
* Peak relative error 4.2e-23
*/
static long double S[9] = {
-1.193945051381510095614E-3L,
7.220599478036909672331E-3L,
-9.622023360406271645744E-3L,
-4.219773360705915470089E-2L,
1.665386113720805206758E-1L,
-4.200263503403344054473E-2L,
-6.558780715202540684668E-1L,
5.772156649015328608253E-1L,
1.000000000000000000000E0L,
};
/* 1/tgamma(-x) = z P(z)
* z(x) = 1/x
* 0 < x < 0.03125
* Peak relative error 5.16e-23
* Relative error spread = 2.5e-24
*/
static long double SN[9] = {
1.133374167243894382010E-3L,
7.220837261893170325704E-3L,
9.621911155035976733706E-3L,
-4.219773343731191721664E-2L,
-1.665386113944413519335E-1L,
-4.200263503402112910504E-2L,
6.558780715202536547116E-1L,
5.772156649015328608727E-1L,
-1.000000000000000000000E0L,
};
static const long double PIL = 3.1415926535897932384626L;
static long double stirf ( long double );
/* Gamma function computed by Stirling's formula.
*/
static long double stirf(long double x)
{
long double y, w, v;
w = 1.0L/x;
/* For large x, use rational coefficients from the analytical expansion. */
if( x > 1024.0L )
w = (((((6.97281375836585777429E-5L * w
+ 7.84039221720066627474E-4L) * w
- 2.29472093621399176955E-4L) * w
- 2.68132716049382716049E-3L) * w
+ 3.47222222222222222222E-3L) * w
+ 8.33333333333333333333E-2L) * w
+ 1.0L;
else
w = 1.0L + w * __polevll( w, STIR, 8 );
y = expl(x);
if( x > MAXSTIR )
{ /* Avoid overflow in pow() */
v = powl( x, 0.5L * x - 0.25L );
y = v * (v / y);
}
else
{
y = powl( x, x - 0.5L ) / y;
}
y = SQTPI * y * w;
return( y );
}
long double
tgammal(long double x)
{
long double p, q, z;
int i;
signgam = 1;
if( isnan(x) )
return(NAN);
if(x == INFINITY)
return(INFINITY);
if(x == -INFINITY)
return(x - x);
if( x == 0.0L )
return( 1.0L / x );
q = fabsl(x);
if( q > 13.0L )
{
if( q > MAXGAML )
goto goverf;
if( x < 0.0L )
{
p = floorl(q);
if( p == q )
return (x - x) / (x - x);
i = p;
if( (i & 1) == 0 )
signgam = -1;
z = q - p;
if( z > 0.5L )
{
p += 1.0L;
z = q - p;
}
z = q * sinl( PIL * z );
z = fabsl(z) * stirf(q);
if( z <= PIL/LDBL_MAX )
{
goverf:
return( signgam * INFINITY);
}
z = PIL/z;
}
else
{
z = stirf(x);
}
return( signgam * z );
}
z = 1.0L;
while( x >= 3.0L )
{
x -= 1.0L;
z *= x;
}
while( x < -0.03125L )
{
z /= x;
x += 1.0L;
}
if( x <= 0.03125L )
goto small;
while( x < 2.0L )
{
z /= x;
x += 1.0L;
}
if( x == 2.0L )
return(z);
x -= 2.0L;
p = __polevll( x, P, 7 );
q = __polevll( x, Q, 8 );
z = z * p / q;
if( z < 0 )
signgam = -1;
return z;
small:
if( x == 0.0L )
return (x - x) / (x - x);
else
{
if( x < 0.0L )
{
x = -x;
q = z / (x * __polevll( x, SN, 8 ));
signgam = -1;
}
else
q = z / (x * __polevll( x, S, 8 ));
}
return q;
}

54
ld80/s_asinhl.c Normal file
View file

@ -0,0 +1,54 @@
/* @(#)s_asinh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* asinhl(x)
* Method :
* Based on
* asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ]
* we have
* asinhl(x) := x if 1+x*x=1,
* := signl(x)*(logl(x)+ln2)) for large |x|, else
* := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else
* := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2)))
*/
#include <math.h>
#include "math_private.h"
static const long double
one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */
ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
huge= 1.000000000000000000e+4900L;
long double
asinhl(long double x)
{
long double t,w;
int32_t hx,ix;
GET_LDOUBLE_EXP(hx,x);
ix = hx&0x7fff;
if(ix==0x7fff) return x+x; /* x is inf or NaN */
if(ix< 0x3fde) { /* |x|<2**-34 */
if(huge+x>one) return x; /* return x inexact except 0 */
}
if(ix>0x4020) { /* |x| > 2**34 */
w = logl(fabsl(x))+ln2;
} else if (ix>0x4000) { /* 2**34 > |x| > 2.0 */
t = fabsl(x);
w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
}
if(hx&0x8000) return -w; else return w;
}

78
ld80/s_ceill.c Normal file
View file

@ -0,0 +1,78 @@
/* @(#)s_ceil.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* ceill(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to ceil(x).
*/
#include <math.h>
#include "math_private.h"
static const long double huge = 1.0e4930L;
long double
ceill(long double x)
{
int32_t i1,jj0;
u_int32_t i,j,se,i0,sx;
GET_LDOUBLE_WORDS(se,i0,i1,x);
sx = (se>>15)&1;
jj0 = (se&0x7fff)-0x3fff;
if(jj0<31) {
if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
if(sx) {se=0x8000;i0=0;i1=0;}
else if((i0|i1)!=0) { se=0x3fff;i0=0;i1=0;}
}
} else {
i = (0x7fffffff)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(sx==0) {
if (jj0>0 && (i0+(0x80000000>>jj0))>i0)
i0+=0x80000000>>jj0;
else
{
i = 0x7fffffff;
++se;
}
}
i0 &= (~i); i1=0;
}
}
} else if (jj0>62) {
if(jj0==0x4000) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = ((u_int32_t)(0xffffffff))>>(jj0-31);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(sx==0) {
if(jj0==31) i0+=1;
else {
j = i1 + (1<<(63-jj0));
if(j<i1) i0+=1; /* got a carry */
i1 = j;
}
}
i1 &= (~i);
}
}
SET_LDOUBLE_WORDS(x,se,i0,i1);
return x;
}

430
ld80/s_erfl.c Normal file
View file

@ -0,0 +1,430 @@
/*
* ====================================================
* 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) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* double erf(double x)
* double erfc(double x)
* x
* 2 |\
* erf(x) = --------- | exp(-t*t)dt
* sqrt(pi) \|
* 0
*
* erfc(x) = 1-erf(x)
* Note that
* erf(-x) = -erf(x)
* erfc(-x) = 2 - erfc(x)
*
* Method:
* 1. For |x| in [0, 0.84375]
* erf(x) = x + x*R(x^2)
* erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
* = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
* Remark. The formula is derived by noting
* erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
* and that
* 2/sqrt(pi) = 1.128379167095512573896158903121545171688
* is close to one. The interval is chosen because the fix
* point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
* near 0.6174), and by some experiment, 0.84375 is chosen to
* guarantee the error is less than one ulp for erf.
*
* 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
* c = 0.84506291151 rounded to single (24 bits)
* erf(x) = sign(x) * (c + P1(s)/Q1(s))
* erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
* 1+(c+P1(s)/Q1(s)) if x < 0
* Remark: here we use the taylor series expansion at x=1.
* erf(1+s) = erf(1) + s*Poly(s)
* = 0.845.. + P1(s)/Q1(s)
* Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
*
* 3. For x in [1.25,1/0.35(~2.857143)],
* erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
* z=1/x^2
* erf(x) = 1 - erfc(x)
*
* 4. For x in [1/0.35,107]
* erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
* = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
* if -6.666<x<0
* = 2.0 - tiny (if x <= -6.666)
* z=1/x^2
* erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
* erf(x) = sign(x)*(1.0 - tiny)
* Note1:
* To compute exp(-x*x-0.5625+R/S), let s be a single
* precision number and s := x; then
* -x*x = -s*s + (s-x)*(s+x)
* exp(-x*x-0.5626+R/S) =
* exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
* Note2:
* Here 4 and 5 make use of the asymptotic series
* exp(-x*x)
* erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
* x*sqrt(pi)
*
* 5. For inf > x >= 107
* erf(x) = sign(x) *(1 - tiny) (raise inexact)
* erfc(x) = tiny*tiny (raise underflow) if x > 0
* = 2 - tiny if x<0
*
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*/
#include <math.h>
#include "math_private.h"
static const long double
tiny = 1e-4931L,
half = 0.5L,
one = 1.0L,
two = 2.0L,
/* c = (float)0.84506291151 */
erx = 0.845062911510467529296875L,
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
/* 2/sqrt(pi) - 1 */
efx = 1.2837916709551257389615890312154517168810E-1L,
/* 8 * (2/sqrt(pi) - 1) */
efx8 = 1.0270333367641005911692712249723613735048E0L,
pp[6] = {
1.122751350964552113068262337278335028553E6L,
-2.808533301997696164408397079650699163276E6L,
-3.314325479115357458197119660818768924100E5L,
-6.848684465326256109712135497895525446398E4L,
-2.657817695110739185591505062971929859314E3L,
-1.655310302737837556654146291646499062882E2L,
},
qq[6] = {
8.745588372054466262548908189000448124232E6L,
3.746038264792471129367533128637019611485E6L,
7.066358783162407559861156173539693900031E5L,
7.448928604824620999413120955705448117056E4L,
4.511583986730994111992253980546131408924E3L,
1.368902937933296323345610240009071254014E2L,
/* 1.000000000000000000000000000000000000000E0 */
},
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
/* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
-0.15625 <= x <= +.25
Peak relative error 8.5e-22 */
pa[8] = {
-1.076952146179812072156734957705102256059E0L,
1.884814957770385593365179835059971587220E2L,
-5.339153975012804282890066622962070115606E1L,
4.435910679869176625928504532109635632618E1L,
1.683219516032328828278557309642929135179E1L,
-2.360236618396952560064259585299045804293E0L,
1.852230047861891953244413872297940938041E0L,
9.394994446747752308256773044667843200719E-2L,
},
qa[7] = {
4.559263722294508998149925774781887811255E2L,
3.289248982200800575749795055149780689738E2L,
2.846070965875643009598627918383314457912E2L,
1.398715859064535039433275722017479994465E2L,
6.060190733759793706299079050985358190726E1L,
2.078695677795422351040502569964299664233E1L,
4.641271134150895940966798357442234498546E0L,
/* 1.000000000000000000000000000000000000000E0 */
},
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
1/2.85711669921875 < 1/x < 1/1.25
Peak relative error 3.1e-21 */
ra[] = {
1.363566591833846324191000679620738857234E-1L,
1.018203167219873573808450274314658434507E1L,
1.862359362334248675526472871224778045594E2L,
1.411622588180721285284945138667933330348E3L,
5.088538459741511988784440103218342840478E3L,
8.928251553922176506858267311750789273656E3L,
7.264436000148052545243018622742770549982E3L,
2.387492459664548651671894725748959751119E3L,
2.220916652813908085449221282808458466556E2L,
},
sa[] = {
-1.382234625202480685182526402169222331847E1L,
-3.315638835627950255832519203687435946482E2L,
-2.949124863912936259747237164260785326692E3L,
-1.246622099070875940506391433635999693661E4L,
-2.673079795851665428695842853070996219632E4L,
-2.880269786660559337358397106518918220991E4L,
-1.450600228493968044773354186390390823713E4L,
-2.874539731125893533960680525192064277816E3L,
-1.402241261419067750237395034116942296027E2L,
/* 1.000000000000000000000000000000000000000E0 */
},
/*
* Coefficients for approximation to erfc in [1/.35,107]
*/
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
1/6.6666259765625 < 1/x < 1/2.85711669921875
Peak relative error 4.2e-22 */
rb[] = {
-4.869587348270494309550558460786501252369E-5L,
-4.030199390527997378549161722412466959403E-3L,
-9.434425866377037610206443566288917589122E-2L,
-9.319032754357658601200655161585539404155E-1L,
-4.273788174307459947350256581445442062291E0L,
-8.842289940696150508373541814064198259278E0L,
-7.069215249419887403187988144752613025255E0L,
-1.401228723639514787920274427443330704764E0L,
},
sb[] = {
4.936254964107175160157544545879293019085E-3L,
1.583457624037795744377163924895349412015E-1L,
1.850647991850328356622940552450636420484E0L,
9.927611557279019463768050710008450625415E0L,
2.531667257649436709617165336779212114570E1L,
2.869752886406743386458304052862814690045E1L,
1.182059497870819562441683560749192539345E1L,
/* 1.000000000000000000000000000000000000000E0 */
},
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
1/107 <= 1/x <= 1/6.6666259765625
Peak relative error 1.1e-21 */
rc[] = {
-8.299617545269701963973537248996670806850E-5L,
-6.243845685115818513578933902532056244108E-3L,
-1.141667210620380223113693474478394397230E-1L,
-7.521343797212024245375240432734425789409E-1L,
-1.765321928311155824664963633786967602934E0L,
-1.029403473103215800456761180695263439188E0L,
},
sc[] = {
8.413244363014929493035952542677768808601E-3L,
2.065114333816877479753334599639158060979E-1L,
1.639064941530797583766364412782135680148E0L,
4.936788463787115555582319302981666347450E0L,
5.005177727208955487404729933261347679090E0L,
/* 1.000000000000000000000000000000000000000E0 */
};
long double
erfl(long double x)
{
long double R, S, P, Q, s, y, z, r;
int32_t ix, i;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff)
{ /* erf(nan)=nan */
i = ((se & 0xffff) >> 15) << 1;
return (long double) (1 - i) + one / x; /* erf(+-inf)=+-1 */
}
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) /* |x|<0.84375 */
{
if (ix < 0x3fde8000) /* |x|<2**-33 */
{
if (ix < 0x00080000)
return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */
return x + efx * x;
}
z = x * x;
r = pp[0] + z * (pp[1]
+ z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
s = qq[0] + z * (qq[1]
+ z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
return x + x * y;
}
if (ix < 0x3fffa000) /* 1.25 */
{ /* 0.84375 <= |x| < 1.25 */
s = fabsl (x) - one;
P = pa[0] + s * (pa[1] + s * (pa[2]
+ s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2]
+ s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0)
return erx + P / Q;
else
return -erx - P / Q;
}
if (ix >= 0x4001d555) /* 6.6666259765625 */
{ /* inf>|x|>=6.666 */
if ((se & 0x8000) == 0)
return one - tiny;
else
return tiny - one;
}
x = fabsl (x);
s = one / (x * x);
if (ix < 0x4000b6db) /* 2.85711669921875 */
{
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
}
else
{ /* |x| >= 1/0.35 */
R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
s * (rb[5] + s * (rb[6] + s * rb[7]))))));
S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
s * (sb[5] + s * (sb[6] + s))))));
}
z = x;
GET_LDOUBLE_WORDS (i, i0, i1, z);
i1 = 0;
SET_LDOUBLE_WORDS (z, i, i0, i1);
r =
expl (-z * z - 0.5625) * expl ((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
return one - r / x;
else
return r / x - one;
}
long double
erfcl(long double x)
{
int32_t hx, ix;
long double R, S, P, Q, s, y, z, r;
u_int32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff)
{ /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (long double) (((se & 0xffff) >> 15) << 1) + one / x;
}
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) /* |x|<0.84375 */
{
if (ix < 0x3fbe0000) /* |x|<2**-65 */
return one - x;
z = x * x;
r = pp[0] + z * (pp[1]
+ z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
s = qq[0] + z * (qq[1]
+ z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
if (ix < 0x3ffd8000) /* x<1/4 */
{
return one - (x + x * y);
}
else
{
r = x * y;
r += (x - half);
return half - r;
}
}
if (ix < 0x3fffa000) /* 1.25 */
{ /* 0.84375 <= |x| < 1.25 */
s = fabsl (x) - one;
P = pa[0] + s * (pa[1] + s * (pa[2]
+ s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2]
+ s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0)
{
z = one - erx;
return z - P / Q;
}
else
{
z = erx + P / Q;
return one + z;
}
}
if (ix < 0x4005d600) /* 107 */
{ /* |x|<107 */
x = fabsl (x);
s = one / (x * x);
if (ix < 0x4000b6db) /* 2.85711669921875 */
{ /* |x| < 1/.35 ~ 2.857143 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
}
else if (ix < 0x4001d555) /* 6.6666259765625 */
{ /* 6.666 > |x| >= 1/.35 ~ 2.857143 */
R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
s * (rb[5] + s * (rb[6] + s * rb[7]))))));
S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
s * (sb[5] + s * (sb[6] + s))))));
}
else
{ /* |x| >= 6.666 */
if (se & 0x8000)
return two - tiny; /* x < -6.666 */
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
s * (sc[4] + s))));
}
z = x;
GET_LDOUBLE_WORDS (hx, i0, i1, z);
i1 = 0;
i0 &= 0xffffff00;
SET_LDOUBLE_WORDS (z, hx, i0, i1);
r = expl (-z * z - 0.5625) *
expl ((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
return r / x;
else
return two - r / x;
}
else
{
if ((se & 0x8000) == 0)
return tiny * tiny;
else
return two - tiny;
}
}

138
ld80/s_expm1l.c Normal file
View file

@ -0,0 +1,138 @@
/* $OpenBSD: s_expm1l.c,v 1.2 2011/07/20 21:02:51 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* expm1l.c
*
* Exponential function, minus 1
* Long double precision
*
*
* SYNOPSIS:
*
* long double x, y, expm1l();
*
* y = expm1l( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power, minus 1.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
* in the basic range [-0.5 ln 2, 0.5 ln 2].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -45,+MAXLOG 200,000 1.2e-19 2.5e-20
*
* ERROR MESSAGES:
*
* message condition value returned
* expm1l overflow x > MAXLOG MAXNUM
*
*/
#include <math.h>
static const long double MAXLOGL = 1.1356523406294143949492E4L;
/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
-.5 ln 2 < x < .5 ln 2
Theoretical peak relative error = 3.4e-22 */
static const long double
P0 = -1.586135578666346600772998894928250240826E4L,
P1 = 2.642771505685952966904660652518429479531E3L,
P2 = -3.423199068835684263987132888286791620673E2L,
P3 = 1.800826371455042224581246202420972737840E1L,
P4 = -5.238523121205561042771939008061958820811E-1L,
Q0 = -9.516813471998079611319047060563358064497E4L,
Q1 = 3.964866271411091674556850458227710004570E4L,
Q2 = -7.207678383830091850230366618190187434796E3L,
Q3 = 7.206038318724600171970199625081491823079E2L,
Q4 = -4.002027679107076077238836622982900945173E1L,
/* Q5 = 1.000000000000000000000000000000000000000E0 */
/* C1 + C2 = ln 2 */
C1 = 6.93145751953125E-1L,
C2 = 1.428606820309417232121458176568075500134E-6L,
/* ln 2^-65 */
minarg = -4.5054566736396445112120088E1L;
static const long double huge = 0x1p10000L;
long double
expm1l(long double x)
{
long double px, qx, xx;
int k;
/* Overflow. */
if (x > MAXLOGL)
return (huge*huge); /* overflow */
if (x == 0.0)
return x;
/* Minimum value. */
if (x < minarg)
return -1.0L;
xx = C1 + C2;
/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
px = floorl (0.5 + x / xx);
k = px;
/* remainder times ln 2 */
x -= px * C1;
x -= px * C2;
/* Approximate exp(remainder ln 2). */
px = (((( P4 * x
+ P3) * x
+ P2) * x
+ P1) * x
+ P0) * x;
qx = (((( x
+ Q4) * x
+ Q3) * x
+ Q2) * x
+ Q1) * x
+ Q0;
xx = x * x;
qx = x + (0.5 * xx + xx * px / qx);
/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
We have qx = exp(remainder ln 2) - 1, so
exp(x) - 1 = 2^k (qx + 1) - 1 = 2^k qx + 2^k - 1. */
px = ldexpl(1.0L, k);
x = px * qx + (px - 1.0);
return x;
}

80
ld80/s_floorl.c Normal file
View file

@ -0,0 +1,80 @@
/* @(#)s_floor.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* floorl(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to floor(x).
*/
#include <math.h>
#include "math_private.h"
static const long double huge = 1.0e4930L;
long double
floorl(long double x)
{
int32_t i1,jj0;
u_int32_t i,j,se,i0,sx;
GET_LDOUBLE_WORDS(se,i0,i1,x);
sx = (se>>15)&1;
jj0 = (se&0x7fff)-0x3fff;
if(jj0<31) {
if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {
if(sx==0)
return 0.0L;
else if(((se&0x7fff)|i0|i1)!=0)
return -1.0L;
}
} else {
i = (0x7fffffff)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(sx) {
if (jj0>0 && (i0+(0x80000000>>jj0))>i0)
i0 += (0x80000000)>>jj0;
else
{
i = 0x7fffffff;
++se;
}
}
i0 &= (~i); i1=0;
}
}
} else if (jj0>62) {
if(jj0==0x4000) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = ((u_int32_t)(0xffffffff))>>(jj0-31);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
if(sx) {
if(jj0==31) i0+=1;
else {
j = i1+(1<<(63-jj0));
if(j<i1) i0 +=1 ; /* got a carry */
i1=j;
}
}
i1 &= (~i);
}
}
SET_LDOUBLE_WORDS(x,se,i0,i1);
return x;
}

191
ld80/s_log1pl.c Normal file
View file

@ -0,0 +1,191 @@
/* $OpenBSD: s_log1pl.c,v 1.3 2013/11/12 20:35:19 martynas Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/* log1pl.c
*
* Relative error logarithm
* Natural logarithm of 1+x, long double precision
*
*
*
* SYNOPSIS:
*
* long double x, y, log1pl();
*
* y = log1pl( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of 1+x.
*
* The argument 1+x is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z^3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20
*
* ERROR MESSAGES:
*
* log singularity: x-1 = 0; returns -INFINITY
* log domain: x-1 < 0; returns NAN
*/
#include <math.h>
#include "math_private.h"
/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 2.32e-20
*/
static long double P[] = {
4.5270000862445199635215E-5L,
4.9854102823193375972212E-1L,
6.5787325942061044846969E0L,
2.9911919328553073277375E1L,
6.0949667980987787057556E1L,
5.7112963590585538103336E1L,
2.0039553499201281259648E1L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
1.5062909083469192043167E1L,
8.3047565967967209469434E1L,
2.2176239823732856465394E2L,
3.0909872225312059774938E2L,
2.1642788614495947685003E2L,
6.0118660497603843919306E1L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
static const long double C1 = 6.9314575195312500000000E-1L;
static const long double C2 = 1.4286068203094172321215E-6L;
#define SQRTH 0.70710678118654752440L
long double
log1pl(long double xm1)
{
long double x, y, z;
int e;
if( isnan(xm1) )
return(xm1);
if( xm1 == INFINITY )
return(xm1);
if(xm1 == 0.0)
return(xm1);
x = xm1 + 1.0L;
/* Test for domain errors. */
if( x <= 0.0L )
{
if( x == 0.0L )
return( -INFINITY );
else
return( NAN );
}
/* Separate mantissa from exponent.
Use frexp so that denormal numbers will be handled properly. */
x = frexpl( x, &e );
/* logarithm using log(x) = z + z^3 P(z)/Q(z),
where z = 2(x-1)/x+1) */
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
z = z + e * C2;
z = z + x;
z = z + e * C1;
return( z );
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
if (e != 0)
x = 2.0 * x - 1.0L;
else
x = xm1;
}
else
{
if (e != 0)
x = x - 1.0L;
else
x = xm1;
}
z = x*x;
y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) );
y = y + e * C2;
z = y - 0.5 * z;
z = z + x;
z = z + e * C1;
return( z );
}

69
ld80/s_modfl.c Normal file
View file

@ -0,0 +1,69 @@
/* @(#)s_modf.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/*
* modfl(long double x, long double *iptr)
* return fraction part of x, and return x's integral part in *iptr.
* Method:
* Bit twiddling.
*
* Exception:
* No exception.
*/
#include <math.h>
#include "math_private.h"
static const long double one = 1.0;
long double
modfl(long double x, long double *iptr)
{
int32_t i0,i1,jj0;
u_int32_t i,se;
GET_LDOUBLE_WORDS(se,i0,i1,x);
jj0 = (se&0x7fff)-0x3fff; /* exponent of x */
if(jj0<32) { /* integer part in high x */
if(jj0<0) { /* |x|<1 */
SET_LDOUBLE_WORDS(*iptr,se&0x8000,0,0); /* *iptr = +-0 */
return x;
} else {
i = (0x7fffffff)>>jj0;
if(((i0&i)|i1)==0) { /* x is integral */
*iptr = x;
SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */
return x;
} else {
SET_LDOUBLE_WORDS(*iptr,se,i0&(~i),0);
return x - *iptr;
}
}
} else if (jj0>63) { /* no fraction part */
*iptr = x*one;
/* We must handle NaNs separately. */
if (jj0 == 0x4000 && ((i0 & 0x7fffffff) | i1))
return x*one;
SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */
return x;
} else { /* fraction part in low x */
i = ((u_int32_t)(0x7fffffff))>>(jj0-32);
if((i1&i)==0) { /* x is integral */
*iptr = x;
SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */
return x;
} else {
SET_LDOUBLE_WORDS(*iptr,se,i0,i1&(~i));
return x - *iptr;
}
}
}

90
ld80/s_nextafterl.c Normal file
View file

@ -0,0 +1,90 @@
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* IEEE functions
* nextafterl(x,y)
* return the next machine floating-point number of x in the
* direction toward y.
* Special cases:
*/
#include <math.h>
#include "math_private.h"
long double
nextafterl(long double x, long double y)
{
int32_t hx,hy,ix,iy;
u_int32_t lx,ly,esx,esy;
GET_LDOUBLE_WORDS(esx,hx,lx,x);
GET_LDOUBLE_WORDS(esy,hy,ly,y);
ix = esx&0x7fff; /* |x| */
iy = esy&0x7fff; /* |y| */
if (((ix==0x7fff)&&((hx&0x7fffffff|lx)!=0)) || /* x is nan */
((iy==0x7fff)&&((hy&0x7fffffff|ly)!=0))) /* y is nan */
return x+y;
if(x==y) return y; /* x=y, return y */
if((ix|hx|lx)==0) { /* x == 0 */
volatile long double u;
SET_LDOUBLE_WORDS(x,esy&0x8000,0,1);/* return +-minsubnormal */
u = x;
u = u * u; /* raise underflow flag */
return x;
}
if(esx<0x8000) { /* x > 0 */
if(ix>iy||((ix==iy) && (hx>hy||((hx==hy)&&(lx>ly))))) {
/* x > y, x -= ulp */
if(lx==0) {
if ((hx&0x7fffffff)==0) esx -= 1;
hx = (hx - 1) | (hx & 0x80000000);
}
lx -= 1;
} else { /* x < y, x += ulp */
lx += 1;
if(lx==0) {
hx = (hx + 1) | (hx & 0x80000000);
if ((hx&0x7fffffff)==0) esx += 1;
}
}
} else { /* x < 0 */
if(esy>=0||(ix>iy||((ix==iy)&&(hx>hy||((hx==hy)&&(lx>ly)))))){
/* x < y, x -= ulp */
if(lx==0) {
if ((hx&0x7fffffff)==0) esx -= 1;
hx = (hx - 1) | (hx & 0x80000000);
}
lx -= 1;
} else { /* x > y, x += ulp */
lx += 1;
if(lx==0) {
hx = (hx + 1) | (hx & 0x80000000);
if ((hx&0x7fffffff)==0) esx += 1;
}
}
}
esy = esx&0x7fff;
if(esy==0x7fff) return x+x; /* overflow */
if(esy==0) {
volatile long double u = x*x; /* underflow */
if(u==x) {
SET_LDOUBLE_WORDS(x,esx,hx,lx);
return x;
}
}
SET_LDOUBLE_WORDS(x,esx,hx,lx);
return x;
}
__strong_alias(nexttowardl, nextafterl);

86
ld80/s_nexttoward.c Normal file
View file

@ -0,0 +1,86 @@
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* IEEE functions
* nexttoward(x,y)
* return the next machine floating-point number of x in the
* direction toward y.
* Special cases:
*/
#include <math.h>
#include <float.h>
#include "math_private.h"
double
nexttoward(double x, long double y)
{
int32_t hx,ix,iy;
u_int32_t lx,hy,ly,esy;
EXTRACT_WORDS(hx,lx,x);
GET_LDOUBLE_WORDS(esy,hy,ly,y);
ix = hx&0x7fffffff; /* |x| */
iy = esy&0x7fff; /* |y| */
if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */
((iy>=0x7fff)&&(hy|ly)!=0)) /* y is nan */
return x+y;
if((long double) x==y) return y; /* x=y, return y */
if((ix|lx)==0) { /* x == 0 */
volatile double u;
INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
u = x;
u = u * u; /* raise underflow flag */
return x;
}
if(hx>=0) { /* x > 0 */
if (esy>=0x8000||((ix>>20)&0x7ff)>iy-0x3c00
|| (((ix>>20)&0x7ff)==iy-0x3c00
&& (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
|| (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
&& (lx<<11)>ly)))) { /* x > y, x -= ulp */
if(lx==0) hx -= 1;
lx -= 1;
} else { /* x < y, x += ulp */
lx += 1;
if(lx==0) hx += 1;
}
} else { /* x < 0 */
if (esy<0x8000||((ix>>20)&0x7ff)>iy-0x3c00
|| (((ix>>20)&0x7ff)==iy-0x3c00
&& (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
|| (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
&& (lx<<11)>ly)))) {/* x < y, x -= ulp */
if(lx==0) hx -= 1;
lx -= 1;
} else { /* x > y, x += ulp */
lx += 1;
if(lx==0) hx += 1;
}
}
hy = hx&0x7ff00000;
if(hy>=0x7ff00000) {
x = x+x; /* overflow */
return x;
}
if(hy<0x00100000) {
volatile double u = x*x; /* underflow */
if(u==x) {
INSERT_WORDS(x,hx,lx);
return x;
}
}
INSERT_WORDS(x,hx,lx);
return x;
}

67
ld80/s_nexttowardf.c Normal file
View file

@ -0,0 +1,67 @@
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
#include <math.h>
#include <float.h>
#include "math_private.h"
float
nexttowardf(float x, long double y)
{
int32_t hx,ix,iy;
u_int32_t hy,ly,esy;
GET_FLOAT_WORD(hx,x);
GET_LDOUBLE_WORDS(esy,hy,ly,y);
ix = hx&0x7fffffff; /* |x| */
iy = esy&0x7fff; /* |y| */
if((ix>0x7f800000) || /* x is nan */
(iy>=0x7fff&&((hy|ly)!=0))) /* y is nan */
return x+y;
if((long double) x==y) return y; /* x=y, return y */
if(ix==0) { /* x == 0 */
volatile float u;
SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
u = x;
u = u * u; /* raise underflow flag */
return x;
}
if(hx>=0) { /* x > 0 */
if(esy>=0x8000||((ix>>23)&0xff)>iy-0x3f80
|| (((ix>>23)&0xff)==iy-0x3f80
&& ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x > y, x -= ulp */
hx -= 1;
} else { /* x < y, x += ulp */
hx += 1;
}
} else { /* x < 0 */
if(esy<0x8000||((ix>>23)&0xff)>iy-0x3f80
|| (((ix>>23)&0xff)==iy-0x3f80
&& ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x < y, x -= ulp */
hx -= 1;
} else { /* x > y, x += ulp */
hx += 1;
}
}
hy = hx&0x7f800000;
if(hy>=0x7f800000) {
x = x+x; /* overflow */
return x;
}
if(hy<0x00800000) {
volatile float u = x*x; /* underflow */
}
SET_FLOAT_WORD(x,hx);
return x;
}

166
ld80/s_remquol.c Normal file
View file

@ -0,0 +1,166 @@
/* @(#)e_fmod.c 1.3 95/01/18 */
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/types.h>
#include <machine/ieee.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include "math_private.h"
#define BIAS (LDBL_MAX_EXP - 1)
/*
* These macros add and remove an explicit integer bit in front of the
* fractional mantissa, if the architecture doesn't have such a bit by
* default already.
*/
#ifdef LDBL_IMPLICIT_NBIT
#define LDBL_NBIT 0
#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
#define HFRAC_BITS EXT_FRACHBITS
#else
#define LDBL_NBIT 0x80000000
#define SET_NBIT(hx) (hx)
#define HFRAC_BITS (EXT_FRACHBITS - 1)
#endif
#define MANL_SHIFT (EXT_FRACLBITS - 1)
static const long double Zero[] = {0.0L, -0.0L};
/*
* Return the IEEE remainder and set *quo to the last n bits of the
* quotient, rounded to the nearest integer. We choose n=31 because
* we wind up computing all the integer bits of the quotient anyway as
* a side-effect of computing the remainder by the shift and subtract
* method. In practice, this is far more bits than are needed to use
* remquo in reduction algorithms.
*
* Assumptions:
* - The low part of the mantissa fits in a manl_t exactly.
* - The high part of the mantissa fits in an int64_t with enough room
* for an explicit integer bit in front of the fractional bits.
*/
long double
remquol(long double x, long double y, int *quo)
{
int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */
uint32_t hy;
uint32_t lx,ly,lz;
uint32_t esx, esy;
int ix,iy,n,q,sx,sxy;
GET_LDOUBLE_WORDS(esx,hx,lx,x);
GET_LDOUBLE_WORDS(esy,hy,ly,y);
sx = esx & 0x8000;
sxy = sx ^ (esy & 0x8000);
esx &= 0x7fff; /* |x| */
esy &= 0x7fff; /* |y| */
SET_LDOUBLE_EXP(x,esx);
SET_LDOUBLE_EXP(y,esy);
/* purge off exception values */
if((esy|hy|ly)==0 || /* y=0 */
(esx == BIAS + LDBL_MAX_EXP) || /* or x not finite */
(esy == BIAS + LDBL_MAX_EXP &&
((hy&~LDBL_NBIT)|ly)!=0)) /* or y is NaN */
return (x*y)/(x*y);
if(esx<=esy) {
if((esx<esy) ||
(hx<=hy &&
(hx<hy ||
lx<ly))) {
q = 0;
goto fixup; /* |x|<|y| return x or x-y */
}
if(hx==hy && lx==ly) {
*quo = 1;
return Zero[sx!=0]; /* |x|=|y| return x*0*/
}
}
/* determine ix = ilogb(x) */
if(esx == 0) { /* subnormal x */
x *= 0x1.0p512;
GET_LDOUBLE_WORDS(esx,hx,lx,x);
ix = esx - (BIAS + 512);
} else {
ix = esx - BIAS;
}
/* determine iy = ilogb(y) */
if(esy == 0) { /* subnormal y */
y *= 0x1.0p512;
GET_LDOUBLE_WORDS(esy,hy,ly,y);
iy = esy - (BIAS + 512);
} else {
iy = esy - BIAS;
}
/* set up {hx,lx}, {hy,ly} and align y to x */
hx = SET_NBIT(hx);
lx = SET_NBIT(lx);
/* fix point fmod */
n = ix - iy;
q = 0;
while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
else {hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; q++;}
q <<= 1;
}
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;q++;}
/* convert back to floating value and restore the sign */
if((hx|lx)==0) { /* return sign(x)*0 */
*quo = (sxy ? -q : q);
return Zero[sx!=0];
}
while(hx<(1ULL<<HFRAC_BITS)) { /* normalize x */
hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
iy -= 1;
}
if (iy < LDBL_MIN_EXP) {
esx = (iy + BIAS + 512) & 0x7fff;
SET_LDOUBLE_WORDS(x,esx,hx,lx);
x *= 0x1p-512;
GET_LDOUBLE_WORDS(esx,hx,lx,x);
} else {
esx = (iy + BIAS) & 0x7fff;
}
SET_LDOUBLE_WORDS(x,esx,hx,lx);
fixup:
y = fabsl(y);
if (y < LDBL_MIN * 2) {
if (x+x>y || (x+x==y && (q & 1))) {
q++;
x-=y;
}
} else if (x>0.5*y || (x==0.5*y && (q & 1))) {
q++;
x-=y;
}
GET_LDOUBLE_EXP(esx,x);
esx ^= sx;
SET_LDOUBLE_EXP(x,esx);
q &= 0x7fffffff;
*quo = (sxy ? -q : q);
return x;
}

79
ld80/s_tanhl.c Normal file
View file

@ -0,0 +1,79 @@
/* @(#)s_tanh.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* tanhl(x)
* Return the Hyperbolic Tangent of x
*
* Method :
* x -x
* e - e
* 0. tanhl(x) is defined to be -----------
* x -x
* e + e
* 1. reduce x to non-negative by tanhl(-x) = -tanhl(x).
* 2. 0 <= x <= 2**-55 : tanhl(x) := x*(one+x)
* -t
* 2**-55 < x <= 1 : tanhl(x) := -----; t = expm1l(-2x)
* t + 2
* 2
* 1 <= x <= 23.0 : tanhl(x) := 1- ----- ; t=expm1l(2x)
* t + 2
* 23.0 < x <= INF : tanhl(x) := 1.
*
* Special cases:
* tanhl(NaN) is NaN;
* only tanhl(0)=0 is exact for finite argument.
*/
#include <math.h>
#include "math_private.h"
static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
long double
tanhl(long double x)
{
long double t,z;
int32_t se;
u_int32_t jj0,jj1,ix;
/* High word of |x|. */
GET_LDOUBLE_WORDS(se,jj0,jj1,x);
ix = se&0x7fff;
/* x is INF or NaN */
if(ix==0x7fff) {
/* for NaN it's not important which branch: tanhl(NaN) = NaN */
if (se&0x8000) return one/x-one; /* tanhl(-inf)= -1; */
else return one/x+one; /* tanhl(+inf)=+1 */
}
/* |x| < 23 */
if (ix < 0x4003 || (ix == 0x4003 && jj0 < 0xb8000000u)) {/* |x|<23 */
if ((ix|jj0|jj1) == 0)
return x; /* x == +- 0 */
if (ix<0x3fc8) /* |x|<2**-55 */
return x*(one+tiny); /* tanh(small) = small */
if (ix>=0x3fff) { /* |x|>=1 */
t = expm1l(two*fabsl(x));
z = one - two/(t+two);
} else {
t = expm1l(-two*fabsl(x));
z= -t/(t+two);
}
/* |x| > 23, return +-1 */
} else {
z = one - tiny; /* raised inexact flag */
}
return (se&0x8000)? -z: z;
}

72
ld80/s_truncl.c Normal file
View file

@ -0,0 +1,72 @@
/*
* ====================================================
* 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.
* ====================================================
*
* From: @(#)s_floor.c 5.1 93/09/24
*/
/*
* truncl(x)
* Return x rounded toward 0 to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to truncl(x).
*/
#include <sys/types.h>
#include <machine/ieee.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include "math_private.h"
#ifdef LDBL_IMPLICIT_NBIT
#define MANH_SIZE (EXT_FRACHBITS + 1)
#else
#define MANH_SIZE EXT_FRACHBITS
#endif
static const long double huge = 1.0e300;
static const float zero[] = { 0.0, -0.0 };
long double
truncl(long double x)
{
int e, es;
uint32_t ix0, ix1;
GET_LDOUBLE_WORDS(es,ix0,ix1,x);
e = (es&0x7fff) - LDBL_MAX_EXP + 1;
if (e < MANH_SIZE - 1) {
if (e < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0)
return (zero[(es&0x8000)!=0]);
} else {
uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
if (((ix0 & m) | ix1) == 0)
return (x); /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
ix0 &= ~m;
ix1 = 0;
}
}
} else if (e < LDBL_MANT_DIG - 1) {
uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
if ((ix1 & m) == 0)
return (x); /* x is integral */
if (huge + x > 0.0) /* raise inexact flag */
ix1 &= ~m;
}
SET_LDOUBLE_WORDS(x,es,ix0,ix1);
return (x);
}