diff options
| author |   martynas <martynas@openbsd.org> | 2008-09-07 20:36:06 +0000 |
|---|---|---|
| committer | martynas <martynas@openbsd.org> | 2008-09-07 20:36:06 +0000 |
| commit | 7b36286a70b46b494e2bca4a889ab49ef62ba86c (patch) | |
| tree | b2b9870e3495ed88715e57b0d4a9888b16264fe0 /lib/libm/src | |
| parent | sparc now requires this bloated library to be -fPIC (diff) | |
| download | wireguard-openbsd-7b36286a70b46b494e2bca4a889ab49ef62ba86c.tar.xz wireguard-openbsd-7b36286a70b46b494e2bca4a889ab49ef62ba86c.zip | |
- replace dtoa w/ David's gdtoa, version 2008-03-15
- provide proper dtoa locks
- use the real strtof implementation
- add strtold, __hdtoa, __hldtoa
- add %a/%A support
- don't lose precision in printf, don't round to double anymore
- implement extended-precision versions of libc functions: fpclassify,
isnan, isinf, signbit, isnormal, isfinite, now that the ieee.h is
fixed
- separate vax versions of strtof, and __hdtoa
- add complex math support. added functions: cacos, casin, catan,
ccos, csin, ctan, cacosh, casinh, catanh, ccosh, csinh, ctanh, cexp,
clog, cabs, cpow, csqrt, carg, cimag, conj, cproj, creal, cacosf,
casinf, catanf, ccosf, csinf, ctanf, cacoshf, casinhf, catanhf,
ccoshf, csinhf, ctanhf, cexpf, clogf, cabsf, cpowf, csqrtf, cargf,
cimagf, conjf, cprojf, crealf
- add fdim, fmax, fmin
- add log2. (adapted implementation e_log.c. could be more acruate
& faster, but it's good enough for now)
- remove wrappers & cruft in libm, supposed to work-around mistakes
in SVID, etc.; use ieee versions. fixes issues in python 2.6 for
djm@
- make _digittoint static
- proper definitions for i386, and amd64 in ieee.h
- sh, powerpc don't really have extended-precision
- add missing definitions for mips64 (quad), m{6,8}k (96-bit) float.h
for LDBL_*
- merge lead to frac for m{6,8}k, for gdtoa to work properly
- add FRAC*BITS & EXT_TO_ARRAY32 definitions in ieee.h, for hdtoa&ldtoa
to use
- add EXT_IMPLICIT_NBIT definition, which indicates implicit
normalization bit
- add regression tests for libc: fpclassify and printf
- arith.h & gd_qnan.h definitions
- update ieee.h: hppa doesn't have quad-precision, hppa64 does
- add missing prototypes to gdtoaimp
- on 64-bit platforms make sure gdtoa doesn't use a long when it
really wants an int
- etc., what i may have forgotten...
- bump libm major, due to removed&changed symbols
- no libc bump, since this is riding on djm's libc major crank from
a day ago
discussed with / requested by / testing theo, sthen@, djm@, jsg@,
merdely@, jsing@, tedu@, brad@, jakemsr@, and others.
looks good to millert@
parts of the diff ok kettenis@
this commit does not include:
- man page changes
Diffstat (limited to 'lib/libm/src')
150 files changed, 3512 insertions, 3180 deletions
diff --git a/lib/libm/src/e_acos.c b/lib/libm/src/e_acos.c index a5301808768..2ae24586cc1 100644 --- a/lib/libm/src/e_acos.c +++ b/lib/libm/src/e_acos.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $"; #endif -/* __ieee754_acos(x) +/* acos(x) * Method : * acos(x) = pi/2 - asin(x) * acos(-x) = pi/2 + asin(x) @@ -35,7 +35,7 @@ static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $"; * if x is NaN, return x itself; * if |x|>1, return NaN with invalid signal. * - * Function needed: __ieee754_sqrt + * Function needed: sqrt */ #include "math.h" @@ -58,7 +58,7 @@ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ double -__ieee754_acos(double x) +acos(double x) { double z,p,q,r,w,s,c,df; int32_t hx,ix; @@ -84,13 +84,13 @@ __ieee754_acos(double x) z = (one+x)*0.5; p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - s = __ieee754_sqrt(z); + s = sqrt(z); r = p/q; w = r*s-pio2_lo; return pi - 2.0*(s+w); } else { /* x > 0.5 */ z = (one-x)*0.5; - s = __ieee754_sqrt(z); + s = sqrt(z); df = s; SET_LOW_WORD(df,0); c = (z-df*df)/(s+df); diff --git a/lib/libm/src/e_acosf.c b/lib/libm/src/e_acosf.c index 5407842186c..57f42b18f8f 100644 --- a/lib/libm/src/e_acosf.c +++ b/lib/libm/src/e_acosf.c @@ -37,7 +37,7 @@ qS3 = -6.8828397989e-01, /* 0xbf303361 */ qS4 = 7.7038154006e-02; /* 0x3d9dc62e */ float -__ieee754_acosf(float x) +acosf(float x) { float z,p,q,r,w,s,c,df; int32_t hx,ix; @@ -60,14 +60,14 @@ __ieee754_acosf(float x) z = (one+x)*(float)0.5; p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); - s = __ieee754_sqrtf(z); + s = sqrtf(z); r = p/q; w = r*s-pio2_lo; return pi - (float)2.0*(s+w); } else { /* x > 0.5 */ int32_t idf; z = (one-x)*(float)0.5; - s = __ieee754_sqrtf(z); + s = sqrtf(z); df = s; GET_FLOAT_WORD(idf,df); SET_FLOAT_WORD(df,idf&0xfffff000); diff --git a/lib/libm/src/e_acosh.c b/lib/libm/src/e_acosh.c index fbfb65e09d0..273b93d02a6 100644 --- a/lib/libm/src/e_acosh.c +++ b/lib/libm/src/e_acosh.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $"; #endif -/* __ieee754_acosh(x) +/* acosh(x) * Method : * Based on * acosh(x) = log [ x + sqrt(x*x-1) ] @@ -36,7 +36,7 @@ one = 1.0, ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ double -__ieee754_acosh(double x) +acosh(double x) { double t; int32_t hx; @@ -48,12 +48,12 @@ __ieee754_acosh(double x) if(hx >=0x7ff00000) { /* x is inf of NaN */ return x+x; } else - return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */ + return log(x)+ln2; /* acosh(huge)=log(2x) */ } else if(((hx-0x3ff00000)|lx)==0) { return 0.0; /* acosh(1) = 0 */ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ t=x*x; - return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one))); + return log(2.0*x-one/(x+sqrt(t-one))); } else { /* 1<x<2 */ t = x-one; return log1p(t+sqrt(2.0*t+t*t)); diff --git a/lib/libm/src/e_acoshf.c b/lib/libm/src/e_acoshf.c index 6e3b1e96870..515d33daf0c 100644 --- a/lib/libm/src/e_acoshf.c +++ b/lib/libm/src/e_acoshf.c @@ -25,7 +25,7 @@ one = 1.0, ln2 = 6.9314718246e-01; /* 0x3f317218 */ float -__ieee754_acoshf(float x) +acoshf(float x) { float t; int32_t hx; @@ -36,12 +36,12 @@ __ieee754_acoshf(float x) if(hx >=0x7f800000) { /* x is inf of NaN */ return x+x; } else - return __ieee754_logf(x)+ln2; /* acosh(huge)=log(2x) */ + return logf(x)+ln2; /* acosh(huge)=log(2x) */ } else if (hx==0x3f800000) { return 0.0; /* acosh(1) = 0 */ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ t=x*x; - return __ieee754_logf((float)2.0*x-one/(x+__ieee754_sqrtf(t-one))); + return logf((float)2.0*x-one/(x+sqrtf(t-one))); } else { /* 1<x<2 */ t = x-one; return log1pf(t+sqrtf((float)2.0*t+t*t)); diff --git a/lib/libm/src/e_asin.c b/lib/libm/src/e_asin.c index 0d358012280..4c732971e40 100644 --- a/lib/libm/src/e_asin.c +++ b/lib/libm/src/e_asin.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_asin.c,v 1.9 1995/05/12 04:57:22 jtc Exp $"; #endif -/* __ieee754_asin(x) +/* asin(x) * Method : * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... * we approximate asin(x) on [0,0.5] by @@ -67,7 +67,7 @@ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ double -__ieee754_asin(double x) +asin(double x) { double t,w,p,q,c,r,s; int32_t hx,ix; @@ -95,7 +95,7 @@ __ieee754_asin(double x) t = w*0.5; p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); - s = __ieee754_sqrt(t); + s = sqrt(t); if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ w = p/q; t = pio2_hi-(2.0*(s+s*w)-pio2_lo); diff --git a/lib/libm/src/e_asinf.c b/lib/libm/src/e_asinf.c index f4e780847ed..3f70199fbc9 100644 --- a/lib/libm/src/e_asinf.c +++ b/lib/libm/src/e_asinf.c @@ -39,7 +39,7 @@ qS3 = -6.8828397989e-01, /* 0xbf303361 */ qS4 = 7.7038154006e-02; /* 0x3d9dc62e */ float -__ieee754_asinf(float x) +asinf(float x) { float t,w,p,q,c,r,s; int32_t hx,ix; @@ -65,7 +65,7 @@ __ieee754_asinf(float x) t = w*(float)0.5; p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); - s = __ieee754_sqrtf(t); + s = sqrtf(t); if(ix>=0x3F79999A) { /* if |x| > 0.975 */ w = p/q; t = pio2_hi-((float)2.0*(s+s*w)-pio2_lo); diff --git a/lib/libm/src/e_atan2.c b/lib/libm/src/e_atan2.c index d469450efe0..f8a0330dcb8 100644 --- a/lib/libm/src/e_atan2.c +++ b/lib/libm/src/e_atan2.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $"; #endif -/* __ieee754_atan2(y,x) +/* atan2(y,x) * Method : * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). * 2. Reduce x to positive by (if x and y are unexceptional): @@ -53,7 +53,7 @@ pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ double -__ieee754_atan2(double y, double x) +atan2(double y, double x) { double z; int32_t k,m,hx,hy,ix,iy; diff --git a/lib/libm/src/e_atan2f.c b/lib/libm/src/e_atan2f.c index ead3b133972..0855bd4ffb5 100644 --- a/lib/libm/src/e_atan2f.c +++ b/lib/libm/src/e_atan2f.c @@ -29,7 +29,7 @@ pi = 3.1415925026e+00, /* 0x40490fda */ pi_lo = 1.5099578832e-07; /* 0x34222168 */ float -__ieee754_atan2f(float y, float x) +atan2f(float y, float x) { float z; int32_t k,m,hx,hy,ix,iy; diff --git a/lib/libm/src/e_atanh.c b/lib/libm/src/e_atanh.c index fd9d952b61f..1cd468b8e56 100644 --- a/lib/libm/src/e_atanh.c +++ b/lib/libm/src/e_atanh.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $"; #endif -/* __ieee754_atanh(x) +/* atanh(x) * Method : * 1.Reduced x to positive by atanh(-x) = -atanh(x) * 2.For x>=0.5 @@ -39,7 +39,7 @@ static const double one = 1.0, huge = 1e300; static const double zero = 0.0; double -__ieee754_atanh(double x) +atanh(double x) { double t; int32_t hx,ix; diff --git a/lib/libm/src/e_atanhf.c b/lib/libm/src/e_atanhf.c index 8d912239224..40322c2cf0c 100644 --- a/lib/libm/src/e_atanhf.c +++ b/lib/libm/src/e_atanhf.c @@ -24,7 +24,7 @@ static const float one = 1.0, huge = 1e30; static const float zero = 0.0; float -__ieee754_atanhf(float x) +atanhf(float x) { float t; int32_t hx,ix; diff --git a/lib/libm/src/e_cosh.c b/lib/libm/src/e_cosh.c index 4d7b61bfdc9..afad33162af 100644 --- a/lib/libm/src/e_cosh.c +++ b/lib/libm/src/e_cosh.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $"; #endif -/* __ieee754_cosh(x) +/* cosh(x) * Method : * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 * 1. Replace x by |x| (cosh(x) = cosh(-x)). @@ -41,7 +41,7 @@ static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $"; static const double one = 1.0, half=0.5, huge = 1.0e300; double -__ieee754_cosh(double x) +cosh(double x) { double t,w; int32_t ix; @@ -64,18 +64,18 @@ __ieee754_cosh(double x) /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ if (ix < 0x40360000) { - t = __ieee754_exp(fabs(x)); + t = exp(fabs(x)); return half*t+half/t; } /* |x| in [22, log(maxdouble)] return half*exp(|x|) */ - if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x)); + if (ix < 0x40862E42) return half*exp(fabs(x)); /* |x| in [log(maxdouble), overflowthresold] */ GET_LOW_WORD(lx,x); if (ix<0x408633CE || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) { - w = __ieee754_exp(half*fabs(x)); + w = exp(half*fabs(x)); t = half*w; return t*w; } diff --git a/lib/libm/src/e_coshf.c b/lib/libm/src/e_coshf.c index ae52f14d7a3..29eb56cfd4e 100644 --- a/lib/libm/src/e_coshf.c +++ b/lib/libm/src/e_coshf.c @@ -24,7 +24,7 @@ static const volatile float huge = 1.0e30; static const float one = 1.0, half=0.5; float -__ieee754_coshf(float x) +coshf(float x) { float t,w; int32_t ix; @@ -45,16 +45,16 @@ __ieee754_coshf(float x) /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ if (ix < 0x41b00000) { - t = __ieee754_expf(fabsf(x)); + t = expf(fabsf(x)); return half*t+half/t; } /* |x| in [22, log(maxdouble)] return half*exp(|x|) */ - if (ix < 0x42b17180) return half*__ieee754_expf(fabsf(x)); + if (ix < 0x42b17180) return half*expf(fabsf(x)); /* |x| in [log(maxdouble), overflowthresold] */ if (ix<=0x42b2d4fc) { - w = __ieee754_expf(half*fabsf(x)); + w = expf(half*fabsf(x)); t = half*w; return t*w; } diff --git a/lib/libm/src/e_exp.c b/lib/libm/src/e_exp.c index c9ecdecdc04..aef00ba006d 100644 --- a/lib/libm/src/e_exp.c +++ b/lib/libm/src/e_exp.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_exp.c,v 1.8 1995/05/10 20:45:03 jtc Exp $"; #endif -/* __ieee754_exp(x) +/* exp(x) * Returns the exponential of x. * * Method @@ -100,7 +100,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ double -__ieee754_exp(double x) /* default IEEE double exp */ +exp(double x) /* default IEEE double exp */ { double y,hi,lo,c,t; int32_t k,xsb; diff --git a/lib/libm/src/e_expf.c b/lib/libm/src/e_expf.c index 2457f286bb9..ac0fc40df09 100644 --- a/lib/libm/src/e_expf.c +++ b/lib/libm/src/e_expf.c @@ -40,7 +40,7 @@ P4 = -1.6533901999e-06, /* 0xb5ddea0e */ P5 = 4.1381369442e-08; /* 0x3331bb4c */ float -__ieee754_expf(float x) /* default IEEE double exp */ +expf(float x) /* default IEEE double exp */ { float y,hi,lo,c,t; int32_t k,xsb; diff --git a/lib/libm/src/e_fmod.c b/lib/libm/src/e_fmod.c index 377dee8e569..1a85f8295c1 100644 --- a/lib/libm/src/e_fmod.c +++ b/lib/libm/src/e_fmod.c @@ -15,7 +15,7 @@ static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $"; #endif /* - * __ieee754_fmod(x,y) + * fmod(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract */ @@ -26,7 +26,7 @@ static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $"; static const double one = 1.0, Zero[] = {0.0, -0.0,}; double -__ieee754_fmod(double x, double y) +fmod(double x, double y) { int32_t n,hx,hy,hz,ix,iy,sx,i; u_int32_t lx,ly,lz; diff --git a/lib/libm/src/e_fmodf.c b/lib/libm/src/e_fmodf.c index 5313a0efdf9..0bb8d172925 100644 --- a/lib/libm/src/e_fmodf.c +++ b/lib/libm/src/e_fmodf.c @@ -18,7 +18,7 @@ static char rcsid[] = "$NetBSD: e_fmodf.c,v 1.4 1995/05/10 20:45:10 jtc Exp $"; #endif /* - * __ieee754_fmodf(x,y) + * fmodf(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract */ @@ -29,7 +29,7 @@ static char rcsid[] = "$NetBSD: e_fmodf.c,v 1.4 1995/05/10 20:45:10 jtc Exp $"; static const float one = 1.0, Zero[] = {0.0, -0.0,}; float -__ieee754_fmodf(float x, float y) +fmodf(float x, float y) { int32_t n,hx,hy,hz,ix,iy,sx,i; diff --git a/lib/libm/src/e_hypot.c b/lib/libm/src/e_hypot.c index 6dc49042b63..23cf67f766d 100644 --- a/lib/libm/src/e_hypot.c +++ b/lib/libm/src/e_hypot.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; #endif -/* __ieee754_hypot(x,y) +/* hypot(x,y) * * Method : * If (assume round-to-nearest) z=x*x+y*y @@ -50,7 +50,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; #include "math_private.h" double -__ieee754_hypot(double x, double y) +hypot(double x, double y) { double a=x,b=y,t1,t2,yy1,y2,w; int32_t j,k,ha,hb; @@ -103,7 +103,7 @@ __ieee754_hypot(double x, double y) t1 = 0; SET_HIGH_WORD(t1,ha); t2 = a-t1; - w = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1))); + w = sqrt(t1*t1-(b*(-b)-t2*(a+t1))); } else { a = a+a; yy1 = 0; @@ -112,7 +112,7 @@ __ieee754_hypot(double x, double y) t1 = 0; SET_HIGH_WORD(t1,ha+0x00100000); t2 = a - t1; - w = __ieee754_sqrt(t1*yy1-(w*(-w)-(t1*y2+t2*b))); + w = sqrt(t1*yy1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { u_int32_t high; diff --git a/lib/libm/src/e_hypotf.c b/lib/libm/src/e_hypotf.c index db781f5dfdc..863899f2ef6 100644 --- a/lib/libm/src/e_hypotf.c +++ b/lib/libm/src/e_hypotf.c @@ -21,7 +21,7 @@ static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $"; #include "math_private.h" float -__ieee754_hypotf(float x, float y) +hypotf(float x, float y) { float a=x,b=y,t1,t2,yy1,y2,w; int32_t j,k,ha,hb; @@ -67,14 +67,14 @@ __ieee754_hypotf(float x, float y) if (w>b) { SET_FLOAT_WORD(t1,ha&0xfffff000); t2 = a-t1; - w = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1))); + w = sqrtf(t1*t1-(b*(-b)-t2*(a+t1))); } else { a = a+a; SET_FLOAT_WORD(yy1,hb&0xfffff000); y2 = b - yy1; SET_FLOAT_WORD(t1,ha+0x00800000); t2 = a - t1; - w = __ieee754_sqrtf(t1*yy1-(w*(-w)-(t1*y2+t2*b))); + w = sqrtf(t1*yy1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { SET_FLOAT_WORD(t1,0x3f800000+(k<<23)); diff --git a/lib/libm/src/e_j0.c b/lib/libm/src/e_j0.c index d49297ca8f3..ef34cb7887b 100644 --- a/lib/libm/src/e_j0.c +++ b/lib/libm/src/e_j0.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $"; #endif -/* __ieee754_j0(x), __ieee754_y0(x) +/* j0(x), y0(x) * Bessel function of the first and second kinds of order zero. * Method -- j0(x): * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ... @@ -82,7 +82,7 @@ S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */ static const double zero = 0.0; double -__ieee754_j0(double x) +j0(double x) { double z, s,c,ss,cc,r,u,v; int32_t hx,ix; @@ -143,7 +143,7 @@ v03 = 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */ v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */ double -__ieee754_y0(double x) +y0(double x) { double z, s,c,ss,cc,u,v; int32_t hx,ix,lx; @@ -187,12 +187,12 @@ __ieee754_y0(double x) return z; } if(ix<=0x3e400000) { /* x < 2**-27 */ - return(u00 + tpi*__ieee754_log(x)); + return(u00 + tpi*log(x)); } z = x*x; u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); v = one+z*(v01+z*(v02+z*(v03+z*v04))); - return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x))); + return(u/v + tpi*(j0(x)*log(x))); } /* The asymptotic expansions of pzero is diff --git a/lib/libm/src/e_j0f.c b/lib/libm/src/e_j0f.c index 02471711017..f5000adf5b6 100644 --- a/lib/libm/src/e_j0f.c +++ b/lib/libm/src/e_j0f.c @@ -40,7 +40,7 @@ S04 = 1.1661400734e-09; /* 0x30a045e8 */ static const float zero = 0.0; float -__ieee754_j0f(float x) +j0f(float x) { float z, s,c,ss,cc,r,u,v; int32_t hx,ix; @@ -101,7 +101,7 @@ v03 = 2.5915085189e-07, /* 0x348b216c */ v04 = 4.4111031494e-10; /* 0x2ff280c2 */ float -__ieee754_y0f(float x) +y0f(float x) { float z, s,c,ss,cc,u,v; int32_t hx,ix; @@ -145,12 +145,12 @@ __ieee754_y0f(float x) return z; } if(ix<=0x32000000) { /* x < 2**-27 */ - return(u00 + tpi*__ieee754_logf(x)); + return(u00 + tpi*logf(x)); } z = x*x; u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); v = one+z*(v01+z*(v02+z*(v03+z*v04))); - return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x))); + return(u/v + tpi*(j0f(x)*logf(x))); } /* The asymptotic expansions of pzero is diff --git a/lib/libm/src/e_j1.c b/lib/libm/src/e_j1.c index 372767ad582..3a783c80f08 100644 --- a/lib/libm/src/e_j1.c +++ b/lib/libm/src/e_j1.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $"; #endif -/* __ieee754_j1(x), __ieee754_y1(x) +/* j1(x), y1(x) * Bessel function of the first and second kinds of order zero. * Method -- j1(x): * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ... @@ -83,7 +83,7 @@ s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */ static const double zero = 0.0; double -__ieee754_j1(double x) +j1(double x) { double z, s,c,ss,cc,r,u,v,y; int32_t hx,ix; @@ -140,7 +140,7 @@ static const double V0[5] = { }; double -__ieee754_y1(double x) +y1(double x) { double z, s,c,ss,cc,u,v; int32_t hx,ix,lx; @@ -185,7 +185,7 @@ __ieee754_y1(double x) z = x*x; u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); - return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x)); + return(x*(u/v) + tpi*(j1(x)*log(x)-one/x)); } /* For x >= 8, the asymptotic expansions of pone is diff --git a/lib/libm/src/e_j1f.c b/lib/libm/src/e_j1f.c index bccc8f68ccc..62b34072ea8 100644 --- a/lib/libm/src/e_j1f.c +++ b/lib/libm/src/e_j1f.c @@ -41,7 +41,7 @@ s05 = 1.2354227016e-11; /* 0x2d59567e */ static const float zero = 0.0; float -__ieee754_j1f(float x) +j1f(float x) { float z, s,c,ss,cc,r,u,v,y; int32_t hx,ix; @@ -98,7 +98,7 @@ static const float V0[5] = { }; float -__ieee754_y1f(float x) +y1f(float x) { float z, s,c,ss,cc,u,v; int32_t hx,ix; @@ -143,7 +143,7 @@ __ieee754_y1f(float x) z = x*x; u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); - return(x*(u/v) + tpi*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x)); + return(x*(u/v) + tpi*(j1f(x)*logf(x)-one/x)); } /* For x >= 8, the asymptotic expansions of pone is diff --git a/lib/libm/src/e_jn.c b/lib/libm/src/e_jn.c index ada373aa6b4..f22231a8e97 100644 --- a/lib/libm/src/e_jn.c +++ b/lib/libm/src/e_jn.c @@ -15,7 +15,7 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $"; #endif /* - * __ieee754_jn(n, x), __ieee754_yn(n, x) + * jn(n, x), yn(n, x) * floating point Bessel's function of the 1st and 2nd kind * of order n * @@ -51,7 +51,7 @@ one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */ static const double zero = 0.00000000000000000000e+00; double -__ieee754_jn(int n, double x) +jn(int n, double x) { int32_t i,hx,ix,lx, sgn; double a, b, temp, di; @@ -69,8 +69,8 @@ __ieee754_jn(int n, double x) x = -x; hx ^= 0x80000000; } - if(n==0) return(__ieee754_j0(x)); - if(n==1) return(__ieee754_j1(x)); + if(n==0) return(j0(x)); + if(n==1) return(j1(x)); sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ x = fabs(x); if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */ @@ -99,8 +99,8 @@ __ieee754_jn(int n, double x) } b = invsqrtpi*temp/sqrt(x); } else { - a = __ieee754_j0(x); - b = __ieee754_j1(x); + a = j0(x); + b = j1(x); for(i=1;i<n;i++){ temp = b; b = b*((double)(i+i)/x) - a; /* avoid underflow */ @@ -176,7 +176,7 @@ __ieee754_jn(int n, double x) */ tmp = n; v = two/x; - tmp = tmp*__ieee754_log(fabs(v*tmp)); + tmp = tmp*log(fabs(v*tmp)); if(tmp<7.09782712893383973096e+02) { for(i=n-1,di=(double)(i+i);i>0;i--){ temp = b; @@ -200,14 +200,14 @@ __ieee754_jn(int n, double x) } } } - b = (t*__ieee754_j0(x)/b); + b = (t*j0(x)/b); } } if(sgn==1) return -b; else return b; } double -__ieee754_yn(int n, double x) +yn(int n, double x) { int32_t i,hx,ix,lx; int32_t sign; @@ -224,8 +224,8 @@ __ieee754_yn(int n, double x) n = -n; sign = 1 - ((n&1)<<1); } - if(n==0) return(__ieee754_y0(x)); - if(n==1) return(sign*__ieee754_y1(x)); + if(n==0) return(y0(x)); + if(n==1) return(sign*y1(x)); if(ix==0x7ff00000) return zero; if(ix>=0x52D00000) { /* x > 2**302 */ /* (x >> n**2) @@ -250,8 +250,8 @@ __ieee754_yn(int n, double x) b = invsqrtpi*temp/sqrt(x); } else { u_int32_t high; - a = __ieee754_y0(x); - b = __ieee754_y1(x); + a = y0(x); + b = y1(x); /* quit if b is -inf */ GET_HIGH_WORD(high,b); for(i=1;i<n&&high!=0xfff00000;i++){ diff --git a/lib/libm/src/e_jnf.c b/lib/libm/src/e_jnf.c index 3512c9f51a2..c96d9670e27 100644 --- a/lib/libm/src/e_jnf.c +++ b/lib/libm/src/e_jnf.c @@ -27,7 +27,7 @@ one = 1.0000000000e+00; /* 0x3F800000 */ static const float zero = 0.0000000000e+00; float -__ieee754_jnf(int n, float x) +jnf(int n, float x) { int32_t i,hx,ix, sgn; float a, b, temp, di; @@ -45,16 +45,16 @@ __ieee754_jnf(int n, float x) x = -x; hx ^= 0x80000000; } - if(n==0) return(__ieee754_j0f(x)); - if(n==1) return(__ieee754_j1f(x)); + if(n==0) return(j0f(x)); + if(n==1) return(j1f(x)); sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ x = fabsf(x); if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */ b = zero; else if((float)n<=x) { /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ - a = __ieee754_j0f(x); - b = __ieee754_j1f(x); + a = j0f(x); + b = j1f(x); for(i=1;i<n;i++){ temp = b; b = b*((float)(i+i)/x) - a; /* avoid underflow */ @@ -129,7 +129,7 @@ __ieee754_jnf(int n, float x) */ tmp = n; v = two/x; - tmp = tmp*__ieee754_logf(fabsf(v*tmp)); + tmp = tmp*logf(fabsf(v*tmp)); if(tmp<(float)8.8721679688e+01) { for(i=n-1,di=(float)(i+i);i>0;i--){ temp = b; @@ -153,14 +153,14 @@ __ieee754_jnf(int n, float x) } } } - b = (t*__ieee754_j0f(x)/b); + b = (t*j0f(x)/b); } } if(sgn==1) return -b; else return b; } float -__ieee754_ynf(int n, float x) +ynf(int n, float x) { int32_t i,hx,ix,ib; int32_t sign; @@ -177,12 +177,12 @@ __ieee754_ynf(int n, float x) n = -n; sign = 1 - ((n&1)<<1); } - if(n==0) return(__ieee754_y0f(x)); - if(n==1) return(sign*__ieee754_y1f(x)); + if(n==0) return(y0f(x)); + if(n==1) return(sign*y1f(x)); if(ix==0x7f800000) return zero; - a = __ieee754_y0f(x); - b = __ieee754_y1f(x); + a = y0f(x); + b = y1f(x); /* quit if b is -inf */ GET_FLOAT_WORD(ib,b); for(i=1;i<n&&ib!=0xff800000;i++){ diff --git a/lib/libm/src/e_lgamma_r.c b/lib/libm/src/e_lgamma_r.c index 1da673bc0c4..dbca95d5ea0 100644 --- a/lib/libm/src/e_lgamma_r.c +++ b/lib/libm/src/e_lgamma_r.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $"; #endif -/* __ieee754_lgamma_r(x, signgamp) +/* lgamma_r(x, signgamp) * Reentrant version of the logarithm of the Gamma function * with user provide pointer for the sign of Gamma(x). * @@ -201,7 +201,7 @@ sin_pi(double x) double -__ieee754_lgamma_r(double x, int *signgamp) +lgamma_r(double x, int *signgamp) { double t,y,z,nadj,p,p1,p2,p3,q,r,w; int i,hx,lx,ix; @@ -216,15 +216,15 @@ __ieee754_lgamma_r(double x, int *signgamp) if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { *signgamp = -1; - return -__ieee754_log(-x); - } else return -__ieee754_log(x); + return - log(-x); + } else return - log(x); } if(hx<0) { if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ return one/zero; t = sin_pi(x); if(t==zero) return one/zero; /* -integer */ - nadj = __ieee754_log(pi/fabs(t*x)); + nadj = log(pi/fabs(t*x)); if(t<zero) *signgamp = -1; x = -x; } @@ -234,7 +234,7 @@ __ieee754_lgamma_r(double x, int *signgamp) /* for x < 2.0 */ else if(ix<0x40000000) { if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */ - r = -__ieee754_log(x); + r = - log(x); if(ix>=0x3FE76944) {y = one-x; i= 0;} else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;} else {y = x; i=2;} @@ -279,18 +279,18 @@ __ieee754_lgamma_r(double x, int *signgamp) case 5: z *= (y+4.0); /* FALLTHRU */ case 4: z *= (y+3.0); /* FALLTHRU */ case 3: z *= (y+2.0); /* FALLTHRU */ - r += __ieee754_log(z); break; + r += log(z); break; } /* 8.0 <= x < 2**58 */ } else if (ix < 0x43900000) { - t = __ieee754_log(x); + t = log(x); z = one/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); r = (x-half)*(t-one)+w; } else /* 2**58 <= x <= inf */ - r = x*(__ieee754_log(x)-one); + r = x*(log(x)-one); if(hx<0) r = nadj - r; return r; } diff --git a/lib/libm/src/e_lgammaf_r.c b/lib/libm/src/e_lgammaf_r.c index 4f67523ecf9..1dfef89f9f9 100644 --- a/lib/libm/src/e_lgammaf_r.c +++ b/lib/libm/src/e_lgammaf_r.c @@ -137,7 +137,7 @@ sin_pif(float x) float -__ieee754_lgammaf_r(float x, int *signgamp) +lgammaf_r(float x, int *signgamp) { float t,y,z,nadj,p,p1,p2,p3,q,r,w; int i,hx,ix; @@ -152,15 +152,15 @@ __ieee754_lgammaf_r(float x, int *signgamp) if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { *signgamp = -1; - return -__ieee754_logf(-x); - } else return -__ieee754_logf(x); + return - logf(-x); + } else return - logf(x); } if(hx<0) { if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ return one/zero; t = sin_pif(x); if(t==zero) return one/zero; /* -integer */ - nadj = __ieee754_logf(pi/fabsf(t*x)); + nadj = logf(pi/fabsf(t*x)); if(t<zero) *signgamp = -1; x = -x; } @@ -170,7 +170,7 @@ __ieee754_lgammaf_r(float x, int *signgamp) /* for x < 2.0 */ else if(ix<0x40000000) { if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ - r = -__ieee754_logf(x); + r = - logf(x); if(ix>=0x3f3b4a20) {y = one-x; i= 0;} else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;} else {y = x; i=2;} @@ -215,18 +215,18 @@ __ieee754_lgammaf_r(float x, int *signgamp) case 5: z *= (y+(float)4.0); /* FALLTHRU */ case 4: z *= (y+(float)3.0); /* FALLTHRU */ case 3: z *= (y+(float)2.0); /* FALLTHRU */ - r += __ieee754_logf(z); break; + r += logf(z); break; } /* 8.0 <= x < 2**58 */ } else if (ix < 0x5c800000) { - t = __ieee754_logf(x); + t = logf(x); z = one/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); r = (x-half)*(t-one)+w; } else /* 2**58 <= x <= inf */ - r = x*(__ieee754_logf(x)-one); + r = x*(logf(x)-one); if(hx<0) r = nadj - r; return r; } diff --git a/lib/libm/src/e_log.c b/lib/libm/src/e_log.c index 968faea8ce5..052deba253c 100644 --- a/lib/libm/src/e_log.c +++ b/lib/libm/src/e_log.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $"; #endif -/* __ieee754_log(x) +/* log(x) * Return the logrithm of x * * Method : @@ -83,7 +83,7 @@ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ static const double zero = 0.0; double -__ieee754_log(double x) +log(double x) { double hfsq,f,s,z,R,w,t1,t2,dk; int32_t k,hx,i,j; diff --git a/lib/libm/src/e_log10.c b/lib/libm/src/e_log10.c index 7b2e09969da..ac8b87cfa03 100644 --- a/lib/libm/src/e_log10.c +++ b/lib/libm/src/e_log10.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_log10.c,v 1.9 1995/05/10 20:45:51 jtc Exp $"; #endif -/* __ieee754_log10(x) +/* log10(x) * Return the base 10 logarithm of x * * Method : @@ -59,7 +59,7 @@ log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ static const double zero = 0.0; double -__ieee754_log10(double x) +log10(double x) { double y,z; int32_t i,k,hx; @@ -81,6 +81,6 @@ __ieee754_log10(double x) hx = (hx&0x000fffff)|((0x3ff-i)<<20); y = (double)(k+i); SET_HIGH_WORD(x,hx); - z = y*log10_2lo + ivln10*__ieee754_log(x); + z = y*log10_2lo + ivln10*log(x); return z+y*log10_2hi; } diff --git a/lib/libm/src/e_log10f.c b/lib/libm/src/e_log10f.c index b491f166e38..d4777bf8409 100644 --- a/lib/libm/src/e_log10f.c +++ b/lib/libm/src/e_log10f.c @@ -29,7 +29,7 @@ log10_2lo = 7.9034151668e-07; /* 0x355427db */ static const float zero = 0.0; float -__ieee754_log10f(float x) +log10f(float x) { float y,z; int32_t i,k,hx; @@ -50,6 +50,6 @@ __ieee754_log10f(float x) hx = (hx&0x007fffff)|((0x7f-i)<<23); y = (float)(k+i); SET_FLOAT_WORD(x,hx); - z = y*log10_2lo + ivln10*__ieee754_logf(x); + z = y*log10_2lo + ivln10*logf(x); return z+y*log10_2hi; } diff --git a/lib/libm/src/e_log2.c b/lib/libm/src/e_log2.c new file mode 100644 index 00000000000..aed716670dc --- /dev/null +++ b/lib/libm/src/e_log2.c @@ -0,0 +1,74 @@ +/* @(#)e_log.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 "math.h" +#include "math_private.h" + +static const double +ln2 = 0.6931471805599452862268, +two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ +Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ +Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ +Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ +Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ +Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ +Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ +Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + +static const double zero = 0.0; + +double +log2(double x) +{ + double hfsq,f,s,z,R,w,t1,t2,dk; + int32_t k,hx,i,j; + u_int32_t lx; + + EXTRACT_WORDS(hx,lx,x); + + k=0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx&0x7fffffff)|lx)==0) + return -two54/zero; /* log(+-0)=-inf */ + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 54; x *= two54; /* subnormal number, scale up x */ + GET_HIGH_WORD(hx,x); + } + if (hx >= 0x7ff00000) return x+x; + k += (hx>>20)-1023; + hx &= 0x000fffff; + i = (hx+0x95f64)&0x100000; + SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ + k += (i>>20); + f = x-1.0; + dk = (double)k; + if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ + if (f==zero) + return (dk); + R = f*f*(0.5-0.33333333333333333*f); + return (dk-(R-f)/ln2); + } + s = f/(2.0+f); + z = s*s; + i = hx-0x6147a; + w = z*z; + j = 0x6b851-hx; + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + i |= j; + R = t2+t1; + if(i>0) { + hfsq=0.5*f*f; + return (dk-(hfsq-s*(hfsq+R)-f)/ln2); + } else + return (dk-((s*(f-R))-f)/ln2); +} diff --git a/lib/libm/src/e_log2f.c b/lib/libm/src/e_log2f.c new file mode 100644 index 00000000000..c454775feff --- /dev/null +++ b/lib/libm/src/e_log2f.c @@ -0,0 +1,76 @@ +/* e_logf.c -- float version of e_log.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * 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" + +static const float +ln2 = 0.6931471805599452862268, +two25 = 3.355443200e+07, /* 0x4c000000 */ +Lg1 = 6.6666668653e-01, /* 3F2AAAAB */ +Lg2 = 4.0000000596e-01, /* 3ECCCCCD */ +Lg3 = 2.8571429849e-01, /* 3E924925 */ +Lg4 = 2.2222198546e-01, /* 3E638E29 */ +Lg5 = 1.8183572590e-01, /* 3E3A3325 */ +Lg6 = 1.5313838422e-01, /* 3E1CD04F */ +Lg7 = 1.4798198640e-01; /* 3E178897 */ + +static const float zero = 0.0; + +float +log2f(float x) +{ + float hfsq,f,s,z,R,w,t1,t2,dk; + int32_t k,ix,i,j; + + GET_FLOAT_WORD(ix,x); + + k=0; + if (ix < 0x00800000) { /* x < 2**-126 */ + if ((ix&0x7fffffff)==0) + return -two25/zero; /* log(+-0)=-inf */ + if (ix<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 25; x *= two25; /* subnormal number, scale up x */ + GET_FLOAT_WORD(ix,x); + } + if (ix >= 0x7f800000) return x+x; + k += (ix>>23)-127; + ix &= 0x007fffff; + i = (ix+(0x95f64<<3))&0x800000; + SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */ + k += (i>>23); + dk = (float)k; + f = x-(float)1.0; + if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */ + if (f==zero) + return (dk); + R = f*f*((float)0.5-(float)0.33333333333333333*f); + return (dk-(R-f)/ln2); + } + s = f/((float)2.0+f); + z = s*s; + i = ix-(0x6147a<<3); + w = z*z; + j = (0x6b851<<3)-ix; + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + i |= j; + R = t2+t1; + if(i>0) { + hfsq=(float)0.5*f*f; + return (dk-(hfsq-s*(hfsq+R)-f)/ln2); + } else + return (dk-((s*(f-R))-f)/ln2); +} diff --git a/lib/libm/src/e_logf.c b/lib/libm/src/e_logf.c index 60d327a089f..9491341e6b2 100644 --- a/lib/libm/src/e_logf.c +++ b/lib/libm/src/e_logf.c @@ -35,7 +35,7 @@ Lg7 = 1.4798198640e-01; /* 3E178897 */ static const float zero = 0.0; float -__ieee754_logf(float x) +logf(float x) { float hfsq,f,s,z,R,w,t1,t2,dk; int32_t k,ix,i,j; diff --git a/lib/libm/src/e_pow.c b/lib/libm/src/e_pow.c index 4c4ad3bd112..2aa035eee83 100644 --- a/lib/libm/src/e_pow.c +++ b/lib/libm/src/e_pow.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $"; #endif -/* __ieee754_pow(x,y) return x**y +/* pow(x,y) return x**y * * n * Method: Let x = 2 * (1+f) @@ -96,7 +96,7 @@ ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ double -__ieee754_pow(double x, double y) +pow(double x, double y) { double z,ax,z_h,z_l,p_h,p_l; double yy1,t1,t2,r,s,t,u,v,w; @@ -152,7 +152,7 @@ __ieee754_pow(double x, double y) if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3fe00000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ - return __ieee754_sqrt(x); + return sqrt(x); } } diff --git a/lib/libm/src/e_powf.c b/lib/libm/src/e_powf.c index ec93660cc3f..9e520f486b9 100644 --- a/lib/libm/src/e_powf.c +++ b/lib/libm/src/e_powf.c @@ -54,7 +54,7 @@ ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ float -__ieee754_powf(float x, float y) +powf(float x, float y) { float z,ax,z_h,z_l,p_h,p_l; float yy1,t1,t2,r,s,t,u,v,w; @@ -103,7 +103,7 @@ __ieee754_powf(float x, float y) if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3f000000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ - return __ieee754_sqrtf(x); + return sqrtf(x); } ax = fabsf(x); diff --git a/lib/libm/src/e_remainder.c b/lib/libm/src/e_remainder.c index cd67c4410ff..623a13b141e 100644 --- a/lib/libm/src/e_remainder.c +++ b/lib/libm/src/e_remainder.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_remainder.c,v 1.8 1995/05/10 20:46:05 jtc Exp $"; #endif -/* __ieee754_remainder(x,p) +/* remainder(x,p) * Return : * returns x REM p = x - [x/p]*p as if in infinite * precise arithmetic, where [x/p] is the (infinite bit) @@ -30,7 +30,7 @@ static const double zero = 0.0; double -__ieee754_remainder(double x, double p) +remainder(double x, double p) { int32_t hx,hp; u_int32_t sx,lx,lp; @@ -50,7 +50,7 @@ __ieee754_remainder(double x, double p) return (x*p)/(x*p); - if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ + if (hp<=0x7fdfffff) x = fmod(x,p+p); /* now x < 2p */ if (((hx-hp)|(lx-lp))==0) return zero*x; x = fabs(x); p = fabs(p); diff --git a/lib/libm/src/e_remainderf.c b/lib/libm/src/e_remainderf.c index 36ec236732e..4405866914d 100644 --- a/lib/libm/src/e_remainderf.c +++ b/lib/libm/src/e_remainderf.c @@ -23,7 +23,7 @@ static char rcsid[] = "$NetBSD: e_remainderf.c,v 1.4 1995/05/10 20:46:08 jtc Exp static const float zero = 0.0; float -__ieee754_remainderf(float x, float p) +remainderf(float x, float p) { int32_t hx,hp; u_int32_t sx; @@ -42,7 +42,7 @@ __ieee754_remainderf(float x, float p) return (x*p)/(x*p); - if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */ + if (hp<=0x7effffff) x = fmodf(x,p+p); /* now x < 2p */ if ((hx-hp)==0) return zero*x; x = fabsf(x); p = fabsf(p); diff --git a/lib/libm/src/e_scalb.c b/lib/libm/src/e_scalb.c index e51f5447d60..d07808896f1 100644 --- a/lib/libm/src/e_scalb.c +++ b/lib/libm/src/e_scalb.c @@ -15,7 +15,7 @@ static char rcsid[] = "$NetBSD: e_scalb.c,v 1.6 1995/05/10 20:46:09 jtc Exp $"; #endif /* - * __ieee754_scalb(x, fn) is provide for + * scalb(x, fn) is provide for * passing various standard test suite. One * should use scalbn() instead. */ @@ -25,7 +25,7 @@ static char rcsid[] = "$NetBSD: e_scalb.c,v 1.6 1995/05/10 20:46:09 jtc Exp $"; #ifdef _SCALB_INT double -__ieee754_scalb(double x, int fn) +scalb(double x, int fn) { return scalbn(x, fn); } @@ -33,7 +33,7 @@ __ieee754_scalb(double x, int fn) #else double -__ieee754_scalb(double x, double fn) +scalb(double x, double fn) { if (isnan(x)||isnan(fn)) return x*fn; if (!finite(fn)) { diff --git a/lib/libm/src/e_scalbf.c b/lib/libm/src/e_scalbf.c index 58587f5d98c..e50179a1835 100644 --- a/lib/libm/src/e_scalbf.c +++ b/lib/libm/src/e_scalbf.c @@ -22,7 +22,7 @@ static char rcsid[] = "$NetBSD: e_scalbf.c,v 1.3 1995/05/10 20:46:12 jtc Exp $"; #ifdef _SCALB_INT float -__ieee754_scalbf(float x, int fn) +scalbf(float x, int fn) { return scalbnf(x,fn); } @@ -30,7 +30,7 @@ __ieee754_scalbf(float x, int fn) #else float -__ieee754_scalbf(float x, float fn) +scalbf(float x, float fn) { if (isnanf(x)||isnanf(fn)) return x*fn; if (!finitef(fn)) { diff --git a/lib/libm/src/e_sinh.c b/lib/libm/src/e_sinh.c index 03cb60b3ef3..68be7b5080f 100644 --- a/lib/libm/src/e_sinh.c +++ b/lib/libm/src/e_sinh.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; #endif -/* __ieee754_sinh(x) +/* sinh(x) * Method : * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 * 1. Replace x by |x| (sinh(-x) = -sinh(x)). @@ -38,7 +38,7 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; static const double one = 1.0, shuge = 1.0e307; double -__ieee754_sinh(double x) +sinh(double x) { double t,w,h; int32_t ix,jx; @@ -63,12 +63,12 @@ __ieee754_sinh(double x) } /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ - if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x)); + if (ix < 0x40862E42) return h*exp(fabs(x)); /* |x| in [log(maxdouble), overflowthresold] */ GET_LOW_WORD(lx,x); if (ix<0x408633CE || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) { - w = __ieee754_exp(0.5*fabs(x)); + w = exp(0.5*fabs(x)); t = h*w; return t*w; } diff --git a/lib/libm/src/e_sinhf.c b/lib/libm/src/e_sinhf.c index a4b27c77ff4..a593c488b95 100644 --- a/lib/libm/src/e_sinhf.c +++ b/lib/libm/src/e_sinhf.c @@ -23,7 +23,7 @@ static char rcsid[] = "$NetBSD: e_sinhf.c,v 1.4 1995/05/10 20:46:15 jtc Exp $"; static const float one = 1.0, shuge = 1.0e37; float -__ieee754_sinhf(float x) +sinhf(float x) { float t,w,h; int32_t ix,jx; @@ -46,11 +46,11 @@ __ieee754_sinhf(float x) } /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ - if (ix < 0x42b17180) return h*__ieee754_expf(fabsf(x)); + if (ix < 0x42b17180) return h*expf(fabsf(x)); /* |x| in [log(maxdouble), overflowthresold] */ if (ix<=0x42b2d4fc) { - w = __ieee754_expf((float)0.5*fabsf(x)); + w = expf((float)0.5*fabsf(x)); t = h*w; return t*w; } diff --git a/lib/libm/src/e_sqrt.c b/lib/libm/src/e_sqrt.c index 4cf344e064a..05355a33291 100644 --- a/lib/libm/src/e_sqrt.c +++ b/lib/libm/src/e_sqrt.c @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $"; #endif -/* __ieee754_sqrt(x) +/* sqrt(x) * Return correctly rounded sqrt. * ------------------------------------------ * | Use the hardware sqrt if you have one | @@ -90,7 +90,7 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $"; static const double one = 1.0, tiny=1.0e-300; double -__ieee754_sqrt(double x) +sqrt(double x) { double z; int32_t sign = (int)0x80000000; diff --git a/lib/libm/src/e_sqrtf.c b/lib/libm/src/e_sqrtf.c index f886fecc583..2dee4848572 100644 --- a/lib/libm/src/e_sqrtf.c +++ b/lib/libm/src/e_sqrtf.c @@ -23,7 +23,7 @@ static char rcsid[] = "$NetBSD: e_sqrtf.c,v 1.4 1995/05/10 20:46:19 jtc Exp $"; static const float one = 1.0, tiny=1.0e-30; float -__ieee754_sqrtf(float x) +sqrtf(float x) { float z; int32_t sign = (int)0x80000000; diff --git a/lib/libm/src/k_standard.c b/lib/libm/src/k_standard.c deleted file mode 100644 index 2ed40b44b87..00000000000 --- a/lib/libm/src/k_standard.c +++ /dev/null @@ -1,773 +0,0 @@ -/* @(#)k_standard.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: k_standard.c,v 1.6 1995/05/10 20:46:35 jtc Exp $"; -#endif - -#include "math.h" -#include "math_private.h" -#include <errno.h> -#include <stdio.h> /* fputs(), stderr */ - -#ifndef _USE_WRITE -#define WRITE2(u,v) fputs(u, stderr) -#else /* !defined(_USE_WRITE) */ -#include <unistd.h> /* write */ -#define WRITE2(u,v) write(2, u, v) -#endif /* !defined(_USE_WRITE) */ - -static const double zero = 0.0; /* used as const */ - -/* - * Standard conformance (non-IEEE) on exception cases. - * Mapping: - * 1 -- acos(|x|>1) - * 2 -- asin(|x|>1) - * 3 -- atan2(+-0,+-0) - * 4 -- hypot overflow - * 5 -- cosh overflow - * 6 -- exp overflow - * 7 -- exp underflow - * 8 -- y0(0) - * 9 -- y0(-ve) - * 10-- y1(0) - * 11-- y1(-ve) - * 12-- yn(0) - * 13-- yn(-ve) - * 14-- lgamma(finite) overflow - * 15-- lgamma(-integer) - * 16-- log(0) - * 17-- log(x<0) - * 18-- log10(0) - * 19-- log10(x<0) - * 20-- pow(0.0,0.0) - * 21-- pow(x,y) overflow - * 22-- pow(x,y) underflow - * 23-- pow(0,negative) - * 24-- pow(neg,non-integral) - * 25-- sinh(finite) overflow - * 26-- sqrt(negative) - * 27-- fmod(x,0) - * 28-- remainder(x,0) - * 29-- acosh(x<1) - * 30-- atanh(|x|>1) - * 31-- atanh(|x|=1) - * 32-- scalb overflow - * 33-- scalb underflow - * 34-- j0(|x|>X_TLOSS) - * 35-- y0(x>X_TLOSS) - * 36-- j1(|x|>X_TLOSS) - * 37-- y1(x>X_TLOSS) - * 38-- jn(|x|>X_TLOSS, n) - * 39-- yn(x>X_TLOSS, n) - * 40-- gamma(finite) overflow - * 41-- gamma(-integer) - * 42-- pow(NaN,0.0) - */ - - -double -__kernel_standard(double x, double y, int type) -{ - struct exception exc; -#ifndef HUGE_VAL /* this is the only routine that uses HUGE_VAL */ -#define HUGE_VAL inf - double inf = 0.0; - - SET_HIGH_WORD(inf,0x7ff00000); /* set inf to infinite */ -#endif - -#ifdef _USE_WRITE - (void) fflush(stdout); -#endif - exc.arg1 = x; - exc.arg2 = y; - switch(type) { - case 1: - case 101: - /* acos(|x|>1) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "acos" : "acosf"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if(_LIB_VERSION == _SVID_) { - (void) WRITE2("acos: DOMAIN error\n", 19); - } - errno = EDOM; - } - break; - case 2: - case 102: - /* asin(|x|>1) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "asin" : "asinf"; - exc.retval = zero; - if(_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if(_LIB_VERSION == _SVID_) { - (void) WRITE2("asin: DOMAIN error\n", 19); - } - errno = EDOM; - } - break; - case 3: - case 103: - /* atan2(+-0,+-0) */ - exc.arg1 = y; - exc.arg2 = x; - exc.type = DOMAIN; - exc.name = type < 100 ? "atan2" : "atan2f"; - exc.retval = zero; - if(_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if(_LIB_VERSION == _SVID_) { - (void) WRITE2("atan2: DOMAIN error\n", 20); - } - errno = EDOM; - } - break; - case 4: - case 104: - /* hypot(finite,finite) overflow */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "hypot" : "hypotf"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 5: - case 105: - /* cosh(finite) overflow */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "cosh" : "coshf"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 6: - case 106: - /* exp(finite) overflow */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "exp" : "expf"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 7: - case 107: - /* exp(finite) underflow */ - exc.type = UNDERFLOW; - exc.name = type < 100 ? "exp" : "expf"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 8: - case 108: - /* y0(0) = -inf */ - exc.type = DOMAIN; /* should be SING for IEEE */ - exc.name = type < 100 ? "y0" : "y0f"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("y0: DOMAIN error\n", 17); - } - errno = EDOM; - } - break; - case 9: - case 109: - /* y0(x<0) = NaN */ - exc.type = DOMAIN; - exc.name = type < 100 ? "y0" : "y0f"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("y0: DOMAIN error\n", 17); - } - errno = EDOM; - } - break; - case 10: - case 110: - /* y1(0) = -inf */ - exc.type = DOMAIN; /* should be SING for IEEE */ - exc.name = type < 100 ? "y1" : "y1f"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("y1: DOMAIN error\n", 17); - } - errno = EDOM; - } - break; - case 11: - case 111: - /* y1(x<0) = NaN */ - exc.type = DOMAIN; - exc.name = type < 100 ? "y1" : "y1f"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("y1: DOMAIN error\n", 17); - } - errno = EDOM; - } - break; - case 12: - case 112: - /* yn(n,0) = -inf */ - exc.type = DOMAIN; /* should be SING for IEEE */ - exc.name = type < 100 ? "yn" : "ynf"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("yn: DOMAIN error\n", 17); - } - errno = EDOM; - } - break; - case 13: - case 113: - /* yn(x<0) = NaN */ - exc.type = DOMAIN; - exc.name = type < 100 ? "yn" : "ynf"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("yn: DOMAIN error\n", 17); - } - errno = EDOM; - } - break; - case 14: - case 114: - /* lgamma(finite) overflow */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "lgamma" : "lgammaf"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 15: - case 115: - /* lgamma(-integer) or lgamma(0) */ - exc.type = SING; - exc.name = type < 100 ? "lgamma" : "lgammaf"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("lgamma: SING error\n", 19); - } - errno = EDOM; - } - break; - case 16: - case 116: - /* log(0) */ - exc.type = SING; - exc.name = type < 100 ? "log" : "logf"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("log: SING error\n", 16); - } - errno = EDOM; - } - break; - case 17: - case 117: - /* log(x<0) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "log" : "logf"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("log: DOMAIN error\n", 18); - } - errno = EDOM; - } - break; - case 18: - case 118: - /* log10(0) */ - exc.type = SING; - exc.name = type < 100 ? "log10" : "log10f"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("log10: SING error\n", 18); - } - errno = EDOM; - } - break; - case 19: - case 119: - /* log10(x<0) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "log10" : "log10f"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("log10: DOMAIN error\n", 20); - } - errno = EDOM; - } - break; - case 20: - case 120: - /* pow(0.0,0.0) */ - /* error only if _LIB_VERSION == _SVID_ */ - exc.type = DOMAIN; - exc.name = type < 100 ? "pow" : "powf"; - exc.retval = zero; - if (_LIB_VERSION != _SVID_) exc.retval = 1.0; - else if (!matherr(&exc)) { - (void) WRITE2("pow(0,0): DOMAIN error\n", 23); - errno = EDOM; - } - break; - case 21: - case 121: - /* pow(x,y) overflow */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "pow" : "powf"; - if (_LIB_VERSION == _SVID_) { - exc.retval = HUGE; - y *= 0.5; - if(x<zero&&rint(y)!=y) exc.retval = -HUGE; - } else { - exc.retval = HUGE_VAL; - y *= 0.5; - if(x<zero&&rint(y)!=y) exc.retval = -HUGE_VAL; - } - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 22: - case 122: - /* pow(x,y) underflow */ - exc.type = UNDERFLOW; - exc.name = type < 100 ? "pow" : "powf"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 23: - case 123: - /* 0**neg */ - exc.type = DOMAIN; - exc.name = type < 100 ? "pow" : "powf"; - if (_LIB_VERSION == _SVID_) - exc.retval = zero; - else - exc.retval = -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("pow(0,neg): DOMAIN error\n", 25); - } - errno = EDOM; - } - break; - case 24: - case 124: - /* neg**non-integral */ - exc.type = DOMAIN; - exc.name = type < 100 ? "pow" : "powf"; - if (_LIB_VERSION == _SVID_) - exc.retval = zero; - else - exc.retval = zero/zero; /* X/Open allow NaN */ - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("neg**non-integral: DOMAIN error\n", 32); - } - errno = EDOM; - } - break; - case 25: - case 125: - /* sinh(finite) overflow */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "sinh" : "sinhf"; - if (_LIB_VERSION == _SVID_) - exc.retval = ( (x>zero) ? HUGE : -HUGE); - else - exc.retval = ( (x>zero) ? HUGE_VAL : -HUGE_VAL); - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 26: - case 126: - /* sqrt(x<0) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "sqrt" : "sqrtf"; - if (_LIB_VERSION == _SVID_) - exc.retval = zero; - else - exc.retval = zero/zero; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("sqrt: DOMAIN error\n", 19); - } - errno = EDOM; - } - break; - case 27: - case 127: - /* fmod(x,0) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "fmod" : "fmodf"; - if (_LIB_VERSION == _SVID_) - exc.retval = x; - else - exc.retval = zero/zero; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("fmod: DOMAIN error\n", 20); - } - errno = EDOM; - } - break; - case 28: - case 128: - /* remainder(x,0) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "remainder" : "remainderf"; - exc.retval = zero/zero; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("remainder: DOMAIN error\n", 24); - } - errno = EDOM; - } - break; - case 29: - case 129: - /* acosh(x<1) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "acosh" : "acoshf"; - exc.retval = zero/zero; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("acosh: DOMAIN error\n", 20); - } - errno = EDOM; - } - break; - case 30: - case 130: - /* atanh(|x|>1) */ - exc.type = DOMAIN; - exc.name = type < 100 ? "atanh" : "atanhf"; - exc.retval = zero/zero; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("atanh: DOMAIN error\n", 20); - } - errno = EDOM; - } - break; - case 31: - case 131: - /* atanh(|x|=1) */ - exc.type = SING; - exc.name = type < 100 ? "atanh" : "atanhf"; - exc.retval = x/zero; /* sign(x)*inf */ - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("atanh: SING error\n", 18); - } - errno = EDOM; - } - break; - case 32: - case 132: - /* scalb overflow; SVID also returns +-HUGE_VAL */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "scalb" : "scalbf"; - exc.retval = x > zero ? HUGE_VAL : -HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 33: - case 133: - /* scalb underflow */ - exc.type = UNDERFLOW; - exc.name = type < 100 ? "scalb" : "scalbf"; - exc.retval = copysign(zero,x); - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 34: - case 134: - /* j0(|x|>X_TLOSS) */ - exc.type = TLOSS; - exc.name = type < 100 ? "j0" : "j0f"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2(exc.name, 2); - (void) WRITE2(": TLOSS error\n", 14); - } - errno = ERANGE; - } - break; - case 35: - case 135: - /* y0(x>X_TLOSS) */ - exc.type = TLOSS; - exc.name = type < 100 ? "y0" : "y0f"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2(exc.name, 2); - (void) WRITE2(": TLOSS error\n", 14); - } - errno = ERANGE; - } - break; - case 36: - case 136: - /* j1(|x|>X_TLOSS) */ - exc.type = TLOSS; - exc.name = type < 100 ? "j1" : "j1f"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2(exc.name, 2); - (void) WRITE2(": TLOSS error\n", 14); - } - errno = ERANGE; - } - break; - case 37: - case 137: - /* y1(x>X_TLOSS) */ - exc.type = TLOSS; - exc.name = type < 100 ? "y1" : "y1f"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2(exc.name, 2); - (void) WRITE2(": TLOSS error\n", 14); - } - errno = ERANGE; - } - break; - case 38: - case 138: - /* jn(|x|>X_TLOSS) */ - exc.type = TLOSS; - exc.name = type < 100 ? "jn" : "jnf"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2(exc.name, 2); - (void) WRITE2(": TLOSS error\n", 14); - } - errno = ERANGE; - } - break; - case 39: - case 139: - /* yn(x>X_TLOSS) */ - exc.type = TLOSS; - exc.name = type < 100 ? "yn" : "ynf"; - exc.retval = zero; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2(exc.name, 2); - (void) WRITE2(": TLOSS error\n", 14); - } - errno = ERANGE; - } - break; - case 40: - case 140: - /* gamma(finite) overflow */ - exc.type = OVERFLOW; - exc.name = type < 100 ? "gamma" : "gammaf"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = ERANGE; - else if (!matherr(&exc)) { - errno = ERANGE; - } - break; - case 41: - case 141: - /* gamma(-integer) or gamma(0) */ - exc.type = SING; - exc.name = type < 100 ? "gamma" : "gammaf"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - if (_LIB_VERSION == _POSIX_) - errno = EDOM; - else if (!matherr(&exc)) { - if (_LIB_VERSION == _SVID_) { - (void) WRITE2("gamma: SING error\n", 18); - } - errno = EDOM; - } - break; - case 42: - case 142: - /* pow(NaN,0.0) */ - /* error only if _LIB_VERSION == _SVID_ & _XOPEN_ */ - exc.type = DOMAIN; - exc.name = type < 100 ? "pow" : "powf"; - exc.retval = x; - if (_LIB_VERSION == _IEEE_ || - _LIB_VERSION == _POSIX_) exc.retval = 1.0; - else if (!matherr(&exc)) { - errno = EDOM; - } - break; - } - return exc.retval; -} diff --git a/lib/libm/src/math_private.h b/lib/libm/src/math_private.h index cdfd07affe1..0fc4f6a78d3 100644 --- a/lib/libm/src/math_private.h +++ b/lib/libm/src/math_private.h @@ -1,4 +1,4 @@ -/* $OpenBSD: math_private.h,v 1.9 2008/07/24 09:40:16 martynas Exp $ */ +/* $OpenBSD: math_private.h,v 1.10 2008/09/07 20:36:09 martynas Exp $ */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -213,73 +213,15 @@ do { \ #endif #endif -/* ieee style elementary functions */ -extern double __ieee754_sqrt(double); -extern double __ieee754_acos(double); -extern double __ieee754_acosh(double); -extern double __ieee754_log(double); -extern double __ieee754_atanh(double); -extern double __ieee754_asin(double); -extern double __ieee754_atan2(double,double); -extern double __ieee754_exp(double); -extern double __ieee754_cosh(double); -extern double __ieee754_fmod(double,double); -extern double __ieee754_pow(double,double); -extern double __ieee754_lgamma_r(double,int *); -extern double __ieee754_gamma_r(double,int *); -extern double __ieee754_lgamma(double); -extern double __ieee754_gamma(double); -extern double __ieee754_log10(double); -extern double __ieee754_sinh(double); -extern double __ieee754_hypot(double,double); -extern double __ieee754_j0(double); -extern double __ieee754_j1(double); -extern double __ieee754_y0(double); -extern double __ieee754_y1(double); -extern double __ieee754_jn(int,double); -extern double __ieee754_yn(int,double); -extern double __ieee754_remainder(double,double); -extern int __ieee754_rem_pio2(double,double*); -extern double __ieee754_scalb(double,double); - /* fdlibm kernel function */ -extern double __kernel_standard(double,double,int); +extern int __ieee754_rem_pio2(double,double*); extern double __kernel_sin(double,double,int); extern double __kernel_cos(double,double); extern double __kernel_tan(double,double,int); extern int __kernel_rem_pio2(double*,double*,int,int,int,const int*); - -/* ieee style elementary float functions */ -extern float __ieee754_sqrtf(float); -extern float __ieee754_acosf(float); -extern float __ieee754_acoshf(float); -extern float __ieee754_logf(float); -extern float __ieee754_atanhf(float); -extern float __ieee754_asinf(float); -extern float __ieee754_atan2f(float,float); -extern float __ieee754_expf(float); -extern float __ieee754_coshf(float); -extern float __ieee754_fmodf(float,float); -extern float __ieee754_powf(float,float); -extern float __ieee754_lgammaf_r(float,int *); -extern float __ieee754_gammaf_r(float,int *); -extern float __ieee754_lgammaf(float); -extern float __ieee754_gammaf(float); -extern float __ieee754_log10f(float); -extern float __ieee754_sinhf(float); -extern float __ieee754_hypotf(float,float); -extern float __ieee754_j0f(float); -extern float __ieee754_j1f(float); -extern float __ieee754_y0f(float); -extern float __ieee754_y1f(float); -extern float __ieee754_jnf(int,float); -extern float __ieee754_ynf(int,float); -extern float __ieee754_remainderf(float,float); -extern int __ieee754_rem_pio2f(float,float*); -extern float __ieee754_scalbf(float,float); - /* float versions of fdlibm kernel functions */ +extern int __ieee754_rem_pio2f(float,float*); extern float __kernel_sinf(float,float,int); extern float __kernel_cosf(float,float); extern float __kernel_tanf(float,float,int); diff --git a/lib/libm/src/s_asinh.c b/lib/libm/src/s_asinh.c index da49239a015..7398aa9d057 100644 --- a/lib/libm/src/s_asinh.c +++ b/lib/libm/src/s_asinh.c @@ -45,13 +45,13 @@ asinh(double x) if(huge+x>one) return x; /* return x inexact except 0 */ } if(ix>0x41b00000) { /* |x| > 2**28 */ - w = __ieee754_log(fabs(x))+ln2; + w = log(fabs(x))+ln2; } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ t = fabs(x); - w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t)); + w = log(2.0*t+one/(sqrt(x*x+one)+t)); } else { /* 2.0 > |x| > 2**-28 */ t = x*x; - w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t))); + w =log1p(fabs(x)+t/(one+sqrt(one+t))); } if(hx>0) return w; else return -w; } diff --git a/lib/libm/src/s_asinhf.c b/lib/libm/src/s_asinhf.c index 50e60b4d29e..e6426142985 100644 --- a/lib/libm/src/s_asinhf.c +++ b/lib/libm/src/s_asinhf.c @@ -37,13 +37,13 @@ asinhf(float x) if(huge+x>one) return x; /* return x inexact except 0 */ } if(ix>0x4d800000) { /* |x| > 2**28 */ - w = __ieee754_logf(fabsf(x))+ln2; + w = logf(fabsf(x))+ln2; } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ t = fabsf(x); - w = __ieee754_logf((float)2.0*t+one/(__ieee754_sqrtf(x*x+one)+t)); + w = logf((float)2.0*t+one/(sqrtf(x*x+one)+t)); } else { /* 2.0 > |x| > 2**-28 */ t = x*x; - w =log1pf(fabsf(x)+t/(one+__ieee754_sqrtf(one+t))); + w =log1pf(fabsf(x)+t/(one+sqrtf(one+t))); } if(hx>0) return w; else return -w; } diff --git a/lib/libm/src/s_cabs.c b/lib/libm/src/s_cabs.c new file mode 100644 index 00000000000..d2f3f942c3f --- /dev/null +++ b/lib/libm/src/s_cabs.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_cabs.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +double +cabs(double complex z) +{ + return hypot(__real__ z, __imag__ z); +} diff --git a/lib/libm/src/s_cabsf.c b/lib/libm/src/s_cabsf.c new file mode 100644 index 00000000000..d4d22f2dbb5 --- /dev/null +++ b/lib/libm/src/s_cabsf.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_cabsf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +float +cabsf(float complex z) +{ + return hypotf(__real__ z, __imag__ z); +} diff --git a/lib/libm/src/s_cacos.c b/lib/libm/src/s_cacos.c new file mode 100644 index 00000000000..96d26626b98 --- /dev/null +++ b/lib/libm/src/s_cacos.c @@ -0,0 +1,60 @@ +/* $OpenBSD: s_cacos.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cacos() + * + * Complex circular arc cosine + * + * + * + * SYNOPSIS: + * + * double complex cacos(); + * double complex z, w; + * + * w = cacos (z); + * + * + * + * DESCRIPTION: + * + * + * w = arccos z = PI/2 - arcsin z. + * + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5200 1.6e-15 2.8e-16 + * IEEE -10,+10 30000 1.8e-14 2.2e-15 + */ + +#include <complex.h> +#include <math.h> + +double complex +cacos(double complex z) +{ + double complex w; + + w = casin (z); + w = (M_PI_2 - creal (w)) - cimag (w) * I; + return (w); +} diff --git a/lib/libm/src/s_cacosf.c b/lib/libm/src/s_cacosf.c new file mode 100644 index 00000000000..5a644c15bd1 --- /dev/null +++ b/lib/libm/src/s_cacosf.c @@ -0,0 +1,60 @@ +/* $OpenBSD: s_cacosf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cacosf() + * + * Complex circular arc cosine + * + * + * + * SYNOPSIS: + * + * void cacosf(); + * cmplxf z, w; + * + * cacosf( &z, &w ); + * + * + * + * DESCRIPTION: + * + * + * w = arccos z = PI/2 - arcsin z. + * + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 9.2e-6 1.2e-6 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +cacosf(float complex z) +{ + float complex w; + + w = casinf( z ); + w = (M_PI_2 - creal (w)) - cimag (w) * I; + return (w); +} diff --git a/lib/libm/src/s_cacosh.c b/lib/libm/src/s_cacosh.c new file mode 100644 index 00000000000..43ba99ad6f9 --- /dev/null +++ b/lib/libm/src/s_cacosh.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_cacosh.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cacosh + * + * Complex inverse hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * double complex cacosh(); + * double complex z, w; + * + * w = cacosh (z); + * + * + * + * DESCRIPTION: + * + * acosh z = i acos z . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.6e-14 2.1e-15 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +cacosh(double complex z) +{ + double complex w; + + w = I * cacos (z); + return (w); +} diff --git a/lib/libm/src/s_cacoshf.c b/lib/libm/src/s_cacoshf.c new file mode 100644 index 00000000000..d81f7843c91 --- /dev/null +++ b/lib/libm/src/s_cacoshf.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_cacoshf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cacoshf + * + * Complex inverse hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * float complex cacoshf(); + * float complex z, w; + * + * w = cacoshf (z); + * + * + * + * DESCRIPTION: + * + * acosh z = i acos z . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.6e-14 2.1e-15 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +cacoshf(float complex z) +{ + float complex w; + + w = I * cacosf (z); + return (w); +} diff --git a/lib/libm/src/s_carg.c b/lib/libm/src/s_carg.c new file mode 100644 index 00000000000..db51662a1d9 --- /dev/null +++ b/lib/libm/src/s_carg.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_carg.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +double +carg(double complex z) +{ + return atan2 (__imag__ z, __real__ z); +} diff --git a/lib/libm/src/s_cargf.c b/lib/libm/src/s_cargf.c new file mode 100644 index 00000000000..12dbcde753d --- /dev/null +++ b/lib/libm/src/s_cargf.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_cargf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +float +cargf(float complex z) +{ + return atan2f (__imag__ z, __real__ z); +} diff --git a/lib/libm/src/s_casin.c b/lib/libm/src/s_casin.c new file mode 100644 index 00000000000..1c1b53541bb --- /dev/null +++ b/lib/libm/src/s_casin.c @@ -0,0 +1,129 @@ +/* $OpenBSD: s_casin.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* casin() + * + * Complex circular arc sine + * + * + * + * SYNOPSIS: + * + * double complex casin(); + * double complex z, w; + * + * w = casin (z); + * + * + * + * DESCRIPTION: + * + * Inverse complex sine: + * + * 2 + * w = -i clog( iz + csqrt( 1 - z ) ). + * + * casin(z) = -i casinh(iz) + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 10100 2.1e-15 3.4e-16 + * IEEE -10,+10 30000 2.2e-14 2.7e-15 + * Larger relative error can be observed for z near zero. + * Also tested by csin(casin(z)) = z. + */ + +#include <complex.h> +#include <math.h> + +double complex +casin(double complex z) +{ + double complex w; + static double complex ca, ct, zz, z2; + double x, y; + + x = creal (z); + y = cimag (z); + + if (y == 0.0) { + if (fabs(x) > 1.0) { + w = M_PI_2 + 0.0 * I; + /*mtherr ("casin", DOMAIN);*/ + } + else { + w = asin (x) + 0.0 * I; + } + return (w); + } + + /* Power series expansion */ + /* + b = cabs(z); + if( b < 0.125 ) { + z2.r = (x - y) * (x + y); + z2.i = 2.0 * x * y; + + cn = 1.0; + n = 1.0; + ca.r = x; + ca.i = y; + sum.r = x; + sum.i = y; + do { + ct.r = z2.r * ca.r - z2.i * ca.i; + ct.i = z2.r * ca.i + z2.i * ca.r; + ca.r = ct.r; + ca.i = ct.i; + + cn *= n; + n += 1.0; + cn /= n; + n += 1.0; + b = cn/n; + + ct.r *= b; + ct.i *= b; + sum.r += ct.r; + sum.i += ct.i; + b = fabs(ct.r) + fabs(ct.i); + } + while( b > MACHEP ); + w->r = sum.r; + w->i = sum.i; + return; + } + */ + + ca = x + y * I; + ct = ca * I; + /* sqrt( 1 - z*z) */ + /* cmul( &ca, &ca, &zz ) */ + /*x * x - y * y */ + zz = (x - y) * (x + y) + (2.0 * x * y) * I; + + zz = 1.0 - creal(zz) - cimag(zz) * I; + z2 = csqrt (zz); + + zz = ct + z2; + zz = clog (zz); + /* multiply by 1/i = -i */ + w = zz * (-1.0 * I); + return (w); +} diff --git a/lib/libm/src/s_casinf.c b/lib/libm/src/s_casinf.c new file mode 100644 index 00000000000..84b2096a4fc --- /dev/null +++ b/lib/libm/src/s_casinf.c @@ -0,0 +1,132 @@ +/* $OpenBSD: s_casinf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* casinf() + * + * Complex circular arc sine + * + * + * + * SYNOPSIS: + * + * void casinf(); + * cmplxf z, w; + * + * casinf( &z, &w ); + * + * + * + * DESCRIPTION: + * + * Inverse complex sine: + * + * 2 + * w = -i clog( iz + csqrt( 1 - z ) ). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.1e-5 1.5e-6 + * Larger relative error can be observed for z near zero. + * + */ + +#include <complex.h> +#include <math.h> + +float complex +casinf(float complex z) +{ + float complex w; + float x, y; + static float complex ca, ct, zz, z2; + /* + float cn, n; + static float a, b, s, t, u, v, y2; + static cmplxf sum; + */ + + x = creal(z); + y = cimag(z); + + if(y == 0.0f) { + if(fabsf(x) > 1.0f) { + w = M_PI_2 + 0.0f * I; + /*mtherr( "casinf", DOMAIN );*/ + } + else { + w = asinf (x) + 0.0f * I; + } + return (w); + } + + /* Power series expansion */ + /* + b = cabsf(z); + if(b < 0.125) { + z2.r = (x - y) * (x + y); + z2.i = 2.0 * x * y; + + cn = 1.0; + n = 1.0; + ca.r = x; + ca.i = y; + sum.r = x; + sum.i = y; + do { + ct.r = z2.r * ca.r - z2.i * ca.i; + ct.i = z2.r * ca.i + z2.i * ca.r; + ca.r = ct.r; + ca.i = ct.i; + + cn *= n; + n += 1.0; + cn /= n; + n += 1.0; + b = cn/n; + + ct.r *= b; + ct.i *= b; + sum.r += ct.r; + sum.i += ct.i; + b = fabsf(ct.r) + fabsf(ct.i); + } + while(b > MACHEPF); + w->r = sum.r; + w->i = sum.i; + return; + } + */ + + + ca = x + y * I; + ct = ca * I; /* iz */ + /* sqrt( 1 - z*z) */ + /* cmul( &ca, &ca, &zz ) */ + /*x * x - y * y */ + zz = (x - y) * (x + y) + (2.0f * x * y) * I; + zz = 1.0f - creal(zz) - cimag(zz) * I; + z2 = csqrtf (zz); + + zz = ct + z2; + zz = clogf (zz); + /* multiply by 1/i = -i */ + w = zz * (-1.0f * I); + return (w); +} diff --git a/lib/libm/src/s_casinh.c b/lib/libm/src/s_casinh.c new file mode 100644 index 00000000000..a249cf15ede --- /dev/null +++ b/lib/libm/src/s_casinh.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_casinh.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* casinh + * + * Complex inverse hyperbolic sine + * + * + * + * SYNOPSIS: + * + * double complex casinh(); + * double complex z, w; + * + * w = casinh (z); + * + * + * + * DESCRIPTION: + * + * casinh z = -i casin iz . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.8e-14 2.6e-15 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +casinh(double complex z) +{ + double complex w; + + w = -1.0 * I * casin (z * I); + return (w); +} diff --git a/lib/libm/src/s_casinhf.c b/lib/libm/src/s_casinhf.c new file mode 100644 index 00000000000..334c29bcab2 --- /dev/null +++ b/lib/libm/src/s_casinhf.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_casinhf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* casinhf + * + * Complex inverse hyperbolic sine + * + * + * + * SYNOPSIS: + * + * float complex casinhf(); + * float complex z, w; + * + * w = casinhf (z); + * + * + * + * DESCRIPTION: + * + * casinh z = -i casin iz . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.8e-14 2.6e-15 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +casinhf(float complex z) +{ + float complex w; + + w = -1.0f * I * casinf (z * I); + return (w); +} diff --git a/lib/libm/src/s_catan.c b/lib/libm/src/s_catan.c new file mode 100644 index 00000000000..ddd7940456e --- /dev/null +++ b/lib/libm/src/s_catan.c @@ -0,0 +1,126 @@ +/* $OpenBSD: s_catan.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* catan() + * + * Complex circular arc tangent + * + * + * + * SYNOPSIS: + * + * double complex catan(); + * double complex z, w; + * + * w = catan (z); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * 1 ( 2x ) + * Re w = - arctan(-----------) + k PI + * 2 ( 2 2) + * (1 - x - y ) + * + * ( 2 2) + * 1 (x + (y+1) ) + * Im w = - log(------------) + * 4 ( 2 2) + * (x + (y-1) ) + * + * Where k is an arbitrary integer. + * + * catan(z) = -i catanh(iz). + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5900 1.3e-16 7.8e-18 + * IEEE -10,+10 30000 2.3e-15 8.5e-17 + * The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, + * had peak relative error 1.5e-16, rms relative error + * 2.9e-17. See also clog(). + */ + +#include <complex.h> +#include <math.h> + +#define MAXNUM 1.0e308 + +static const double DP1 = 3.14159265160560607910E0; +static const double DP2 = 1.98418714791870343106E-9; +static const double DP3 = 1.14423774522196636802E-17; + +static double +_redupi(double x) +{ + double t; + long i; + + t = x/M_PI; + if(t >= 0.0) + t += 0.5; + else + t -= 0.5; + + i = t; /* the multiple */ + t = i; + t = ((x - t * DP1) - t * DP2) - t * DP3; + return (t); +} + +double complex +catan(double complex z) +{ + double complex w; + double a, t, x, x2, y; + + x = creal (z); + y = cimag (z); + + if ((x == 0.0) && (y > 1.0)) + goto ovrf; + + x2 = x * x; + a = 1.0 - x2 - (y * y); + if (a == 0.0) + goto ovrf; + + t = 0.5 * atan2 (2.0 * x, a); + w = _redupi (t); + + t = y - 1.0; + a = x2 + (t * t); + if (a == 0.0) + goto ovrf; + + t = y + 1.0; + a = (x2 + (t * t))/a; + w = w + (0.25 * log (a)) * I; + return (w); + +ovrf: + /*mtherr ("catan", OVERFLOW);*/ + w = MAXNUM + MAXNUM * I; + return (w); +} diff --git a/lib/libm/src/s_catanf.c b/lib/libm/src/s_catanf.c new file mode 100644 index 00000000000..e112e7f090a --- /dev/null +++ b/lib/libm/src/s_catanf.c @@ -0,0 +1,124 @@ +/* $OpenBSD: s_catanf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* catanf() + * + * Complex circular arc tangent + * + * + * + * SYNOPSIS: + * + * float complex catanf(); + * float complex z, w; + * + * w = catanf( z ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * 1 ( 2x ) + * Re w = - arctan(-----------) + k PI + * 2 ( 2 2) + * (1 - x - y ) + * + * ( 2 2) + * 1 (x + (y+1) ) + * Im w = - log(------------) + * 4 ( 2 2) + * (x + (y-1) ) + * + * Where k is an arbitrary integer. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.3e-6 5.2e-8 + * + */ + +#include <complex.h> +#include <math.h> + +#define MAXNUMF 1.0e38F + +static const double DP1 = 3.140625; +static const double DP2 = 9.67502593994140625E-4; +static const double DP3 = 1.509957990978376432E-7; + +static float +_redupif(float xx) +{ + float x, t; + long i; + + x = xx; + t = x/(float)M_PI; + if(t >= 0.0) + t += 0.5; + else + t -= 0.5; + + i = t; /* the multiple */ + t = i; + t = ((x - t * DP1) - t * DP2) - t * DP3; + return(t); +} + +float complex +catanf(float complex z) +{ + float complex w; + float a, t, x, x2, y; + + x = creal (z); + y = cimag (z); + + if((x == 0.0f) && (y > 1.0f)) + goto ovrf; + + x2 = x * x; + a = 1.0f - x2 - (y * y); + if (a == 0.0f) + goto ovrf; + + t = 0.5f * atan2f(2.0f * x, a); + w = _redupif(t); + + t = y - 1.0f; + a = x2 + (t * t); + if(a == 0.0f) + goto ovrf; + + t = y + 1.0f; + a = (x2 + (t * t))/a; + w = w + (0.25f * logf (a)) * I; + return (w); + +ovrf: + /*mtherr( "catanf", OVERFLOW );*/ + w = MAXNUMF + MAXNUMF * I; + return (w); +} diff --git a/lib/libm/src/s_catanh.c b/lib/libm/src/s_catanh.c new file mode 100644 index 00000000000..7a6d6f99646 --- /dev/null +++ b/lib/libm/src/s_catanh.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_catanh.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* catanh + * + * Complex inverse hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * double complex catanh(); + * double complex z, w; + * + * w = catanh (z); + * + * + * + * DESCRIPTION: + * + * Inverse tanh, equal to -i catan (iz); + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.3e-16 6.2e-17 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +catanh(double complex z) +{ + double complex w; + + w = -1.0 * I * catan (z * I); + return (w); +} diff --git a/lib/libm/src/s_catanhf.c b/lib/libm/src/s_catanhf.c new file mode 100644 index 00000000000..35fc8f7541d --- /dev/null +++ b/lib/libm/src/s_catanhf.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_catanhf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* catanhf + * + * Complex inverse hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * float complex catanhf(); + * float complex z, w; + * + * w = catanhf (z); + * + * + * + * DESCRIPTION: + * + * Inverse tanh, equal to -i catan (iz); + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.3e-16 6.2e-17 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +catanhf(float complex z) +{ + float complex w; + + w = -1.0f * I * catanf (z * I); + return (w); +} diff --git a/lib/libm/src/s_ccos.c b/lib/libm/src/s_ccos.c new file mode 100644 index 00000000000..264143ab24c --- /dev/null +++ b/lib/libm/src/s_ccos.c @@ -0,0 +1,84 @@ +/* $OpenBSD: s_ccos.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ccos() + * + * Complex circular cosine + * + * + * + * SYNOPSIS: + * + * double complex ccos(); + * double complex z, w; + * + * w = ccos (z); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = cos x cosh y - i sin x sinh y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8400 4.5e-17 1.3e-17 + * IEEE -10,+10 30000 3.8e-16 1.0e-16 + */ + +#include <complex.h> +#include <math.h> + +/* calculate cosh and sinh */ + +static void +_cchsh(double x, double *c, double *s) +{ + double e, ei; + + if (fabs(x) <= 0.5) { + *c = cosh(x); + *s = sinh(x); + } + else { + e = exp(x); + ei = 0.5/e; + e = 0.5 * e; + *s = e - ei; + *c = e + ei; + } +} + +double complex +ccos(double complex z) +{ + double complex w; + double ch, sh; + + _cchsh( cimag(z), &ch, &sh ); + w = cos(creal (z)) * ch - (sin (creal (z)) * sh) * I; + return (w); +} diff --git a/lib/libm/src/s_ccosf.c b/lib/libm/src/s_ccosf.c new file mode 100644 index 00000000000..299e06b1276 --- /dev/null +++ b/lib/libm/src/s_ccosf.c @@ -0,0 +1,84 @@ +/* $OpenBSD: s_ccosf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ccosf() + * + * Complex circular cosine + * + * + * + * SYNOPSIS: + * + * void ccosf(); + * cmplxf z, w; + * + * ccosf( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = cos x cosh y - i sin x sinh y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.8e-7 5.5e-8 + */ + +#include <complex.h> +#include <math.h> + +/* calculate cosh and sinh */ + +static void +_cchshf(float xx, float *c, float *s) +{ + float x, e, ei; + + x = xx; + if(fabsf(x) <= 0.5f) { + *c = coshf(x); + *s = sinhf(x); + } + else { + e = expf(x); + ei = 0.5f/e; + e = 0.5f * e; + *s = e - ei; + *c = e + ei; + } +} + +float complex +ccosf(float complex z) +{ + float complex w; + float ch, sh; + + _cchshf( cimag(z), &ch, &sh ); + w = cosf( creal(z) ) * ch + ( -sinf( creal(z) ) * sh) * I; + return (w); +} diff --git a/lib/libm/src/s_ccosh.c b/lib/libm/src/s_ccosh.c new file mode 100644 index 00000000000..0267b6f2a15 --- /dev/null +++ b/lib/libm/src/s_ccosh.c @@ -0,0 +1,58 @@ +/* $OpenBSD: s_ccosh.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ccosh + * + * Complex hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * double complex ccosh(); + * double complex z, w; + * + * w = ccosh (z); + * + * + * + * DESCRIPTION: + * + * ccosh(z) = cosh x cos y + i sinh x sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.9e-16 8.1e-17 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +ccosh(double complex z) +{ + double complex w; + double x, y; + + x = creal(z); + y = cimag(z); + w = cosh (x) * cos (y) + (sinh (x) * sin (y)) * I; + return (w); +} diff --git a/lib/libm/src/s_ccoshf.c b/lib/libm/src/s_ccoshf.c new file mode 100644 index 00000000000..527a47afb98 --- /dev/null +++ b/lib/libm/src/s_ccoshf.c @@ -0,0 +1,58 @@ +/* $OpenBSD: s_ccoshf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ccoshf + * + * Complex hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * float complex ccoshf(); + * float complex z, w; + * + * w = ccoshf (z); + * + * + * + * DESCRIPTION: + * + * ccosh(z) = cosh x cos y + i sinh x sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.9e-16 8.1e-17 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +ccoshf(float complex z) +{ + float complex w; + float x, y; + + x = creal(z); + y = cimag(z); + w = coshf (x) * cosf (y) + (sinhf (x) * sinf (y)) * I; + return (w); +} diff --git a/lib/libm/src/s_cexp.c b/lib/libm/src/s_cexp.c new file mode 100644 index 00000000000..52241a8eb5e --- /dev/null +++ b/lib/libm/src/s_cexp.c @@ -0,0 +1,70 @@ +/* $OpenBSD: s_cexp.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cexp() + * + * Complex exponential function + * + * + * + * SYNOPSIS: + * + * double complex cexp (); + * double complex z, w; + * + * w = cexp (z); + * + * + * + * DESCRIPTION: + * + * Returns the exponential of the complex argument z + * into the complex result w. + * + * If + * z = x + iy, + * r = exp(x), + * + * then + * + * w = r cos y + i r sin y. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8700 3.7e-17 1.1e-17 + * IEEE -10,+10 30000 3.0e-16 8.7e-17 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +cexp(double complex z) +{ + double complex w; + double r, x, y; + + x = creal (z); + y = cimag (z); + r = exp (x); + w = r * cos (y) + r * sin (y) * I; + return (w); +} diff --git a/lib/libm/src/s_cexpf.c b/lib/libm/src/s_cexpf.c new file mode 100644 index 00000000000..041bf84318d --- /dev/null +++ b/lib/libm/src/s_cexpf.c @@ -0,0 +1,67 @@ +/* $OpenBSD: s_cexpf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cexpf() + * + * Complex exponential function + * + * + * + * SYNOPSIS: + * + * void cexpf(); + * cmplxf z, w; + * + * cexpf( &z, &w ); + * + * + * + * DESCRIPTION: + * + * Returns the exponential of the complex argument z + * into the complex result w. + * + * If + * z = x + iy, + * r = exp(x), + * + * then + * + * w = r cos y + i r sin y. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.4e-7 4.5e-8 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +cexpf(float complex z) +{ + float complex w; + float r; + + r = expf( creal(z) ); + w = r * cosf( cimag(z) ) + r * sinf( cimag(z) ) * I; + return (w); +} diff --git a/lib/libm/src/s_cimag.c b/lib/libm/src/s_cimag.c new file mode 100644 index 00000000000..85aaa2975d1 --- /dev/null +++ b/lib/libm/src/s_cimag.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_cimag.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +double +cimag(double complex z) +{ + return __imag__ z; +} diff --git a/lib/libm/src/s_cimagf.c b/lib/libm/src/s_cimagf.c new file mode 100644 index 00000000000..609851fbbda --- /dev/null +++ b/lib/libm/src/s_cimagf.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_cimagf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +float +cimagf(float complex z) +{ + return __imag__ z; +} diff --git a/lib/libm/src/s_clog.c b/lib/libm/src/s_clog.c new file mode 100644 index 00000000000..31163700607 --- /dev/null +++ b/lib/libm/src/s_clog.c @@ -0,0 +1,72 @@ +/* $OpenBSD: s_clog.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* clog.c + * + * Complex natural logarithm + * + * + * + * SYNOPSIS: + * + * double complex clog(); + * double complex z, w; + * + * w = clog (z); + * + * + * + * DESCRIPTION: + * + * Returns complex logarithm to the base e (2.718...) of + * the complex argument x. + * + * If z = x + iy, r = sqrt( x**2 + y**2 ), + * then + * w = log(r) + i arctan(y/x). + * + * The arctangent ranges from -PI to +PI. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 7000 8.5e-17 1.9e-17 + * IEEE -10,+10 30000 5.0e-15 1.1e-16 + * + * Larger relative error can be observed for z near 1 +i0. + * In IEEE arithmetic the peak absolute error is 5.2e-16, rms + * absolute error 1.0e-16. + */ + +#include <complex.h> +#include <math.h> + +double complex +clog(double complex z) +{ + double complex w; + double p, rr; + + /*rr = sqrt( z->r * z->r + z->i * z->i );*/ + rr = cabs(z); + p = log(rr); + rr = atan2 (cimag (z), creal (z)); + w = p + rr * I; + return (w); +} diff --git a/lib/libm/src/s_clogf.c b/lib/libm/src/s_clogf.c new file mode 100644 index 00000000000..2cea0f218be --- /dev/null +++ b/lib/libm/src/s_clogf.c @@ -0,0 +1,72 @@ +/* $OpenBSD: s_clogf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* clogf.c + * + * Complex natural logarithm + * + * + * + * SYNOPSIS: + * + * void clogf(); + * cmplxf z, w; + * + * clogf( &z, &w ); + * + * + * + * DESCRIPTION: + * + * Returns complex logarithm to the base e (2.718...) of + * the complex argument x. + * + * If z = x + iy, r = sqrt( x**2 + y**2 ), + * then + * w = log(r) + i arctan(y/x). + * + * The arctangent ranges from -PI to +PI. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.9e-6 6.2e-8 + * + * Larger relative error can be observed for z near 1 +i0. + * In IEEE arithmetic the peak absolute error is 3.1e-7. + * + */ + +#include <complex.h> +#include <math.h> + +float complex +clogf(float complex z) +{ + float complex w; + float p, rr, x, y; + + x = creal(z); + y = cimag(z); + rr = atan2f(y, x); + p = cabsf(z); + p = logf(p); + w = p + rr * I; + return (w); +} diff --git a/lib/libm/src/s_conj.c b/lib/libm/src/s_conj.c new file mode 100644 index 00000000000..fd38a9d034c --- /dev/null +++ b/lib/libm/src/s_conj.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_conj.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +double complex +conj(double complex z) +{ + return ~z; +} diff --git a/lib/libm/src/s_conjf.c b/lib/libm/src/s_conjf.c new file mode 100644 index 00000000000..d5b765790c6 --- /dev/null +++ b/lib/libm/src/s_conjf.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_conjf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +float complex +conjf(float complex z) +{ + return ~z; +} diff --git a/lib/libm/src/s_cpow.c b/lib/libm/src/s_cpow.c new file mode 100644 index 00000000000..6b0537393ce --- /dev/null +++ b/lib/libm/src/s_cpow.c @@ -0,0 +1,71 @@ +/* $OpenBSD: s_cpow.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cpow + * + * Complex power function + * + * + * + * SYNOPSIS: + * + * double complex cpow(); + * double complex a, z, w; + * + * w = cpow (a, z); + * + * + * + * DESCRIPTION: + * + * Raises complex A to the complex Zth power. + * Definition is per AMS55 # 4.2.8, + * analytically equivalent to cpow(a,z) = cexp(z clog(a)). + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 9.4e-15 1.5e-15 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +cpow(double complex a, double complex z) +{ + double complex w; + double x, y, r, theta, absa, arga; + + x = creal (z); + y = cimag (z); + absa = cabs (a); + if (absa == 0.0) { + return (0.0 + 0.0 * I); + } + arga = carg (a); + r = pow (absa, x); + theta = x * arga; + if (y != 0.0) { + r = r * exp (-y * arga); + theta = theta + y * log (absa); + } + w = r * cos (theta) + (r * sin (theta)) * I; + return (w); +} diff --git a/lib/libm/src/s_cpowf.c b/lib/libm/src/s_cpowf.c new file mode 100644 index 00000000000..8a24596d338 --- /dev/null +++ b/lib/libm/src/s_cpowf.c @@ -0,0 +1,71 @@ +/* $OpenBSD: s_cpowf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* cpowf + * + * Complex power function + * + * + * + * SYNOPSIS: + * + * float complex cpowf(); + * float complex a, z, w; + * + * w = cpowf (a, z); + * + * + * + * DESCRIPTION: + * + * Raises complex A to the complex Zth power. + * Definition is per AMS55 # 4.2.8, + * analytically equivalent to cpow(a,z) = cexp(z clog(a)). + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 9.4e-15 1.5e-15 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +cpowf(float complex a, float complex z) +{ + float complex w; + float x, y, r, theta, absa, arga; + + x = creal (z); + y = cimag (z); + absa = cabsf (a); + if (absa == 0.0f) { + return (0.0f + 0.0f * I); + } + arga = cargf (a); + r = powf (absa, x); + theta = x * arga; + if (y != 0.0f) { + r = r * expf (-y * arga); + theta = theta + y * logf (absa); + } + w = r * cosf (theta) + (r * sinf (theta)) * I; + return (w); +} diff --git a/lib/libm/src/s_cproj.c b/lib/libm/src/s_cproj.c new file mode 100644 index 00000000000..7064e9285c3 --- /dev/null +++ b/lib/libm/src/s_cproj.c @@ -0,0 +1,32 @@ +/* $OpenBSD: s_cproj.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +double complex +cproj(double complex z) +{ + double complex res; + + if (isinf(__real__ z) || isinf(__imag__ z)) { + __real__ res = INFINITY; + __imag__ res = copysign(0.0, __imag__ z); + } + + return res; +} diff --git a/lib/libm/src/s_cprojf.c b/lib/libm/src/s_cprojf.c new file mode 100644 index 00000000000..af754ed3d70 --- /dev/null +++ b/lib/libm/src/s_cprojf.c @@ -0,0 +1,32 @@ +/* $OpenBSD: s_cprojf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +float complex +cprojf(float complex z) +{ + float complex res; + + if (isinf(__real__ z) || isinf(__imag__ z)) { + __real__ res = INFINITY; + __imag__ res = copysign(0.0, __imag__ z); + } + + return res; +} diff --git a/lib/libm/src/s_creal.c b/lib/libm/src/s_creal.c new file mode 100644 index 00000000000..55c36e8d839 --- /dev/null +++ b/lib/libm/src/s_creal.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_creal.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +double +creal(double complex z) +{ + return __real__ z; +} diff --git a/lib/libm/src/s_crealf.c b/lib/libm/src/s_crealf.c new file mode 100644 index 00000000000..12c54ea615b --- /dev/null +++ b/lib/libm/src/s_crealf.c @@ -0,0 +1,25 @@ +/* $OpenBSD: s_crealf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* + * Copyright (c) 2008 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 <complex.h> +#include <math.h> + +float +crealf(float complex z) +{ + return __real__ z; +} diff --git a/lib/libm/src/s_csin.c b/lib/libm/src/s_csin.c new file mode 100644 index 00000000000..ad249de7ac0 --- /dev/null +++ b/lib/libm/src/s_csin.c @@ -0,0 +1,86 @@ +/* $OpenBSD: s_csin.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* csin() + * + * Complex circular sine + * + * + * + * SYNOPSIS: + * + * double complex csin(); + * double complex z, w; + * + * w = csin (z); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = sin x cosh y + i cos x sinh y. + * + * csin(z) = -i csinh(iz). + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8400 5.3e-17 1.3e-17 + * IEEE -10,+10 30000 3.8e-16 1.0e-16 + * Also tested by csin(casin(z)) = z. + * + */ + +#include <complex.h> +#include <math.h> + +/* calculate cosh and sinh */ + +static void +cchsh(double x, double *c, double *s) +{ + double e, ei; + + if (fabs(x) <= 0.5) { + *c = cosh(x); + *s = sinh(x); + } + else { + e = exp(x); + ei = 0.5/e; + e = 0.5 * e; + *s = e - ei; + *c = e + ei; + } +} + +double complex +csin(double complex z) +{ + double complex w; + double ch, sh; + + cchsh( cimag (z), &ch, &sh ); + w = sin (creal(z)) * ch + (cos (creal(z)) * sh) * I; + return (w); +} diff --git a/lib/libm/src/s_csinf.c b/lib/libm/src/s_csinf.c new file mode 100644 index 00000000000..aacc59f9dc6 --- /dev/null +++ b/lib/libm/src/s_csinf.c @@ -0,0 +1,85 @@ +/* $OpenBSD: s_csinf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* csinf() + * + * Complex circular sine + * + * + * + * SYNOPSIS: + * + * void csinf(); + * cmplxf z, w; + * + * csinf( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = sin x cosh y + i cos x sinh y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.9e-7 5.5e-8 + * + */ + +#include <complex.h> +#include <math.h> + +/* calculate cosh and sinh */ + +static void +cchshf(float xx, float *c, float *s) +{ + float x, e, ei; + + x = xx; + if(fabsf(x) <= 0.5f) { + *c = coshf(x); + *s = sinhf(x); + } + else { + e = expf(x); + ei = 0.5f/e; + e = 0.5f * e; + *s = e - ei; + *c = e + ei; + } +} + +float complex +csinf(float complex z) +{ + float complex w; + float ch, sh; + + cchshf((float) cimag(z), &ch, &sh); + w = sinf(creal(z)) * ch + (cosf(creal(z)) * sh) * I; + return (w); +} diff --git a/lib/libm/src/s_csinh.c b/lib/libm/src/s_csinh.c new file mode 100644 index 00000000000..ece633cda8f --- /dev/null +++ b/lib/libm/src/s_csinh.c @@ -0,0 +1,57 @@ +/* $OpenBSD: s_csinh.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* csinh + * + * Complex hyperbolic sine + * + * + * + * SYNOPSIS: + * + * double complex csinh(); + * double complex z, w; + * + * w = csinh (z); + * + * DESCRIPTION: + * + * csinh z = (cexp(z) - cexp(-z))/2 + * = sinh x * cos y + i cosh x * sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 3.1e-16 8.2e-17 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +csinh(double complex z) +{ + double complex w; + double x, y; + + x = creal(z); + y = cimag(z); + w = sinh (x) * cos (y) + (cosh (x) * sin (y)) * I; + return (w); +} diff --git a/lib/libm/src/s_csinhf.c b/lib/libm/src/s_csinhf.c new file mode 100644 index 00000000000..170b4405640 --- /dev/null +++ b/lib/libm/src/s_csinhf.c @@ -0,0 +1,57 @@ +/* $OpenBSD: s_csinhf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* csinhf + * + * Complex hyperbolic sine + * + * + * + * SYNOPSIS: + * + * float complex csinhf(); + * float complex z, w; + * + * w = csinhf (z); + * + * DESCRIPTION: + * + * csinh z = (cexp(z) - cexp(-z))/2 + * = sinh x * cos y + i cosh x * sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 3.1e-16 8.2e-17 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +csinhf(float complex z) +{ + float complex w; + float x, y; + + x = creal(z); + y = cimag(z); + w = sinhf (x) * cosf (y) + (coshf (x) * sinf (y)) * I; + return (w); +} diff --git a/lib/libm/src/s_csqrt.c b/lib/libm/src/s_csqrt.c new file mode 100644 index 00000000000..3ed4aa56076 --- /dev/null +++ b/lib/libm/src/s_csqrt.c @@ -0,0 +1,131 @@ +/* $OpenBSD: s_csqrt.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* csqrt() + * + * Complex square root + * + * + * + * SYNOPSIS: + * + * double complex csqrt(); + * double complex z, w; + * + * w = csqrt (z); + * + * + * + * DESCRIPTION: + * + * + * If z = x + iy, r = |z|, then + * + * 1/2 + * Re w = [ (r + x)/2 ] , + * + * 1/2 + * Im w = [ (r - x)/2 ] . + * + * Cancellation error in r-x or r+x is avoided by using the + * identity 2 Re w Im w = y. + * + * Note that -w is also a square root of z. The root chosen + * is always in the right half plane and Im w has the same sign as y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 25000 3.2e-17 9.6e-18 + * IEEE -10,+10 1,000,000 2.9e-16 6.1e-17 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +csqrt(double complex z) +{ + double complex w; + double x, y, r, t, scale; + + x = creal (z); + y = cimag (z); + + if (y == 0.0) { + if (x == 0.0) { + w = 0.0 + y * I; + } + else { + r = fabs (x); + r = sqrt (r); + if (x < 0.0) { + w = 0.0 + r * I; + } + else { + w = r + y * I; + } + } + return (w); + } + if (x == 0.0) { + r = fabs (y); + r = sqrt (0.5*r); + if (y > 0) + w = r + r * I; + else + w = r - r * I; + return (w); + } + /* Rescale to avoid internal overflow or underflow. */ + if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) { + x *= 0.25; + y *= 0.25; + scale = 2.0; + } + else { + x *= 1.8014398509481984e16; /* 2^54 */ + y *= 1.8014398509481984e16; + scale = 7.450580596923828125e-9; /* 2^-27 */ +#if 0 + x *= 4.0; + y *= 4.0; + scale = 0.5; +#endif + } + w = x + y * I; + r = cabs(w); + if (x > 0) { + t = sqrt(0.5 * r + 0.5 * x); + r = scale * fabs((0.5 * y) / t); + t *= scale; + } + else { + r = sqrt( 0.5 * r - 0.5 * x ); + t = scale * fabs( (0.5 * y) / r ); + r *= scale; + } + if (y < 0) + w = t - r * I; + else + w = t + r * I; + return (w); +} diff --git a/lib/libm/src/s_csqrtf.c b/lib/libm/src/s_csqrtf.c new file mode 100644 index 00000000000..9f9e8982b0e --- /dev/null +++ b/lib/libm/src/s_csqrtf.c @@ -0,0 +1,131 @@ +/* $OpenBSD: s_csqrtf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* csqrtf() + * + * Complex square root + * + * + * + * SYNOPSIS: + * + * float complex csqrtf(); + * float complex z, w; + * + * w = csqrtf( z ); + * + * + * + * DESCRIPTION: + * + * + * If z = x + iy, r = |z|, then + * + * 1/2 + * Re w = [ (r + x)/2 ] , + * + * 1/2 + * Im w = [ (r - x)/2 ] . + * + * Cancellation error in r-x or r+x is avoided by using the + * identity 2 Re w Im w = y. + * + * Note that -w is also a square root of z. The root chosen + * is always in the right half plane and Im w has the same sign as y. + * + * + * + * ACCURACY: + * + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 1,000,000 1.8e-7 3.5e-8 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +csqrtf(float complex z) +{ + float complex w; + float x, y, r, t, scale; + + x = creal(z); + y = cimag(z); + + if(y == 0.0f) { + if (x < 0.0f) { + w = 0.0f + sqrtf(-x) * I; + return (w); + } + else if (x == 0.0f) { + return (0.0f + y * I); + } + else { + w = sqrtf(x) + y * I; + return (w); + } + } + + if (x == 0.0f) { + r = fabsf(y); + r = sqrtf(0.5f*r); + if(y > 0) + w = r + r * I; + else + w = r - r * I; + return (w); + } + + /* Rescale to avoid internal overflow or underflow. */ + if ((fabsf(x) > 4.0f) || (fabsf(y) > 4.0f)) { + x *= 0.25f; + y *= 0.25f; + scale = 2.0f; + } + else { + x *= 6.7108864e7f; /* 2^26 */ + y *= 6.7108864e7f; + scale = 1.220703125e-4f; /* 2^-13 */ +#if 0 + x *= 4.0f; + y *= 4.0f; + scale = 0.5f; +#endif + } + w = x + y * I; + r = cabsf(w); + if (x > 0) { + t = sqrtf( 0.5f * r + 0.5f * x ); + r = scale * fabsf((0.5f * y) / t); + t *= scale; + } + else { + r = sqrtf(0.5f * r - 0.5f * x); + t = scale * fabsf((0.5f * y) / r); + r *= scale; + } + + if (y < 0) + w = t - r * I; + else + w = t + r * I; + return (w); +} diff --git a/lib/libm/src/s_ctan.c b/lib/libm/src/s_ctan.c new file mode 100644 index 00000000000..4db6fb18e83 --- /dev/null +++ b/lib/libm/src/s_ctan.c @@ -0,0 +1,152 @@ +/* $OpenBSD: s_ctan.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ctan() + * + * Complex circular tangent + * + * + * + * SYNOPSIS: + * + * double complex ctan(); + * double complex z, w; + * + * w = ctan (z); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * sin 2x + i sinh 2y + * w = --------------------. + * cos 2x + cosh 2y + * + * On the real axis the denominator is zero at odd multiples + * of PI/2. The denominator is evaluated by its Taylor + * series near these points. + * + * ctan(z) = -i ctanh(iz). + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5200 7.1e-17 1.6e-17 + * IEEE -10,+10 30000 7.2e-16 1.2e-16 + * Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. + */ + +#include <complex.h> +#include <math.h> + +#define MACHEP 1.1e-16 +#define MAXNUM 1.0e308 + +static const double DP1 = 3.14159265160560607910E0; +static const double DP2 = 1.98418714791870343106E-9; +static const double DP3 = 1.14423774522196636802E-17; + +static double +_redupi(double x) +{ + double t; + long i; + + t = x/M_PI; + if (t >= 0.0) + t += 0.5; + else + t -= 0.5; + + i = t; /* the multiple */ + t = i; + t = ((x - t * DP1) - t * DP2) - t * DP3; + return (t); +} + +/* Taylor series expansion for cosh(2y) - cos(2x) */ + +static double +_ctans(double complex z) +{ + double f, x, x2, y, y2, rn, t; + double d; + + x = fabs (2.0 * creal (z)); + y = fabs (2.0 * cimag(z)); + + x = _redupi(x); + + x = x * x; + y = y * y; + x2 = 1.0; + y2 = 1.0; + f = 1.0; + rn = 0.0; + d = 0.0; + do { + rn += 1.0; + f *= rn; + rn += 1.0; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 + x2; + t /= f; + d += t; + + rn += 1.0; + f *= rn; + rn += 1.0; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 - x2; + t /= f; + d += t; + } + while (fabs(t/d) > MACHEP) + ; + return (d); +} + +double complex +ctan(double complex z) +{ + double complex w; + double d; + + d = cos (2.0 * creal (z)) + cosh (2.0 * cimag (z)); + + if (fabs(d) < 0.25) + d = _ctans (z); + + if (d == 0.0) { + /*mtherr ("ctan", OVERFLOW);*/ + w = MAXNUM + MAXNUM * I; + return (w); + } + + w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I; + return (w); +} diff --git a/lib/libm/src/s_ctanf.c b/lib/libm/src/s_ctanf.c new file mode 100644 index 00000000000..ca67d602e55 --- /dev/null +++ b/lib/libm/src/s_ctanf.c @@ -0,0 +1,148 @@ +/* $OpenBSD: s_ctanf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ctanf() + * + * Complex circular tangent + * + * + * + * SYNOPSIS: + * + * void ctanf(); + * cmplxf z, w; + * + * ctanf( &z, &w ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * sin 2x + i sinh 2y + * w = --------------------. + * cos 2x + cosh 2y + * + * On the real axis the denominator is zero at odd multiples + * of PI/2. The denominator is evaluated by its Taylor + * series near these points. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 3.3e-7 5.1e-8 + */ + +#include <complex.h> +#include <math.h> + +#define MACHEPF 3.0e-8 +#define MAXNUMF 1.0e38f + +static const double DP1 = 3.140625; +static const double DP2 = 9.67502593994140625E-4; +static const double DP3 = 1.509957990978376432E-7; + +static float +_redupif(float xx) +{ + float x, t; + long i; + + x = xx; + t = x/(float)M_PI; + if(t >= 0.0) + t += 0.5; + else + t -= 0.5; + + i = t; /* the multiple */ + t = i; + t = ((x - t * DP1) - t * DP2) - t * DP3; + return(t); +} + +/* Taylor series expansion for cosh(2y) - cos(2x) */ + +static float +_ctansf(float complex z) +{ + float f, x, x2, y, y2, rn, t, d; + + x = fabsf(2.0f * creal(z)); + y = fabsf(2.0f * cimag(z)); + + x = _redupif(x); + + x = x * x; + y = y * y; + x2 = 1.0f; + y2 = 1.0f; + f = 1.0f; + rn = 0.0f; + d = 0.0f; + do { + rn += 1.0f; + f *= rn; + rn += 1.0f; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 + x2; + t /= f; + d += t; + + rn += 1.0f; + f *= rn; + rn += 1.0f; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 - x2; + t /= f; + d += t; + } + while (fabsf(t/d) > MACHEPF) + ; + return(d); +} + +float complex +ctanf(float complex z) +{ + float complex w; + float d; + + d = cosf( 2.0f * creal(z) ) + coshf( 2.0f * cimag(z) ); + + if(fabsf(d) < 0.25f) + d = _ctansf(z); + + if (d == 0.0f) { + /*mtherr( "ctanf", OVERFLOW );*/ + w = MAXNUMF + MAXNUMF * I; + return (w); + } + w = sinf (2.0f * creal(z)) / d + (sinhf (2.0f * cimag(z)) / d) * I; + return (w); +} diff --git a/lib/libm/src/s_ctanh.c b/lib/libm/src/s_ctanh.c new file mode 100644 index 00000000000..30c20155672 --- /dev/null +++ b/lib/libm/src/s_ctanh.c @@ -0,0 +1,59 @@ +/* $OpenBSD: s_ctanh.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ctanh + * + * Complex hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * double complex ctanh(); + * double complex z, w; + * + * w = ctanh (z); + * + * + * + * DESCRIPTION: + * + * tanh z = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y) . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.7e-14 2.4e-16 + * + */ + +#include <complex.h> +#include <math.h> + +double complex +ctanh(double complex z) +{ + double complex w; + double x, y, d; + + x = creal(z); + y = cimag(z); + d = cosh (2.0 * x) + cos (2.0 * y); + w = sinh (2.0 * x) / d + (sin (2.0 * y) / d) * I; + return (w); +} diff --git a/lib/libm/src/s_ctanhf.c b/lib/libm/src/s_ctanhf.c new file mode 100644 index 00000000000..a0f61614452 --- /dev/null +++ b/lib/libm/src/s_ctanhf.c @@ -0,0 +1,59 @@ +/* $OpenBSD: s_ctanhf.c,v 1.1 2008/09/07 20:36:09 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. + */ + +/* ctanhf + * + * Complex hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * float complex ctanhf(); + * float complex z, w; + * + * w = ctanhf (z); + * + * + * + * DESCRIPTION: + * + * tanh z = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y) . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.7e-14 2.4e-16 + * + */ + +#include <complex.h> +#include <math.h> + +float complex +ctanhf(float complex z) +{ + float complex w; + float x, y, d; + + x = creal(z); + y = cimag(z); + d = coshf (2.0f * x) + cosf (2.0f * y); + w = sinhf (2.0f * x) / d + (sinf (2.0f * y) / d) * I; + return (w); +} diff --git a/lib/libm/src/s_erf.c b/lib/libm/src/s_erf.c index 844e2f5565a..22944c1bc83 100644 --- a/lib/libm/src/s_erf.c +++ b/lib/libm/src/s_erf.c @@ -235,7 +235,7 @@ erf(double x) } z = x; SET_LOW_WORD(z,0); - r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S); + r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S); if(hx>=0) return one-r/x; else return r/x-one; } @@ -293,8 +293,7 @@ erfc(double x) } z = x; SET_LOW_WORD(z,0); - r = __ieee754_exp(-z*z-0.5625)* - __ieee754_exp((z-x)*(z+x)+R/S); + r = exp(-z*z-0.5625) * exp((z-x)*(z+x)+R/S); if(hx>0) return r/x; else return two-r/x; } else { if(hx>0) return tiny*tiny; else return two-tiny; diff --git a/lib/libm/src/s_erff.c b/lib/libm/src/s_erff.c index faff192c21e..7eadb30dfa3 100644 --- a/lib/libm/src/s_erff.c +++ b/lib/libm/src/s_erff.c @@ -144,7 +144,7 @@ erff(float x) } GET_FLOAT_WORD(ix,x); SET_FLOAT_WORD(z,ix&0xfffff000); - r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S); + r = expf(-z*z-(float)0.5625)*expf((z-x)*(z+x)+R/S); if(hx>=0) return one-r/x; else return r/x-one; } @@ -202,8 +202,7 @@ erfcf(float x) } GET_FLOAT_WORD(ix,x); SET_FLOAT_WORD(z,ix&0xfffff000); - r = __ieee754_expf(-z*z-(float)0.5625)* - __ieee754_expf((z-x)*(z+x)+R/S); + r = expf(-z*z-(float)0.5625) * expf((z-x)*(z+x)+R/S); if(hx>0) return r/x; else return two-r/x; } else { if(hx>0) return tiny*tiny; else return two-tiny; diff --git a/lib/libm/src/s_fdim.c b/lib/libm/src/s_fdim.c new file mode 100644 index 00000000000..6af4ba038b9 --- /dev/null +++ b/lib/libm/src/s_fdim.c @@ -0,0 +1,44 @@ +/* $OpenBSD: s_fdim.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/*- + * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include <math.h> + +#define DECL(type, fn) \ +type \ +fn(type x, type y) \ +{ \ + \ + if (isnan(x)) \ + return (x); \ + if (isnan(y)) \ + return (y); \ + return (x > y ? x - y : 0.0); \ +} + +DECL(double, fdim) +DECL(float, fdimf) +DECL(long double, fdiml) diff --git a/lib/libm/src/s_fmax.c b/lib/libm/src/s_fmax.c new file mode 100644 index 00000000000..76f91272f2c --- /dev/null +++ b/lib/libm/src/s_fmax.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_fmax.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/*- + * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include <sys/types.h> +#include <machine/ieee.h> +#include <math.h> + +double +fmax(double x, double y) +{ + struct ieee_double *px = (struct ieee_double *)&x; + struct ieee_double *py = (struct ieee_double *)&y; + + /* Check for NaNs to avoid raising spurious exceptions. */ + if (px->dbl_exp == DBL_EXP_INFNAN && + (px->dbl_frach | px->dbl_fracl) != 0) + return (y); + if (py->dbl_exp == DBL_EXP_INFNAN && + (py->dbl_frach | py->dbl_fracl) != 0) + return (x); + + /* Handle comparisons of signed zeroes. */ + if (px->dbl_sign != py->dbl_sign && + px->dbl_sign == 1) + return (y); + if (px->dbl_sign != py->dbl_sign && + px->dbl_sign == 0) + return (x); + + return (x > y ? x : y); +} diff --git a/lib/libm/src/s_fmaxf.c b/lib/libm/src/s_fmaxf.c new file mode 100644 index 00000000000..472372950d8 --- /dev/null +++ b/lib/libm/src/s_fmaxf.c @@ -0,0 +1,51 @@ +/* $OpenBSD: s_fmaxf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/*- + * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include <sys/types.h> +#include <machine/ieee.h> +#include <math.h> + +float +fmaxf(float x, float y) +{ + struct ieee_single *px = (struct ieee_single *)&x; + struct ieee_single *py = (struct ieee_single *)&y; + + /* Check for NaNs to avoid raising spurious exceptions. */ + if (px->sng_exp == SNG_EXP_INFNAN && px->sng_frac != 0) + return (y); + if (py->sng_exp == SNG_EXP_INFNAN && py->sng_frac != 0) + return (x); + + /* Handle comparisons of signed zeroes. */ + if (px->sng_sign != py->sng_sign && px->sng_sign == 1) + return (y); + if (px->sng_sign != py->sng_sign && px->sng_sign == 0) + return (x); + + return (x > y ? x : y); +} diff --git a/lib/libm/src/s_fmin.c b/lib/libm/src/s_fmin.c new file mode 100644 index 00000000000..47c5d9a67a9 --- /dev/null +++ b/lib/libm/src/s_fmin.c @@ -0,0 +1,55 @@ +/* $OpenBSD: s_fmin.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/*- + * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include <sys/types.h> +#include <machine/ieee.h> +#include <math.h> + +double +fmin(double x, double y) +{ + struct ieee_double *px = (struct ieee_double *)&x; + struct ieee_double *py = (struct ieee_double *)&y; + + /* Check for NaNs to avoid raising spurious exceptions. */ + if (px->dbl_exp == DBL_EXP_INFNAN && + (px->dbl_frach | px->dbl_fracl) != 0) + return (y); + if (py->dbl_exp == DBL_EXP_INFNAN && + (py->dbl_frach | py->dbl_fracl) != 0) + return (x); + + /* Handle comparisons of signed zeroes. */ + if (px->dbl_sign != py->dbl_sign && + py->dbl_sign == 1) + return (y); + if (px->dbl_sign != py->dbl_sign && + py->dbl_sign == 0) + return (x); + + return (x < y ? x : y); +} diff --git a/lib/libm/src/s_fminf.c b/lib/libm/src/s_fminf.c new file mode 100644 index 00000000000..37de5c63e83 --- /dev/null +++ b/lib/libm/src/s_fminf.c @@ -0,0 +1,51 @@ +/* $OpenBSD: s_fminf.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/*- + * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include <sys/types.h> +#include <machine/ieee.h> +#include <math.h> + +float +fminf(float x, float y) +{ + struct ieee_single *px = (struct ieee_single *)&x; + struct ieee_single *py = (struct ieee_single *)&y; + + /* Check for NaNs to avoid raising spurious exceptions. */ + if (px->sng_exp == SNG_EXP_INFNAN && px->sng_frac != 0) + return (y); + if (py->sng_exp == SNG_EXP_INFNAN && py->sng_frac != 0) + return (x); + + /* Handle comparisons of signed zeroes. */ + if (px->sng_sign != py->sng_sign && py->sng_sign == 1) + return (y); + if (px->sng_sign != py->sng_sign && py->sng_sign == 0) + return (x); + + return (x < y ? x : y); +} diff --git a/lib/libm/src/s_lib_version.c b/lib/libm/src/s_lib_version.c deleted file mode 100644 index c4cfae37adb..00000000000 --- a/lib/libm/src/s_lib_version.c +++ /dev/null @@ -1,39 +0,0 @@ -/* @(#)s_lib_ver.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_lib_version.c,v 1.6 1995/05/10 20:47:44 jtc Exp $"; -#endif - -/* - * MACRO for standards - */ - -#include "math.h" -#include "math_private.h" - -/* - * define and initialize _LIB_VERSION - */ -#ifdef _POSIX_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _POSIX_; -#else -#ifdef _XOPEN_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _XOPEN_; -#else -#ifdef _SVID3_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _SVID_; -#else /* default _IEEE_MODE */ -_LIB_VERSION_TYPE _LIB_VERSION = _IEEE_; -#endif -#endif -#endif diff --git a/lib/libm/src/s_matherr.c b/lib/libm/src/s_matherr.c deleted file mode 100644 index 9dfa8c7dd85..00000000000 --- a/lib/libm/src/s_matherr.c +++ /dev/null @@ -1,26 +0,0 @@ -/* @(#)s_matherr.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_matherr.c,v 1.6 1995/05/10 20:47:53 jtc Exp $"; -#endif - -#include "math.h" -#include "math_private.h" - -int -matherr(struct exception *x) -{ - int n=0; - if(x->arg1!=x->arg1) return 0; - return n; -} diff --git a/lib/libm/src/s_nan.c b/lib/libm/src/s_nan.c index f847545722c..3e99cfda129 100644 --- a/lib/libm/src/s_nan.c +++ b/lib/libm/src/s_nan.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_nan.c,v 1.3 2008/08/08 00:41:17 martynas Exp $ */ +/* $OpenBSD: s_nan.c,v 1.4 2008/09/07 20:36:09 martynas Exp $ */ /*- * Copyright (c) 2007 David Schultz * All rights reserved. @@ -40,7 +40,7 @@ /* * OpenBSD's ctype doesn't have digittoint, so we define it here. */ -int +static int _digittoint(int c) { if (!isxdigit(c)) diff --git a/lib/libm/src/s_significand.c b/lib/libm/src/s_significand.c index 1751d59e72b..cf4faa24d04 100644 --- a/lib/libm/src/s_significand.c +++ b/lib/libm/src/s_significand.c @@ -26,5 +26,5 @@ static char rcsid[] = "$NetBSD: s_significand.c,v 1.6 1995/05/10 20:48:11 jtc Ex double significand(double x) { - return __ieee754_scalb(x,(double) -ilogb(x)); + return scalb(x,(double) -ilogb(x)); } diff --git a/lib/libm/src/s_significandf.c b/lib/libm/src/s_significandf.c index 5a7c0965475..085fabc93d2 100644 --- a/lib/libm/src/s_significandf.c +++ b/lib/libm/src/s_significandf.c @@ -23,5 +23,5 @@ static char rcsid[] = "$NetBSD: s_significandf.c,v 1.3 1995/05/10 20:48:13 jtc E float significandf(float x) { - return __ieee754_scalbf(x,(float) -ilogbf(x)); + return scalbf(x,(float) -ilogbf(x)); } diff --git a/lib/libm/src/w_acos.c b/lib/libm/src/w_acos.c deleted file mode 100644 index 142a6f9c076..00000000000 --- a/lib/libm/src/w_acos.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_acos.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_acos.c,v 1.6 1995/05/10 20:48:26 jtc Exp $"; -#endif - -/* - * wrap_acos(x) - */ - -#include "math.h" -#include "math_private.h" - -double -acos(double x) /* wrapper acos */ -{ -#ifdef _IEEE_LIBM - return __ieee754_acos(x); -#else - double z; - z = __ieee754_acos(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(fabs(x)>1.0) { - return __kernel_standard(x,x,1); /* acos(|x|>1) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_acosf.c b/lib/libm/src/w_acosf.c deleted file mode 100644 index ee56e6df6e4..00000000000 --- a/lib/libm/src/w_acosf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_acosf.c -- float version of w_acos.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_acosf.c,v 1.3 1995/05/10 20:48:29 jtc Exp $"; -#endif - -/* - * wrap_acosf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -acosf(float x) /* wrapper acosf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_acosf(x); -#else - float z; - z = __ieee754_acosf(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - if(fabsf(x)>(float)1.0) { - /* acosf(|x|>1) */ - return (float)__kernel_standard((double)x,(double)x,101); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_acosh.c b/lib/libm/src/w_acosh.c deleted file mode 100644 index 545b77bd352..00000000000 --- a/lib/libm/src/w_acosh.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_acosh.c,v 1.6 1995/05/10 20:48:31 jtc Exp $"; -#endif - -/* - * wrapper acosh(x) - */ - -#include "math.h" -#include "math_private.h" - -double -acosh(double x) /* wrapper acosh */ -{ -#ifdef _IEEE_LIBM - return __ieee754_acosh(x); -#else - double z; - z = __ieee754_acosh(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(x<1.0) { - return __kernel_standard(x,x,29); /* acosh(x<1) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_acoshf.c b/lib/libm/src/w_acoshf.c deleted file mode 100644 index 62e48b0f15b..00000000000 --- a/lib/libm/src/w_acoshf.c +++ /dev/null @@ -1,43 +0,0 @@ -/* w_acoshf.c -- float version of w_acosh.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_acoshf.c,v 1.3 1995/05/10 20:48:33 jtc Exp $"; -#endif - -/* - * wrapper acoshf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -acoshf(float x) /* wrapper acoshf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_acoshf(x); -#else - float z; - z = __ieee754_acoshf(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - if(x<(float)1.0) { - /* acosh(x<1) */ - return (float)__kernel_standard((double)x,(double)x,129); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_asin.c b/lib/libm/src/w_asin.c deleted file mode 100644 index f039e729f16..00000000000 --- a/lib/libm/src/w_asin.c +++ /dev/null @@ -1,39 +0,0 @@ -/* @(#)w_asin.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_asin.c,v 1.6 1995/05/10 20:48:35 jtc Exp $"; -#endif - -/* - * wrapper asin(x) - */ - - -#include "math.h" -#include "math_private.h" - -double -asin(double x) /* wrapper asin */ -{ -#ifdef _IEEE_LIBM - return __ieee754_asin(x); -#else - double z; - z = __ieee754_asin(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(fabs(x)>1.0) { - return __kernel_standard(x,x,2); /* asin(|x|>1) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_asinf.c b/lib/libm/src/w_asinf.c deleted file mode 100644 index 3d2c94d3c6c..00000000000 --- a/lib/libm/src/w_asinf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_asinf.c -- float version of w_asin.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_asinf.c,v 1.3 1995/05/10 20:48:37 jtc Exp $"; -#endif - -/* - * wrapper asinf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -asinf(float x) /* wrapper asinf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_asinf(x); -#else - float z; - z = __ieee754_asinf(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - if(fabsf(x)>(float)1.0) { - /* asinf(|x|>1) */ - return (float)__kernel_standard((double)x,(double)x,102); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_atan2.c b/lib/libm/src/w_atan2.c deleted file mode 100644 index 135836af5b2..00000000000 --- a/lib/libm/src/w_atan2.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_atan2.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_atan2.c,v 1.6 1995/05/10 20:48:39 jtc Exp $"; -#endif - -/* - * wrapper atan2(y,x) - */ - -#include "math.h" -#include "math_private.h" - -double -atan2(double y, double x) /* wrapper atan2 */ -{ -#ifdef _IEEE_LIBM - return __ieee754_atan2(y,x); -#else - double z; - z = __ieee754_atan2(y,x); - if(_LIB_VERSION == _IEEE_||isnan(x)||isnan(y)) return z; - if(x==0.0&&y==0.0) { - return __kernel_standard(y,x,3); /* atan2(+-0,+-0) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_atan2f.c b/lib/libm/src/w_atan2f.c deleted file mode 100644 index a1c9763b0ad..00000000000 --- a/lib/libm/src/w_atan2f.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_atan2f.c -- float version of w_atan2.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_atan2f.c,v 1.3 1995/05/10 20:48:42 jtc Exp $"; -#endif - -/* - * wrapper atan2f(y,x) - */ - -#include "math.h" -#include "math_private.h" - -float -atan2f(float y, float x) /* wrapper atan2f */ -{ -#ifdef _IEEE_LIBM - return __ieee754_atan2f(y,x); -#else - float z; - z = __ieee754_atan2f(y,x); - if(_LIB_VERSION == _IEEE_||isnanf(x)||isnanf(y)) return z; - if(x==(float)0.0&&y==(float)0.0) { - /* atan2f(+-0,+-0) */ - return (float)__kernel_standard((double)y,(double)x,103); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_atanh.c b/lib/libm/src/w_atanh.c deleted file mode 100644 index 29cd574c018..00000000000 --- a/lib/libm/src/w_atanh.c +++ /dev/null @@ -1,42 +0,0 @@ -/* @(#)w_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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_atanh.c,v 1.6 1995/05/10 20:48:43 jtc Exp $"; -#endif - -/* - * wrapper atanh(x) - */ - -#include "math.h" -#include "math_private.h" - -double -atanh(double x) /* wrapper atanh */ -{ -#ifdef _IEEE_LIBM - return __ieee754_atanh(x); -#else - double z,y; - z = __ieee754_atanh(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - y = fabs(x); - if(y>=1.0) { - if(y>1.0) - return __kernel_standard(x,x,30); /* atanh(|x|>1) */ - else - return __kernel_standard(x,x,31); /* atanh(|x|==1) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_atanhf.c b/lib/libm/src/w_atanhf.c deleted file mode 100644 index f33d64ebe62..00000000000 --- a/lib/libm/src/w_atanhf.c +++ /dev/null @@ -1,47 +0,0 @@ -/* w_atanhf.c -- float version of w_atanh.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_atanhf.c,v 1.3 1995/05/10 20:48:45 jtc Exp $"; -#endif - -/* - * wrapper atanhf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -atanhf(float x) /* wrapper atanhf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_atanhf(x); -#else - float z,y; - z = __ieee754_atanhf(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - y = fabsf(x); - if(y>=(float)1.0) { - if(y>(float)1.0) - /* atanhf(|x|>1) */ - return (float)__kernel_standard((double)x,(double)x,130); - else - /* atanhf(|x|==1) */ - return (float)__kernel_standard((double)x,(double)x,131); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_cabs.c b/lib/libm/src/w_cabs.c deleted file mode 100644 index d851d339e58..00000000000 --- a/lib/libm/src/w_cabs.c +++ /dev/null @@ -1,19 +0,0 @@ -/* - * cabs() wrapper for hypot(). - * - * Written by J.T. Conklin, <jtc@wimsey.com> - * Placed into the Public Domain, 1994. - */ - -#include <math.h> - -struct complex { - double x; - double y; -}; - -double -cabs(struct complex z) -{ - return hypot(z.x, z.y); -} diff --git a/lib/libm/src/w_cabsf.c b/lib/libm/src/w_cabsf.c deleted file mode 100644 index d98d0422f28..00000000000 --- a/lib/libm/src/w_cabsf.c +++ /dev/null @@ -1,20 +0,0 @@ -/* - * cabsf() wrapper for hypotf(). - * - * Written by J.T. Conklin, <jtc@wimsey.com> - * Placed into the Public Domain, 1994. - */ - -#include "math.h" -#include "math_private.h" - -struct complex { - float x; - float y; -}; - -float -cabsf(struct complex z) -{ - return hypotf(z.x, z.y); -} diff --git a/lib/libm/src/w_cosh.c b/lib/libm/src/w_cosh.c deleted file mode 100644 index 4eb375f6ec1..00000000000 --- a/lib/libm/src/w_cosh.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_cosh.c,v 1.6 1995/05/10 20:48:47 jtc Exp $"; -#endif - -/* - * wrapper cosh(x) - */ - -#include "math.h" -#include "math_private.h" - -double -cosh(double x) /* wrapper cosh */ -{ -#ifdef _IEEE_LIBM - return __ieee754_cosh(x); -#else - double z; - z = __ieee754_cosh(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(fabs(x)>7.10475860073943863426e+02) { - return __kernel_standard(x,x,5); /* cosh overflow */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_coshf.c b/lib/libm/src/w_coshf.c deleted file mode 100644 index df7915d1e72..00000000000 --- a/lib/libm/src/w_coshf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_coshf.c -- float version of w_cosh.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_coshf.c,v 1.3 1995/05/10 20:48:49 jtc Exp $"; -#endif - -/* - * wrapper coshf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -coshf(float x) /* wrapper coshf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_coshf(x); -#else - float z; - z = __ieee754_coshf(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - if(fabsf(x)>(float)8.9415985107e+01) { - /* cosh overflow */ - return (float)__kernel_standard((double)x,(double)x,105); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_exp.c b/lib/libm/src/w_exp.c deleted file mode 100644 index 655b438ca96..00000000000 --- a/lib/libm/src/w_exp.c +++ /dev/null @@ -1,45 +0,0 @@ -/* @(#)w_exp.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_exp.c,v 1.6 1995/05/10 20:48:51 jtc Exp $"; -#endif - -/* - * wrapper exp(x) - */ - -#include "math.h" -#include "math_private.h" - -static const double -o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ -u_threshold= -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */ - -double -exp(double x) /* wrapper exp */ -{ -#ifdef _IEEE_LIBM - return __ieee754_exp(x); -#else - double z; - z = __ieee754_exp(x); - if(_LIB_VERSION == _IEEE_) return z; - if(finite(x)) { - if(x>o_threshold) - return __kernel_standard(x,x,6); /* exp overflow */ - else if(x<u_threshold) - return __kernel_standard(x,x,7); /* exp underflow */ - } - return z; -#endif -} diff --git a/lib/libm/src/w_expf.c b/lib/libm/src/w_expf.c deleted file mode 100644 index f11948b4e74..00000000000 --- a/lib/libm/src/w_expf.c +++ /dev/null @@ -1,50 +0,0 @@ -/* w_expf.c -- float version of w_exp.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_expf.c,v 1.3 1995/05/10 20:48:53 jtc Exp $"; -#endif - -/* - * wrapper expf(x) - */ - -#include "math.h" -#include "math_private.h" - -static const float -o_threshold= 8.8721679688e+01, /* 0x42b17180 */ -u_threshold= -1.0397208405e+02; /* 0xc2cff1b5 */ - -float -expf(float x) /* wrapper expf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_expf(x); -#else - float z; - z = __ieee754_expf(x); - if(_LIB_VERSION == _IEEE_) return z; - if(finitef(x)) { - if(x>o_threshold) - /* exp overflow */ - return (float)__kernel_standard((double)x,(double)x,106); - else if(x<u_threshold) - /* exp underflow */ - return (float)__kernel_standard((double)x,(double)x,107); - } - return z; -#endif -} diff --git a/lib/libm/src/w_fmod.c b/lib/libm/src/w_fmod.c deleted file mode 100644 index 1f1504c13d5..00000000000 --- a/lib/libm/src/w_fmod.c +++ /dev/null @@ -1,39 +0,0 @@ -/* @(#)w_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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_fmod.c,v 1.6 1995/05/10 20:48:55 jtc Exp $"; -#endif - -/* - * wrapper fmod(x,y) - */ - -#include "math.h" -#include "math_private.h" - - -double -fmod(double x, double y) /* wrapper fmod */ -{ -#ifdef _IEEE_LIBM - return __ieee754_fmod(x,y); -#else - double z; - z = __ieee754_fmod(x,y); - if(_LIB_VERSION == _IEEE_ ||isnan(y)||isnan(x)) return z; - if(y==0.0) { - return __kernel_standard(x,y,27); /* fmod(x,0) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_fmodf.c b/lib/libm/src/w_fmodf.c deleted file mode 100644 index c0deb804aba..00000000000 --- a/lib/libm/src/w_fmodf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_fmodf.c -- float version of w_fmod.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_fmodf.c,v 1.3 1995/05/10 20:48:57 jtc Exp $"; -#endif - -/* - * wrapper fmodf(x,y) - */ - -#include "math.h" -#include "math_private.h" - -float -fmodf(float x, float y) /* wrapper fmodf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_fmodf(x,y); -#else - float z; - z = __ieee754_fmodf(x,y); - if(_LIB_VERSION == _IEEE_ ||isnanf(y)||isnanf(x)) return z; - if(y==(float)0.0) { - /* fmodf(x,0) */ - return (float)__kernel_standard((double)x,(double)y,127); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_gamma.c b/lib/libm/src/w_gamma.c index baa7c711248..92112add553 100644 --- a/lib/libm/src/w_gamma.c +++ b/lib/libm/src/w_gamma.c @@ -28,18 +28,5 @@ extern int signgam; double gamma(double x) { -#ifdef _IEEE_LIBM - return __ieee754_lgamma_r(x,&signgam); -#else - double y; - y = __ieee754_lgamma_r(x,&signgam); - if(_LIB_VERSION == _IEEE_) return y; - if(!finite(y)&&finite(x)) { - if(floor(x)==x&&x<=0.0) - return __kernel_standard(x,x,41); /* gamma pole */ - else - return __kernel_standard(x,x,40); /* gamma overflow */ - } else - return y; -#endif + return lgamma_r(x,&signgam); } diff --git a/lib/libm/src/w_gamma_r.c b/lib/libm/src/w_gamma_r.c index 8157923f203..dcefb1f4832 100644 --- a/lib/libm/src/w_gamma_r.c +++ b/lib/libm/src/w_gamma_r.c @@ -24,18 +24,5 @@ static char rcsid[] = "$NetBSD: w_gamma_r.c,v 1.7 1995/11/20 22:06:45 jtc Exp $" double gamma_r(double x, int *signgamp) /* wrapper lgamma_r */ { -#ifdef _IEEE_LIBM - return __ieee754_lgamma_r(x,signgamp); -#else - double y; - y = __ieee754_lgamma_r(x,signgamp); - if(_LIB_VERSION == _IEEE_) return y; - if(!finite(y)&&finite(x)) { - if(floor(x)==x&&x<=0.0) - return __kernel_standard(x,x,41); /* gamma pole */ - else - return __kernel_standard(x,x,40); /* gamma overflow */ - } else - return y; -#endif + return lgamma_r(x,signgamp); } diff --git a/lib/libm/src/w_gammaf.c b/lib/libm/src/w_gammaf.c deleted file mode 100644 index cf668d2bf10..00000000000 --- a/lib/libm/src/w_gammaf.c +++ /dev/null @@ -1,44 +0,0 @@ -/* w_gammaf.c -- float version of w_gamma.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_gammaf.c,v 1.4 1995/11/20 22:06:48 jtc Exp $"; -#endif - -#include "math.h" -#include "math_private.h" - -extern int signgam; - -float -gammaf(float x) -{ -#ifdef _IEEE_LIBM - return __ieee754_lgammaf_r(x,&signgam); -#else - float y; - y = __ieee754_lgammaf_r(x,&signgam); - if(_LIB_VERSION == _IEEE_) return y; - if(!finitef(y)&&finitef(x)) { - if(floorf(x)==x&&x<=(float)0.0) - /* gammaf pole */ - return (float)__kernel_standard((double)x,(double)x,141); - else - /* gammaf overflow */ - return (float)__kernel_standard((double)x,(double)x,140); - } else - return y; -#endif -} diff --git a/lib/libm/src/w_gammaf_r.c b/lib/libm/src/w_gammaf_r.c deleted file mode 100644 index 16f4b43b17f..00000000000 --- a/lib/libm/src/w_gammaf_r.c +++ /dev/null @@ -1,46 +0,0 @@ -/* w_gammaf_r.c -- float version of w_gamma_r.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_gammaf_r.c,v 1.4 1995/11/20 22:06:50 jtc Exp $"; -#endif - -/* - * wrapper float gammaf_r(float x, int *signgamp) - */ - -#include "math.h" -#include "math_private.h" - -float -gammaf_r(float x, int *signgamp) /* wrapper lgammaf_r */ -{ -#ifdef _IEEE_LIBM - return __ieee754_lgammaf_r(x,signgamp); -#else - float y; - y = __ieee754_lgammaf_r(x,signgamp); - if(_LIB_VERSION == _IEEE_) return y; - if(!finitef(y)&&finitef(x)) { - if(floorf(x)==x&&x<=(float)0.0) - /* gammaf pole */ - return (float)__kernel_standard((double)x,(double)x,141); - else - /* gamma overflow */ - return (float)__kernel_standard((double)x,(double)x,140); - } else - return y; -#endif -} diff --git a/lib/libm/src/w_hypot.c b/lib/libm/src/w_hypot.c deleted file mode 100644 index 57bdf938dbe..00000000000 --- a/lib/libm/src/w_hypot.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_hypot.c,v 1.6 1995/05/10 20:49:07 jtc Exp $"; -#endif - -/* - * wrapper hypot(x,y) - */ - -#include "math.h" -#include "math_private.h" - -double -hypot(double x, double y)/* wrapper hypot */ -{ -#ifdef _IEEE_LIBM - return __ieee754_hypot(x,y); -#else - double z; - z = __ieee754_hypot(x,y); - if(_LIB_VERSION == _IEEE_) return z; - if((!finite(z))&&finite(x)&&finite(y)) - return __kernel_standard(x,y,4); /* hypot overflow */ - else - return z; -#endif -} diff --git a/lib/libm/src/w_hypotf.c b/lib/libm/src/w_hypotf.c deleted file mode 100644 index d0f83fc7115..00000000000 --- a/lib/libm/src/w_hypotf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_hypotf.c -- float version of w_hypot.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_hypotf.c,v 1.3 1995/05/10 20:49:09 jtc Exp $"; -#endif - -/* - * wrapper hypotf(x,y) - */ - -#include "math.h" -#include "math_private.h" - -float -hypotf(float x, float y) /* wrapper hypotf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_hypotf(x,y); -#else - float z; - z = __ieee754_hypotf(x,y); - if(_LIB_VERSION == _IEEE_) return z; - if((!finitef(z))&&finitef(x)&&finitef(y)) - /* hypot overflow */ - return (float)__kernel_standard((double)x,(double)y,104); - else - return z; -#endif -} diff --git a/lib/libm/src/w_j0.c b/lib/libm/src/w_j0.c deleted file mode 100644 index 6e617427920..00000000000 --- a/lib/libm/src/w_j0.c +++ /dev/null @@ -1,61 +0,0 @@ -/* @(#)w_j0.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_j0.c,v 1.6 1995/05/10 20:49:11 jtc Exp $"; -#endif - -/* - * wrapper j0(double x), y0(double x) - */ - -#include "math.h" -#include "math_private.h" - -double -j0(double x) /* wrapper j0 */ -{ -#ifdef _IEEE_LIBM - return __ieee754_j0(x); -#else - double z = __ieee754_j0(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(fabs(x)>X_TLOSS) { - return __kernel_standard(x,x,34); /* j0(|x|>X_TLOSS) */ - } else - return z; -#endif -} - -double -y0(double x) /* wrapper y0 */ -{ -#ifdef _IEEE_LIBM - return __ieee754_y0(x); -#else - double z; - z = __ieee754_y0(x); - if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; - if(x <= 0.0){ - if(x==0.0) - /* d= -one/(x-x); */ - return __kernel_standard(x,x,8); - else - /* d = zero/(x-x); */ - return __kernel_standard(x,x,9); - } - if(x>X_TLOSS) { - return __kernel_standard(x,x,35); /* y0(x>X_TLOSS) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_j0f.c b/lib/libm/src/w_j0f.c deleted file mode 100644 index 93ec42b366a..00000000000 --- a/lib/libm/src/w_j0f.c +++ /dev/null @@ -1,66 +0,0 @@ -/* w_j0f.c -- float version of w_j0.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_j0f.c,v 1.3 1995/05/10 20:49:13 jtc Exp $"; -#endif - -/* - * wrapper j0f(float x), y0f(float x) - */ - -#include "math.h" -#include "math_private.h" - -float -j0f(float x) /* wrapper j0f */ -{ -#ifdef _IEEE_LIBM - return __ieee754_j0f(x); -#else - float z = __ieee754_j0f(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - if(fabsf(x)>(float)X_TLOSS) { - /* j0f(|x|>X_TLOSS) */ - return (float)__kernel_standard((double)x,(double)x,134); - } else - return z; -#endif -} - -float -y0f(float x) /* wrapper y0f */ -{ -#ifdef _IEEE_LIBM - return __ieee754_y0f(x); -#else - float z; - z = __ieee754_y0f(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x) ) return z; - if(x <= (float)0.0){ - if(x==(float)0.0) - /* d= -one/(x-x); */ - return (float)__kernel_standard((double)x,(double)x,108); - else - /* d = zero/(x-x); */ - return (float)__kernel_standard((double)x,(double)x,109); - } - if(x>(float)X_TLOSS) { - /* y0(x>X_TLOSS) */ - return (float)__kernel_standard((double)x,(double)x,135); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_j1.c b/lib/libm/src/w_j1.c deleted file mode 100644 index 23ca8befaf6..00000000000 --- a/lib/libm/src/w_j1.c +++ /dev/null @@ -1,62 +0,0 @@ -/* @(#)w_j1.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_j1.c,v 1.6 1995/05/10 20:49:15 jtc Exp $"; -#endif - -/* - * wrapper of j1,y1 - */ - -#include "math.h" -#include "math_private.h" - -double -j1(double x) /* wrapper j1 */ -{ -#ifdef _IEEE_LIBM - return __ieee754_j1(x); -#else - double z; - z = __ieee754_j1(x); - if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; - if(fabs(x)>X_TLOSS) { - return __kernel_standard(x,x,36); /* j1(|x|>X_TLOSS) */ - } else - return z; -#endif -} - -double -y1(double x) /* wrapper y1 */ -{ -#ifdef _IEEE_LIBM - return __ieee754_y1(x); -#else - double z; - z = __ieee754_y1(x); - if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; - if(x <= 0.0){ - if(x==0.0) - /* d= -one/(x-x); */ - return __kernel_standard(x,x,10); - else - /* d = zero/(x-x); */ - return __kernel_standard(x,x,11); - } - if(x>X_TLOSS) { - return __kernel_standard(x,x,37); /* y1(x>X_TLOSS) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_j1f.c b/lib/libm/src/w_j1f.c deleted file mode 100644 index 770c9ceab61..00000000000 --- a/lib/libm/src/w_j1f.c +++ /dev/null @@ -1,67 +0,0 @@ -/* w_j1f.c -- float version of w_j1.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_j1f.c,v 1.3 1995/05/10 20:49:17 jtc Exp $"; -#endif - -/* - * wrapper of j1f,y1f - */ - -#include "math.h" -#include "math_private.h" - -float -j1f(float x) /* wrapper j1f */ -{ -#ifdef _IEEE_LIBM - return __ieee754_j1f(x); -#else - float z; - z = __ieee754_j1f(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x) ) return z; - if(fabsf(x)>(float)X_TLOSS) { - /* j1(|x|>X_TLOSS) */ - return (float)__kernel_standard((double)x,(double)x,136); - } else - return z; -#endif -} - -float -y1f(float x) /* wrapper y1f */ -{ -#ifdef _IEEE_LIBM - return __ieee754_y1f(x); -#else - float z; - z = __ieee754_y1f(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x) ) return z; - if(x <= (float)0.0){ - if(x==(float)0.0) - /* d= -one/(x-x); */ - return (float)__kernel_standard((double)x,(double)x,110); - else - /* d = zero/(x-x); */ - return (float)__kernel_standard((double)x,(double)x,111); - } - if(x>(float)X_TLOSS) { - /* y1(x>X_TLOSS) */ - return (float)__kernel_standard((double)x,(double)x,137); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_jn.c b/lib/libm/src/w_jn.c deleted file mode 100644 index 2ea879e0bf3..00000000000 --- a/lib/libm/src/w_jn.c +++ /dev/null @@ -1,84 +0,0 @@ -/* @(#)w_jn.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_jn.c,v 1.6 1995/05/10 20:49:19 jtc Exp $"; -#endif - -/* - * wrapper jn(int n, double x), yn(int n, double x) - * floating point Bessel's function of the 1st and 2nd kind - * of order n - * - * Special cases: - * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; - * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. - * Note 2. About jn(n,x), yn(n,x) - * For n=0, j0(x) is called, - * for n=1, j1(x) is called, - * for n<x, forward recursion us used starting - * from values of j0(x) and j1(x). - * for n>x, a continued fraction approximation to - * j(n,x)/j(n-1,x) is evaluated and then backward - * recursion is used starting from a supposed value - * for j(n,x). The resulting value of j(0,x) is - * compared with the actual value to correct the - * supposed value of j(n,x). - * - * yn(n,x) is similar in all respects, except - * that forward recursion is used for all - * values of n>1. - * - */ - -#include "math.h" -#include "math_private.h" - -double -jn(int n, double x) /* wrapper jn */ -{ -#ifdef _IEEE_LIBM - return __ieee754_jn(n,x); -#else - double z; - z = __ieee754_jn(n,x); - if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; - if(fabs(x)>X_TLOSS) { - return __kernel_standard((double)n,x,38); /* jn(|x|>X_TLOSS,n) */ - } else - return z; -#endif -} - -double -yn(int n, double x) /* wrapper yn */ -{ -#ifdef _IEEE_LIBM - return __ieee754_yn(n,x); -#else - double z; - z = __ieee754_yn(n,x); - if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; - if(x <= 0.0){ - if(x==0.0) - /* d= -one/(x-x); */ - return __kernel_standard((double)n,x,12); - else - /* d = zero/(x-x); */ - return __kernel_standard((double)n,x,13); - } - if(x>X_TLOSS) { - return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_jnf.c b/lib/libm/src/w_jnf.c deleted file mode 100644 index a45f3b878f5..00000000000 --- a/lib/libm/src/w_jnf.c +++ /dev/null @@ -1,63 +0,0 @@ -/* w_jnf.c -- float version of w_jn.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_jnf.c,v 1.3 1995/05/10 20:49:21 jtc Exp $"; -#endif - -#include "math.h" -#include "math_private.h" - -float -jnf(int n, float x) /* wrapper jnf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_jnf(n,x); -#else - float z; - z = __ieee754_jnf(n,x); - if(_LIB_VERSION == _IEEE_ || isnanf(x) ) return z; - if(fabsf(x)>(float)X_TLOSS) { - /* jn(|x|>X_TLOSS,n) */ - return (float)__kernel_standard((double)n,(double)x,138); - } else - return z; -#endif -} - -float -ynf(int n, float x) /* wrapper ynf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_ynf(n,x); -#else - float z; - z = __ieee754_ynf(n,x); - if(_LIB_VERSION == _IEEE_ || isnanf(x) ) return z; - if(x <= (float)0.0){ - if(x==(float)0.0) - /* d= -one/(x-x); */ - return (float)__kernel_standard((double)n,(double)x,112); - else - /* d = zero/(x-x); */ - return (float)__kernel_standard((double)n,(double)x,113); - } - if(x>(float)X_TLOSS) { - /* yn(x>X_TLOSS,n) */ - return (float)__kernel_standard((double)n,(double)x,139); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_lgamma.c b/lib/libm/src/w_lgamma.c deleted file mode 100644 index 055fdb774fa..00000000000 --- a/lib/libm/src/w_lgamma.c +++ /dev/null @@ -1,45 +0,0 @@ -/* @(#)w_lgamma.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_lgamma.c,v 1.6 1995/05/10 20:49:24 jtc Exp $"; -#endif - -/* double lgamma(double x) - * Return the logarithm of the Gamma function of x. - * - * Method: call __ieee754_lgamma_r - */ - -#include "math.h" -#include "math_private.h" - -extern int signgam; - -double -lgamma(double x) -{ -#ifdef _IEEE_LIBM - return __ieee754_lgamma_r(x,&signgam); -#else - double y; - y = __ieee754_lgamma_r(x,&signgam); - if(_LIB_VERSION == _IEEE_) return y; - if(!finite(y)&&finite(x)) { - if(floor(x)==x&&x<=0.0) - return __kernel_standard(x,x,15); /* lgamma pole */ - else - return __kernel_standard(x,x,14); /* lgamma overflow */ - } else - return y; -#endif -} diff --git a/lib/libm/src/w_lgamma_r.c b/lib/libm/src/w_lgamma_r.c deleted file mode 100644 index 5e9d69ce239..00000000000 --- a/lib/libm/src/w_lgamma_r.c +++ /dev/null @@ -1,42 +0,0 @@ -/* @(#)wr_lgamma.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_lgamma_r.c,v 1.6 1995/05/10 20:49:27 jtc Exp $"; -#endif - -/* - * wrapper double lgamma_r(double x, int *signgamp) - */ - -#include "math.h" -#include "math_private.h" - - -double -lgamma_r(double x, int *signgamp) /* wrapper lgamma_r */ -{ -#ifdef _IEEE_LIBM - return __ieee754_lgamma_r(x,signgamp); -#else - double y; - y = __ieee754_lgamma_r(x,signgamp); - if(_LIB_VERSION == _IEEE_) return y; - if(!finite(y)&&finite(x)) { - if(floor(x)==x&&x<=0.0) - return __kernel_standard(x,x,15); /* lgamma pole */ - else - return __kernel_standard(x,x,14); /* lgamma overflow */ - } else - return y; -#endif -} diff --git a/lib/libm/src/w_lgammaf.c b/lib/libm/src/w_lgammaf.c deleted file mode 100644 index 0ec7716761f..00000000000 --- a/lib/libm/src/w_lgammaf.c +++ /dev/null @@ -1,44 +0,0 @@ -/* w_lgammaf.c -- float version of w_lgamma.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_lgammaf.c,v 1.3 1995/05/10 20:49:30 jtc Exp $"; -#endif - -#include "math.h" -#include "math_private.h" - -extern int signgam; - -float -lgammaf(float x) -{ -#ifdef _IEEE_LIBM - return __ieee754_lgammaf_r(x,&signgam); -#else - float y; - y = __ieee754_lgammaf_r(x,&signgam); - if(_LIB_VERSION == _IEEE_) return y; - if(!finitef(y)&&finitef(x)) { - if(floorf(x)==x&&x<=(float)0.0) - /* lgamma pole */ - return (float)__kernel_standard((double)x,(double)x,115); - else - /* lgamma overflow */ - return (float)__kernel_standard((double)x,(double)x,114); - } else - return y; -#endif -} diff --git a/lib/libm/src/w_lgammaf_r.c b/lib/libm/src/w_lgammaf_r.c deleted file mode 100644 index 49973d3637d..00000000000 --- a/lib/libm/src/w_lgammaf_r.c +++ /dev/null @@ -1,47 +0,0 @@ -/* w_lgammaf_r.c -- float version of w_lgamma_r.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_lgammaf_r.c,v 1.3 1995/05/10 20:49:32 jtc Exp $"; -#endif - -/* - * wrapper float lgammaf_r(float x, int *signgamp) - */ - -#include "math.h" -#include "math_private.h" - - -float -lgammaf_r(float x, int *signgamp) /* wrapper lgammaf_r */ -{ -#ifdef _IEEE_LIBM - return __ieee754_lgammaf_r(x,signgamp); -#else - float y; - y = __ieee754_lgammaf_r(x,signgamp); - if(_LIB_VERSION == _IEEE_) return y; - if(!finitef(y)&&finitef(x)) { - if(floorf(x)==x&&x<=(float)0.0) - /* lgamma pole */ - return (float)__kernel_standard((double)x,(double)x,115); - else - /* lgamma overflow */ - return (float)__kernel_standard((double)x,(double)x,114); - } else - return y; -#endif -} diff --git a/lib/libm/src/w_log.c b/lib/libm/src/w_log.c deleted file mode 100644 index b4dcdef1f05..00000000000 --- a/lib/libm/src/w_log.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_log.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_log.c,v 1.6 1995/05/10 20:49:33 jtc Exp $"; -#endif - -/* - * wrapper log(x) - */ - -#include "math.h" -#include "math_private.h" - -double -log(double x) /* wrapper log */ -{ -#ifdef _IEEE_LIBM - return __ieee754_log(x); -#else - double z; - z = __ieee754_log(x); - if(_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0) return z; - if(x==0.0) - return __kernel_standard(x,x,16); /* log(0) */ - else - return __kernel_standard(x,x,17); /* log(x<0) */ -#endif -} diff --git a/lib/libm/src/w_log10.c b/lib/libm/src/w_log10.c deleted file mode 100644 index f0b4e7c1740..00000000000 --- a/lib/libm/src/w_log10.c +++ /dev/null @@ -1,41 +0,0 @@ -/* @(#)w_log10.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_log10.c,v 1.6 1995/05/10 20:49:35 jtc Exp $"; -#endif - -/* - * wrapper log10(X) - */ - -#include "math.h" -#include "math_private.h" - -double -log10(double x) /* wrapper log10 */ -{ -#ifdef _IEEE_LIBM - return __ieee754_log10(x); -#else - double z; - z = __ieee754_log10(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(x<=0.0) { - if(x==0.0) - return __kernel_standard(x,x,18); /* log10(0) */ - else - return __kernel_standard(x,x,19); /* log10(x<0) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_log10f.c b/lib/libm/src/w_log10f.c deleted file mode 100644 index 4907438d820..00000000000 --- a/lib/libm/src/w_log10f.c +++ /dev/null @@ -1,46 +0,0 @@ -/* w_log10f.c -- float version of w_log10.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_log10f.c,v 1.3 1995/05/10 20:49:37 jtc Exp $"; -#endif - -/* - * wrapper log10f(X) - */ - -#include "math.h" -#include "math_private.h" - -float -log10f(float x) /* wrapper log10f */ -{ -#ifdef _IEEE_LIBM - return __ieee754_log10f(x); -#else - float z; - z = __ieee754_log10f(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - if(x<=(float)0.0) { - if(x==(float)0.0) - /* log10(0) */ - return (float)__kernel_standard((double)x,(double)x,118); - else - /* log10(x<0) */ - return (float)__kernel_standard((double)x,(double)x,119); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_logf.c b/lib/libm/src/w_logf.c deleted file mode 100644 index 5db12a6beaa..00000000000 --- a/lib/libm/src/w_logf.c +++ /dev/null @@ -1,43 +0,0 @@ -/* w_logf.c -- float version of w_log.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_logf.c,v 1.3 1995/05/10 20:49:40 jtc Exp $"; -#endif - -/* - * wrapper logf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -logf(float x) /* wrapper logf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_logf(x); -#else - float z; - z = __ieee754_logf(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x) || x > (float)0.0) return z; - if(x==(float)0.0) - /* logf(0) */ - return (float)__kernel_standard((double)x,(double)x,116); - else - /* logf(x<0) */ - return (float)__kernel_standard((double)x,(double)x,117); -#endif -} diff --git a/lib/libm/src/w_pow.c b/lib/libm/src/w_pow.c deleted file mode 100644 index d1bf6374b9c..00000000000 --- a/lib/libm/src/w_pow.c +++ /dev/null @@ -1,56 +0,0 @@ - - -/* @(#)w_pow.c 5.2 93/10/01 */ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -/* - * wrapper pow(x,y) return x**y - */ - -#include "math.h" -#include "math_private.h" - -double -pow(double x, double y) /* wrapper pow */ -{ -#ifdef _IEEE_LIBM - return __ieee754_pow(x,y); -#else - double z; - z=__ieee754_pow(x,y); - if(_LIB_VERSION == _IEEE_|| isnan(y)) return z; - if(isnan(x)) { - if(y==0.0) - return __kernel_standard(x,y,42); /* pow(NaN,0.0) */ - else - return z; - } - if(x==0.0){ - if(y==0.0) - return __kernel_standard(x,y,20); /* pow(0.0,0.0) */ - if(finite(y)&&y<0.0) - return __kernel_standard(x,y,23); /* pow(0.0,negative) */ - return z; - } - if(!finite(z)) { - if(finite(x)&&finite(y)) { - if(isnan(z)) - return __kernel_standard(x,y,24); /* pow neg**non-int */ - else - return __kernel_standard(x,y,21); /* pow overflow */ - } - } - if(z==0.0&&finite(x)&&finite(y)) - return __kernel_standard(x,y,22); /* pow underflow */ - return z; -#endif -} diff --git a/lib/libm/src/w_powf.c b/lib/libm/src/w_powf.c deleted file mode 100644 index 6d4772a3830..00000000000 --- a/lib/libm/src/w_powf.c +++ /dev/null @@ -1,67 +0,0 @@ -/* w_powf.c -- float version of w_pow.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_powf.c,v 1.3 1995/05/10 20:49:41 jtc Exp $"; -#endif - -/* - * wrapper powf(x,y) return x**y - */ - -#include "math.h" -#include "math_private.h" - -float -powf(float x, float y) /* wrapper powf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_powf(x,y); -#else - float z; - z=__ieee754_powf(x,y); - if(_LIB_VERSION == _IEEE_|| isnanf(y)) return z; - if(isnanf(x)) { - if(y==(float)0.0) - /* powf(NaN,0.0) */ - return (float)__kernel_standard((double)x,(double)y,142); - else - return z; - } - if(x==(float)0.0){ - if(y==(float)0.0) - /* powf(0.0,0.0) */ - return (float)__kernel_standard((double)x,(double)y,120); - if(finitef(y)&&y<(float)0.0) - /* powf(0.0,negative) */ - return (float)__kernel_standard((double)x,(double)y,123); - return z; - } - if(!finitef(z)) { - if(finitef(x)&&finitef(y)) { - if(isnanf(z)) - /* powf neg**non-int */ - return (float)__kernel_standard((double)x,(double)y,124); - else - /* powf overflow */ - return (float)__kernel_standard((double)x,(double)y,121); - } - } - if(z==(float)0.0&&finitef(x)&&finitef(y)) - /* powf underflow */ - return (float)__kernel_standard((double)x,(double)y,122); - return z; -#endif -} diff --git a/lib/libm/src/w_remainder.c b/lib/libm/src/w_remainder.c deleted file mode 100644 index 50a41e77c4b..00000000000 --- a/lib/libm/src/w_remainder.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_remainder.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_remainder.c,v 1.6 1995/05/10 20:49:44 jtc Exp $"; -#endif - -/* - * wrapper remainder(x,p) - */ - -#include "math.h" -#include "math_private.h" - -double -remainder(double x, double y) /* wrapper remainder */ -{ -#ifdef _IEEE_LIBM - return __ieee754_remainder(x,y); -#else - double z; - z = __ieee754_remainder(x,y); - if(_LIB_VERSION == _IEEE_ || isnan(y)) return z; - if(y==0.0) - return __kernel_standard(x,y,28); /* remainder(x,0) */ - else - return z; -#endif -} diff --git a/lib/libm/src/w_remainderf.c b/lib/libm/src/w_remainderf.c deleted file mode 100644 index 0339c546a94..00000000000 --- a/lib/libm/src/w_remainderf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_remainderf.c -- float version of w_remainder.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_remainderf.c,v 1.3 1995/05/10 20:49:46 jtc Exp $"; -#endif - -/* - * wrapper remainderf(x,p) - */ - -#include "math.h" -#include "math_private.h" - -float -remainderf(float x, float y) /* wrapper remainder */ -{ -#ifdef _IEEE_LIBM - return __ieee754_remainderf(x,y); -#else - float z; - z = __ieee754_remainderf(x,y); - if(_LIB_VERSION == _IEEE_ || isnanf(y)) return z; - if(y==(float)0.0) - /* remainder(x,0) */ - return (float)__kernel_standard((double)x,(double)y,128); - else - return z; -#endif -} diff --git a/lib/libm/src/w_scalb.c b/lib/libm/src/w_scalb.c deleted file mode 100644 index d956ab75ee8..00000000000 --- a/lib/libm/src/w_scalb.c +++ /dev/null @@ -1,52 +0,0 @@ -/* @(#)w_scalb.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_scalb.c,v 1.6 1995/05/10 20:49:48 jtc Exp $"; -#endif - -/* - * wrapper scalb(double x, double fn) is provide for - * passing various standard test suite. One - * should use scalbn() instead. - */ - -#include "math.h" -#include "math_private.h" - -#include <errno.h> - -double -#ifdef _SCALB_INT -scalb(double x, int fn) /* wrapper scalb */ -#else -scalb(double x, double fn) /* wrapper scalb */ -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_scalb(x,fn); -#else - double z; - z = __ieee754_scalb(x,fn); - if(_LIB_VERSION == _IEEE_) return z; - if(!(finite(z)||isnan(z))&&finite(x)) { - return __kernel_standard(x,(double)fn,32); /* scalb overflow */ - } - if(z==0.0&&z!=x) { - return __kernel_standard(x,(double)fn,33); /* scalb underflow */ - } -#ifndef _SCALB_INT - if(!finite(fn)) errno = ERANGE; -#endif - return z; -#endif -} diff --git a/lib/libm/src/w_scalbf.c b/lib/libm/src/w_scalbf.c deleted file mode 100644 index 9542611f1b6..00000000000 --- a/lib/libm/src/w_scalbf.c +++ /dev/null @@ -1,57 +0,0 @@ -/* w_scalbf.c -- float version of w_scalb.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_scalbf.c,v 1.3 1995/05/10 20:49:50 jtc Exp $"; -#endif - -/* - * wrapper scalbf(float x, float fn) is provide for - * passing various standard test suite. One - * should use scalbn() instead. - */ - -#include "math.h" -#include "math_private.h" - -#include <errno.h> - -float -#ifdef _SCALB_INT -scalbf(float x, int fn) /* wrapper scalbf */ -#else -scalbf(float x, float fn) /* wrapper scalbf */ -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_scalbf(x,fn); -#else - float z; - z = __ieee754_scalbf(x,fn); - if(_LIB_VERSION == _IEEE_) return z; - if(!(finitef(z)||isnanf(z))&&finitef(x)) { - /* scalbf overflow */ - return (float)__kernel_standard((double)x,(double)fn,132); - } - if(z==(float)0.0&&z!=x) { - /* scalbf underflow */ - return (float)__kernel_standard((double)x,(double)fn,133); - } -#ifndef _SCALB_INT - if(!finitef(fn)) errno = ERANGE; -#endif - return z; -#endif -} diff --git a/lib/libm/src/w_sinh.c b/lib/libm/src/w_sinh.c deleted file mode 100644 index c02a0ae0531..00000000000 --- a/lib/libm/src/w_sinh.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_sinh.c,v 1.6 1995/05/10 20:49:51 jtc Exp $"; -#endif - -/* - * wrapper sinh(x) - */ - -#include "math.h" -#include "math_private.h" - -double -sinh(double x) /* wrapper sinh */ -{ -#ifdef _IEEE_LIBM - return __ieee754_sinh(x); -#else - double z; - z = __ieee754_sinh(x); - if(_LIB_VERSION == _IEEE_) return z; - if(!finite(z)&&finite(x)) { - return __kernel_standard(x,x,25); /* sinh overflow */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_sinhf.c b/lib/libm/src/w_sinhf.c deleted file mode 100644 index 58df8cd515a..00000000000 --- a/lib/libm/src/w_sinhf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_sinhf.c -- float version of w_sinh.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_sinhf.c,v 1.3 1995/05/10 20:49:54 jtc Exp $"; -#endif - -/* - * wrapper sinhf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -sinhf(float x) /* wrapper sinhf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_sinhf(x); -#else - float z; - z = __ieee754_sinhf(x); - if(_LIB_VERSION == _IEEE_) return z; - if(!finitef(z)&&finitef(x)) { - /* sinhf overflow */ - return (float)__kernel_standard((double)x,(double)x,125); - } else - return z; -#endif -} diff --git a/lib/libm/src/w_sqrt.c b/lib/libm/src/w_sqrt.c deleted file mode 100644 index fec804b282c..00000000000 --- a/lib/libm/src/w_sqrt.c +++ /dev/null @@ -1,38 +0,0 @@ -/* @(#)w_sqrt.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_sqrt.c,v 1.6 1995/05/10 20:49:55 jtc Exp $"; -#endif - -/* - * wrapper sqrt(x) - */ - -#include "math.h" -#include "math_private.h" - -double -sqrt(double x) /* wrapper sqrt */ -{ -#ifdef _IEEE_LIBM - return __ieee754_sqrt(x); -#else - double z; - z = __ieee754_sqrt(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(x<0.0) { - return __kernel_standard(x,x,26); /* sqrt(negative) */ - } else - return z; -#endif -} diff --git a/lib/libm/src/w_sqrtf.c b/lib/libm/src/w_sqrtf.c deleted file mode 100644 index fb21a63529c..00000000000 --- a/lib/libm/src/w_sqrtf.c +++ /dev/null @@ -1,42 +0,0 @@ -/* w_sqrtf.c -- float version of w_sqrt.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_sqrtf.c,v 1.3 1995/05/10 20:49:59 jtc Exp $"; -#endif - -/* - * wrapper sqrtf(x) - */ - -#include "math.h" -#include "math_private.h" - -float -sqrtf(float x) /* wrapper sqrtf */ -{ -#ifdef _IEEE_LIBM - return __ieee754_sqrtf(x); -#else - float z; - z = __ieee754_sqrtf(x); - if(_LIB_VERSION == _IEEE_ || isnanf(x)) return z; - if(x<(float)0.0) { - /* sqrtf(negative) */ - return (float)__kernel_standard((double)x,(double)x,126); - } else - return z; -#endif -} |