mirror of
https://github.com/adtools/clib2.git
synced 2025-12-08 14:59:05 +00:00
- Ported acosh(), acoshf(), asinh(), asinhf(), lgamma(), lgammaf(),
remainder() and remainderf() for C99. git-svn-id: file:///Users/olsen/Code/migration-svn-zu-git/logical-line-staging/clib2/trunk@14965 87f5fb63-7c3d-0410-a384-fd976d0f7a62
This commit is contained in:
Olaf Barthel
@ -1,5 +1,5 @@
|
||||
#
|
||||
# $Id: GNUmakefile.68k,v 1.57 2005-05-30 08:10:36 obarthel Exp $
|
||||
# $Id: GNUmakefile.68k,v 1.58 2005-05-30 08:47:25 obarthel Exp $
|
||||
#
|
||||
# :ts=8
|
||||
#
|
||||
@ -519,13 +519,14 @@ MATH_LIB = \
|
||||
complex_creal.o \
|
||||
complex_crealf.o \
|
||||
complex_creall.o \
|
||||
math_coshf.o \
|
||||
math_sinhf.o \
|
||||
math_tanhf.o \
|
||||
math_acos.o \
|
||||
math_acosf.o \
|
||||
math_acosh.o \
|
||||
math_acoshf.o \
|
||||
math_asin.o \
|
||||
math_asinf.o \
|
||||
math_asinh.o \
|
||||
math_asinhf.o \
|
||||
math_atan.o \
|
||||
math_atan2.o \
|
||||
math_atan2f.o \
|
||||
@ -539,6 +540,7 @@ MATH_LIB = \
|
||||
math_cos.o \
|
||||
math_cosf.o \
|
||||
math_cosh.o \
|
||||
math_coshf.o \
|
||||
math_erf.o \
|
||||
math_erfc.o \
|
||||
math_erfcf.o \
|
||||
@ -572,6 +574,8 @@ MATH_LIB = \
|
||||
math_kernel_tanf.o \
|
||||
math_ldexp.o \
|
||||
math_ldexpf.o \
|
||||
math_lgamma.o \
|
||||
math_lgammaf.o \
|
||||
math_log.o \
|
||||
math_log10.o \
|
||||
math_log10f.o \
|
||||
@ -588,6 +592,8 @@ MATH_LIB = \
|
||||
math_nextafterf.o \
|
||||
math_pow.o \
|
||||
math_powf.o \
|
||||
math_remainder.o \
|
||||
math_remainderf.o \
|
||||
math_rem_pio2f.o \
|
||||
math_rint.o \
|
||||
math_rintf.o \
|
||||
@ -597,11 +603,13 @@ MATH_LIB = \
|
||||
math_sin.o \
|
||||
math_sinf.o \
|
||||
math_sinh.o \
|
||||
math_sinhf.o \
|
||||
math_sqrt.o \
|
||||
math_sqrtf.o \
|
||||
math_tan.o \
|
||||
math_tanf.o \
|
||||
math_tanh.o \
|
||||
math_tanhf.o \
|
||||
stdio_asprintf.o \
|
||||
stdio_flush.o \
|
||||
stdio_flush_all_files.o \
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
#
|
||||
# $Id: GNUmakefile.os4,v 1.61 2005-05-30 08:10:37 obarthel Exp $
|
||||
# $Id: GNUmakefile.os4,v 1.62 2005-05-30 08:47:26 obarthel Exp $
|
||||
#
|
||||
# :ts=8
|
||||
#
|
||||
@ -519,13 +519,14 @@ MATH_LIB = \
|
||||
complex_creal.o \
|
||||
complex_crealf.o \
|
||||
complex_creall.o \
|
||||
math_coshf.o \
|
||||
math_sinhf.o \
|
||||
math_tanhf.o \
|
||||
math_acos.o \
|
||||
math_acosf.o \
|
||||
math_acosh.o \
|
||||
math_acoshf.o \
|
||||
math_asin.o \
|
||||
math_asinf.o \
|
||||
math_asinh.o \
|
||||
math_asinhf.o \
|
||||
math_atan.o \
|
||||
math_atan2.o \
|
||||
math_atan2f.o \
|
||||
@ -539,6 +540,7 @@ MATH_LIB = \
|
||||
math_cos.o \
|
||||
math_cosf.o \
|
||||
math_cosh.o \
|
||||
math_coshf.o \
|
||||
math_erf.o \
|
||||
math_erfc.o \
|
||||
math_erfcf.o \
|
||||
@ -578,6 +580,8 @@ MATH_LIB = \
|
||||
math_kernel_tanf.o \
|
||||
math_ldexp.o \
|
||||
math_ldexpf.o \
|
||||
math_lgamma.o \
|
||||
math_lgammaf.o \
|
||||
math_log.o \
|
||||
math_log10.o \
|
||||
math_log10f.o \
|
||||
@ -594,6 +598,8 @@ MATH_LIB = \
|
||||
math_nextafterf.o \
|
||||
math_pow.o \
|
||||
math_powf.o \
|
||||
math_remainder.o \
|
||||
math_remainderf.o \
|
||||
math_rem_pio2f.o \
|
||||
math_rint.o \
|
||||
math_rintf.o \
|
||||
@ -603,11 +609,13 @@ MATH_LIB = \
|
||||
math_sin.o \
|
||||
math_sinf.o \
|
||||
math_sinh.o \
|
||||
math_sinhf.o \
|
||||
math_sqrt.o \
|
||||
math_sqrtf.o \
|
||||
math_tan.o \
|
||||
math_tanf.o \
|
||||
math_tanh.o \
|
||||
math_tanhf.o \
|
||||
stdio_asprintf.o \
|
||||
stdio_flush.o \
|
||||
stdio_flush_all_files.o \
|
||||
|
||||
10
library/TODO
10
library/TODO
@ -1,16 +1,10 @@
|
||||
C99 math functions:
|
||||
|
||||
(functions generally missing, including their "float" counterparts)
|
||||
acosh
|
||||
acoshf
|
||||
asinh
|
||||
asinhf
|
||||
exp2
|
||||
exp2f
|
||||
fma
|
||||
fmaf
|
||||
lgamma
|
||||
lgammaf
|
||||
log2
|
||||
log2f
|
||||
lrint
|
||||
@ -19,10 +13,6 @@ C99 math functions:
|
||||
lroundf
|
||||
nearbyint
|
||||
nearbyintf
|
||||
nexttoward
|
||||
nexttowardf
|
||||
remainder
|
||||
remainderf
|
||||
remquo
|
||||
remquof
|
||||
round
|
||||
|
||||
@ -73,6 +73,9 @@
|
||||
- Ported cosf(), coshf(), sinf(), sinhf(), tanf(), tanhf()
|
||||
and hypotf() for C99.
|
||||
|
||||
- Ported acosh(), acoshf(), asinh(), asinhf(), lgamma(), lgammaf(),
|
||||
remainder() and remainderf() for C99.
|
||||
|
||||
|
||||
c.lib 1.192 (12.5.2005)
|
||||
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math.h,v 1.14 2005-05-30 08:10:41 obarthel Exp $
|
||||
* $Id: math.h,v 1.15 2005-05-30 08:47:41 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -242,53 +242,49 @@ extern float tanhf(float x);
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
extern float logbf(float x);
|
||||
extern double logb(double x);
|
||||
|
||||
extern float hypotf(float x, float y);
|
||||
extern double hypot(double x,double y);
|
||||
|
||||
extern float nanf(const char *tagp);
|
||||
extern double nan(const char *tagp);
|
||||
|
||||
extern float nextafterf(float x,float y);
|
||||
extern double nextafter(double x,double y);
|
||||
|
||||
extern float copysignf(float x, float y);
|
||||
extern double copysign(double x, double y);
|
||||
|
||||
extern float fdimf(float x,float y);
|
||||
extern double fdim(double x,double y);
|
||||
|
||||
extern float fminf(float x,float y);
|
||||
extern double fmin(double x,double y);
|
||||
|
||||
extern float fmaxf(float x,float y);
|
||||
extern double fmax(double x,double y);
|
||||
|
||||
extern float scalbnf (float x, int n);
|
||||
extern double scalbn (double x, int n);
|
||||
|
||||
extern float acoshf(float x);
|
||||
extern float asinhf(float x);
|
||||
extern float cbrtf(float x);
|
||||
extern double cbrt(double x);
|
||||
|
||||
extern float erff(float x);
|
||||
extern double erf(double x);
|
||||
|
||||
extern float copysignf(float x, float y);
|
||||
extern float erfcf(float x);
|
||||
extern double erfc(double x);
|
||||
|
||||
extern float erff(float x);
|
||||
extern float expm1f(float x);
|
||||
extern double expm1(double x);
|
||||
|
||||
extern int ilogbf(float x);
|
||||
extern int ilogb(double x);
|
||||
|
||||
extern float fdimf(float x,float y);
|
||||
extern float fmaxf(float x,float y);
|
||||
extern float fminf(float x,float y);
|
||||
extern float hypotf(float x, float y);
|
||||
extern float lgammaf(float x);
|
||||
extern float log1pf(float x);
|
||||
extern double log1p(double x);
|
||||
|
||||
extern float logbf(float x);
|
||||
extern float nanf(const char *tagp);
|
||||
extern float nextafterf(float x,float y);
|
||||
extern float remainderf(float x, float p);
|
||||
extern float rintf(float x);
|
||||
extern float scalbnf (float x, int n);
|
||||
extern int ilogbf(float x);
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
extern double acosh(double x);
|
||||
extern double asinh(double x);
|
||||
extern double cbrt(double x);
|
||||
extern double copysign(double x, double y);
|
||||
extern double erf(double x);
|
||||
extern double erfc(double x);
|
||||
extern double expm1(double x);
|
||||
extern double fdim(double x,double y);
|
||||
extern double fmax(double x,double y);
|
||||
extern double fmin(double x,double y);
|
||||
extern double hypot(double x,double y);
|
||||
extern double lgamma(double x);
|
||||
extern double log1p(double x);
|
||||
extern double logb(double x);
|
||||
extern double nan(const char *tagp);
|
||||
extern double nextafter(double x,double y);
|
||||
extern double remainder(double x, double p);
|
||||
extern double rint(double x);
|
||||
extern double scalbn (double x, int n);
|
||||
extern int ilogb(double x);
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_acosh.c,v 1.1 2005-05-29 11:19:00 obarthel Exp $
|
||||
* $Id: math_acosh.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -29,6 +29,15 @@
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* PowerPC math library based in part on work by Sun Microsystems
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef _MATH_HEADERS_H
|
||||
@ -41,11 +50,33 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const double
|
||||
one = 1.0,
|
||||
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
|
||||
|
||||
double
|
||||
acosh(double x)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
double t;
|
||||
LONG hx;
|
||||
ULONG lx;
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
if(hx<0x3ff00000) { /* x < 1 */
|
||||
return (x-x)/(x-x);
|
||||
} else if(hx >=0x41b00000) { /* x > 2**28 */
|
||||
if(hx >=0x7ff00000) { /* x is inf of NaN */
|
||||
return x+x;
|
||||
} else
|
||||
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 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));
|
||||
}
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_acoshf.c,v 1.1 2005-05-29 11:19:00 obarthel Exp $
|
||||
* $Id: math_acoshf.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -29,6 +29,18 @@
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* PowerPC math library based in part on work by Sun Microsystems
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
#ifndef _MATH_HEADERS_H
|
||||
@ -41,11 +53,32 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const float
|
||||
one = 1.0,
|
||||
ln2 = 6.9314718246e-01; /* 0x3f317218 */
|
||||
|
||||
float
|
||||
acoshf(float x)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
float t;
|
||||
LONG hx;
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
if(hx<0x3f800000) { /* x < 1 */
|
||||
return (x-x)/(x-x);
|
||||
} else if(hx >=0x4d800000) { /* x > 2**28 */
|
||||
if(hx >=0x7f800000) { /* x is inf of NaN */
|
||||
return x+x;
|
||||
} else
|
||||
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 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));
|
||||
}
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_asinh.c,v 1.1 2005-05-29 11:19:00 obarthel Exp $
|
||||
* $Id: math_asinh.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -29,6 +29,15 @@
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* PowerPC math library based in part on work by Sun Microsystems
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef _MATH_HEADERS_H
|
||||
@ -41,11 +50,32 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const double
|
||||
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
|
||||
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
|
||||
huge= 1.00000000000000000000e+300;
|
||||
|
||||
double
|
||||
asinh(double x)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
double t,w;
|
||||
LONG hx,ix;
|
||||
GET_HIGH_WORD(hx,x);
|
||||
ix = hx&0x7fffffff;
|
||||
if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */
|
||||
if(ix< 0x3e300000) { /* |x|<2**-28 */
|
||||
if(huge+x>one) return x; /* return x inexact except 0 */
|
||||
}
|
||||
if(ix>0x41b00000) { /* |x| > 2**28 */
|
||||
w = log(fabs(x))+ln2;
|
||||
} else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
|
||||
t = fabs(x);
|
||||
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+sqrt(one+t)));
|
||||
}
|
||||
if(hx>0) return w; else return -w;
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_asinhf.c,v 1.1 2005-05-29 11:19:00 obarthel Exp $
|
||||
* $Id: math_asinhf.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -29,6 +29,18 @@
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* PowerPC math library based in part on work by Sun Microsystems
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
#ifndef _MATH_HEADERS_H
|
||||
@ -41,11 +53,32 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const float
|
||||
one = 1.0000000000e+00, /* 0x3F800000 */
|
||||
ln2 = 6.9314718246e-01, /* 0x3f317218 */
|
||||
huge= 1.0000000000e+30;
|
||||
|
||||
float
|
||||
asinhf(float x)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
float t,w;
|
||||
LONG hx,ix;
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
ix = hx&0x7fffffff;
|
||||
if(ix>=0x7f800000) return x+x; /* x is inf or NaN */
|
||||
if(ix< 0x31800000) { /* |x|<2**-28 */
|
||||
if(huge+x>one) return x; /* return x inexact except 0 */
|
||||
}
|
||||
if(ix>0x4d800000) { /* |x| > 2**28 */
|
||||
w = logf(fabsf(x))+ln2;
|
||||
} else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
|
||||
t = fabsf(x);
|
||||
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+sqrtf(one+t)));
|
||||
}
|
||||
if(hx>0) return w; else return -w;
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_lgamma.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $
|
||||
* $Id: math_lgamma.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -41,11 +41,211 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const double
|
||||
two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
|
||||
half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
|
||||
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
|
||||
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
|
||||
a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
|
||||
a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
|
||||
a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
|
||||
a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
|
||||
a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
|
||||
a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
|
||||
a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
|
||||
a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
|
||||
a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
|
||||
a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
|
||||
a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
|
||||
a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
|
||||
tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
|
||||
tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
|
||||
/* tt = -(tail of tf) */
|
||||
tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
|
||||
t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
|
||||
t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
|
||||
t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
|
||||
t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
|
||||
t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
|
||||
t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
|
||||
t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
|
||||
t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
|
||||
t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
|
||||
t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
|
||||
t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
|
||||
t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
|
||||
t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
|
||||
t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
|
||||
t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
|
||||
u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
|
||||
u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
|
||||
u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
|
||||
u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
|
||||
u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
|
||||
u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
|
||||
v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
|
||||
v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
|
||||
v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
|
||||
v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
|
||||
v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
|
||||
s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
|
||||
s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
|
||||
s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
|
||||
s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
|
||||
s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
|
||||
s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
|
||||
s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
|
||||
r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
|
||||
r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
|
||||
r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
|
||||
r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
|
||||
r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
|
||||
r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
|
||||
w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
|
||||
w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
|
||||
w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
|
||||
w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
|
||||
w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
|
||||
w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
|
||||
w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
|
||||
|
||||
static const double zero= 0.00000000000000000000e+00;
|
||||
|
||||
static double
|
||||
sin_pi(double x)
|
||||
{
|
||||
double y,z;
|
||||
LONG n,ix;
|
||||
|
||||
GET_HIGH_WORD(ix,x);
|
||||
ix &= 0x7fffffff;
|
||||
|
||||
if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
|
||||
y = -x; /* x is assume negative */
|
||||
|
||||
/*
|
||||
* argument reduction, make sure inexact flag not raised if input
|
||||
* is an integer
|
||||
*/
|
||||
z = floor(y);
|
||||
if(z!=y) { /* inexact anyway */
|
||||
y *= 0.5;
|
||||
y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */
|
||||
n = (LONG) (y*4.0);
|
||||
} else {
|
||||
if(ix>=0x43400000) {
|
||||
y = zero; n = 0; /* y must be even */
|
||||
} else {
|
||||
if(ix<0x43300000) z = y+two52; /* exact */
|
||||
GET_LOW_WORD(n,z);
|
||||
n &= 1;
|
||||
y = n;
|
||||
n<<= 2;
|
||||
}
|
||||
}
|
||||
switch (n) {
|
||||
case 0: y = __kernel_sin(pi*y,zero,0); break;
|
||||
case 1:
|
||||
case 2: y = __kernel_cos(pi*(0.5-y),zero); break;
|
||||
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;
|
||||
default: y = __kernel_sin(pi*(y-2.0),zero,0); break;
|
||||
}
|
||||
return -y;
|
||||
}
|
||||
|
||||
double
|
||||
lgamma(double x)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
double t,y,z,nadj=0.0,p,p1,p2,p3,q,r,w;
|
||||
LONG i,hx,lx,ix;
|
||||
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
|
||||
/* purge off +-inf, NaN, +-0, and negative arguments */
|
||||
ix = hx&0x7fffffff;
|
||||
if(ix>=0x7ff00000) return x*x;
|
||||
if((ix|lx)==0) return one/zero;
|
||||
if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */
|
||||
if(hx<0) {
|
||||
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 = log(pi/fabs(t*x));
|
||||
x = -x;
|
||||
}
|
||||
|
||||
/* purge off 1 and 2 */
|
||||
if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
|
||||
/* for x < 2.0 */
|
||||
else if(ix<0x40000000) {
|
||||
if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-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;}
|
||||
} else {
|
||||
r = zero;
|
||||
if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */
|
||||
else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */
|
||||
else {y=x-one;i=2;}
|
||||
}
|
||||
switch(i) {
|
||||
case 0:
|
||||
z = y*y;
|
||||
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-0.5*y); break;
|
||||
case 1:
|
||||
z = y*y;
|
||||
w = z*y;
|
||||
p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
|
||||
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:
|
||||
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);
|
||||
}
|
||||
}
|
||||
else if(ix<0x40200000) { /* x < 8.0 */
|
||||
i = (LONG)x;
|
||||
t = zero;
|
||||
y = x-(double)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;
|
||||
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
|
||||
switch(i) {
|
||||
case 7: z *= (y+6.0); /* FALLTHRU */
|
||||
case 6: z *= (y+5.0); /* FALLTHRU */
|
||||
case 5: z *= (y+4.0); /* FALLTHRU */
|
||||
case 4: z *= (y+3.0); /* FALLTHRU */
|
||||
case 3: z *= (y+2.0); /* FALLTHRU */
|
||||
r += log(z); break;
|
||||
}
|
||||
/* 8.0 <= x < 2**58 */
|
||||
} else if (ix < 0x43900000) {
|
||||
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*(log(x)-one);
|
||||
if(hx<0) r = nadj - r;
|
||||
return r;
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_lgammaf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $
|
||||
* $Id: math_lgammaf.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -29,6 +29,18 @@
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* PowerPC math library based in part on work by Sun Microsystems
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
#ifndef _MATH_HEADERS_H
|
||||
@ -41,11 +53,212 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const float
|
||||
two23= 8.3886080000e+06, /* 0x4b000000 */
|
||||
half= 5.0000000000e-01, /* 0x3f000000 */
|
||||
one = 1.0000000000e+00, /* 0x3f800000 */
|
||||
pi = 3.1415927410e+00, /* 0x40490fdb */
|
||||
a0 = 7.7215664089e-02, /* 0x3d9e233f */
|
||||
a1 = 3.2246702909e-01, /* 0x3ea51a66 */
|
||||
a2 = 6.7352302372e-02, /* 0x3d89f001 */
|
||||
a3 = 2.0580807701e-02, /* 0x3ca89915 */
|
||||
a4 = 7.3855509982e-03, /* 0x3bf2027e */
|
||||
a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */
|
||||
a6 = 1.1927076848e-03, /* 0x3a9c54a1 */
|
||||
a7 = 5.1006977446e-04, /* 0x3a05b634 */
|
||||
a8 = 2.2086278477e-04, /* 0x39679767 */
|
||||
a9 = 1.0801156895e-04, /* 0x38e28445 */
|
||||
a10 = 2.5214456400e-05, /* 0x37d383a2 */
|
||||
a11 = 4.4864096708e-05, /* 0x383c2c75 */
|
||||
tc = 1.4616321325e+00, /* 0x3fbb16c3 */
|
||||
tf = -1.2148628384e-01, /* 0xbdf8cdcd */
|
||||
/* tt = -(tail of tf) */
|
||||
tt = 6.6971006518e-09, /* 0x31e61c52 */
|
||||
t0 = 4.8383611441e-01, /* 0x3ef7b95e */
|
||||
t1 = -1.4758771658e-01, /* 0xbe17213c */
|
||||
t2 = 6.4624942839e-02, /* 0x3d845a15 */
|
||||
t3 = -3.2788541168e-02, /* 0xbd064d47 */
|
||||
t4 = 1.7970675603e-02, /* 0x3c93373d */
|
||||
t5 = -1.0314224288e-02, /* 0xbc28fcfe */
|
||||
t6 = 6.1005386524e-03, /* 0x3bc7e707 */
|
||||
t7 = -3.6845202558e-03, /* 0xbb7177fe */
|
||||
t8 = 2.2596477065e-03, /* 0x3b141699 */
|
||||
t9 = -1.4034647029e-03, /* 0xbab7f476 */
|
||||
t10 = 8.8108185446e-04, /* 0x3a66f867 */
|
||||
t11 = -5.3859531181e-04, /* 0xba0d3085 */
|
||||
t12 = 3.1563205994e-04, /* 0x39a57b6b */
|
||||
t13 = -3.1275415677e-04, /* 0xb9a3f927 */
|
||||
t14 = 3.3552918467e-04, /* 0x39afe9f7 */
|
||||
u0 = -7.7215664089e-02, /* 0xbd9e233f */
|
||||
u1 = 6.3282704353e-01, /* 0x3f2200f4 */
|
||||
u2 = 1.4549225569e+00, /* 0x3fba3ae7 */
|
||||
u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */
|
||||
u4 = 2.2896373272e-01, /* 0x3e6a7578 */
|
||||
u5 = 1.3381091878e-02, /* 0x3c5b3c5e */
|
||||
v1 = 2.4559779167e+00, /* 0x401d2ebe */
|
||||
v2 = 2.1284897327e+00, /* 0x4008392d */
|
||||
v3 = 7.6928514242e-01, /* 0x3f44efdf */
|
||||
v4 = 1.0422264785e-01, /* 0x3dd572af */
|
||||
v5 = 3.2170924824e-03, /* 0x3b52d5db */
|
||||
s0 = -7.7215664089e-02, /* 0xbd9e233f */
|
||||
s1 = 2.1498242021e-01, /* 0x3e5c245a */
|
||||
s2 = 3.2577878237e-01, /* 0x3ea6cc7a */
|
||||
s3 = 1.4635047317e-01, /* 0x3e15dce6 */
|
||||
s4 = 2.6642270386e-02, /* 0x3cda40e4 */
|
||||
s5 = 1.8402845599e-03, /* 0x3af135b4 */
|
||||
s6 = 3.1947532989e-05, /* 0x3805ff67 */
|
||||
r1 = 1.3920053244e+00, /* 0x3fb22d3b */
|
||||
r2 = 7.2193557024e-01, /* 0x3f38d0c5 */
|
||||
r3 = 1.7193385959e-01, /* 0x3e300f6e */
|
||||
r4 = 1.8645919859e-02, /* 0x3c98bf54 */
|
||||
r5 = 7.7794247773e-04, /* 0x3a4beed6 */
|
||||
r6 = 7.3266842264e-06, /* 0x36f5d7bd */
|
||||
w0 = 4.1893854737e-01, /* 0x3ed67f1d */
|
||||
w1 = 8.3333335817e-02, /* 0x3daaaaab */
|
||||
w2 = -2.7777778450e-03, /* 0xbb360b61 */
|
||||
w3 = 7.9365057172e-04, /* 0x3a500cfd */
|
||||
w4 = -5.9518753551e-04, /* 0xba1c065c */
|
||||
w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
|
||||
w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
|
||||
|
||||
static const float zero= 0.0000000000e+00;
|
||||
|
||||
static float
|
||||
sin_pif(float x)
|
||||
{
|
||||
float y,z;
|
||||
LONG n,ix;
|
||||
|
||||
GET_FLOAT_WORD(ix,x);
|
||||
ix &= 0x7fffffff;
|
||||
|
||||
if(ix<0x3e800000) return __kernel_sinf(pi*x,zero,0);
|
||||
y = -x; /* x is assume negative */
|
||||
|
||||
/*
|
||||
* argument reduction, make sure inexact flag not raised if input
|
||||
* is an integer
|
||||
*/
|
||||
z = floorf(y);
|
||||
if(z!=y) { /* inexact anyway */
|
||||
y *= (float)0.5;
|
||||
y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */
|
||||
n = (LONG) (y*(float)4.0);
|
||||
} else {
|
||||
if(ix>=0x4b800000) {
|
||||
y = zero; n = 0; /* y must be even */
|
||||
} else {
|
||||
if(ix<0x4b000000) z = y+two23; /* exact */
|
||||
GET_FLOAT_WORD(n,z);
|
||||
n &= 1;
|
||||
y = n;
|
||||
n<<= 2;
|
||||
}
|
||||
}
|
||||
switch (n) {
|
||||
case 0: y = __kernel_sinf(pi*y,zero,0); break;
|
||||
case 1:
|
||||
case 2: y = __kernel_cosf(pi*((float)0.5-y),zero); break;
|
||||
case 3:
|
||||
case 4: y = __kernel_sinf(pi*(one-y),zero,0); break;
|
||||
case 5:
|
||||
case 6: y = -__kernel_cosf(pi*(y-(float)1.5),zero); break;
|
||||
default: y = __kernel_sinf(pi*(y-(float)2.0),zero,0); break;
|
||||
}
|
||||
return -y;
|
||||
}
|
||||
|
||||
|
||||
float
|
||||
lgammaf(float x)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
float t,y,z,nadj=0.0,p,p1,p2,p3,q,r,w;
|
||||
LONG i,hx,ix;
|
||||
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
|
||||
/* purge off +-inf, NaN, +-0, and negative arguments */
|
||||
ix = hx&0x7fffffff;
|
||||
if(ix>=0x7f800000) return x*x;
|
||||
if(ix==0) return one/zero;
|
||||
if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */
|
||||
if(hx<0) {
|
||||
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 = logf(pi/fabsf(t*x));
|
||||
x = -x;
|
||||
}
|
||||
|
||||
/* purge off 1 and 2 */
|
||||
if (ix==0x3f800000||ix==0x40000000) r = 0;
|
||||
/* for x < 2.0 */
|
||||
else if(ix<0x40000000) {
|
||||
if(ix<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(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;}
|
||||
} else {
|
||||
r = zero;
|
||||
if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */
|
||||
else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */
|
||||
else {y=x-one;i=2;}
|
||||
}
|
||||
switch(i) {
|
||||
case 0:
|
||||
z = y*y;
|
||||
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-(float)0.5*y); break;
|
||||
case 1:
|
||||
z = y*y;
|
||||
w = z*y;
|
||||
p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
|
||||
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:
|
||||
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 += (-(float)0.5*y + p1/p2);
|
||||
}
|
||||
}
|
||||
else if(ix<0x41000000) { /* x < 8.0 */
|
||||
i = (LONG)x;
|
||||
t = zero;
|
||||
y = x-(float)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;
|
||||
z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
|
||||
switch(i) {
|
||||
case 7: z *= (y+(float)6.0); /* FALLTHRU */
|
||||
case 6: z *= (y+(float)5.0); /* FALLTHRU */
|
||||
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 += logf(z); break;
|
||||
}
|
||||
/* 8.0 <= x < 2**58 */
|
||||
} else if (ix < 0x5c800000) {
|
||||
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*(logf(x)-one);
|
||||
if(hx<0) r = nadj - r;
|
||||
return r;
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_remainder.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $
|
||||
* $Id: math_remainder.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -29,6 +29,15 @@
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* PowerPC math library based in part on work by Sun Microsystems
|
||||
* 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.
|
||||
*/
|
||||
|
||||
#ifndef _MATH_HEADERS_H
|
||||
@ -41,11 +50,48 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const double zero = 0.0;
|
||||
|
||||
double
|
||||
remainder(double x,double y)
|
||||
remainder(double x, double p)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
LONG hx,hp;
|
||||
ULONG sx,lx,lp;
|
||||
double p_half;
|
||||
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
EXTRACT_WORDS(hp,lp,p);
|
||||
sx = hx&0x80000000U;
|
||||
hp &= 0x7fffffff;
|
||||
hx &= 0x7fffffff;
|
||||
|
||||
/* purge off exception values */
|
||||
if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
|
||||
if((hx>=0x7ff00000)|| /* x not finite */
|
||||
((hp>=0x7ff00000)&& /* p is NaN */
|
||||
(((hp-0x7ff00000)|lp)!=0)))
|
||||
return (x*p)/(x*p);
|
||||
|
||||
|
||||
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);
|
||||
if (hp<0x00200000) {
|
||||
if(x+x>p) {
|
||||
x-=p;
|
||||
if(x+x>=p) x -= p;
|
||||
}
|
||||
} else {
|
||||
p_half = 0.5*p;
|
||||
if(x>p_half) {
|
||||
x-=p;
|
||||
if(x>=p_half) x -= p;
|
||||
}
|
||||
}
|
||||
GET_HIGH_WORD(hx,x);
|
||||
SET_HIGH_WORD(x,hx^sx);
|
||||
return x;
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* $Id: math_remainderf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $
|
||||
* $Id: math_remainderf.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
|
||||
*
|
||||
* :ts=4
|
||||
*
|
||||
@ -29,6 +29,18 @@
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* PowerPC math library based in part on work by Sun Microsystems
|
||||
* 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.
|
||||
*
|
||||
*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
#ifndef _MATH_HEADERS_H
|
||||
@ -41,11 +53,47 @@
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
static const float zero = 0.0;
|
||||
|
||||
float
|
||||
remainderf(float x,float y)
|
||||
remainderf(float x, float p)
|
||||
{
|
||||
/* ZZZ unimplemented */
|
||||
return(0);
|
||||
LONG hx,hp;
|
||||
ULONG sx;
|
||||
float p_half;
|
||||
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
GET_FLOAT_WORD(hp,p);
|
||||
sx = hx&0x80000000U;
|
||||
hp &= 0x7fffffff;
|
||||
hx &= 0x7fffffff;
|
||||
|
||||
/* purge off exception values */
|
||||
if(hp==0) return (x*p)/(x*p); /* p = 0 */
|
||||
if((hx>=0x7f800000)|| /* x not finite */
|
||||
((hp>0x7f800000))) /* p is NaN */
|
||||
return (x*p)/(x*p);
|
||||
|
||||
|
||||
if (hp<=0x7effffff) x = fmodf(x,p+p); /* now x < 2p */
|
||||
if ((hx-hp)==0) return zero*x;
|
||||
x = fabsf(x);
|
||||
p = fabsf(p);
|
||||
if (hp<0x01000000) {
|
||||
if(x+x>p) {
|
||||
x-=p;
|
||||
if(x+x>=p) x -= p;
|
||||
}
|
||||
} else {
|
||||
p_half = (float)0.5*p;
|
||||
if(x>p_half) {
|
||||
x-=p;
|
||||
if(x>=p_half) x -= p;
|
||||
}
|
||||
}
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
SET_FLOAT_WORD(x,hx^sx);
|
||||
return x;
|
||||
}
|
||||
|
||||
/****************************************************************************/
|
||||
|
||||
Reference in New Issue
Block a user