From 359d84fb9e0a22ae9d6def2bc1a002542354c50c Mon Sep 17 00:00:00 2001 From: Raymond Toy <toy.raymond@gmail.com> Date: Sat, 2 Aug 2014 08:35:39 -0700 Subject: [PATCH] Import inverse trig functions from fdlibm, as is. --- src/lisp/e_acos.c | 105 ++++++++++++++++++++++++++++++++++++ src/lisp/e_asin.c | 114 +++++++++++++++++++++++++++++++++++++++ src/lisp/s_atan.c | 134 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 353 insertions(+) create mode 100644 src/lisp/e_acos.c create mode 100644 src/lisp/e_asin.c create mode 100644 src/lisp/s_atan.c diff --git a/src/lisp/e_acos.c b/src/lisp/e_acos.c new file mode 100644 index 000000000..d7c9ed225 --- /dev/null +++ b/src/lisp/e_acos.c @@ -0,0 +1,105 @@ + +/* @(#)e_acos.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_acos(x) + * Method : + * acos(x) = pi/2 - asin(x) + * acos(-x) = pi/2 + asin(x) + * For |x|<=0.5 + * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) + * For x>0.5 + * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) + * = 2asin(sqrt((1-x)/2)) + * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) + * = 2f + (2c + 2s*z*R(z)) + * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term + * for f so that f+c ~ sqrt(z). + * For x<-0.5 + * acos(x) = pi - 2asin(sqrt((1-|x|)/2)) + * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + * Function needed: sqrt + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ +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 */ +pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ +pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ +pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ +qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ +qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ +qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ +qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + +#ifdef __STDC__ + double __ieee754_acos(double x) +#else + double __ieee754_acos(x) + double x; +#endif +{ + double z,p,q,r,w,s,c,df; + int hx,ix; + hx = __HI(x); + ix = hx&0x7fffffff; + if(ix>=0x3ff00000) { /* |x| >= 1 */ + if(((ix-0x3ff00000)|__LO(x))==0) { /* |x|==1 */ + if(hx>0) return 0.0; /* acos(1) = 0 */ + else return pi+2.0*pio2_lo; /* acos(-1)= pi */ + } + return (x-x)/(x-x); /* acos(|x|>1) is NaN */ + } + if(ix<0x3fe00000) { /* |x| < 0.5 */ + if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/ + z = x*x; + p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); + q = one+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 */ + z = (one+x)*0.5; + p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); + q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); + s = sqrt(z); + r = p/q; + w = r*s-pio2_lo; + return pi - 2.0*(s+w); + } else { /* x > 0.5 */ + z = (one-x)*0.5; + s = sqrt(z); + df = s; + __LO(df) = 0; + c = (z-df*df)/(s+df); + p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5))))); + q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4))); + r = p/q; + w = r*s+c; + return 2.0*(df+w); + } +} diff --git a/src/lisp/e_asin.c b/src/lisp/e_asin.c new file mode 100644 index 000000000..8e37e22ad --- /dev/null +++ b/src/lisp/e_asin.c @@ -0,0 +1,114 @@ + +/* @(#)e_asin.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* __ieee754_asin(x) + * Method : + * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... + * we approximate asin(x) on [0,0.5] by + * asin(x) = x + x*x^2*R(x^2) + * where + * R(x^2) is a rational approximation of (asin(x)-x)/x^3 + * and its remez error is bounded by + * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) + * + * For x in [0.5,1] + * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) + * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2; + * then for x>0.98 + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo) + * For x<=0.98, let pio4_hi = pio2_hi/2, then + * f = hi part of s; + * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z) + * and + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo) + * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + */ + + +#include "fdlibm.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +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 */ + /* coefficient for R(x^2) */ +pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ +pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ +pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ +pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ +pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ +pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ +qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ +qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ +qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ +qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + +#ifdef __STDC__ + double __ieee754_asin(double x) +#else + double __ieee754_asin(x) + double x; +#endif +{ + double t,w,p,q,c,r,s; + int hx,ix; + hx = __HI(x); + ix = hx&0x7fffffff; + if(ix>= 0x3ff00000) { /* |x|>= 1 */ + if(((ix-0x3ff00000)|__LO(x))==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<0x3e400000) { /* if |x| < 2**-27 */ + if(huge+x>one) return x;/* return x with inexact if x!=0*/ + } else + t = x*x; + p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); + q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); + w = p/q; + return x+x*w; + } + /* 1> |x|>= 0.5 */ + w = one-fabs(x); + t = w*0.5; + p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); + q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); + s = sqrt(t); + if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ + w = p/q; + t = pio2_hi-(2.0*(s+s*w)-pio2_lo); + } else { + w = s; + __LO(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); + } + if(hx>0) return t; else return -t; +} diff --git a/src/lisp/s_atan.c b/src/lisp/s_atan.c new file mode 100644 index 000000000..0093eafdf --- /dev/null +++ b/src/lisp/s_atan.c @@ -0,0 +1,134 @@ + +/* @(#)s_atan.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* atan(x) + * Method + * 1. Reduce x to positive by atan(x) = -atan(-x). + * 2. According to the integer k=4t+0.25 chopped, t=x, the argument + * is further reduced to one of the following intervals and the + * arctangent of t is evaluated by the corresponding formula: + * + * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) + * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) + * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) + * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) + * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const double atanhi[] = { +#else +static double atanhi[] = { +#endif + 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ + 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ + 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ + 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ +}; + +#ifdef __STDC__ +static const double atanlo[] = { +#else +static double atanlo[] = { +#endif + 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ + 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ + 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ + 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ +}; + +#ifdef __STDC__ +static const double aT[] = { +#else +static double aT[] = { +#endif + 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ + -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ + 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ + -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ + 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ + -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ + 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ + -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ + 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ + -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ + 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ +}; + +#ifdef __STDC__ + static const double +#else + static double +#endif +one = 1.0, +huge = 1.0e300; + +#ifdef __STDC__ + double atan(double x) +#else + double atan(x) + double x; +#endif +{ + double w,s1,s2,z; + int ix,hx,id; + + hx = __HI(x); + ix = hx&0x7fffffff; + if(ix>=0x44100000) { /* if |x| >= 2^66 */ + if(ix>0x7ff00000|| + (ix==0x7ff00000&&(__LO(x)!=0))) + return x+x; /* NaN */ + if(hx>0) return atanhi[3]+atanlo[3]; + else return -atanhi[3]-atanlo[3]; + } if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ + if (ix < 0x3e200000) { /* |x| < 2^-29 */ + if(huge+x>one) return x; /* raise inexact */ + } + id = -1; + } else { + x = fabs(x); + if (ix < 0x3ff30000) { /* |x| < 1.1875 */ + if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */ + id = 0; x = (2.0*x-one)/(2.0+x); + } else { /* 11/16<=|x|< 19/16 */ + id = 1; x = (x-one)/(x+one); + } + } else { + if (ix < 0x40038000) { /* |x| < 2.4375 */ + id = 2; x = (x-1.5)/(one+1.5*x); + } else { /* 2.4375 <= |x| < 2^66 */ + id = 3; x = -1.0/x; + } + }} + /* end of argument reduction */ + z = x*x; + w = z*z; + /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ + s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); + s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); + if (id<0) return x - x*(s1+s2); + else { + z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); + return (hx<0)? -z:z; + } +} -- GitLab