1 /* 2 * Copyright (c) 1985, 1993 3 * The Regents of the University of California. All rights reserved. 4 * 5 * %sccs.include.redist.c% 6 */ 7 8 #ifndef lint 9 static char sccsid[] = "@(#)asincos.c 8.1 (Berkeley) 06/04/93"; 10 #endif /* not lint */ 11 12 /* ASIN(X) 13 * RETURNS ARC SINE OF X 14 * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) 15 * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. 16 * 17 * Required system supported functions: 18 * copysign(x,y) 19 * sqrt(x) 20 * 21 * Required kernel function: 22 * atan2(y,x) 23 * 24 * Method : 25 * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is 26 * computed as follows 27 * 1-x*x if x < 0.5, 28 * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5. 29 * 30 * Special cases: 31 * if x is NaN, return x itself; 32 * if |x|>1, return NaN. 33 * 34 * Accuracy: 35 * 1) If atan2() uses machine PI, then 36 * 37 * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded; 38 * and PI is the exact pi rounded to machine precision (see atan2 for 39 * details): 40 * 41 * in decimal: 42 * pi = 3.141592653589793 23846264338327 ..... 43 * 53 bits PI = 3.141592653589793 115997963 ..... , 44 * 56 bits PI = 3.141592653589793 227020265 ..... , 45 * 46 * in hexadecimal: 47 * pi = 3.243F6A8885A308D313198A2E.... 48 * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps 49 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps 50 * 51 * In a test run with more than 200,000 random arguments on a VAX, the 52 * maximum observed error in ulps (units in the last place) was 53 * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x))); 54 * 55 * 2) If atan2() uses true pi, then 56 * 57 * asin(x) returns the exact asin(x) with error below about 2 ulps. 58 * 59 * In a test run with more than 1,024,000 random arguments on a VAX, the 60 * maximum observed error in ulps (units in the last place) was 61 * 1.99 ulps. 62 */ 63 64 double asin(x) 65 double x; 66 { 67 double s,t,copysign(),atan2(),sqrt(),one=1.0; 68 #if !defined(vax)&&!defined(tahoe) 69 if(x!=x) return(x); /* x is NaN */ 70 #endif /* !defined(vax)&&!defined(tahoe) */ 71 s=copysign(x,one); 72 if(s <= 0.5) 73 return(atan2(x,sqrt(one-x*x))); 74 else 75 { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); } 76 77 } 78 79 /* ACOS(X) 80 * RETURNS ARC COS OF X 81 * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) 82 * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. 83 * 84 * Required system supported functions: 85 * copysign(x,y) 86 * sqrt(x) 87 * 88 * Required kernel function: 89 * atan2(y,x) 90 * 91 * Method : 92 * ________ 93 * / 1 - x 94 * acos(x) = 2*atan2( / -------- , 1 ) . 95 * \/ 1 + x 96 * 97 * Special cases: 98 * if x is NaN, return x itself; 99 * if |x|>1, return NaN. 100 * 101 * Accuracy: 102 * 1) If atan2() uses machine PI, then 103 * 104 * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded; 105 * and PI is the exact pi rounded to machine precision (see atan2 for 106 * details): 107 * 108 * in decimal: 109 * pi = 3.141592653589793 23846264338327 ..... 110 * 53 bits PI = 3.141592653589793 115997963 ..... , 111 * 56 bits PI = 3.141592653589793 227020265 ..... , 112 * 113 * in hexadecimal: 114 * pi = 3.243F6A8885A308D313198A2E.... 115 * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps 116 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps 117 * 118 * In a test run with more than 200,000 random arguments on a VAX, the 119 * maximum observed error in ulps (units in the last place) was 120 * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x))); 121 * 122 * 2) If atan2() uses true pi, then 123 * 124 * acos(x) returns the exact acos(x) with error below about 2 ulps. 125 * 126 * In a test run with more than 1,024,000 random arguments on a VAX, the 127 * maximum observed error in ulps (units in the last place) was 128 * 2.15 ulps. 129 */ 130 131 double acos(x) 132 double x; 133 { 134 double t,copysign(),atan2(),sqrt(),one=1.0; 135 #if !defined(vax)&&!defined(tahoe) 136 if(x!=x) return(x); 137 #endif /* !defined(vax)&&!defined(tahoe) */ 138 if( x != -1.0) 139 t=atan2(sqrt((one-x)/(one+x)),one); 140 else 141 t=atan2(one,0.0); /* t = PI/2 */ 142 return(t+t); 143 } 144