projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| 3cd1de16 ZAL |
1 | /* |
| 2 | * Copyright (c) 1985 Regents of the University of California. | |
| 3 | * | |
| 4 | * Use and reproduction of this software are granted in accordance with | |
| 5 | * the terms and conditions specified in the Berkeley Software License | |
| 6 | * Agreement (in particular, this entails acknowledgement of the programs' | |
| 7 | * source, and inclusion of this notice) with the additional understanding | |
| 8 | * that all recipients should regard themselves as participants in an | |
| 9 | * ongoing research project and hence should feel obligated to report | |
| 10 | * their experiences (good or bad) with these elementary function codes, | |
| 11 | * using "sendbug 4bsd-bugs@BERKELEY", to the authors. | |
| 12 | */ | |
| 13 | ||
| 14 | #ifndef lint | |
| a62df508 GK |
15 | static char sccsid[] = |
| 16 | "@(#)asincos.c 1.1 (Berkeley) 8/21/85; 1.2 (ucb.elefunt) %G%"; | |
| 3cd1de16 ZAL |
17 | #endif not lint |
| 18 | ||
| 19 | /* ASIN(X) | |
| 20 | * RETURNS ARC SINE OF X | |
| 21 | * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) | |
| 22 | * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. | |
| 23 | * | |
| 24 | * Required system supported functions: | |
| 25 | * copysign(x,y) | |
| 26 | * sqrt(x) | |
| 27 | * | |
| 28 | * Required kernel function: | |
| 29 | * atan2(y,x) | |
| 30 | * | |
| 31 | * Method : | |
| 32 | * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is | |
| 33 | * computed as follows | |
| 34 | * 1-x*x if x < 0.5, | |
| 35 | * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5. | |
| 36 | * | |
| 37 | * Special cases: | |
| 38 | * if x is NaN, return x itself; | |
| 39 | * if |x|>1, return NaN. | |
| 40 | * | |
| 41 | * Accuracy: | |
| 42 | * 1) If atan2() uses machine PI, then | |
| 43 | * | |
| 44 | * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded; | |
| 45 | * and PI is the exact pi rounded to machine precision (see atan2 for | |
| 46 | * details): | |
| 47 | * | |
| 48 | * in decimal: | |
| 49 | * pi = 3.141592653589793 23846264338327 ..... | |
| 50 | * 53 bits PI = 3.141592653589793 115997963 ..... , | |
| 51 | * 56 bits PI = 3.141592653589793 227020265 ..... , | |
| 52 | * | |
| 53 | * in hexadecimal: | |
| 54 | * pi = 3.243F6A8885A308D313198A2E.... | |
| 55 | * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps | |
| 56 | * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps | |
| 57 | * | |
| 58 | * In a test run with more than 200,000 random arguments on a VAX, the | |
| 59 | * maximum observed error in ulps (units in the last place) was | |
| 60 | * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x))); | |
| 61 | * | |
| 62 | * 2) If atan2() uses true pi, then | |
| 63 | * | |
| 64 | * asin(x) returns the exact asin(x) with error below about 2 ulps. | |
| 65 | * | |
| 66 | * In a test run with more than 1,024,000 random arguments on a VAX, the | |
| 67 | * maximum observed error in ulps (units in the last place) was | |
| 68 | * 1.99 ulps. | |
| 69 | */ | |
| 70 | ||
| 71 | double asin(x) | |
| 72 | double x; | |
| 73 | { | |
| 74 | double s,t,copysign(),atan2(),sqrt(),one=1.0; | |
| 75 | #ifndef VAX | |
| 76 | if(x!=x) return(x); /* x is NaN */ | |
| 77 | #endif | |
| 78 | s=copysign(x,one); | |
| 79 | if(s <= 0.5) | |
| 80 | return(atan2(x,sqrt(one-x*x))); | |
| 81 | else | |
| 82 | { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); } | |
| 83 | ||
| 84 | } | |
| 85 | ||
| 86 | /* ACOS(X) | |
| 87 | * RETURNS ARC COS OF X | |
| 88 | * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) | |
| 89 | * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. | |
| 90 | * | |
| 91 | * Required system supported functions: | |
| 92 | * copysign(x,y) | |
| 93 | * sqrt(x) | |
| 94 | * | |
| 95 | * Required kernel function: | |
| 96 | * atan2(y,x) | |
| 97 | * | |
| 98 | * Method : | |
| 99 | * ________ | |
| 100 | * / 1 - x | |
| 101 | * acos(x) = 2*atan2( / -------- , 1 ) . | |
| 102 | * \/ 1 + x | |
| 103 | * | |
| 104 | * Special cases: | |
| 105 | * if x is NaN, return x itself; | |
| 106 | * if |x|>1, return NaN. | |
| 107 | * | |
| 108 | * Accuracy: | |
| 109 | * 1) If atan2() uses machine PI, then | |
| 110 | * | |
| 111 | * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded; | |
| 112 | * and PI is the exact pi rounded to machine precision (see atan2 for | |
| 113 | * details): | |
| 114 | * | |
| 115 | * in decimal: | |
| 116 | * pi = 3.141592653589793 23846264338327 ..... | |
| 117 | * 53 bits PI = 3.141592653589793 115997963 ..... , | |
| 118 | * 56 bits PI = 3.141592653589793 227020265 ..... , | |
| 119 | * | |
| 120 | * in hexadecimal: | |
| 121 | * pi = 3.243F6A8885A308D313198A2E.... | |
| 122 | * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps | |
| 123 | * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps | |
| 124 | * | |
| 125 | * In a test run with more than 200,000 random arguments on a VAX, the | |
| 126 | * maximum observed error in ulps (units in the last place) was | |
| 127 | * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x))); | |
| 128 | * | |
| 129 | * 2) If atan2() uses true pi, then | |
| 130 | * | |
| 131 | * acos(x) returns the exact acos(x) with error below about 2 ulps. | |
| 132 | * | |
| 133 | * In a test run with more than 1,024,000 random arguments on a VAX, the | |
| 134 | * maximum observed error in ulps (units in the last place) was | |
| 135 | * 2.15 ulps. | |
| 136 | */ | |
| 137 | ||
| 138 | double acos(x) | |
| 139 | double x; | |
| 140 | { | |
| 141 | double t,copysign(),atan2(),sqrt(),one=1.0; | |
| 142 | #ifndef VAX | |
| 143 | if(x!=x) return(x); | |
| 144 | #endif | |
| 145 | if( x != -1.0) | |
| 146 | t=atan2(sqrt((one-x)/(one+x)),one); | |
| 147 | else | |
| 148 | t=atan2(one,0.0); /* t = PI/2 */ | |
| 149 | return(t+t); | |
| 150 | } |