|
|
1.1 ! root 1: /* ! 2: * Copyright (c) 1985 Regents of the University of California. ! 3: * All rights reserved. ! 4: * ! 5: * Redistribution and use in source and binary forms are permitted ! 6: * provided that the above copyright notice and this paragraph are ! 7: * duplicated in all such forms and that any documentation, ! 8: * advertising materials, and other materials related to such ! 9: * distribution and use acknowledge that the software was developed ! 10: * by the University of California, Berkeley. The name of the ! 11: * University may not be used to endorse or promote products derived ! 12: * from this software without specific prior written permission. ! 13: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR ! 14: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED ! 15: * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. ! 16: * ! 17: * All recipients should regard themselves as participants in an ongoing ! 18: * research project and hence should feel obligated to report their ! 19: * experiences (good or bad) with these elementary function codes, using ! 20: * the sendbug(8) program, to the authors. ! 21: */ ! 22: ! 23: #ifndef lint ! 24: static char sccsid[] = "@(#)asincos.c 5.3 (Berkeley) 6/30/88"; ! 25: #endif /* not lint */ ! 26: ! 27: /* ASIN(X) ! 28: * RETURNS ARC SINE OF X ! 29: * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) ! 30: * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. ! 31: * ! 32: * Required system supported functions: ! 33: * copysign(x,y) ! 34: * sqrt(x) ! 35: * ! 36: * Required kernel function: ! 37: * atan2(y,x) ! 38: * ! 39: * Method : ! 40: * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is ! 41: * computed as follows ! 42: * 1-x*x if x < 0.5, ! 43: * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5. ! 44: * ! 45: * Special cases: ! 46: * if x is NaN, return x itself; ! 47: * if |x|>1, return NaN. ! 48: * ! 49: * Accuracy: ! 50: * 1) If atan2() uses machine PI, then ! 51: * ! 52: * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded; ! 53: * and PI is the exact pi rounded to machine precision (see atan2 for ! 54: * details): ! 55: * ! 56: * in decimal: ! 57: * pi = 3.141592653589793 23846264338327 ..... ! 58: * 53 bits PI = 3.141592653589793 115997963 ..... , ! 59: * 56 bits PI = 3.141592653589793 227020265 ..... , ! 60: * ! 61: * in hexadecimal: ! 62: * pi = 3.243F6A8885A308D313198A2E.... ! 63: * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps ! 64: * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps ! 65: * ! 66: * In a test run with more than 200,000 random arguments on a VAX, the ! 67: * maximum observed error in ulps (units in the last place) was ! 68: * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x))); ! 69: * ! 70: * 2) If atan2() uses true pi, then ! 71: * ! 72: * asin(x) returns the exact asin(x) with error below about 2 ulps. ! 73: * ! 74: * In a test run with more than 1,024,000 random arguments on a VAX, the ! 75: * maximum observed error in ulps (units in the last place) was ! 76: * 1.99 ulps. ! 77: */ ! 78: ! 79: double asin(x) ! 80: double x; ! 81: { ! 82: double s,t,copysign(),atan2(),sqrt(),one=1.0; ! 83: #if !defined(vax)&&!defined(tahoe) ! 84: if(x!=x) return(x); /* x is NaN */ ! 85: #endif /* !defined(vax)&&!defined(tahoe) */ ! 86: s=copysign(x,one); ! 87: if(s <= 0.5) ! 88: return(atan2(x,sqrt(one-x*x))); ! 89: else ! 90: { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); } ! 91: ! 92: } ! 93: ! 94: /* ACOS(X) ! 95: * RETURNS ARC COS OF X ! 96: * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) ! 97: * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. ! 98: * ! 99: * Required system supported functions: ! 100: * copysign(x,y) ! 101: * sqrt(x) ! 102: * ! 103: * Required kernel function: ! 104: * atan2(y,x) ! 105: * ! 106: * Method : ! 107: * ________ ! 108: * / 1 - x ! 109: * acos(x) = 2*atan2( / -------- , 1 ) . ! 110: * \/ 1 + x ! 111: * ! 112: * Special cases: ! 113: * if x is NaN, return x itself; ! 114: * if |x|>1, return NaN. ! 115: * ! 116: * Accuracy: ! 117: * 1) If atan2() uses machine PI, then ! 118: * ! 119: * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded; ! 120: * and PI is the exact pi rounded to machine precision (see atan2 for ! 121: * details): ! 122: * ! 123: * in decimal: ! 124: * pi = 3.141592653589793 23846264338327 ..... ! 125: * 53 bits PI = 3.141592653589793 115997963 ..... , ! 126: * 56 bits PI = 3.141592653589793 227020265 ..... , ! 127: * ! 128: * in hexadecimal: ! 129: * pi = 3.243F6A8885A308D313198A2E.... ! 130: * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps ! 131: * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps ! 132: * ! 133: * In a test run with more than 200,000 random arguments on a VAX, the ! 134: * maximum observed error in ulps (units in the last place) was ! 135: * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x))); ! 136: * ! 137: * 2) If atan2() uses true pi, then ! 138: * ! 139: * acos(x) returns the exact acos(x) with error below about 2 ulps. ! 140: * ! 141: * In a test run with more than 1,024,000 random arguments on a VAX, the ! 142: * maximum observed error in ulps (units in the last place) was ! 143: * 2.15 ulps. ! 144: */ ! 145: ! 146: double acos(x) ! 147: double x; ! 148: { ! 149: double t,copysign(),atan2(),sqrt(),one=1.0; ! 150: #if !defined(vax)&&!defined(tahoe) ! 151: if(x!=x) return(x); ! 152: #endif /* !defined(vax)&&!defined(tahoe) */ ! 153: if( x != -1.0) ! 154: t=atan2(sqrt((one-x)/(one+x)),one); ! 155: else ! 156: t=atan2(one,0.0); /* t = PI/2 */ ! 157: return(t+t); ! 158: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.