/* NFS 20-Nov-85 14:34:56 from trig.c */ /* sin(), cos(), tan(), atan(), atan2(), asin(), acos(), sinh(), cosh() */ /* tanh(), asinh(), acosh(), atanh() */ asm (" export Libm"); mesa double LibmSupport←copysign(), LibmSupport←scalb(), LibmSupport←logb(); mesa double LibmSupport←drem(), Libm←sqrt(), Libm←expm1(), LibmSupport←expE(); mesa double Libm←exp(), Libm←log1p(); mesa int LibmSupport←finite(); #define copysign(x,y) (LibmSupport←copysign(x,y)) #define scalb(x,n) (LibmSupport←scalb(x,n)) #define logb(x) (LibmSupport←logb(x)) #define finite(x) (LibmSupport←finite(x)) #define drem(x,p) (LibmSupport←drem(x,p)) #define sqrt(x) (Libm←sqrt(x)) #define expm1(x) (Libm←expm1(x)) #define exp←←E(x,c) (LibmSupport←expE(x,c)) #define exp(x) (Libm←exp(x)) #define log1p(x) (Libm←log1p(x)) /* SIN(X), COS(X), TAN(X) * RETURN THE SINE, COSINE, AND TANGENT OF X RESPECTIVELY * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY W. Kahan and K.C. NG, 8/17/85. * * Required system supported functions: * copysign(x,y) * finite(x) * drem(x,p) * * Static kernel functions: * sin←←S(z) ....sin←←S(x*x) return (sin(x)-x)/x * cos←←C(z) ....cos←←C(x*x) return cos(x)-1-x*x/2 * * Method. * Let S and C denote the polynomial approximations to sin and cos * respectively on [-PI/4, +PI/4]. * * SIN and COS: * 1. Reduce the argument into [-PI , +PI] by the remainder function. * 2. For x in (-PI,+PI), there are three cases: * case 1: |x| < PI/4 * case 2: PI/4 <= |x| < 3PI/4 * case 3: 3PI/4 <= |x|. * SIN and COS of x are computed by: * * sin(x) cos(x) remark * ---------------------------------------------------------- * case 1 S(x) C(x) * case 2 sign(x)*C(y) S(y) y=PI/2-|x| * case 3 S(y) -C(y) y=sign(x)*(PI-|x|) * ---------------------------------------------------------- * * TAN: * 1. Reduce the argument into [-PI/2 , +PI/2] by the remainder function. * 2. For x in (-PI/2,+PI/2), there are two cases: * case 1: |x| < PI/4 * case 2: PI/4 <= |x| < PI/2 * TAN of x is computed by: * * tan (x) remark * ---------------------------------------------------------- * case 1 S(x)/C(x) * case 2 C(y)/S(y) y=sign(x)*(PI/2-|x|) * ---------------------------------------------------------- * * Notes: * 1. S(y) and C(y) were computed by: * S(y) = y+y*sin←←S(y*y) * C(y) = 1-(y*y/2-cos←←C(x*x)) ... if y*y/2 < thresh, * = 0.5-((y*y/2-0.5)-cos←←C(x*x)) ... if y*y/2 >= thresh. * where * thresh = 0.5*(acos(3/4)**2) * * 2. For better accuracy, we use the following formula for S/C for tan * (k=0): let ss=sin←←S(y*y), and cc=cos←←C(y*y), then * * y+y*ss (y*y/2-cc)+ss * S(y)/C(y) = -------- = y + y * ---------------. * C C * * * Special cases: * Let trig be any of sin, cos, or tan. * trig(+-INF) is NaN, with signals; * trig(NaN) is that NaN; * trig(n*PI/2) is exact for any integer n, provided n*PI is * representable; otherwise, trig(x) is inexact. * * Accuracy: * trig(x) returns the exact trig(x*pi/PI) nearly rounded, where * * Decimal: * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , * 56 bits PI = 3.141592653589793 227020265 ..... , * * Hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps * * In a test run with 1,024,000 random arguments on a VAX, the maximum * observed errors (compared with the exact trig(x*pi/PI)) were * tan(x) : 2.09 ulps (around 4.716340404662354) * sin(x) : .861 ulps * cos(x) : .857 ulps * * 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. */ static double thresh = 2.6117239648121182150E-1 , /*Hex 2↑ -2 * 1.0B70C6D604DD4 */ PIo4 = 7.8539816339744827900E-1 , /*Hex 2↑ -1 * 1.921FB54442D18 */ PIo2 = 1.5707963267948965580E0 , /*Hex 2↑ 0 * 1.921FB54442D18 */ PI3o4 = 2.3561944901923448370E0 , /*Hex 2↑ 1 * 1.2D97C7F3321D2 */ PI = 3.1415926535897931160E0 , /*Hex 2↑ 1 * 1.921FB54442D18 */ PI2 = 6.2831853071795862320E0 ; /*Hex 2↑ 2 * 1.921FB54442D18 */ static double zero=0, one=1, negone= -1, half=1.0/2.0, small=1E-10, /* 1+small**2==1; better values for small: small = 1.5E-9 for VAX D = 1.2E-8 for IEEE Double = 2.8E-10 for IEEE Extended */ big=1E20; /* big = 1/(small**2) */ asm (" exportproc ←tan, Libm"); double tan(x) double x; { double LibmSupport←copysign(), DoubleReal←drem(), cos←←C(),sin←←S(),a,z,ss,cc,c; int finite(),k; /* tan(NaN) and tan(INF) must be NaN */ if(!finite(x)) return(x-x); x=drem(x,PI); /* reduce x into [-PI/2, PI/2] */ a=copysign(x,one); /* ... = abs(x) */ if ( a >= PIo4 ) {k=1; x = copysign( PIo2 - a , x ); } else { k=0; if(a < small ) { big + a; return(x); }} z = x*x; cc = cos←←C(z); ss = sin←←S(z); z = z*half ; /* Next get c = cos(x) accurately */ c = (z >= thresh )? half-((z-half)-cc) : one-(z-cc); if (k==0) return ( x + (x*(z-(cc-ss)))/c ); /* sin/cos */ return( c/(x+x*ss) ); /* ... cos/sin */ } asm (" exportproc ←sin, Libm"); double sin(x) double x; { double LibmSupport←copysign(), DoubleReal←drem(), sin←←S(),cos←←C(),a,c,z; int finite(); /* sin(NaN) and sin(INF) must be NaN */ if(!finite(x)) return(x-x); x=drem(x,PI2); /* reduce x into [-PI, PI] */ a=copysign(x,one); if( a >= PIo4 ) { if( a >= PI3o4 ) /* .. in [3PI/4, PI ] */ x=copysign((a=PI-a),x); else { /* .. in [PI/4, 3PI/4] */ a=PIo2-a; /* return sign(x)*C(PI/2-|x|) */ z=a*a; c=cos←←C(z); z=z*half; a=(z>=thresh)?half-((z-half)-c):one-(z-c); return(copysign(a,x)); } } /* return S(x) */ if( a < small) { big + a; return(x);} return(x+x*sin←←S(x*x)); } asm (" exportproc ←cos, Libm"); double cos(x) double x; { double LibmSupport←copysign(), DoubleReal←drem(), sin←←S(),cos←←C(),a,c,z,s=1.0; int finite(); /* cos(NaN) and cos(INF) must be NaN */ if(!finite(x)) return(x-x); x=drem(x,PI2); /* reduce x into [-PI, PI] */ a=copysign(x,one); if ( a >= PIo4 ) { if ( a >= PI3o4 ) /* .. in [3PI/4, PI ] */ { a=PI-a; s= negone; } else /* .. in [PI/4, 3PI/4] */ /* return S(PI/2-|x|) */ { a=PIo2-a; return(a+a*sin←←S(a*a));} } /* return s*C(a) */ if( a < small) { big + a; return(s);} z=a*a; c=cos←←C(z); z=z*half; a=(z>=thresh)?half-((z-half)-c):one-(z-c); return(copysign(a,s)); } /* sin←←S(x*x) * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) * CODED IN C BY K.C. NG, 1/21/85; * REVISED BY K.C. NG on 8/13/85. * * sin(x*k) - x * RETURN --------------- on [-PI/4,PI/4] , where k=pi/PI, PI is the rounded * x * value of pi in machine precision: * * Decimal: * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , * 56 bits PI = 3.141592653589793 227020265 ..... , * * Hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 * * Method: * 1. Let z=x*x. Create a polynomial approximation to * (sin(k*x)-x)/x = z*(S0 + S1*z↑1 + ... + S5*z↑5). * Then * sin←←S(x*x) = z*(S0 + S1*z↑1 + ... + S5*z↑5) * * The coefficient S's are obtained by a special Remez algorithm. * * Accuracy: * In the absence of rounding error, the approximation has absolute error * less than 2**(-61.11) for VAX D FORMAT, 2**(-57.45) for IEEE DOUBLE. * * 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. * */ static double S0 = -1.6666666666666463126E-1 , /*Hex 2↑ -3 * -1.555555555550C */ S1 = 8.3333333332992771264E-3 , /*Hex 2↑ -7 * 1.111111110C461 */ S2 = -1.9841269816180999116E-4 , /*Hex 2↑-13 * -1.A01A019746345 */ S3 = 2.7557309793219876880E-6 , /*Hex 2↑-19 * 1.71DE3209CDCD9 */ S4 = -2.5050225177523807003E-8 , /*Hex 2↑-26 * -1.AE5C0E319A4EF */ S5 = 1.5868926979889205164E-10 ; /*Hex 2↑-33 * 1.5CF61DF672B13 */ static double sin←←S(z) double z; { return(z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5)))))); } /* cos←←C(x*x) * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) * CODED IN C BY K.C. NG, 1/21/85; * REVISED BY K.C. NG on 8/13/85. * * x*x * RETURN cos(k*x) - 1 + ----- on [-PI/4,PI/4], where k = pi/PI, * 2 * PI is the rounded value of pi in machine precision : * * Decimal: * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , * 56 bits PI = 3.141592653589793 227020265 ..... , * * Hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 * * * Method: * 1. Let z=x*x. Create a polynomial approximation to * cos(k*x)-1+z/2 = z*z*(C0 + C1*z↑1 + ... + C5*z↑5) * then * cos←←C(z) = z*z*(C0 + C1*z↑1 + ... + C5*z↑5) * * The coefficient C's are obtained by a special Remez algorithm. * * Accuracy: * In the absence of rounding error, the approximation has absolute error * less than 2**(-64) for VAX D FORMAT, 2**(-58.3) for IEEE DOUBLE. * * * 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. * */ static double C0 = 4.1666666666666504759E-2 , /*Hex 2↑ -5 * 1.555555555553E */ C1 = -1.3888888888865301516E-3 , /*Hex 2↑-10 * -1.6C16C16C14199 */ C2 = 2.4801587269650015769E-5 , /*Hex 2↑-16 * 1.A01A01971CAEB */ C3 = -2.7557304623183959811E-7 , /*Hex 2↑-22 * -1.27E4F1314AD1A */ C4 = 2.0873958177697780076E-9 , /*Hex 2↑-29 * 1.1EE3B60DDDC8C */ C5 = -1.1250289076471311557E-11 ; /*Hex 2↑-37 * -1.8BD5986B2A52E */ static double cos←←C(z) double z; { return(z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5)))))); } /* ATAN(X) * RETURNS ARC TANGENT OF X * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. * * Required kernel function: * atan2(y,x) * * Method: * atan(x) = atan2(x,1.0). * * Special case: * if x is NaN, return x itself. * * Accuracy: * 1) If atan2() uses machine PI, then * * atan(x) returns (PI/pi) * (the exact arc tangent of x) nearly rounded; * and PI is the exact pi rounded to machine precision (see atan2 for * details): * * in decimal: * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , * 56 bits PI = 3.141592653589793 227020265 ..... , * * in hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps * * In a test run with more than 200,000 random arguments on a VAX, the * maximum observed error in ulps (units in the last place) was * 0.86 ulps. (comparing against (PI/pi)*(exact atan(x))). * * 2) If atan2() uses true pi, then * * atan(x) returns the exact atan(x) with error below about 2 ulps. * * In a test run with more than 1,024,000 random arguments on a VAX, the * maximum observed error in ulps (units in the last place) was * 0.85 ulps. */ asm (" exportproc ←atan, Libm"); double atan(x) double x; { double atan2(),one=1.0; return(atan2(x,one)); } /* ATAN2(Y,X) * RETURN ARG (X+iY) * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85. * * Required system supported functions : * copysign(x,y) * scalb(x,y) * logb(x) * * Method : * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). * 2. Reduce x to positive by (if x and y are unexceptional): * ARG (x+iy) = arctan(y/x) ... if x > 0, * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument * is further reduced to one of the following intervals and the * arctangent of y/x is evaluated by the corresponding formula: * * [0,7/16] atan(y/x) = t - t↑3*(a1+t↑2*(a2+...(a10+t↑2*a11)...) * [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) ) * [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) ) * [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) ) * [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y ) * * Special cases: * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). * * ARG( NAN , (anything) ) is NaN; * ARG( (anything), NaN ) is NaN; * ARG(+(anything but NaN), +-0) is +-0 ; * ARG(-(anything but NaN), +-0) is +-PI ; * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; * ARG( -INF,+-(anything but INF and NaN) ) is +-PI; * ARG( +INF,+-INF ) is +-PI/4 ; * ARG( -INF,+-INF ) is +-3PI/4; * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; * * Accuracy: * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, * where * * in decimal: * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , * 56 bits PI = 3.141592653589793 227020265 ..... , * * in hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps * * In a test run with 356,000 random argument on [-1,1] * [-1,1] on a * VAX, the maximum observed error was 1.41 ulps (units of the last place) * compared with (PI/pi)*(the exact ARG(x+iy)). * * Note: * We use machine PI (the true pi rounded) in place of the actual * value of pi for all the trig and inverse trig functions. In general, * if trig is one of sin, cos, tan, then computed trig(y) returns the * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the * trig functions have period PI, and trig(arctrig(x)) returns x for * all critical values x. * * 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. */ static double athfhi = 4.6364760900080609352E-1 , /*Hex 2↑ -2 * 1.DAC670561BB4F */ athflo = 4.6249969567426939759E-18 , /*Hex 2↑-58 * 1.5543B8F253271 */ at1fhi = 9.8279372324732905408E-1 , /*Hex 2↑ -1 * 1.F730BD281F69B */ at1flo = -2.4407677060164810007E-17 , /*Hex 2↑-56 * -1.C23DFEFEAE6B5 */ a1 = 3.3333333333333942106E-1 , /*Hex 2↑ -2 * 1.55555555555C3 */ a2 = -1.9999999999979536924E-1 , /*Hex 2↑ -3 * -1.9999999997CCD */ a3 = 1.4285714278004377209E-1 , /*Hex 2↑ -3 * 1.24924921EC1D7 */ a4 = -1.1111110579344973814E-1 , /*Hex 2↑ -4 * -1.C71C7059AF280 */ a5 = 9.0908906105474668324E-2 , /*Hex 2↑ -4 * 1.745CE5AA35DB2 */ a6 = -7.6919217767468239799E-2 , /*Hex 2↑ -4 * -1.3B0FA54BEC400 */ a7 = 6.6614695906082474486E-2 , /*Hex 2↑ -4 * 1.10DA924597FFF */ a8 = -5.8358371008508623523E-2 , /*Hex 2↑ -5 * -1.DE125FDDBD793 */ a9 = 4.9850617156082015213E-2 , /*Hex 2↑ -5 * 1.9860524BDD807 */ a10 = -3.6700606902093604877E-2 , /*Hex 2↑ -5 * -1.2CA6C04C6937A */ a11 = 1.6438029044759730479E-2 ; /*Hex 2↑ -6 * 1.0D52174A1BB54 */ asm (" exportproc ←atan2, Libm"); double atan2(y,x) double y,x; { static double zero=0, one=1, small=1.0E-9, big=1.0E18; double LibmSupport←copysign(), LibmSupport←logb(), LibmSupport←scalb(),t,z,signy,signx,hi,lo; int LibmSupport←finite(), k,m; /* if x or y is NAN */ if(x!=x) return(x); if(y!=y) return(y); /* copy down the sign of y and x */ signy = copysign(one,y) ; signx = copysign(one,x) ; /* if x is 1.0, goto begin */ if(x==1) { y=copysign(y,one); t=y; if(finite(t)) goto begin;} /* when y = 0 */ if(y==zero) return((signx==one)?y:copysign(PI,signy)); /* when x = 0 */ if(x==zero) return(copysign(PIo2,signy)); /* when x is INF */ if(!finite(x)) if(!finite(y)) return(copysign((signx==one)?PIo4:3*PIo4,signy)); else return(copysign((signx==one)?zero:PI,signy)); /* when y is INF */ if(!finite(y)) return(copysign(PIo2,signy)); /* compute y/x */ x=copysign(x,one); y=copysign(y,one); if((m=(k=logb(y))-logb(x)) > 60) t=big+big; else if(m < -80 ) t=y/x; else { t = y/x ; y = scalb(y,-k); x=scalb(x,-k); } /* begin argument reduction */ begin: if (t < 2.4375) { /* truncate 4(t+1/16) to integer for branching */ k = 4 * (t+0.0625); switch (k) { /* t is in [0,7/16] */ case 0: case 1: if (t < small) { big + small ; /* raise inexact flag */ return (copysign((signx>zero)?t:PI-t,signy)); } hi = zero; lo = zero; break; /* t is in [7/16,11/16] */ case 2: hi = athfhi; lo = athflo; z = x+x; t = ( (y+y) - x ) / ( z + y ); break; /* t is in [11/16,19/16] */ case 3: case 4: hi = PIo4; lo = zero; t = ( y - x ) / ( x + y ); break; /* t is in [19/16,39/16] */ default: hi = at1fhi; lo = at1flo; z = y-x; y=y+y+y; t = x+x; t = ( (z+z)-x ) / ( t + y ); break; } } /* end of if (t < 2.4375) */ else { hi = PIo2; lo = zero; /* t is in [2.4375, big] */ if (t <= big) t = - x / y; /* t is in [big, INF] */ else { big+small; /* raise inexact flag */ t = zero; } } /* end of argument reduction */ /* compute atan(t) for t in [-.4375, .4375] */ z = t*t; z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ z*(a9+z*(a10+z*a11))))))))))); z = lo - z; z += t; z += hi; return(copysign((signx>zero)?z:PI-z,signy)); } /* ASIN(X) * RETURNS ARC SINE OF X * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. * * Required system supported functions: * copysign(x,y) * sqrt(x) * * Required kernel function: * atan2(y,x) * * Method : * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is * computed as follows * 1-x*x if x < 0.5, * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5. * * Special cases: * if x is NaN, return x itself; * if |x|>1, return NaN. * * Accuracy: * 1) If atan2() uses machine PI, then * * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded; * and PI is the exact pi rounded to machine precision (see atan2 for * details): * * in decimal: * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , * 56 bits PI = 3.141592653589793 227020265 ..... , * * in hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps * * In a test run with more than 200,000 random arguments on a VAX, the * maximum observed error in ulps (units in the last place) was * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x))); * * 2) If atan2() uses true pi, then * * asin(x) returns the exact asin(x) with error below about 2 ulps. * * In a test run with more than 1,024,000 random arguments on a VAX, the * maximum observed error in ulps (units in the last place) was * 1.99 ulps. */ asm (" exportproc ←asin, Libm"); double asin(x) double x; { double s,t, LibmSupport←copysign(), Libm←atan2(), sqrt(), one=1.0; if(x!=x) return(x); /* x is NaN */ s=copysign(x,one); if(s <= 0.5) return(atan2(x,sqrt(one-x*x))); else { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); } } /* ACOS(X) * RETURNS ARC COS OF X * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. * * Required system supported functions: * copysign(x,y) * sqrt(x) * * Required kernel function: * atan2(y,x) * * Method : * ←←←←←←←← * / 1 - x * acos(x) = 2*atan2( / -------- , 1 ) . * \/ 1 + x * * Special cases: * if x is NaN, return x itself; * if |x|>1, return NaN. * * Accuracy: * 1) If atan2() uses machine PI, then * * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded; * and PI is the exact pi rounded to machine precision (see atan2 for * details): * * in decimal: * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , * 56 bits PI = 3.141592653589793 227020265 ..... , * * in hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps * * In a test run with more than 200,000 random arguments on a VAX, the * maximum observed error in ulps (units in the last place) was * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x))); * * 2) If atan2() uses true pi, then * * acos(x) returns the exact acos(x) with error below about 2 ulps. * * In a test run with more than 1,024,000 random arguments on a VAX, the * maximum observed error in ulps (units in the last place) was * 2.15 ulps. */ asm (" exportproc ←acos, Libm"); double acos(x) double x; { double t, LibmSupport←copysign(), Libm←atan2(), sqrt(), one=1.0; if(x!=x) return(x); if( x != -1.0) t=atan2(sqrt((one-x)/(one+x)),one); else t=atan2(one,0.0); /* t = PI/2 */ return(t+t); } /* SINH(X) * RETURN THE HYPERBOLIC SINE OF X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85. * * Required system supported functions : * copysign(x,y) * scalb(x,N) * * Required kernel functions: * expm1(x) ...return exp(x)-1 * * Method : * 1. reduce x to non-negative by sinh(-x) = - sinh(x). * 2. * * expm1(x) + expm1(x)/(expm1(x)+1) * 0 <= x <= lnovfl : sinh(x) := -------------------------------- * 2 * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow) * lnovfl+ln2 < x < INF : overflow to INF * * * Special cases: * sinh(x) is x if x is +INF, -INF, or NaN. * only sinh(0)=0 is exact for finite argument. * * Accuracy: * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In * a test run with 1,024,000 random arguments on a VAX, the maximum * observed error was 1.93 ulps (units in the last place). * * 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. */ double static mln2hi = 7.0978271289338397310E2 , /*Hex 2↑ 10 * 1.62E42FEFA39EF */ mln2lo = 2.3747039373786107478E-14 , /*Hex 2↑-45 * 1.ABC9E3B39803F */ lnovfl = 7.0978271289338397310E2 ; /*Hex 2↑ 9 * 1.62E42FEFA39EF */ static max = 1023 ; asm (" exportproc ←sinh, Libm"); double sinh(x) double x; { static double one=1.0, half=1.0/2.0 ; double expm1(), t, LibmSupport←scalb(), LibmSupport←copysign(), sign; if(x!=x) return(x); /* x is NaN */ sign=copysign(one,x); x=copysign(x,one); if(x<lnovfl) {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));} else if(x <= lnovfl+0.7) /* subtract x by ln(2↑(max+1)) and return 2↑max*exp(x) to avoid unnecessary overflow */ return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign)); else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */ return( expm1(x)*sign ); } /* TANH(X) * RETURN THE HYPERBOLIC TANGENT OF X * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85. * * Required system supported functions : * copysign(x,y) * finite(x) * * Required kernel function: * expm1(x) ...exp(x)-1 * * Method : * 1. reduce x to non-negative by tanh(-x) = - tanh(x). * 2. * 0 < x <= 1.e-10 : tanh(x) := x * -expm1(-2x) * 1.e-10 < x <= 1 : tanh(x) := -------------- * expm1(-2x) + 2 * 2 * 1 <= x <= 22.0 : tanh(x) := 1 - --------------- * expm1(2x) + 2 * 22.0 < x <= INF : tanh(x) := 1. * * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1. * * Special cases: * tanh(NaN) is NaN; * only tanh(0)=0 is exact for finite argument. * * Accuracy: * tanh(x) returns the exact hyperbolic tangent of x nealy rounded. * In a test run with 1,024,000 random arguments on a VAX, the maximum * observed error was 2.22 ulps (units in the last place). */ asm (" exportproc ←tanh, Libm"); double tanh(x) double x; { static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10; double expm1(), t, LibmSupport←copysign(), sign; int finite(); if(x!=x) return(x); /* x is NaN */ sign=copysign(one,x); x=copysign(x,one); if(x < 22.0) if( x > one ) return(copysign(one-two/(expm1(x+x)+two),sign)); else if ( x > small ) {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));} else /* raise the INEXACT flag for non-zero x */ {big+x; return(copysign(x,sign));} else if(finite(x)) return (sign+1.0E-37); /* raise the INEXACT flag */ else return(sign); /* x is +- INF */ } /* COSH(X) * RETURN THE HYPERBOLIC COSINE OF X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85. * * Required system supported functions : * copysign(x,y) * scalb(x,N) * * Required kernel function: * exp(x) * exp←←E(x,c) ...return exp(x+c)-1-x for |x|<0.3465 * * Method : * 1. Replace x by |x|. * 2. * [ exp(x) - 1 ]↑2 * 0 <= x <= 0.3465 : cosh(x) := 1 + ------------------- * 2*exp(x) * * exp(x) + 1/exp(x) * 0.3465 <= x <= 22 : cosh(x) := ------------------- * 2 * 22 <= x <= lnovfl : cosh(x) := exp(x)/2 * lnovfl <= x <= lnovfl+log(2) * : cosh(x) := exp(x)/2 (avoid overflow) * log(2)+lnovfl < x < INF: overflow to INF * * Note: .3465 is a number near one half of ln2. * * Special cases: * cosh(x) is x if x is +INF, -INF, or NaN. * only cosh(0)=1 is exact for finite x. * * Accuracy: * cosh(x) returns the exact hyperbolic cosine of x nearly rounded. * In a test run with 768,000 random arguments on a VAX, the maximum * observed error was 1.23 ulps (units in the last place). * * 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. */ asm (" exportproc ←cosh, Libm"); double cosh(x) double x; { static double half=1.0/2.0,one=1.0, small=1.0E-18; /* fl(1+small)==1 */ double Libmsupport←scalb(),LibmSupport←copysign(),exp(), LibmSupport←exp←←E(),t; if(x!=x) return(x); /* x is NaN */ if((x=copysign(x,one)) <= 22) if(x<0.3465) if(x<small) return(one+x); else {t=x+exp←←E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); } else /* for x lies in [0.3465,22] */ { t=exp(x); return((t+one/t)*half); } if( lnovfl <= x && x <= (lnovfl+0.7)) /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2↑(max+1)) * and return 2↑max*exp(x) to avoid unnecessary overflow */ return(scalb(exp((x-mln2hi)-mln2lo), max)); else return(exp(x)*half); /* for large x, cosh(x)=exp(x)/2 */ } /* ASINH(X) * RETURN THE INVERSE HYPERBOLIC SINE OF X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 2/16/85; * REVISED BY K.C. NG on 3/7/85, 3/24/85, 4/16/85. * * Required system supported functions : * copysign(x,y) * sqrt(x) * * Required kernel function: * log1p(x) ...return log(1+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)*(log1p(x)+ln2)) if sqrt(1+x*x)=x, else * := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)↑2)) ) * * Accuracy: * asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded. * In a test run with 52,000 random arguments on a VAX, the maximum * observed error was 1.58 ulps (units in the last place). * * 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. */ static double ln2hi = 6.9314718036912381649E-1 , /*Hex 2↑ -1 * 1.62E42FEE00000 */ ln2lo = 1.9082149292705877000E-10 ; /*Hex 2↑-33 * 1.A39EF35793C76 */ asm (" exportproc ←asinh, Libm"); double asinh(x) double x; { double LibmSupport←copysign(),log1p(),Libm←sqrt(),t,s; static double small=1.0E-10, /* fl(1+small*small) == 1 */ big =1.0E20, /* fl(1+big) == big */ one =1.0 ; if(x!=x) return(x); /* x is NaN */ if((t=copysign(x,one))>small) if(t<big) { s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); } else /* if |x| > big */ {s=log1p(t)+ln2lo; return(copysign(s+ln2hi,x));} else /* if |x| < small */ return(x); } /* ACOSH(X) * RETURN THE INVERSE HYPERBOLIC COSINE OF X * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 2/16/85; * REVISED BY K.C. NG on 3/6/85, 3/24/85, 4/16/85, 8/17/85. * * Required system supported functions : * sqrt(x) * * Required kernel function: * log1p(x) ...return log(1+x) * * Method : * Based on * acosh(x) = log [ x + sqrt(x*x-1) ] * we have * acosh(x) := log1p(x)+ln2, if (x > 1.0E20); else * acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) . * These formulae avoid the over/underflow complication. * * Special cases: * acosh(x) is NaN with signal if x<1. * acosh(NaN) is NaN without signal. * * Accuracy: * acosh(x) returns the exact inverse hyperbolic cosine of x nearly * rounded. In a test run with 512,000 random arguments on a VAX, the * maximum observed error was 3.30 ulps (units of the last place) at * x=1.0070493753568216 . * * 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. */ asm (" exportproc ←acosh, Libm"); double acosh(x) double x; { double log1p(),sqrt(),t,big=1.E20; /* big+1==big */ if(x!=x) return(x); /* x is NaN */ /* return log1p(x) + log(2) if x is large */ if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);} t=sqrt(x-1.0); return(log1p(t*(t+sqrt(x+1.0)))); } /* ATANH(X) * RETURN THE HYPERBOLIC ARC TANGENT OF X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85. * * Required kernel function: * log1p(x) ...return log(1+x) * * Method : * Return * 1 2x x * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) * 2 1 - 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. * * Accuracy: * atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded. * In a test run with 512,000 random arguments on a VAX, the maximum * observed error was 1.87 ulps (units in the last place) at * x= -3.8962076028810414000e-03. */ asm (" exportproc ←atanh, Libm"); double atanh(x) double x; { double LibmSupport←copysign(),log1p(),z; z = copysign(0.5,x); x = copysign(x,1.0); x = x/(1.0-x); return( z*log1p(x+x) ); }