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