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