Skip to content

T
Third Party Musl

OpenHarmony / Third Party Musl
1 年多 前同步成功

通知 37
Star 125
Fork 0
T
Third Party Musl
提交 969ddbc4 编写于 作者: R Rich Felker
浏览文件
操作

Merge remote-tracking branch 'nsz/math'

上级 9cb58993 a8f73bb1
隐藏空白更改
内联 并排
Showing with 803 addition and 907 deletion
include/complex.h
浏览文件 @ 969ddbc4
......@@ -115,11 +115,9 @@ long double creall(long double complex);
#define __CMPLX(x, y, t) \
((union { _Complex t __z; t __xy[2]; }){.__xy = {(x),(y)}}.__z)
#if __STDC_VERSION__ >= 201112L
#define CMPLX(x, y) __CMPLX(x, y, double)
#define CMPLXF(x, y) __CMPLX(x, y, float)
#define CMPLXL(x, y) __CMPLX(x, y, long double)
#endif
#ifdef __cplusplus
}
......
include/math.h
浏览文件 @ 969ddbc4
......@@ -399,14 +399,6 @@ float ynf(int, float);
#ifdef _GNU_SOURCE
long double lgammal_r(long double, int*);
long double j0l(long double);
long double j1l(long double);
long double jnl(int, long double);
long double y0l(long double);
long double y1l(long double);
long double ynl(int, long double);
void sincos(double, double*, double*);
void sincosf(float, float*, float*);
void sincosl(long double, long double*, long double*);
......
include/tgmath.h
浏览文件 @ 969ddbc4
......@@ -59,10 +59,12 @@ sizeof(double) == sizeof(long double)
/* function selection */
#define __tg_real(fun, x) (__RETCAST(x)( \
#define __tg_real_nocast(fun, x) ( \
__FLT(x) ? fun ## f (x) : \
__LDBL(x) ? fun ## l (x) : \
fun(x) ))
fun(x) )
#define __tg_real(fun, x) (__RETCAST(x)__tg_real_nocast(fun, x))
#define __tg_real_2_1(fun, x, y) (__RETCAST(x)( \
__FLT(x) ? fun ## f (x, y) : \
......@@ -217,18 +219,18 @@ sizeof(double) == sizeof(long double)
#define fmod(x,y) __tg_real_2(fmod, (x), (y))
#define frexp(x,y) __tg_real_2_1(frexp, (x), (y))
#define hypot(x,y) __tg_real_2(hypot, (x), (y))
#define ilogb(x) __tg_real(ilogb, (x))
#define ilogb(x) __tg_real_nocast(ilogb, (x))
#define ldexp(x,y) __tg_real_2_1(ldexp, (x), (y))
#define lgamma(x) __tg_real(lgamma, (x))
#define llrint(x) __tg_real(llrint, (x))
#define llround(x) __tg_real(llround, (x))
#define llrint(x) __tg_real_nocast(llrint, (x))
#define llround(x) __tg_real_nocast(llround, (x))
#define log(x) __tg_real_complex(log, (x))
#define log10(x) __tg_real(log10, (x))
#define log1p(x) __tg_real(log1p, (x))
#define log2(x) __tg_real(log2, (x))
#define logb(x) __tg_real(logb, (x))
#define lrint(x) __tg_real(lrint, (x))
#define lround(x) __tg_real(lround, (x))
#define lrint(x) __tg_real_nocast(lrint, (x))
#define lround(x) __tg_real_nocast(lround, (x))
#define nearbyint(x) __tg_real(nearbyint, (x))
#define nextafter(x,y) __tg_real_2(nextafter, (x), (y))
#define nexttoward(x,y) __tg_real_2(nexttoward, (x), (y))
......
src/internal/libm.h
浏览文件 @ 969ddbc4
......@@ -168,12 +168,4 @@ long double __p1evll(long double, const long double *, int);
#define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval))
#endif
/* complex */
#ifndef CMPLX
#define CMPLX(x, y) __CMPLX(x, y, double)
#define CMPLXF(x, y) __CMPLX(x, y, float)
#define CMPLXL(x, y) __CMPLX(x, y, long double)
#endif
#endif
src/math/__invtrigl.c
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/ld80/invtrig.c */
/*-
* Copyright (c) 2008 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 "__invtrigl.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* coefficients used in asinl() and acosl() */
const long double
static const long double
pS0 = 1.66666666666666666631e-01L,
pS1 = -4.16313987993683104320e-01L,
pS2 = 3.69068046323246813704e-01L,
......@@ -44,9 +16,16 @@ qS3 = -1.68285799854822427013e+00L,
qS4 = 3.90699412641738801874e-01L,
qS5 = -3.14365703596053263322e-02L;
const long double pi_hi = 3.1415926535897932384626433832795L;
const long double pi_lo = -5.01655761266833202345e-20L;
const long double pio2_hi = 1.57079632679489661926L;
const long double pio2_lo = -2.50827880633416601173e-20L;
/* used in asinl() and acosl() */
/* R(x^2) is a rational approximation of (asin(x)-x)/x^3 with Remez algorithm */
long double __invtrigl_R(long double z)
{
long double p, q;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*(pS5+z*pS6))))));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*(qS4+z*qS5))));
return p/q;
}
#endif
src/math/__invtrigl.h
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/ld80/invtrig.h */
/*-
* Copyright (c) 2008 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 "libm.h"
#include <float.h>
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
#define BIAS (LDBL_MAX_EXP - 1)
#define MANH_SIZE LDBL_MANH_SIZE
/* Constants shared by the long double inverse trig functions. */
#define pS0 __pS0
#define pS1 __pS1
#define pS2 __pS2
#define pS3 __pS3
#define pS4 __pS4
#define pS5 __pS5
#define pS6 __pS6
#define qS1 __qS1
#define qS2 __qS2
#define qS3 __qS3
#define qS4 __qS4
#define qS5 __qS5
#define pi_hi __pi_hi
#define pi_lo __pi_lo
/* shared by acosl, asinl and atan2l */
#define pio2_hi __pio2_hi
#define pio2_lo __pio2_lo
extern const long double pio2_hi, pio2_lo;
extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6;
extern const long double qS1, qS2, qS3, qS4, qS5;
extern const long double pi_hi, pi_lo, pio2_hi, pio2_lo;
static long double P(long double x)
{
return x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 +
x * (pS4 + x * (pS5 + x * pS6))))));
}
static long double Q(long double x)
{
return 1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5))));
}
long double __invtrigl_R(long double z);
#endif
src/math/__rem_pio2f.c
浏览文件 @ 969ddbc4
......@@ -24,7 +24,7 @@
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1: first 25 bits of pi/2
* pio2_1t: pi/2 - pio2_1
*/
static const double
......@@ -41,7 +41,7 @@ int __rem_pio2f(float x, double *y)
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
/* 33+53 bit pi is good enough for medium size */
/* 25+53 bit pi is good enough for medium size */
if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */
/* Use a specialized rint() to get fn. Assume round-to-nearest. */
STRICT_ASSIGN(double, fn, x*invpio2 + 0x1.8p52);
......
src/math/acos.c
浏览文件 @ 969ddbc4
......@@ -36,12 +36,8 @@
#include "libm.h"
static const double
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
// FIXME
static const volatile double
pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
static const double
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
......@@ -53,49 +49,53 @@ qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
static double R(double z)
{
double p, q;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
return p/q;
}
double acos(double x)
{
double z,p,q,r,w,s,c,df;
int32_t hx,ix;
double z,w,s,c,df;
uint32_t hx,ix;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x3ff00000) { /* |x| >= 1 */
/* |x| >= 1 or nan */
if (ix >= 0x3ff00000) {
uint32_t lx;
GET_LOW_WORD(lx,x);
if ((ix-0x3ff00000 | lx) == 0) { /* |x|==1 */
if (hx > 0) return 0.0; /* acos(1) = 0 */
return pi + 2.0*pio2_lo; /* acos(-1)= pi */
if ((ix-0x3ff00000 | lx) == 0) {
/* acos(1)=0, acos(-1)=pi */
if (hx >> 31)
return 2*pio2_hi + 0x1p-1000;
return 0;
}
return (x-x)/(x-x); /* acos(|x|>1) is NaN */
return 0/(x-x);
}
if (ix < 0x3fe00000) { /* |x| < 0.5 */
/* |x| < 0.5 */
if (ix < 0x3fe00000) {
if (ix <= 0x3c600000) /* |x| < 2**-57 */
return pio2_hi + pio2_lo;
z = x*x;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */
return pio2_hi + 0x1p-1000;
return pio2_hi - (x - (pio2_lo-x*R(x*x)));
}
/* x < -0.5 */
if (hx >> 31) {
z = (1.0+x)*0.5;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
s = sqrt(z);
r = p/q;
w = r*s-pio2_lo;
return pi - 2.0*(s+w);
} else { /* x > 0.5 */
z = (1.0-x)*0.5;
s = sqrt(z);
df = s;
SET_LOW_WORD(df,0);
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
w = r*s+c;
return 2.0*(df+w);
w = R(z)*s-pio2_lo;
return 2*(pio2_hi - (s+w));
}
/* x > 0.5 */
z = (1.0-x)*0.5;
s = sqrt(z);
df = s;
SET_LOW_WORD(df,0);
c = (z-df*df)/(s+df);
w = R(z)*s+c;
return 2*(df+w);
}
src/math/acosf.c
浏览文件 @ 969ddbc4
......@@ -16,59 +16,56 @@
#include "libm.h"
static const float
pi = 3.1415925026e+00, /* 0x40490fda */
pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */
static const volatile float
pio2_lo = 7.5497894159e-08; /* 0x33a22168 */
static const float
pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
pio2_lo = 7.5497894159e-08, /* 0x33a22168 */
pS0 = 1.6666586697e-01,
pS1 = -4.2743422091e-02,
pS2 = -8.6563630030e-03,
qS1 = -7.0662963390e-01;
static float R(float z)
{
float p, q;
p = z*(pS0+z*(pS1+z*pS2));
q = 1.0f+z*qS1;
return p/q;
}
float acosf(float x)
{
float z,p,q,r,w,s,c,df;
int32_t hx,ix;
float z,w,s,c,df;
uint32_t hx,ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x3f800000) { /* |x| >= 1 */
if (ix == 0x3f800000) { /* |x| == 1 */
if (hx > 0) return 0.0f; /* acos(1) = 0 */
return pi + 2.0f*pio2_lo; /* acos(-1)= pi */
/* |x| >= 1 or nan */
if (ix >= 0x3f800000) {
if (ix == 0x3f800000) {
if (hx >> 31)
return 2*pio2_hi + 0x1p-120f;
return 0;
}
return (x-x)/(x-x); /* acos(|x|>1) is NaN */
return 0/(x-x);
}
if (ix < 0x3f000000) { /* |x| < 0.5 */
/* |x| < 0.5 */
if (ix < 0x3f000000) {
if (ix <= 0x32800000) /* |x| < 2**-26 */
return pio2_hi + pio2_lo;
z = x*x;
p = z*(pS0+z*(pS1+z*pS2));
q = 1.0f+z*qS1;
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */
z = (1.0f+x)*0.5f;
p = z*(pS0+z*(pS1+z*pS2));
q = 1.0f+z*qS1;
s = sqrtf(z);
r = p/q;
w = r*s-pio2_lo;
return pi - 2.0f*(s+w);
} else { /* x > 0.5 */
int32_t idf;
z = (1.0f-x)*0.5f;
return pio2_hi + 0x1p-120f;
return pio2_hi - (x - (pio2_lo-x*R(x*x)));
}
/* x < -0.5 */
if (hx >> 31) {
z = (1+x)*0.5f;
s = sqrtf(z);
df = s;
GET_FLOAT_WORD(idf,df);
SET_FLOAT_WORD(df,idf&0xfffff000);
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*pS2));
q = 1.0f+z*qS1;
r = p/q;
w = r*s+c;
return 2.0f*(df+w);
w = R(z)*s-pio2_lo;
return 2*(pio2_hi - (s+w));
}
/* x > 0.5 */
z = (1-x)*0.5f;
s = sqrtf(z);
GET_FLOAT_WORD(hx,s);
SET_FLOAT_WORD(df,hx&0xfffff000);
c = (z-df*df)/(s+df);
w = R(z)*s+c;
return 2*(df+w);
}
src/math/acosh.c
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/e_acosh.c */
/*
* ====================================================
* 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.
* ====================================================
*
*/
/* acosh(x)
* Method :
* Based on
* acosh(x) = log [ x + sqrt(x*x-1) ]
* we have
* acosh(x) := log(x)+ln2, if x is large; else
* acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
* acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
*
* Special cases:
* acosh(x) is NaN with signal if x<1.
* acosh(NaN) is NaN without signal.
*/
#include "libm.h"
static const double
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
/* acosh(x) = log(x + sqrt(x*x-1)) */
double acosh(double x)
{
double t;
int32_t hx;
uint32_t lx;
union {double f; uint64_t i;} u = {.f = x};
unsigned e = u.i >> 52 & 0x7ff;
/* x < 1 domain error is handled in the called functions */
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;
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 - 1.0/(x+sqrt(t-1.0)));
} else { /* 1 < x < 2 */
t = x-1.0;
return log1p(t + sqrt(2.0*t+t*t));
}
if (e < 0x3ff + 1)
/* |x| < 2, up to 2ulp error in [1,1.125] */
return log1p(x-1 + sqrt((x-1)*(x-1)+2*(x-1)));
if (e < 0x3ff + 26)
/* |x| < 0x1p26 */
return log(2*x - 1/(x+sqrt(x*x-1)));
/* |x| >= 0x1p26 or nan */
return log(x) + 0.693147180559945309417232121458176568;
}
src/math/acoshf.c
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/e_acoshf.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 "libm.h"
static const float
ln2 = 6.9314718246e-01; /* 0x3f317218 */
/* acosh(x) = log(x + sqrt(x*x-1)) */
float acoshf(float x)
{
float t;
int32_t hx;
union {float f; int32_t i;} u = {.f = x};
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;
return logf(x) + ln2; /* acosh(huge)=log(2x) */
} else if (hx == 0x3f800000) {
return 0.0f; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x;
return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f)));
} else { /* 1 < x < 2 */
t = x-1.0f;
return log1pf(t + sqrtf(2.0f*t+t*t));
}
if (u.i < 0x3f800000+(1<<23))
/* x < 2, invalid if x < 1 or nan */
/* up to 2ulp error in [1,1.125] */
return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1)));
if (u.i < 0x3f800000+(12<<23))
/* x < 0x1p12 */
return logf(2*x - 1/(x+sqrtf(x*x-1)));
/* x >= 0x1p12 */
return logf(x) + 0.693147180559945309417232121458176568f;
}
src/math/acoshl.c
浏览文件 @ 969ddbc4
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_acoshl.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* acoshl(x)
* Method :
* Based on
* acoshl(x) = logl [ x + sqrtl(x*x-1) ]
* we have
* acoshl(x) := logl(x)+ln2, if x is large; else
* acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
* acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
*
* Special cases:
* acoshl(x) is NaN with signal if x<1.
* acoshl(NaN) is NaN without signal.
*/
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
......@@ -31,29 +6,20 @@ long double acoshl(long double x)
return acosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
/* acosh(x) = log(x + sqrt(x*x-1)) */
long double acoshl(long double x)
{
long double t;
uint32_t se,i0,i1;
union {
long double f;
struct{uint64_t m; int16_t se; uint16_t pad;} i;
} u = {.f = x};
GET_LDOUBLE_WORDS(se, i0, i1, x);
if (se < 0x3fff || se & 0x8000) { /* x < 1 */
return (x-x)/(x-x);
} else if (se >= 0x401d) { /* x > 2**30 */
if (se >= 0x7fff) /* x is inf or NaN */
return x+x;
return logl(x) + ln2; /* acoshl(huge) = logl(2x) */
} else if (((se-0x3fff)|i0|i1) == 0) {
return 0.0; /* acosh(1) = 0 */
} else if (se > 0x4000) { /* x > 2 */
t = x*x;
return logl(2.0*x - 1.0/(x + sqrtl(t - 1.0)));
}
/* 1 < x <= 2 */
t = x - 1.0;
return log1pl(t + sqrtl(2.0*t + t*t));
if (u.i.se < 0x3fff + 1)
/* x < 2, invalid if x < 1 or nan */
return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1)));
if (u.i.se < 0x3fff + 32)
/* x < 0x1p32 */
return logl(2*x - 1/(x+sqrtl(x*x-1)));
return logl(x) + 0.693147180559945309417232121458176568L;
}
#endif
src/math/acosl.c
浏览文件 @ 969ddbc4
......@@ -23,55 +23,46 @@ long double acosl(long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
#define ACOS_CONST (BIAS - 65) /* 2**-65 */
long double acosl(long double x)
{
union IEEEl2bits u;
long double z, p, q, r, w, s, c, df;
long double z, w, s, c, df;
int16_t expsign, expt;
u.e = x;
expsign = u.xbits.expsign;
expt = expsign & 0x7fff;
if (expt >= BIAS) { /* |x| >= 1 */
if (expt == BIAS &&
/* |x| >= 1 or nan */
if (expt >= 0x3fff) {
if (expt == 0x3fff &&
((u.bits.manh & ~LDBL_NBIT) | u.bits.manl) == 0) {
if (expsign > 0)
return 0.0; /* acos(1) = 0 */
else
// FIXME
return pi_hi + 2.0 * pio2_lo; /* acos(-1)= pi */
return 0; /* acos(1) = 0 */
return 2*pio2_hi + 0x1p-1000; /* acos(-1)= pi */
}
return (x - x) / (x - x); /* acos(|x|>1) is NaN */
return 0/(x-x); /* acos(|x|>1) is NaN */
}
if (expt < BIAS - 1) { /* |x| < 0.5 */
if (expt < ACOS_CONST)
return pio2_hi + pio2_lo; /* x tiny: acosl=pi/2 */
z = x * x;
p = P(z);
q = Q(z);
r = p / q;
return pio2_hi - (x - (pio2_lo - x * r));
} else if (expsign < 0) { /* x < -0.5 */
/* |x| < 0.5 */
if (expt < 0x3fff - 1) {
if (expt < 0x3fff - 65)
return pio2_hi + 0x1p-1000; /* x < 0x1p-65: acosl(x)=pi/2 */
return pio2_hi - (x - (pio2_lo - x * __invtrigl_R(x*x)));
}
/* x < -0.5 */
if (expsign < 0) {
z = (1.0 + x) * 0.5;
p = P(z);
q = Q(z);
s = sqrtl(z);
r = p / q;
w = r * s - pio2_lo;
return pi_hi - 2.0 * (s + w);
} else { /* x > 0.5 */
z = (1.0 - x) * 0.5;
s = sqrtl(z);
u.e = s;
u.bits.manl = 0;
df = u.e;
c = (z - df * df) / (s + df);
p = P(z);
q = Q(z);
r = p / q;
w = r * s + c;
return 2.0 * (df + w);
w = __invtrigl_R(z) * s - pio2_lo;
return 2*(pio2_hi - (s + w));
}
/* x > 0.5 */
z = (1.0 - x) * 0.5;
s = sqrtl(z);
u.e = s;
u.bits.manl = 0;
df = u.e;
c = (z - df * df) / (s + df);
w = __invtrigl_R(z) * s + c;
return 2*(df + w);
}
#endif
src/math/asin.c
浏览文件 @ 969ddbc4
......@@ -42,10 +42,8 @@
#include "libm.h"
static const double
huge = 1.000e+300,
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
/* coefficients for R(x^2) */
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
......@@ -58,51 +56,54 @@ qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
static double R(double z)
{
double p, q;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
return p/q;
}
double asin(double x)
{
double t=0.0,w,p,q,c,r,s;
int32_t hx,ix;
double z,r,s;
uint32_t hx,ix;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x3ff00000) { /* |x|>= 1 */
/* |x| >= 1 or nan */
if (ix >= 0x3ff00000) {
uint32_t lx;
GET_LOW_WORD(lx, x);
if ((ix-0x3ff00000 | lx) == 0)
/* asin(1) = +-pi/2 with inexact */
return x*pio2_hi + x*pio2_lo;
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3fe00000) { /* |x|<0.5 */
if (ix < 0x3e500000) { /* if |x| < 2**-26 */
if (huge+x > 1.0)
return x; /* return x with inexact if x!=0*/
return x*pio2_hi + 0x1p-1000;
return 0/(x-x);
}
/* |x| < 0.5 */
if (ix < 0x3fe00000) {
if (ix < 0x3e500000) {
/* |x|<0x1p-26, return x with inexact if x!=0*/
FORCE_EVAL(x + 0x1p1000);
return x;
}
t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q;
return x + x*w;
return x + x*R(x*x);
}
/* 1 > |x| >= 0.5 */
w = 1.0 - fabs(x);
t = w*0.5;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
s = sqrt(t);
if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */
w = p/q;
t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
z = (1 - fabs(x))*0.5;
s = sqrt(z);
r = R(z);
if (ix >= 0x3fef3333) { /* if |x| > 0.975 */
x = pio2_hi-(2*(s+s*r)-pio2_lo);
} else {
w = s;
SET_LOW_WORD(w,0);
c = (t-w*w)/(s+w);
r = p/q;
p = 2.0*s*r-(pio2_lo-2.0*c);
q = pio4_hi - 2.0*w;
t = pio4_hi - (p-q);
double f,c;
/* f+c = sqrt(z) */
f = s;
SET_LOW_WORD(f,0);
c = (z-f*f)/(s+f);
x = 0.5*pio2_hi - (2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
}
if (hx > 0)
return t;
return -t;
if (hx >> 31)
return -x;
return x;
}
src/math/asinf.c
浏览文件 @ 969ddbc4
......@@ -12,52 +12,51 @@
* is preserved.
* ====================================================
*/
#include "libm.h"
static const double
pio2 = 1.570796326794896558e+00;
static const float
huge = 1.000e+30,
/* coefficients for R(x^2) */
pS0 = 1.6666586697e-01,
pS1 = -4.2743422091e-02,
pS2 = -8.6563630030e-03,
qS1 = -7.0662963390e-01;
static const double
pio2 = 1.570796326794896558e+00;
static float R(float z)
{
float p, q;
p = z*(pS0+z*(pS1+z*pS2));
q = 1.0f+z*qS1;
return p/q;
}
float asinf(float x)
{
double s;
float t,w,p,q;
int32_t hx,ix;
float z;
uint32_t hx,ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x3f800000) { /* |x| >= 1 */
if (ix == 0x3f800000) /* |x| == 1 */
return x*pio2; /* asin(+-1) = +-pi/2 with inexact */
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3f000000) { /* |x|<0.5 */
return x*pio2 + 0x1p-120f; /* asin(+-1) = +-pi/2 with inexact */
return 0/(x-x); /* asin(|x|>1) is NaN */
}
if (ix < 0x3f000000) { /* |x| < 0.5 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
if (huge+x > 1.0f)
return x; /* return x with inexact if x!=0 */
FORCE_EVAL(x + 0x1p120f);
return x; /* return x with inexact if x!=0 */
}
t = x*x;
p = t*(pS0+t*(pS1+t*pS2));
q = 1.0f+t*qS1;
w = p/q;
return x + x*w;
return x + x*R(x*x);
}
/* 1 > |x| >= 0.5 */
w = 1.0f - fabsf(x);
t = w*0.5f;
p = t*(pS0+t*(pS1+t*pS2));
q = 1.0f+t*qS1;
s = sqrt(t);
w = p/q;
t = pio2-2.0*(s+s*w);
if (hx > 0)
return t;
return -t;
z = (1 - fabsf(x))*0.5f;
s = sqrt(z);
x = pio2 - 2*(s+s*R(z));
if (hx >> 31)
return -x;
return x;
}
src/math/asinh.c
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/s_asinh.c */
/*
* ====================================================
* 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.
* ====================================================
*/
/* asinh(x)
* Method :
* Based on
* asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
* we have
* asinh(x) := x if 1+x*x=1,
* := sign(x)*(log(x)+ln2)) for large |x|, else
* := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
*/
#include "libm.h"
static const double
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge= 1.00000000000000000000e+300;
/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
double asinh(double x)
{
double t,w;
int32_t hx,ix;
union {double f; uint64_t i;} u = {.f = x};
unsigned e = u.i >> 52 & 0x7ff;
unsigned s = u.i >> 63;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) /* x is inf or NaN */
return x+x;
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
if (huge+x > 1.0)
return x;
}
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 + 1.0/(sqrt(x*x+1.0)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t)));
/* |x| */
u.i &= (uint64_t)-1/2;
x = u.f;
if (e >= 0x3ff + 26) {
/* |x| >= 0x1p26 or inf or nan */
x = log(x) + 0.693147180559945309417232121458176568;
} else if (e >= 0x3ff + 1) {
/* |x| >= 2 */
x = log(2*x + 1/(sqrt(x*x+1)+x));
} else if (e >= 0x3ff - 26) {
/* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */
x = log1p(x + x*x/(sqrt(x*x+1)+1));
} else {
/* |x| < 0x1p-26, raise inexact if x != 0 */
FORCE_EVAL(x + 0x1p1000);
}
if (hx > 0)
return w;
return -w;
return s ? -x : x;
}
src/math/asinhf.c
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/s_asinhf.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 "libm.h"
static const float
ln2 = 6.9314718246e-01, /* 0x3f317218 */
huge= 1.0000000000e+30;
/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
float asinhf(float x)
{
float t,w;
int32_t hx,ix;
union {float f; uint32_t i;} u = {.f = x};
uint32_t i = u.i & 0x7fffffff;
unsigned s = u.i >> 31;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) /* x is inf or NaN */
return x+x;
if (ix < 0x31800000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
if (huge+x > 1.0f)
return x;
}
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(2.0f*t + 1.0f/(sqrtf(x*x+1.0f)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t)));
/* |x| */
u.i = i;
x = u.f;
if (i >= 0x3f800000 + (12<<23)) {
/* |x| >= 0x1p12 or inf or nan */
x = logf(x) + 0.693147180559945309417232121458176568f;
} else if (i >= 0x3f800000 + (1<<23)) {
/* |x| >= 2 */
x = logf(2*x + 1/(sqrtf(x*x+1)+x));
} else if (i >= 0x3f800000 - (12<<23)) {
/* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */
x = log1pf(x + x*x/(sqrtf(x*x+1)+1));
} else {
/* |x| < 0x1p-12, raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p120f);
}
if (hx > 0)
return w;
return -w;
return s ? -x : x;
}
src/math/asinhl.c
浏览文件 @ 969ddbc4
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_asinhl.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* asinhl(x)
* Method :
* Based on
* asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ]
* we have
* asinhl(x) := x if 1+x*x=1,
* := signl(x)*(logl(x)+ln2)) for large |x|, else
* := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else
* := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2)))
*/
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
......@@ -28,35 +6,33 @@ long double asinhl(long double x)
return asinh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
huge = 1.000000000000000000e+4900L;
/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
long double asinhl(long double x)
{
long double t,w;
int32_t hx,ix;
union {
long double f;
struct{uint64_t m; uint16_t se; uint16_t pad;} i;
} u = {.f = x};
unsigned e = u.i.se & 0x7fff;
unsigned s = u.i.se >> 15;
GET_LDOUBLE_EXP(hx, x);
ix = hx & 0x7fff;
if (ix == 0x7fff)
return x + x; /* x is inf or NaN */
if (ix < 0x3fde) { /* |x| < 2**-34 */
/* return x, raise inexact if x != 0 */
if (huge+x > 1.0)
return x;
}
if (ix > 0x4020) { /* |x| > 2**34 */
w = logl(fabsl(x)) + ln2;
} else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */
t = fabsl(x);
w = logl(2.0*t + 1.0/(sqrtl(x*x + 1.0) + t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1pl(fabsl(x) + t/(1.0 + sqrtl(1.0 + t)));
/* |x| */
u.i.se = e;
x = u.f;
if (e >= 0x3fff + 32) {
/* |x| >= 0x1p32 or inf or nan */
x = logl(x) + 0.693147180559945309417232121458176568L;
} else if (e >= 0x3fff + 1) {
/* |x| >= 2 */
x = logl(2*x + 1/(sqrtl(x*x+1)+x));
} else if (e >= 0x3fff - 32) {
/* |x| >= 0x1p-32 */
x = log1pl(x + x*x/(sqrtl(x*x+1)+1));
} else {
/* |x| < 0x1p-32, raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p1000);
}
if (hx & 0x8000)
return -w;
return w;
return s ? -x : x;
}
#endif
src/math/asinl.c
浏览文件 @ 969ddbc4
......@@ -23,60 +23,49 @@ long double asinl(long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
static const long double huge = 1.000e+300;
static const long double pio4_hi = 7.85398163397448309628e-01L;
#define ASIN_LINEAR (BIAS - 32) /* 2**-32 */
/* 0.95 */
#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
#define THRESH ((0xe666666666666666ULL>>(64-(LDBL_MANH_SIZE-1)))|LDBL_NBIT)
long double asinl(long double x)
{
union IEEEl2bits u;
long double t=0.0,w,p,q,c,r,s;
int16_t expsign, expt;
long double z,r,s;
uint16_t expsign, expt;
u.e = x;
expsign = u.xbits.expsign;
expt = expsign & 0x7fff;
if (expt >= BIAS) { /* |x|>= 1 */
if (expt == BIAS &&
if (expt >= 0x3fff) { /* |x| >= 1 or nan */
if (expt == 0x3fff &&
((u.bits.manh&~LDBL_NBIT)|u.bits.manl) == 0)
/* asin(1)=+-pi/2 with inexact */
return x*pio2_hi + x*pio2_lo;
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (expt < BIAS-1) { /* |x|<0.5 */
if (expt < ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */
/* asin(+-1)=+-pi/2 with inexact */
return x*pio2_hi + 0x1p-1000;
return 0/(x-x);
}
if (expt < 0x3fff - 1) { /* |x| < 0.5 */
if (expt < 0x3fff - 32) { /* |x|<0x1p-32, asinl(x)=x */
/* return x with inexact if x!=0 */
if (huge+x > 1.0)
return x;
FORCE_EVAL(x + 0x1p1000);
return x;
}
t = x*x;
p = P(t);
q = Q(t);
w = p/q;
return x + x*w;
return x + x*__invtrigl_R(x*x);
}
/* 1 > |x| >= 0.5 */
w = 1.0 - fabsl(x);
t = w*0.5;
p = P(t);
q = Q(t);
s = sqrtl(t);
z = (1.0 - fabsl(x))*0.5;
s = sqrtl(z);
r = __invtrigl_R(z);
if (u.bits.manh >= THRESH) { /* if |x| is close to 1 */
w = p/q;
t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
x = pio2_hi - (2*(s+s*r)-pio2_lo);
} else {
long double f, c;
u.e = s;
u.bits.manl = 0;
w = u.e;
c = (t-w*w)/(s+w);
r = p/q;
p = 2.0*s*r-(pio2_lo-2.0*c);
q = pio4_hi-2.0*w;
t = pio4_hi-(p-q);
f = u.e;
c = (z-f*f)/(s+f);
x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
}
if (expsign > 0)
return t;
return -t;
if (expsign>>15)
return -x;
return x;
}
#endif
src/math/atan.c
浏览文件 @ 969ddbc4
......@@ -60,32 +60,26 @@ static const double aT[] = {
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
};
static const double huge = 1.0e300;
double atan(double x)
{
double w,s1,s2,z;
int32_t ix,hx,id;
uint32_t ix,sign;
int id;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
GET_HIGH_WORD(ix, x);
sign = ix >> 31;
ix &= 0x7fffffff;
if (ix >= 0x44100000) { /* if |x| >= 2^66 */
uint32_t low;
GET_LOW_WORD(low, x);
if (ix > 0x7ff00000 ||
(ix == 0x7ff00000 && low != 0)) /* NaN */
return x+x;
if (hx > 0)
return atanhi[3] + *(volatile double *)&atanlo[3];
else
return -atanhi[3] - *(volatile double *)&atanlo[3];
if (isnan(x))
return x;
z = atanhi[3] + 0x1p-1000;
return sign ? -z : z;
}
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e400000) { /* |x| < 2^-27 */
/* raise inexact */
if (huge+x > 1.0)
return x;
/* raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p1000);
return x;
}
id = -1;
} else {
......@@ -117,5 +111,5 @@ double atan(double x)
if (id < 0)
return x - x*(s1+s2);
z = atanhi[id] - (x*(s1+s2) - atanlo[id] - x);
return hx < 0 ? -z : z;
return sign ? -z : z;
}
src/math/atan2l.c
浏览文件 @ 969ddbc4
......@@ -24,8 +24,6 @@ long double atan2l(long double y, long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
// FIXME:
static const volatile long double tiny = 1.0e-300;
long double atan2l(long double y, long double x)
{
......@@ -40,12 +38,12 @@ long double atan2l(long double y, long double x)
ux.e = x;
expsignx = ux.xbits.expsign;
exptx = expsignx & 0x7fff;
if ((exptx==BIAS+LDBL_MAX_EXP &&
if ((exptx==0x7fff &&
((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) || /* x is NaN */
(expty==BIAS+LDBL_MAX_EXP &&
(expty==0x7fff &&
((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* y is NaN */
return x+y;
if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) /* x=1.0 */
if (expsignx==0x3fff && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) /* x=1.0 */
return atanl(y);
m = ((expsigny>>15)&1) | ((expsignx>>14)&2); /* 2*sign(x)+sign(y) */
......@@ -54,39 +52,39 @@ long double atan2l(long double y, long double x)
switch(m) {
case 0:
case 1: return y; /* atan(+-0,+anything)=+-0 */
case 2: return pi_hi+tiny; /* atan(+0,-anything) = pi */
case 3: return -pi_hi-tiny; /* atan(-0,-anything) =-pi */
case 2: return 2*pio2_hi+0x1p-1000; /* atan(+0,-anything) = pi */
case 3: return -2*pio2_hi-0x1p-1000; /* atan(-0,-anything) =-pi */
}
}
/* when x = 0 */
if (exptx==0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
return expsigny < 0 ? -pio2_hi-tiny : pio2_hi+tiny;
return expsigny < 0 ? -pio2_hi-0x1p-1000 : pio2_hi+0x1p-1000;
/* when x is INF */
if (exptx == BIAS+LDBL_MAX_EXP) {
if (expty == BIAS+LDBL_MAX_EXP) {
if (exptx == 0x7fff) {
if (expty == 0x7fff) {
switch(m) {
case 0: return pio2_hi*0.5+tiny; /* atan(+INF,+INF) */
case 1: return -pio2_hi*0.5-tiny; /* atan(-INF,+INF) */
case 2: return 1.5*pio2_hi+tiny; /* atan(+INF,-INF) */
case 3: return -1.5*pio2_hi-tiny; /* atan(-INF,-INF) */
case 0: return pio2_hi*0.5+0x1p-1000; /* atan(+INF,+INF) */
case 1: return -pio2_hi*0.5-0x1p-1000; /* atan(-INF,+INF) */
case 2: return 1.5*pio2_hi+0x1p-1000; /* atan(+INF,-INF) */
case 3: return -1.5*pio2_hi-0x1p-1000; /* atan(-INF,-INF) */
}
} else {
switch(m) {
case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi_hi+tiny; /* atan(+...,-INF) */
case 3: return -pi_hi-tiny; /* atan(-...,-INF) */
case 2: return 2*pio2_hi+0x1p-1000; /* atan(+...,-INF) */
case 3: return -2*pio2_hi-0x1p-1000; /* atan(-...,-INF) */
}
}
}
/* when y is INF */
if (expty == BIAS+LDBL_MAX_EXP)
return expsigny < 0 ? -pio2_hi-tiny : pio2_hi+tiny;
if (expty == 0x7fff)
return expsigny < 0 ? -pio2_hi-0x1p-1000 : pio2_hi+0x1p-1000;
/* compute y/x */
k = expty-exptx;
if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */
z = pio2_hi+pio2_lo;
z = pio2_hi+0x1p-1000;
m &= 1;
} else if (expsignx < 0 && k < -LDBL_MANT_DIG-2) /* |y/x| tiny, x<0 */
z = 0.0;
......@@ -95,9 +93,9 @@ long double atan2l(long double y, long double x)
switch (m) {
case 0: return z; /* atan(+,+) */
case 1: return -z; /* atan(-,+) */
case 2: return pi_hi-(z-pi_lo); /* atan(+,-) */
case 2: return 2*pio2_hi-(z-2*pio2_lo); /* atan(+,-) */
default: /* case 3 */
return (z-pi_lo)-pi_hi; /* atan(-,-) */
return (z-2*pio2_lo)-2*pio2_hi; /* atan(-,-) */
}
}
#endif
src/math/atanf.c
浏览文件 @ 969ddbc4
......@@ -38,28 +38,26 @@ static const float aT[] = {
6.1687607318e-02,
};
static const float huge = 1.0e30;
float atanf(float x)
{
float w,s1,s2,z;
int32_t ix,hx,id;
uint32_t ix,sign;
int id;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
GET_FLOAT_WORD(ix, x);
sign = ix>>31;
ix &= 0x7fffffff;
if (ix >= 0x4c800000) { /* if |x| >= 2**26 */
if (ix > 0x7f800000) /* NaN */
return x+x;
if (hx > 0)
return atanhi[3] + *(volatile float *)&atanlo[3];
else
return -atanhi[3] - *(volatile float *)&atanlo[3];
if (isnan(x))
return x;
z = atanhi[3] + 0x1p-120f;
return sign ? -z : z;
}
if (ix < 0x3ee00000) { /* |x| < 0.4375 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* raise inexact */
if(huge+x>1.0f)
return x;
/* raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p120f);
return x;
}
id = -1;
} else {
......@@ -91,5 +89,5 @@ float atanf(float x)
if (id < 0)
return x - x*(s1+s2);
z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
return hx < 0 ? -z : z;
return sign ? -z : z;
}
src/math/atanh.c
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/e_atanh.c */
/*
* ====================================================
* 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.
* ====================================================
*
*/
/* atanh(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
* 1 2x x
* atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x
*
* For x<0.5
* atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
*
* Special cases:
* atanh(x) is NaN if |x| > 1 with signal;
* atanh(NaN) is that NaN with no signal;
* atanh(+-1) is +-INF with signal.
*
*/
#include "libm.h"
static const double huge = 1e300;
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
double atanh(double x)
{
double t;
int32_t hx,ix;
uint32_t lx;
union {double f; uint64_t i;} u = {.f = x};
unsigned e = u.i >> 52 & 0x7ff;
unsigned s = u.i >> 63;
/* |x| */
u.i &= (uint64_t)-1/2;
x = u.f;
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3ff00000)
return x/0.0;
if (ix < 0x3e300000 && (huge+x) > 0.0) /* x < 2**-28 */
return x;
SET_HIGH_WORD(x, ix);
if (ix < 0x3fe00000) { /* x < 0.5 */
t = x+x;
t = 0.5*log1p(t + t*x/(1.0-x));
} else
t = 0.5*log1p((x+x)/(1.0-x));
if (hx >= 0)
return t;
return -t;
if (e < 0x3ff - 1) {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5*log1p(2*x + 2*x*x/(1-x));
} else {
x = 0.5*log1p(2*x/(1-x));
}
return s ? -x : x;
}
src/math/atanhf.c
浏览文件 @ 969ddbc4
/* origin: FreeBSD /usr/src/lib/msun/src/e_atanhf.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 "libm.h"
static const float huge = 1e30;
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
float atanhf(float x)
{
float t;
int32_t hx,ix;
union {float f; uint32_t i;} u = {.f = x};
unsigned s = u.i >> 31;
/* |x| */
u.i &= 0x7fffffff;
x = u.f;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix > 0x3f800000) /* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3f800000)
return x/0.0f;
if (ix < 0x31800000 && huge+x > 0.0f) /* x < 2**-28 */
return x;
SET_FLOAT_WORD(x, ix);
if (ix < 0x3f000000) { /* x < 0.5 */
t = x+x;
t = 0.5f*log1pf(t + t*x/(1.0f-x));
} else
t = 0.5f*log1pf((x+x)/(1.0f-x));
if (hx >= 0)
return t;
return -t;
if (u.i < 0x3f800000 - (1<<23)) {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5f*log1pf(2*x + 2*x*x/(1-x));
} else {
x = 0.5f*log1pf(2*x/(1-x));
}
return s ? -x : x;
}
src/math/atanhl.c
浏览文件 @ 969ddbc4
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_atanh.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* atanhl(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
* 1 2x x
* atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x
*
* For x<0.5
* atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
*
* Special cases:
* atanhl(x) is NaN if |x| > 1 with signal;
* atanhl(NaN) is that NaN with no signal;
* atanhl(+-1) is +-INF with signal.
*/
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
......@@ -34,31 +6,26 @@ long double atanhl(long double x)
return atanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double huge = 1e4900L;
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
long double atanhl(long double x)
{
long double t;
int32_t ix;
uint32_t se,i0,i1;
union {
long double f;
struct{uint64_t m; uint16_t se; uint16_t pad;} i;
} u = {.f = x};
unsigned e = u.i.se & 0x7fff;
unsigned s = u.i.se >> 15;
/* |x| */
u.i.se = e;
x = u.f;
GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff;
if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31)) > 0x3fff)
/* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3fff)
return x/0.0;
if (ix < 0x3fe3 && huge+x > 0.0) /* x < 2**-28 */
return x;
SET_LDOUBLE_EXP(x, ix);
if (ix < 0x3ffe) { /* x < 0.5 */
t = x + x;
t = 0.5*log1pl(t + t*x/(1.0 - x));
} else
t = 0.5*log1pl((x + x)/(1.0 - x));
if (se <= 0x7fff)
return t;
return -t;
if (e < 0x3fff - 1) {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5*log1pl(2*x + 2*x*x/(1-x));
} else {
x = 0.5*log1pl(2*x/(1-x));
}
return s ? -x : x;
}
#endif
src/math/atanl.c
浏览文件 @ 969ddbc4
......@@ -22,11 +22,6 @@ long double atanl(long double x)
return atan(x);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
#define ATAN_CONST (BIAS + 65) /* 2**65 */
#define ATAN_LINEAR (BIAS - 32) /* 2**-32 */
static const long double huge = 1.0e300;
static const long double atanhi[] = {
4.63647609000806116202e-01L,
......@@ -81,29 +76,27 @@ long double atanl(long double x)
u.e = x;
expsign = u.xbits.expsign;
expt = expsign & 0x7fff;
if (expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */
if (expt == BIAS + LDBL_MAX_EXP &&
if (expt >= 0x3fff + 65) { /* if |x| is large, atan(x)~=pi/2 */
if (expt == 0x7fff &&
((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0) /* NaN */
return x+x;
if (expsign > 0)
return atanhi[3]+atanlo[3];
else
return -atanhi[3]-atanlo[3];
z = atanhi[3] + 0x1p-1000;
return expsign < 0 ? -z : z;
}
/* Extract the exponent and the first few bits of the mantissa. */
/* XXX There should be a more convenient way to do this. */
expman = (expt << 8) | ((u.bits.manh >> (MANH_SIZE - 9)) & 0xff);
if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */
/* raise inexact */
if (huge+x > 1.0)
return x;
expman = (expt << 8) | ((u.bits.manh >> (LDBL_MANH_SIZE - 9)) & 0xff);
if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < 0x3fff - 32) { /* if |x| is small, atanl(x)~=x */
/* raise inexact if x!=0 */
FORCE_EVAL(x + 0x1p1000);
return x;
}
id = -1;
} else {
x = fabsl(x);
if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */
if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
if (expman < (0x3fff << 8) + 0x30) { /* |x| < 1.1875 */
if (expman < ((0x3fff - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
id = 0;
x = (2.0*x-1.0)/(2.0+x);
} else { /* 11/16 <= |x| < 19/16 */
......@@ -111,7 +104,7 @@ long double atanl(long double x)
x = (x-1.0)/(x+1.0);
}
} else {
if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
if (expman < ((0x3fff + 1) << 8) + 0x38) { /* |x| < 2.4375 */
id = 2;
x = (x-1.5)/(1.0+1.5*x);
} else { /* 2.4375 <= |x| < 2^ATAN_CONST */
......
src/math/cosh.c
浏览文件 @ 969ddbc4
......@@ -23,7 +23,7 @@
* 2
* 22 <= x <= lnovft : cosh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : cosh(x) := huge*huge (overflow)
* ln2ovft < x : cosh(x) := inf (overflow)
*
* Special cases:
* cosh(x) is |x| if x is +INF, -INF, or NaN.
......@@ -32,43 +32,40 @@
#include "libm.h"
static const double huge = 1.0e300;
double cosh(double x)
{
double t, w;
int32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
union {double f; uint64_t i;} u = {.f = x};
uint32_t ix;
double t;
/* x is INF or NaN */
if (ix >= 0x7ff00000)
return x*x;
/* |x| */
u.i &= (uint64_t)-1/2;
x = u.f;
ix = u.i >> 32;
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3fd62e43) {
t = expm1(fabs(x));
w = 1.0+t;
t = expm1(x);
if (ix < 0x3c800000)
return w; /* cosh(tiny) = 1 */
return 1.0 + (t*t)/(w+w);
return 1;
return 1 + t*t/(2*(1+t));
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x40360000) {
t = exp(fabs(x));
t = exp(x);
return 0.5*t + 0.5/t;
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862E42)
return 0.5*exp(fabs(x));
if (ix < 0x40862e42)
return 0.5*exp(x);
/* |x| in [log(maxdouble), overflowthresold] */
if (ix <= 0x408633CE)
return __expo2(fabs(x));
if (ix <= 0x408633ce)
return __expo2(x);
/* |x| > overflowthresold, cosh(x) overflow */
return huge*huge;
/* overflow (or nan) */
x *= 0x1p1023;
return x;
}
src/math/coshf.c
浏览文件 @ 969ddbc4
......@@ -15,43 +15,40 @@
#include "libm.h"
static const float huge = 1.0e30;
float coshf(float x)
{
float t, w;
int32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
union {float f; uint32_t i;} u = {.f = x};
uint32_t ix;
float t;
/* x is INF or NaN */
if (ix >= 0x7f800000)
return x*x;
/* |x| */
u.i &= 0x7fffffff;
x = u.f;
ix = u.i;
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3eb17218) {
t = expm1f(fabsf(x));
w = 1.0f+t;
if (ix<0x39800000)
return 1.0f; /* cosh(tiny) = 1 */
return 1.0f + (t*t)/(w+w);
t = expm1f(x);
if (ix < 0x39800000)
return 1;
return 1 + t*t/(2*(1+t));
}
/* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x41100000) {
t = expf(fabsf(x));
t = expf(x);
return 0.5f*t + 0.5f/t;
}
/* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */
if (ix < 0x42b17217)
return 0.5f*expf(fabsf(x));
return 0.5f*expf(x);
/* |x| in [log(maxfloat), overflowthresold] */
if (ix <= 0x42b2d4fc)
return __expo2f(fabsf(x));
return __expo2f(x);
/* |x| > overflowthresold, cosh(x) overflow */
return huge*huge;
/* |x| > overflowthresold or nan */
x *= 0x1p127f;
return x;
}
src/math/coshl.c
浏览文件 @ 969ddbc4
......@@ -23,7 +23,7 @@
* 2
* 22 <= x <= lnovft : coshl(x) := expl(x)/2
* lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2)
* ln2ovft < x : coshl(x) := huge*huge (overflow)
* ln2ovft < x : coshl(x) := inf (overflow)
*
* Special cases:
* coshl(x) is |x| if x is +INF, -INF, or NaN.
......@@ -38,49 +38,48 @@ long double coshl(long double x)
return cosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double huge = 1.0e4900L;
long double coshl(long double x)
{
long double t,w;
int32_t ex;
union {
long double f;
struct{uint64_t m; uint16_t se; uint16_t pad;} i;
} u = {.f = x};
unsigned ex = u.i.se & 0x7fff;
long double t;
uint32_t mx,lx;
/* High word of |x|. */
GET_LDOUBLE_WORDS(ex, mx, lx, x);
ex &= 0x7fff;
/* x is INF or NaN */
if (ex == 0x7fff) return x*x;
/* |x| */
u.i.se = ex;
x = u.f;
mx = u.i.m >> 32;
lx = u.i.m;
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
if (ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
t = expm1l(fabsl(x));
w = 1.0 + t;
if (ex < 0x3fbc) return w; /* cosh(tiny) = 1 */
return 1.0+(t*t)/(w+w);
if (ex < 0x3fff-2 || (ex == 0x3fff-2 && mx < 0xb17217f7)) {
t = expm1l(x);
if (ex < 0x3fff-64)
return 1;
return 1 + t*t/(2*(1+t));
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
t = expl(fabsl(x));
if (ex < 0x3fff+4 || (ex == 0x3fff+4 && mx < 0xb0000000)) {
t = expl(x);
return 0.5*t + 0.5/t;
}
/* |x| in [22, ln(maxdouble)] return 0.5*exp(|x|) */
if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
return 0.5*expl(fabsl(x));
if (ex < 0x3fff+13 || (ex == 0x3fff+13 && mx < 0xb1700000))
return 0.5*expl(x);
/* |x| in [log(maxdouble), log(2*maxdouble)) */
if (ex == 0x400c && (mx < 0xb174ddc0u ||
(mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
{
w = expl(0.5*fabsl(x));
t = 0.5*w;
return t*w;
if (ex == 0x3fff+13 && (mx < 0xb174ddc0 ||
(mx == 0xb174ddc0 && lx < 0x31aec0eb))) {
t = expl(0.5*x);
return 0.5*t*t;
}
/* |x| >= log(2*maxdouble), cosh(x) overflow */
return huge*huge;
/* |x| >= log(2*maxdouble) or nan */
return x*0x1p16383L;
}
#endif
src/math/i386/exp.s
浏览文件 @ 969ddbc4
......@@ -50,12 +50,6 @@ expf:
flds 4(%esp)
jmp 2f
.global expl
.type expl,@function
expl:
fldt 4(%esp)
jmp 2f
.global exp
.type exp,@function
exp:
......
src/math/i386/expl.s
浏览文件 @ 969ddbc4
# see exp.s
# exp(x) = 2^hi + 2^hi (2^lo - 1)
# where hi+lo = log2e*x with 128bit precision
# exact log2e*x calculation depends on nearest rounding mode
.global expl
.type expl,@function
expl:
fldt 4(%esp)
# special cases: 2*x is +-inf, nan or |x| < 0x1p-32
# check (exponent|0x8000)+2 < 0xbfff+2-32
movw 12(%esp), %ax
movw %ax, %dx
orw $0x8000, %dx
addw $2, %dx
cmpw $0xbfff-30, %dx
jnb 3f
cmpw $1, %dx
jbe 1f
# if |x|<0x1p-32 return 1+x
fld1
jmp 2f
1: testw %ax, %ax
jns 1f
# if 2*x == -inf,-nan return -0/x
fldz
fchs
fdivp
ret
# if 2*x == inf,nan return 2*x
1: fld %st(0)
2: faddp
ret
# should be 0x1.71547652b82fe178p0 == 0x3fff b8aa3b29 5c17f0bc
# it will be wrong on non-nearest rounding mode
3: fldl2e
# subl $32, %esp
subl $44, %esp
# hi = log2e_hi*x
# 2^hi = exp2l(hi)
fmul %st(1),%st
fld %st(0)
fstpt (%esp)
fstpt 16(%esp)
fstpt 32(%esp)
call exp2l
# if 2^hi == inf return 2^hi
fld %st(0)
fstpt (%esp)
cmpw $0x7fff, 8(%esp)
je 1f
fldt 32(%esp)
fldt 16(%esp)
# fpu stack: 2^hi x hi
# exact mult: x*log2e
fld %st(1) # x
# c = 0x1p32+1
pushl $0x41f00000
pushl $0x00100000
fldl (%esp)
# xh = x - c*x + c*x
# xl = x - xh
fmulp
fld %st(2)
fsub %st(1), %st
faddp
fld %st(2)
fsub %st(1), %st
# yh = log2e_hi - c*log2e_hi + c*log2e_hi
pushl $0x3ff71547
pushl $0x65200000
fldl (%esp)
# fpu stack: 2^hi x hi xh xl yh
# lo = hi - xh*yh + xl*yh
fld %st(2)
fmul %st(1), %st
fsubp %st, %st(4)
fmul %st(1), %st
faddp %st, %st(3)
# yl = log2e_hi - yh
pushl $0x3de705fc
pushl $0x2f000000
fldl (%esp)
# fpu stack: 2^hi x lo xh xl yl
# lo += xh*yl + xl*yl
fmul %st, %st(2)
fmulp %st, %st(1)
fxch %st(2)
faddp
faddp
# log2e_lo
pushl $0xbfbe
pushl $0x82f0025f
pushl $0x2dc582ee
fldt (%esp)
addl $36,%esp
# fpu stack: 2^hi x lo log2e_lo
# lo += log2e_lo*x
# return 2^hi + 2^hi (2^lo - 1)
fmulp %st, %st(2)
faddp
f2xm1
fmul %st(1), %st
faddp
1: addl $44, %esp
ret
src/math/tgamma.c
浏览文件 @ 969ddbc4
#include <math.h>
/*
"A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
"Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
"An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
// FIXME: use lanczos approximation
approximation method:
double __lgamma_r(double, int *);
(x - 0.5) S(x)
Gamma(x) = (x + g - 0.5) * ----------------
exp(x + g - 0.5)
with
a1 a2 a3 aN
S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
x + 1 x + 2 x + 3 x + N
with a0, a1, a2, a3,.. aN constants which depend on g.
for x < 0 the following reflection formula is used:
Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
most ideas and constants are from boost and python
*/
#include "libm.h"
static const double pi = 3.141592653589793238462643383279502884;
/* sin(pi x) with x > 0 && isnormal(x) assumption */
static double sinpi(double x)
{
int n;
/* argument reduction: x = |x| mod 2 */
/* spurious inexact when x is odd int */
x = x * 0.5;
x = 2 * (x - floor(x));
/* reduce x into [-.25,.25] */
n = 4 * x;
n = (n+1)/2;
x -= n * 0.5;
x *= pi;
switch (n) {
default: /* case 4 */
case 0:
return __sin(x, 0, 0);
case 1:
return __cos(x, 0);
case 2:
/* sin(0-x) and -sin(x) have different sign at 0 */
return __sin(0-x, 0, 0);
case 3:
return -__cos(x, 0);
}
}
#define N 12
//static const double g = 6.024680040776729583740234375;
static const double gmhalf = 5.524680040776729583740234375;
static const double Snum[N+1] = {
23531376880.410759688572007674451636754734846804940,
42919803642.649098768957899047001988850926355848959,
35711959237.355668049440185451547166705960488635843,
17921034426.037209699919755754458931112671403265390,
6039542586.3520280050642916443072979210699388420708,
1439720407.3117216736632230727949123939715485786772,
248874557.86205415651146038641322942321632125127801,
31426415.585400194380614231628318205362874684987640,
2876370.6289353724412254090516208496135991145378768,
186056.26539522349504029498971604569928220784236328,
8071.6720023658162106380029022722506138218516325024,
210.82427775157934587250973392071336271166969580291,
2.5066282746310002701649081771338373386264310793408,
};
static const double Sden[N+1] = {
0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535,
2637558, 357423, 32670, 1925, 66, 1,
};
/* n! for small integer n */
static const double fact[] = {
1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0,
479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0,
355687428096000.0, 6402373705728000.0, 121645100408832000.0,
2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0,
};
/* S(x) rational function for positive x */
static double S(double x)
{
double num = 0, den = 0;
int i;
/* to avoid overflow handle large x differently */
if (x < 8)
for (i = N; i >= 0; i--) {
num = num * x + Snum[i];
den = den * x + Sden[i];
}
else
for (i = 0; i <= N; i++) {
num = num / x + Snum[i];
den = den / x + Sden[i];
}
return num/den;
}
double tgamma(double x)
{
int sign;
double y;
double absx, y, dy, z, r;
y = exp(__lgamma_r(x, &sign));
if (sign < 0)
y = -y;
return y;
/* special cases */
if (!isfinite(x))
/* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
return x + INFINITY;
/* integer arguments */
/* raise inexact when non-integer */
if (x == floor(x)) {
if (x == 0)
/* tgamma(+-0)=+-inf with divide-by-zero */
return 1/x;
if (x < 0)
return 0/0.0;
if (x <= sizeof fact/sizeof *fact)
return fact[(int)x - 1];
}
absx = fabs(x);
/* x ~ 0: tgamma(x) ~ 1/x */
if (absx < 0x1p-54)
return 1/x;
/* x >= 172: tgamma(x)=inf with overflow */
/* x =< -184: tgamma(x)=+-0 with underflow */
if (absx >= 184) {
if (x < 0) {
if (floor(x) * 0.5 == floor(x * 0.5))
return 0;
return -0.0;
}
x *= 0x1p1023;
return x;
}
/* handle the error of x + g - 0.5 */
y = absx + gmhalf;
if (absx > gmhalf) {
dy = y - absx;
dy -= gmhalf;
} else {
dy = y - gmhalf;
dy -= absx;
}
z = absx - 0.5;
r = S(absx) * exp(-y);
if (x < 0) {
/* reflection formula for negative x */
r = -pi / (sinpi(absx) * absx * r);
dy = -dy;
z = -z;
}
r += dy * (gmhalf+0.5) * r / y;
z = pow(y, 0.5*z);
r = r * z * z;
return r;
}
#if 0
double __lgamma_r(double x, int *sign)
{
double r, absx, z, zz, w;
*sign = 1;
/* special cases */
if (!isfinite(x))
/* lgamma(nan)=nan, lgamma(+-inf)=inf */
return x*x;
/* integer arguments */
if (x == floor(x) && x <= 2) {
/* n <= 0: lgamma(n)=inf with divbyzero */
/* n == 1,2: lgamma(n)=0 */
if (x <= 0)
return 1/0.0;
return 0;
}
absx = fabs(x);
/* lgamma(x) ~ -log(|x|) for tiny |x| */
if (absx < 0x1p-54) {
*sign = 1 - 2*!!signbit(x);
return -log(absx);
}
/* use tgamma for smaller |x| */
if (absx < 128) {
x = tgamma(x);
*sign = 1 - 2*!!signbit(x);
return log(fabs(x));
}
/* second term (log(S)-g) could be more precise here.. */
/* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5));
if (x < 0) {
/* reflection formula for negative x */
x = sinpi(absx);
*sign = 2*!!signbit(x) - 1;
r = log(pi/(fabs(x)*absx)) - r;
}
return r;
}
weak_alias(__lgamma_r, lgamma_r);
#endif
src/math/tgammaf.c
浏览文件 @ 969ddbc4
#include <math.h>
// FIXME: use lanczos approximation
float __lgammaf_r(float, int *);
float tgammaf(float x)
{
int sign;
float y;
y = exp(__lgammaf_r(x, &sign));
if (sign < 0)
y = -y;
return y;
return tgamma(x);
}
Markdown is supported
0% 或 .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑！
想要评论请 注册 或