mirror of
https://github.com/UpsilonNumworks/Upsilon.git
synced 2026-03-18 21:30:38 +01:00
[poincare] Use OpenBSD function for cosh, sinh, tanh and arccosh,
arcsinh and arctanh Change-Id: I64f6718ebdd042512ce9b9db78dffa3c943471ff
This commit is contained in:
@@ -20,28 +20,36 @@ objs += $(addprefix liba/src/, \
|
||||
objs += $(addprefix liba/src/external/openbsd/, \
|
||||
a_atanf.o \
|
||||
e_acosf.o \
|
||||
e_acoshf.o \
|
||||
e_asinf.o \
|
||||
e_atanhf.o \
|
||||
e_coshf.o \
|
||||
e_expf.o \
|
||||
e_lgammaf_r.o \
|
||||
e_log10f.o \
|
||||
e_logf.o \
|
||||
e_powf.o \
|
||||
e_rem_pio2f.o \
|
||||
e_sinhf.o \
|
||||
e_sqrtf.o \
|
||||
k_cosf.o \
|
||||
k_rem_pio2f.o \
|
||||
k_sinf.o \
|
||||
k_tanf.o \
|
||||
s_asinhf.o\
|
||||
s_ceilf.o \
|
||||
s_copysignf.o \
|
||||
s_cosf.o \
|
||||
s_expm1f.o\
|
||||
s_fabsf.o \
|
||||
s_floorf.o \
|
||||
s_log1pf.o \
|
||||
s_roundf.o \
|
||||
s_scalbnf.o \
|
||||
s_signgam.o \
|
||||
s_sinf.o \
|
||||
s_tanf.o \
|
||||
s_tanhf.o \
|
||||
w_lgammaf.o \
|
||||
)
|
||||
|
||||
@@ -52,8 +60,10 @@ ifeq ($(DEBUG),1)
|
||||
# code. That's annoying because Clang rightfully decides to emit double-to-float
|
||||
# aeabi conversions when building in -O0 mode, and we really don't want to code
|
||||
# such functions. A simple workaround is to always build this file -Os.
|
||||
liba/src/external/openbsd/s_roundf.o: CFLAGS += -Os
|
||||
liba/src/external/openbsd/e_expf.o: CFLAGS += -Os
|
||||
liba/src/external/openbsd/s_expm1f.o: CFLAGS += -Os
|
||||
liba/src/external/openbsd/s_log1pf.o: CFLAGS += -Os
|
||||
liba/src/external/openbsd/s_roundf.o: CFLAGS += -Os
|
||||
endif
|
||||
|
||||
tests += $(addprefix liba/test/, \
|
||||
|
||||
@@ -22,24 +22,32 @@ int isnanf(float x);
|
||||
#define isnan(x) isnanf(x)
|
||||
|
||||
float acosf(float x);
|
||||
float acoshf(float x);
|
||||
float asinf(float x);
|
||||
float asinhf(float x);
|
||||
float atanf(float x);
|
||||
float atanhf(float x);
|
||||
float ceilf(float x);
|
||||
float copysignf(float x, float y);
|
||||
float cosf(float x);
|
||||
float coshf(float x);
|
||||
float expf(float x);
|
||||
float expm1f(float x);
|
||||
float fabsf(float x);
|
||||
float floorf(float x);
|
||||
float lgammaf(float x);
|
||||
float lgammaf_r(float x, int *signgamp);
|
||||
float log1pf(float x);
|
||||
float log10f(float x);
|
||||
float logf(float x);
|
||||
float powf(float x, float y);
|
||||
float roundf(float x);
|
||||
float scalbnf(float x, int n);
|
||||
float sinf(float x);
|
||||
float sinhf(float x);
|
||||
float sqrtf(float x);
|
||||
float tanf(float x);
|
||||
float tanhf(float x);
|
||||
|
||||
LIBA_END_DECLS
|
||||
|
||||
|
||||
45
liba/src/external/openbsd/e_acoshf.c
vendored
Normal file
45
liba/src/external/openbsd/e_acoshf.c
vendored
Normal file
@@ -0,0 +1,45 @@
|
||||
/* e_acoshf.c -- float version of e_acosh.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const float
|
||||
one = 1.0,
|
||||
ln2 = 6.9314718246e-01; /* 0x3f317218 */
|
||||
|
||||
float
|
||||
acoshf(float x)
|
||||
{
|
||||
float t;
|
||||
int32_t 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));
|
||||
}
|
||||
}
|
||||
41
liba/src/external/openbsd/e_atanhf.c
vendored
Normal file
41
liba/src/external/openbsd/e_atanhf.c
vendored
Normal file
@@ -0,0 +1,41 @@
|
||||
/* e_atanhf.c -- float version of e_atanh.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const float one = 1.0, huge = 1e30;
|
||||
static const float zero = 0.0;
|
||||
|
||||
float
|
||||
atanhf(float x)
|
||||
{
|
||||
float t;
|
||||
int32_t hx,ix;
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
ix = hx&0x7fffffff;
|
||||
if (ix>0x3f800000) /* |x|>1 */
|
||||
return (x-x)/(x-x);
|
||||
if(ix==0x3f800000)
|
||||
return x/zero;
|
||||
if(ix<0x31800000&&(huge+x)>zero) return x; /* x<2**-28 */
|
||||
SET_FLOAT_WORD(x,ix);
|
||||
if(ix<0x3f000000) { /* x < 0.5 */
|
||||
t = x+x;
|
||||
t = (float)0.5*log1pf(t+t*x/(one-x));
|
||||
} else
|
||||
t = (float)0.5*log1pf((x+x)/(one-x));
|
||||
if(hx>=0) return t; else return -t;
|
||||
}
|
||||
60
liba/src/external/openbsd/e_coshf.c
vendored
Normal file
60
liba/src/external/openbsd/e_coshf.c
vendored
Normal file
@@ -0,0 +1,60 @@
|
||||
/* e_coshf.c -- float version of e_cosh.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const volatile float huge = 1.0e30;
|
||||
static const float one = 1.0, half=0.5;
|
||||
|
||||
float
|
||||
coshf(float x)
|
||||
{
|
||||
float t,w;
|
||||
int32_t ix;
|
||||
|
||||
GET_FLOAT_WORD(ix,x);
|
||||
ix &= 0x7fffffff;
|
||||
|
||||
/* x is INF or NaN */
|
||||
if(ix>=0x7f800000) return x*x;
|
||||
|
||||
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
|
||||
if(ix<0x3eb17218) {
|
||||
t = expm1f(fabsf(x));
|
||||
w = one+t;
|
||||
if (ix<0x24000000) return w; /* cosh(tiny) = 1 */
|
||||
return one+(t*t)/(w+w);
|
||||
}
|
||||
|
||||
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
|
||||
if (ix < 0x41b00000) {
|
||||
t = expf(fabsf(x));
|
||||
return half*t+half/t;
|
||||
}
|
||||
|
||||
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
|
||||
if (ix < 0x42b17180) return half*expf(fabsf(x));
|
||||
|
||||
/* |x| in [log(maxdouble), overflowthresold] */
|
||||
if (ix<=0x42b2d4fc) {
|
||||
w = expf(half*fabsf(x));
|
||||
t = half*w;
|
||||
return t*w;
|
||||
}
|
||||
|
||||
/* |x| > overflowthresold, cosh(x) overflow */
|
||||
return huge*huge;
|
||||
}
|
||||
56
liba/src/external/openbsd/e_sinhf.c
vendored
Normal file
56
liba/src/external/openbsd/e_sinhf.c
vendored
Normal file
@@ -0,0 +1,56 @@
|
||||
/* e_sinhf.c -- float version of e_sinh.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const float one = 1.0, shuge = 1.0e37;
|
||||
|
||||
float
|
||||
sinhf(float x)
|
||||
{
|
||||
float t,w,h;
|
||||
int32_t ix,jx;
|
||||
|
||||
GET_FLOAT_WORD(jx,x);
|
||||
ix = jx&0x7fffffff;
|
||||
|
||||
/* x is INF or NaN */
|
||||
if(ix>=0x7f800000) return x+x;
|
||||
|
||||
h = 0.5;
|
||||
if (jx<0) h = -h;
|
||||
/* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
|
||||
if (ix < 0x41b00000) { /* |x|<22 */
|
||||
if (ix<0x31800000) /* |x|<2**-28 */
|
||||
if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
|
||||
t = expm1f(fabsf(x));
|
||||
if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one));
|
||||
return h*(t+t/(t+one));
|
||||
}
|
||||
|
||||
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
|
||||
if (ix < 0x42b17180) return h*expf(fabsf(x));
|
||||
|
||||
/* |x| in [log(maxdouble), overflowthresold] */
|
||||
if (ix<=0x42b2d4fc) {
|
||||
w = expf((float)0.5*fabsf(x));
|
||||
t = h*w;
|
||||
return t*w;
|
||||
}
|
||||
|
||||
/* |x| > overflowthresold, sinh(x) overflow */
|
||||
return x*shuge;
|
||||
}
|
||||
45
liba/src/external/openbsd/s_asinhf.c
vendored
Normal file
45
liba/src/external/openbsd/s_asinhf.c
vendored
Normal file
@@ -0,0 +1,45 @@
|
||||
/* s_asinhf.c -- float version of s_asinh.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const float
|
||||
one = 1.0000000000e+00, /* 0x3F800000 */
|
||||
ln2 = 6.9314718246e-01, /* 0x3f317218 */
|
||||
huge= 1.0000000000e+30;
|
||||
|
||||
float
|
||||
asinhf(float x)
|
||||
{
|
||||
float t,w;
|
||||
int32_t 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;
|
||||
}
|
||||
122
liba/src/external/openbsd/s_expm1f.c
vendored
Normal file
122
liba/src/external/openbsd/s_expm1f.c
vendored
Normal file
@@ -0,0 +1,122 @@
|
||||
/* s_expm1f.c -- float version of s_expm1.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const volatile float huge = 1.0e+30, tiny = 1.0e-30;
|
||||
|
||||
static const float
|
||||
one = 1.0,
|
||||
o_threshold = 8.8721679688e+01,/* 0x42b17180 */
|
||||
ln2_hi = 6.9313812256e-01,/* 0x3f317180 */
|
||||
ln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */
|
||||
invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */
|
||||
/* scaled coefficients related to expm1 */
|
||||
Q1 = -3.3333335072e-02, /* 0xbd088889 */
|
||||
Q2 = 1.5873016091e-03, /* 0x3ad00d01 */
|
||||
Q3 = -7.9365076090e-05, /* 0xb8a670cd */
|
||||
Q4 = 4.0082177293e-06, /* 0x36867e54 */
|
||||
Q5 = -2.0109921195e-07; /* 0xb457edbb */
|
||||
|
||||
float
|
||||
expm1f(float x)
|
||||
{
|
||||
float y,hi,lo,c,t,e,hxs,hfx,r1;
|
||||
int32_t k,xsb;
|
||||
u_int32_t hx;
|
||||
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
xsb = hx&0x80000000; /* sign bit of x */
|
||||
if(xsb==0) y=x; else y= -x; /* y = |x| */
|
||||
hx &= 0x7fffffff; /* high word of |x| */
|
||||
|
||||
/* filter out huge and non-finite argument */
|
||||
if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */
|
||||
if(hx >= 0x42b17218) { /* if |x|>=88.721... */
|
||||
if(hx>0x7f800000)
|
||||
return x+x; /* NaN */
|
||||
if(hx==0x7f800000)
|
||||
return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
|
||||
if(x > o_threshold) return huge*huge; /* overflow */
|
||||
}
|
||||
if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */
|
||||
if(x+tiny<(float)0.0) /* raise inexact */
|
||||
return tiny-one; /* return -1 */
|
||||
}
|
||||
}
|
||||
|
||||
/* argument reduction */
|
||||
if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
|
||||
if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
|
||||
if(xsb==0)
|
||||
{hi = x - ln2_hi; lo = ln2_lo; k = 1;}
|
||||
else
|
||||
{hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
|
||||
} else {
|
||||
k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5);
|
||||
t = k;
|
||||
hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
|
||||
lo = t*ln2_lo;
|
||||
}
|
||||
x = hi - lo;
|
||||
c = (hi-x)-lo;
|
||||
}
|
||||
else if(hx < 0x33000000) { /* when |x|<2**-25, return x */
|
||||
t = huge+x; /* return x with inexact flags when x!=0 */
|
||||
return x - (t-(huge+x));
|
||||
}
|
||||
else k = 0;
|
||||
|
||||
/* x is now in primary range */
|
||||
hfx = (float)0.5*x;
|
||||
hxs = x*hfx;
|
||||
r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
|
||||
t = (float)3.0-r1*hfx;
|
||||
e = hxs*((r1-t)/((float)6.0 - x*t));
|
||||
if(k==0) return x - (x*e-hxs); /* c is 0 */
|
||||
else {
|
||||
e = (x*(e-c)-c);
|
||||
e -= hxs;
|
||||
if(k== -1) return (float)0.5*(x-e)-(float)0.5;
|
||||
if(k==1) {
|
||||
if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5));
|
||||
else return one+(float)2.0*(x-e);
|
||||
}
|
||||
if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
|
||||
int32_t i;
|
||||
y = one-(e-x);
|
||||
GET_FLOAT_WORD(i,y);
|
||||
SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
|
||||
return y-one;
|
||||
}
|
||||
t = one;
|
||||
if(k<23) {
|
||||
int32_t i;
|
||||
SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
|
||||
y = t-(e-x);
|
||||
GET_FLOAT_WORD(i,y);
|
||||
SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
|
||||
} else {
|
||||
int32_t i;
|
||||
SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */
|
||||
y = x-(e+t);
|
||||
y += one;
|
||||
GET_FLOAT_WORD(i,y);
|
||||
SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
|
||||
}
|
||||
}
|
||||
return y;
|
||||
}
|
||||
96
liba/src/external/openbsd/s_log1pf.c
vendored
Normal file
96
liba/src/external/openbsd/s_log1pf.c
vendored
Normal file
@@ -0,0 +1,96 @@
|
||||
/* s_log1pf.c -- float version of s_log1p.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const float
|
||||
ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
|
||||
ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
|
||||
two25 = 3.355443200e+07, /* 0x4c000000 */
|
||||
Lp1 = 6.6666668653e-01, /* 3F2AAAAB */
|
||||
Lp2 = 4.0000000596e-01, /* 3ECCCCCD */
|
||||
Lp3 = 2.8571429849e-01, /* 3E924925 */
|
||||
Lp4 = 2.2222198546e-01, /* 3E638E29 */
|
||||
Lp5 = 1.8183572590e-01, /* 3E3A3325 */
|
||||
Lp6 = 1.5313838422e-01, /* 3E1CD04F */
|
||||
Lp7 = 1.4798198640e-01; /* 3E178897 */
|
||||
|
||||
static const float zero = 0.0;
|
||||
|
||||
float
|
||||
log1pf(float x)
|
||||
{
|
||||
float hfsq,f,c,s,z,R,u;
|
||||
int32_t k,hx,hu,ax;
|
||||
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
ax = hx&0x7fffffff;
|
||||
|
||||
k = 1;
|
||||
if (hx < 0x3ed413d7) { /* x < 0.41422 */
|
||||
if(ax>=0x3f800000) { /* x <= -1.0 */
|
||||
if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
|
||||
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
|
||||
}
|
||||
if(ax<0x31000000) { /* |x| < 2**-29 */
|
||||
if(two25+x>zero /* raise inexact */
|
||||
&&ax<0x24800000) /* |x| < 2**-54 */
|
||||
return x;
|
||||
else
|
||||
return x - x*x*(float)0.5;
|
||||
}
|
||||
if(hx>0||hx<=((int32_t)0xbe95f61f)) {
|
||||
k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
|
||||
}
|
||||
if (hx >= 0x7f800000) return x+x;
|
||||
if(k!=0) {
|
||||
if(hx<0x5a000000) {
|
||||
u = (float)1.0+x;
|
||||
GET_FLOAT_WORD(hu,u);
|
||||
k = (hu>>23)-127;
|
||||
/* correction term */
|
||||
c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0);
|
||||
c /= u;
|
||||
} else {
|
||||
u = x;
|
||||
GET_FLOAT_WORD(hu,u);
|
||||
k = (hu>>23)-127;
|
||||
c = 0;
|
||||
}
|
||||
hu &= 0x007fffff;
|
||||
if(hu<0x3504f7) {
|
||||
SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */
|
||||
} else {
|
||||
k += 1;
|
||||
SET_FLOAT_WORD(u,hu|0x3f000000); /* normalize u/2 */
|
||||
hu = (0x00800000-hu)>>2;
|
||||
}
|
||||
f = u-(float)1.0;
|
||||
}
|
||||
hfsq=(float)0.5*f*f;
|
||||
if(hu==0) { /* |f| < 2**-20 */
|
||||
if(f==zero) if(k==0) return zero;
|
||||
else {c += k*ln2_lo; return k*ln2_hi+c;}
|
||||
R = hfsq*((float)1.0-(float)0.66666666666666666*f);
|
||||
if(k==0) return f-R; else
|
||||
return k*ln2_hi-((R-(k*ln2_lo+c))-f);
|
||||
}
|
||||
s = f/((float)2.0+f);
|
||||
z = s*s;
|
||||
R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
|
||||
if(k==0) return f-(hfsq-s*(hfsq+R)); else
|
||||
return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
|
||||
}
|
||||
54
liba/src/external/openbsd/s_tanhf.c
vendored
Normal file
54
liba/src/external/openbsd/s_tanhf.c
vendored
Normal file
@@ -0,0 +1,54 @@
|
||||
/* s_tanhf.c -- float version of s_tanh.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
static const float one=1.0, two=2.0, tiny = 1.0e-30;
|
||||
|
||||
float
|
||||
tanhf(float x)
|
||||
{
|
||||
float t,z;
|
||||
int32_t jx,ix;
|
||||
|
||||
GET_FLOAT_WORD(jx,x);
|
||||
ix = jx&0x7fffffff;
|
||||
|
||||
/* x is INF or NaN */
|
||||
if(ix>=0x7f800000) {
|
||||
if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */
|
||||
else return one/x-one; /* tanh(NaN) = NaN */
|
||||
}
|
||||
|
||||
/* |x| < 22 */
|
||||
if (ix < 0x41b00000) { /* |x|<22 */
|
||||
if (ix == 0)
|
||||
return x; /* x == +-0 */
|
||||
if (ix<0x24000000) /* |x|<2**-55 */
|
||||
return x*(one+x); /* tanh(small) = small */
|
||||
if (ix>=0x3f800000) { /* |x|>=1 */
|
||||
t = expm1f(two*fabsf(x));
|
||||
z = one - two/(t+two);
|
||||
} else {
|
||||
t = expm1f(-two*fabsf(x));
|
||||
z= -t/(t+two);
|
||||
}
|
||||
/* |x| > 22, return +-1 */
|
||||
} else {
|
||||
z = one - tiny; /* raised inexact flag */
|
||||
}
|
||||
return (jx>=0)? z: -z;
|
||||
}
|
||||
@@ -122,6 +122,9 @@ float Complex::b() {
|
||||
}
|
||||
|
||||
float Complex::r() const {
|
||||
if (m_b == 0) {
|
||||
return fabsf(m_a);
|
||||
}
|
||||
return sqrtf(m_a*m_a + m_b*m_b);
|
||||
}
|
||||
|
||||
|
||||
@@ -27,7 +27,7 @@ Expression * HyperbolicArcCosine::cloneWithDifferentOperands(Expression** newOpe
|
||||
float HyperbolicArcCosine::privateApproximate(Context& context, AngleUnit angleUnit) const {
|
||||
assert(angleUnit != AngleUnit::Default);
|
||||
float x = m_args[0]->approximate(context, angleUnit);
|
||||
return logf(x+sqrtf(x*x-1.0f));
|
||||
return acoshf(x);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@@ -27,7 +27,7 @@ Expression * HyperbolicArcSine::cloneWithDifferentOperands(Expression** newOpera
|
||||
float HyperbolicArcSine::privateApproximate(Context& context, AngleUnit angleUnit) const {
|
||||
assert(angleUnit != AngleUnit::Default);
|
||||
float x = m_args[0]->approximate(context, angleUnit);
|
||||
return logf(x+sqrtf(x*x+1.0f));
|
||||
return asinhf(x);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@@ -27,7 +27,7 @@ Expression * HyperbolicArcTangent::cloneWithDifferentOperands(Expression** newOp
|
||||
float HyperbolicArcTangent::privateApproximate(Context& context, AngleUnit angleUnit) const {
|
||||
assert(angleUnit != AngleUnit::Default);
|
||||
float x = m_args[0]->approximate(context, angleUnit);
|
||||
return 0.5f*logf((1.0f+x)/(1.0f-x));
|
||||
return atanhf(x);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@@ -31,7 +31,7 @@ Expression * HyperbolicCosine::cloneWithDifferentOperands(Expression** newOperan
|
||||
|
||||
float HyperbolicCosine::privateApproximate(Context& context, AngleUnit angleUnit) const {
|
||||
assert(angleUnit != AngleUnit::Default);
|
||||
return (expf(m_args[0]->approximate(context, angleUnit))+expf(-m_args[0]->approximate(context, angleUnit)))/2.0f;
|
||||
return coshf(m_args[0]->approximate(context, angleUnit));
|
||||
}
|
||||
|
||||
Expression * HyperbolicCosine::privateEvaluate(Context& context, AngleUnit angleUnit) const {
|
||||
@@ -42,6 +42,12 @@ Expression * HyperbolicCosine::privateEvaluate(Context& context, AngleUnit angle
|
||||
delete evaluation;
|
||||
return new Complex(Complex::Float(NAN));
|
||||
}
|
||||
/* Float case */
|
||||
if (((Complex *)evaluation)->b() == 0) {
|
||||
delete evaluation;
|
||||
return Function::privateEvaluate(context, angleUnit);
|
||||
}
|
||||
/* Complex case */
|
||||
Expression * arguments[2];
|
||||
arguments[0] = new Complex(Complex::Float(M_E));
|
||||
arguments[1] = evaluation;
|
||||
|
||||
@@ -31,7 +31,7 @@ Expression * HyperbolicSine::cloneWithDifferentOperands(Expression** newOperands
|
||||
|
||||
float HyperbolicSine::privateApproximate(Context& context, AngleUnit angleUnit) const {
|
||||
assert(angleUnit != AngleUnit::Default);
|
||||
return (expf(m_args[0]->approximate(context, angleUnit))-expf(-m_args[0]->approximate(context, angleUnit)))/2.0f;
|
||||
return sinhf(m_args[0]->approximate(context, angleUnit));
|
||||
}
|
||||
|
||||
Expression * HyperbolicSine::privateEvaluate(Context& context, AngleUnit angleUnit) const {
|
||||
@@ -42,6 +42,12 @@ Expression * HyperbolicSine::privateEvaluate(Context& context, AngleUnit angleUn
|
||||
delete evaluation;
|
||||
return new Complex(Complex::Float(NAN));
|
||||
}
|
||||
/* Float case */
|
||||
if (((Complex *)evaluation)->b() == 0) {
|
||||
delete evaluation;
|
||||
return Function::privateEvaluate(context, angleUnit);
|
||||
}
|
||||
/* Complex case */
|
||||
Expression * arguments[2];
|
||||
arguments[0] = new Complex(Complex::Float(M_E));
|
||||
arguments[1] = evaluation;
|
||||
|
||||
@@ -30,8 +30,7 @@ Expression * HyperbolicTangent::cloneWithDifferentOperands(Expression** newOpera
|
||||
|
||||
float HyperbolicTangent::privateApproximate(Context& context, AngleUnit angleUnit) const {
|
||||
assert(angleUnit != AngleUnit::Default);
|
||||
return (expf(m_args[0]->approximate(context, angleUnit))-expf(-m_args[0]->approximate(context, angleUnit)))/
|
||||
(expf(m_args[0]->approximate(context, angleUnit))+expf(-m_args[0]->approximate(context, angleUnit)));
|
||||
return tanhf(m_args[0]->approximate(context, angleUnit));
|
||||
}
|
||||
|
||||
Expression * HyperbolicTangent::privateEvaluate(Context& context, AngleUnit angleUnit) const {
|
||||
@@ -42,6 +41,12 @@ Expression * HyperbolicTangent::privateEvaluate(Context& context, AngleUnit angl
|
||||
delete evaluation;
|
||||
return new Complex(Complex::Float(NAN));
|
||||
}
|
||||
/* Float case */
|
||||
if (((Complex *)evaluation)->b() == 0) {
|
||||
delete evaluation;
|
||||
return Function::privateEvaluate(context, angleUnit);
|
||||
}
|
||||
/* Complex case */
|
||||
Expression * arguments[2];
|
||||
arguments[0] = new HyperbolicSine();
|
||||
((Function *)arguments[0])->setArgument(&evaluation, 1, true);
|
||||
|
||||
Reference in New Issue
Block a user