n--;
}
n &= 7;
- y = y - z + n * 0.25L;
+ y = y - z + n * 0.25;
switch (n) {
case 0: y = __kernel_sinl(pi*y,zero,0); break;
case 1:
case 2: y = __kernel_cosl(pi*(0.5-y),zero); break;
- case 3:
+ case 3:
case 4: y = __kernel_sinl(pi*(one-y),zero,0); break;
case 5:
case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break;
return -y;
}
-
long double
lgammal_r(long double x, int *signgamp)
{
long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
uint64_t llx,lx;
int i;
- uint16_t hx;
+ uint16_t hx,ix;
- EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
-
- if((hx & 0x7fff) == 0x7fff) { /* erfl(nan)=nan */
- i = (hx>>15)<<1;
- return (1-i)+one/x; /* erfl(+-inf)=+-1 */
- }
+ EXTRACT_LDBL128_WORDS(hx,lx,llx,x);
- /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+ /* purge +-Inf and NaNs */
*signgamp = 1;
- if((hx & 0x7fff) == 0x7fff) /* x is +-Inf or NaN */
- return x*x;
- if((hx==0||hx==0x8000)&&lx==0) {
- if (hx&0x8000)
- *signgamp = -1;
- return one/vzero;
- }
+ ix = hx&0x7fff;
+ if(ix==0x7fff) return x*x;
- /* purge off tiny and negative arguments */
- if(fabsl(x)<0x1p-119L) {
- if(hx&0x8000) {
- *signgamp = -1;
- return -logl(-x);
- } else return -logl(x);
+ /* purge +-0 and tiny arguments */
+ *signgamp = 1-2*(hx>>15);
+ if(ix<0x3fff-116) { /* |x|<2**-(p+3), return -log(|x|) */
+ if((ix|lx|llx)==0)
+ return one/vzero;
+ return -logl(fabsl(x));
}
+
+ /* purge negative integers and start evaluation for other x < 0 */
if(hx&0x8000) {
- if(fabsl(x)>=0x1p112)
+ *signgamp = 1;
+ if(ix>=0x3fff+112) /* |x|>=2**(p-1), must be -integer */
return one/vzero;
t = sin_pil(x);
if(t==zero) return one/vzero;
x = -x;
}
- if(x == 1 || x ==2) r = 0;
- else if(x<2) {
- if(x<=0.8999996185302734) {
+ /* purge 1 and 2 */
+ if((ix==0x3fff || ix==0x4000) && (lx|llx)==0) r = 0;
+ /* for x < 2.0 */
+ else if(ix<0x4000) {
+ if(x<=8.9999961853027344e-01) {
r = -logl(x);
- if(x>=0.7315998077392578) {y = 1-x; i= 0;}
- else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;}
+ if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;}
+ else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;}
else {y = x; i=2;}
} else {
- r = 0;
- if(x>=1.7316312789916992) {y=2-x;i=0;}
- else if(x>=1.2316322326660156) {y=x-tc;i=1;}
+ r = 0;
+ if(x>=1.7316312789916992e+00) {y=2-x;i=0;}
+ else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;}
else {y=x-1;i=2;}
}
switch(i) {
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+
z*(a17+z*(a19+z*(a21+z*a23)))))))))));
p = y*p1+p2;
- r += (p-y/2); break;
+ r += p-y/2; break;
case 1:
p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+
y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+
y*(t31+y*t32))))))))))))))))))))))))))))));
- r += (tf + p); break;
+ r += tf + p; break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+
y*(u8+y*(u9+y*u10))))))))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+
y*(v8+y*(v9+y*(v10+y*v11))))))))));
- r += (-y/2 + p1/p2);
+ r += p1/p2-y/2;
}
}
- else if(x<8) {
+ /* x < 8.0 */
+ else if(ix<0x4002) {
i = x;
y = x-i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+
case 3: z *= (y+2); /* FALLTHRU */
r += logl(z); break;
}
- } else if (x < 0x1p119L) {
+ /* 8.0 <= x < 2**(p+3) */
+ } else if (ix<0x3fff+116) {
t = logl(x);
z = one/x;
y = z*z;
y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+
y*(w17+y*w18)))))))))))))))));
r = (x-half)*(t-one)+w;
+ /* 2**(p+3) <= x <= inf */
} else
r = x*(logl(x)-1);
if(hx&0x8000) r = nadj - r;
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
} else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */
if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */
RETURNI(x+x);
- } else
+ } else
RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */
} else if (hx == 0x3fff && x == 1) {
RETURNI(0.0); /* acosh(1) = 0 */
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
if (ix < 0x3ffe) { /* |x| < 0.5, or misnormal */
t = x+x;
t = 0.5*log1pl(t+t*x/(one-x));
- } else
+ } else
t = 0.5*log1pl((x+x)/(one-x));
RETURNI((hx & 0x8000) == 0 ? t : -t);
}
-
/* @(#)e_lgamma_r.c 1.3 95/01/18 */
/*
* ====================================================
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
- *
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/* __ieee754_lgamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function
- * with user provide pointer for the sign of Gamma(x).
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
*
* Method:
* 1. Argument Reduction for 0 < x <= 8
- * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* for example,
* by
* 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
- * where
+ * where
* |w - f(z)| < 2**-58.74
- *
+ *
* 4. For negative x, since (G is gamma function)
* -x*G(-x)*G(x) = pi/sin(pi*x),
* we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
- * Hence, for x<0, signgam = sign(sin(pi*x)) and
+ * Hence, for x<0, signgam = sign(sin(pi*x)) and
* lgamma(x) = log(|Gamma(x)|)
* = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
- * Note: one should avoid compute pi*(-x) directly in the
+ * Note: one should avoid compute pi*(-x) directly in the
* computation of sin(pi*(-x)).
- *
+ *
* 5. Special Cases
* lgamma(2+s) ~ s*(1-Euler) for tiny s
* lgamma(1) = lgamma(2) = 0
* lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
* lgamma(inf) = inf
* lgamma(-inf) = inf (bug for bug compatible with C99!?)
- *
*/
#include <float.h>
switch (n) {
case 0: y = __kernel_sin(pi*y,zero,0); break;
- case 1:
+ case 1:
case 2: y = __kernel_cos(pi*(0.5-y),zero); break;
- case 3:
+ case 3:
case 4: y = __kernel_sin(pi*(one-y),zero,0); break;
case 5:
case 6: y = -__kernel_cos(pi*(y-1.5),zero); break;
double
__ieee754_lgamma_r(double x, int *signgamp)
{
- double t,y,z,nadj,p,p1,p2,p3,q,r,w;
+ double nadj,p,p1,p2,p3,q,r,t,w,y,z;
int32_t hx;
int i,ix,lx;
EXTRACT_WORDS(hx,lx,x);
- /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+ /* purge +-Inf and NaNs */
*signgamp = 1;
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return x*x;
- if((ix|lx)==0) {
- if(hx<0)
- *signgamp = -1;
- return one/vzero;
- }
- if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */
- if(hx<0) {
- *signgamp = -1;
- return -__ieee754_log(-x);
- } else return -__ieee754_log(x);
+
+ /* purge +-0 and tiny arguments */
+ *signgamp = 1-2*((uint32_t)hx>>31);
+ if(ix<0x3c700000) { /* |x|<2**-56, return -log(|x|) */
+ if((ix|lx)==0)
+ return one/vzero;
+ return -__ieee754_log(fabs(x));
}
+
+ /* purge negative integers and start evaluation for other x < 0 */
if(hx<0) {
+ *signgamp = 1;
if(ix>=0x43300000) /* |x|>=2**52, must be -integer */
return one/vzero;
t = sin_pi(x);
x = -x;
}
- /* purge off 1 and 2 */
+ /* purge 1 and 2 */
if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
/* for x < 2.0 */
else if(ix<0x40000000) {
p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
p = y*p1+p2;
- r += (p-y/2); break;
+ r += p-y/2; break;
case 1:
z = y*y;
w = z*y;
p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
p = z*p1-(tt-w*(p2+y*p3));
- r += (tf + p); break;
- case 2:
+ r += tf + p; break;
+ case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
- r += (-0.5*y + p1/p2);
+ r += p1/p2-y/2;
}
}
- else if(ix<0x40200000) { /* x < 8.0 */
- i = (int)x;
- y = x-(double)i;
+ /* x < 8.0 */
+ else if(ix<0x40200000) {
+ i = x;
+ y = x-i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
- r = half*y+p/q;
+ r = y/2+p/q;
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
switch(i) {
case 7: z *= (y+6); /* FALLTHRU */
case 3: z *= (y+2); /* FALLTHRU */
r += __ieee754_log(z); break;
}
- /* 8.0 <= x < 2**58 */
- } else if (ix < 0x43900000) {
+ /* 8.0 <= x < 2**56 */
+ } else if (ix < 0x43700000) {
t = __ieee754_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 */
+ } else
+ /* 2**56 <= x <= inf */
r = x*(__ieee754_log(x)-one);
if(hx<0) r = nadj - r;
return r;
#if (LDBL_MANT_DIG == 53)
__weak_reference(lgamma_r, lgammal_r);
#endif
-
float
__ieee754_lgammaf_r(float x, int *signgamp)
{
- float t,y,z,nadj,p,p1,p2,p3,q,r,w;
+ float nadj,p,p1,p2,p3,q,r,t,w,y,z;
int32_t hx;
int i,ix;
GET_FLOAT_WORD(hx,x);
- /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+ /* purge +-Inf and NaNs */
*signgamp = 1;
ix = hx&0x7fffffff;
if(ix>=0x7f800000) return x*x;
- if(ix==0) {
- if(hx<0)
- *signgamp = -1;
- return one/vzero;
- }
- if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */
- if(hx<0) {
- *signgamp = -1;
- return -__ieee754_logf(-x);
- } else return -__ieee754_logf(x);
+
+ /* purge +-0 and tiny arguments */
+ *signgamp = 1-2*((uint32_t)hx>>31);
+ if(ix<0x32000000) { /* |x|<2**-27, return -log(|x|) */
+ if(ix==0)
+ return one/vzero;
+ return -__ieee754_logf(fabsf(x));
}
+
+ /* purge negative integers and start evaluation for other x < 0 */
if(hx<0) {
- if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */
+ *signgamp = 1;
+ if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */
return one/vzero;
t = sin_pif(x);
if(t==zero) return one/vzero; /* -integer */
x = -x;
}
- /* purge off 1 and 2 */
+ /* purge 1 and 2 */
if (ix==0x3f800000||ix==0x40000000) r = 0;
/* for x < 2.0 */
else if(ix<0x40000000) {
p1 = a0+z*(a2+z*a4);
p2 = z*(a1+z*(a3+z*a5));
p = y*p1+p2;
- r += (p-y/2); break;
+ r += p-y/2; break;
case 1:
p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7)))));
- r += (tf + p); break;
+ r += tf + p; break;
case 2:
p1 = y*(u0+y*(u1+y*u2));
p2 = one+y*(v1+y*(v2+y*v3));
- r += (p1/p2-y/2);
+ r += p1/p2-y/2;
}
}
- else if(ix<0x41000000) { /* x < 8.0 */
+ /* x < 8.0 */
+ else if(ix<0x41000000) {
i = x;
y = x-i;
p = y*(s0+y*(s1+y*(s2+y*s3)));
case 3: z *= (y+2); /* FALLTHRU */
r += __ieee754_logf(z); break;
}
- /* 8.0 <= x < 2**24 */
- } else if (ix < 0x4b800000) {
+ /* 8.0 <= x < 2**27 */
+ } else if (ix < 0x4d000000) {
t = __ieee754_logf(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*w2);
r = (x-half)*(t-one)+w;
} else
- /* 2**24 <= x <= inf */
+ /* 2**27 <= x <= inf */
r = x*(__ieee754_logf(x)-one);
if(hx<0) r = nadj - r;
return r;
+++ /dev/null
-/*-
- * Copyright (c) 2004-2011 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.
- *
- * $FreeBSD$
- */
-
-#ifndef _FENV_H_
-#error "This file is meant to be included only by <fenv.h>."
-#endif
-
-/*
- * This file implements the functionality of <fenv.h> on platforms that
- * lack an FPU and use softfloat in libc for floating point. To use it,
- * you must write an <fenv.h> that provides the following:
- *
- * - a typedef for fenv_t, which may be an integer or struct type
- * - a typedef for fexcept_t (XXX This file assumes fexcept_t is a
- * simple integer type containing the exception mask.)
- * - definitions of FE_* constants for the five exceptions and four
- * rounding modes in IEEE 754, as described in fenv(3)
- * - a definition, and the corresponding external symbol, for FE_DFL_ENV
- * - a macro __set_env(env, flags, mask, rnd), which sets the given fenv_t
- * from the exception flags, mask, and rounding mode
- * - macros __env_flags(env), __env_mask(env), and __env_round(env), which
- * extract fields from an fenv_t
- * - a definition of __fenv_static
- *
- * If the architecture supports an optional FPU, it's recommended that you
- * define fenv_t and fexcept_t to match the hardware ABI. Otherwise, it
- * doesn't matter how you define them.
- */
-
-extern int __softfloat_float_exception_flags;
-extern int __softfloat_float_exception_mask;
-extern int __softfloat_float_rounding_mode;
-void __softfloat_float_raise(int);
-
-__fenv_static inline int
-feclearexcept(int __excepts)
-{
-
- __softfloat_float_exception_flags &= ~__excepts;
- return (0);
-}
-
-__fenv_static inline int
-fegetexceptflag(fexcept_t *__flagp, int __excepts)
-{
-
- *__flagp = __softfloat_float_exception_flags & __excepts;
- return (0);
-}
-
-__fenv_static inline int
-fesetexceptflag(const fexcept_t *__flagp, int __excepts)
-{
-
- __softfloat_float_exception_flags &= ~__excepts;
- __softfloat_float_exception_flags |= *__flagp & __excepts;
- return (0);
-}
-
-__fenv_static inline int
-feraiseexcept(int __excepts)
-{
-
- __softfloat_float_raise(__excepts);
- return (0);
-}
-
-__fenv_static inline int
-fetestexcept(int __excepts)
-{
-
- return (__softfloat_float_exception_flags & __excepts);
-}
-
-__fenv_static inline int
-fegetround(void)
-{
-
- return (__softfloat_float_rounding_mode);
-}
-
-__fenv_static inline int
-fesetround(int __round)
-{
-
- __softfloat_float_rounding_mode = __round;
- return (0);
-}
-
-__fenv_static inline int
-fegetenv(fenv_t *__envp)
-{
-
- __set_env(*__envp, __softfloat_float_exception_flags,
- __softfloat_float_exception_mask, __softfloat_float_rounding_mode);
- return (0);
-}
-
-__fenv_static inline int
-feholdexcept(fenv_t *__envp)
-{
- fenv_t __env;
-
- fegetenv(__envp);
- __softfloat_float_exception_flags = 0;
- __softfloat_float_exception_mask = 0;
- return (0);
-}
-
-__fenv_static inline int
-fesetenv(const fenv_t *__envp)
-{
-
- __softfloat_float_exception_flags = __env_flags(*__envp);
- __softfloat_float_exception_mask = __env_mask(*__envp);
- __softfloat_float_rounding_mode = __env_round(*__envp);
- return (0);
-}
-
-__fenv_static inline int
-feupdateenv(const fenv_t *__envp)
-{
- int __oflags = __softfloat_float_exception_flags;
-
- fesetenv(__envp);
- feraiseexcept(__oflags);
- return (0);
-}
-
-#if __BSD_VISIBLE
-
-/* We currently provide no external definitions of the functions below. */
-
-static inline int
-feenableexcept(int __mask)
-{
- int __omask = __softfloat_float_exception_mask;
-
- __softfloat_float_exception_mask |= __mask;
- return (__omask);
-}
-
-static inline int
-fedisableexcept(int __mask)
-{
- int __omask = __softfloat_float_exception_mask;
-
- __softfloat_float_exception_mask &= ~__mask;
- return (__omask);
-}
-
-static inline int
-fegetexcept(void)
-{
-
- return (__softfloat_float_exception_mask);
-}
-
-#endif /* __BSD_VISIBLE */
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
OSDN Git Service
