projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| 7fde0e2e KB |
1 | /*- |
| 2 | * Copyright (c) 1990 The Regents of the University of California. | |
| 3 | * All rights reserved. | |
| 4 | * | |
| 5 | * This code is derived from software contributed to Berkeley by | |
| 6 | * the Systems Programming Group of the University of Utah Computer | |
| 7 | * Science Department. | |
| 8 | * | |
| 9 | * %sccs.include.redist.c% | |
| 10 | */ | |
| 11 | ||
| 12 | #ifndef lint | |
| 13 | static char sccsid[] = "@(#)atan2.c 5.1 (Berkeley) %G%"; | |
| 14 | #endif /* not lint */ | |
| 15 | ||
| 16 | /* | |
| 17 | * ATAN2(Y,X) | |
| 18 | * RETURN ARG (X+iY) | |
| 19 | * DOUBLE PRECISION (IEEE DOUBLE 53 BITS) | |
| 20 | * | |
| 21 | * Scaled down version to weed out special cases. "Normal" cases are | |
| 22 | * handled by calling atan2__A(), an assembly coded support routine in | |
| 23 | * support.s. | |
| 24 | * | |
| 25 | * Required system supported functions : | |
| 26 | * copysign(x,y) | |
| 27 | * atan2__A(y,x) | |
| 28 | * | |
| 29 | * Method : | |
| 30 | * 1. Deal with special cases | |
| 31 | * 2. Call atan2__A() to do the others | |
| 32 | * | |
| 33 | * Special cases: | |
| 34 | * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). | |
| 35 | * | |
| 36 | * ARG( NAN , (anything) ) is NaN; | |
| 37 | * ARG( (anything), NaN ) is NaN; | |
| 38 | * ARG(+(anything but NaN), +-0) is +-0 ; | |
| 39 | * ARG(-(anything but NaN), +-0) is +-PI ; | |
| 40 | * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; | |
| 41 | * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; | |
| 42 | * ARG( -INF,+-(anything but INF and NaN) ) is +-PI; | |
| 43 | * ARG( +INF,+-INF ) is +-PI/4 ; | |
| 44 | * ARG( -INF,+-INF ) is +-3PI/4; | |
| 45 | * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; | |
| 46 | * | |
| 47 | * Accuracy: | |
| 48 | * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, | |
| 49 | * where | |
| 50 | * | |
| 51 | * in decimal: | |
| 52 | * pi = 3.141592653589793 23846264338327 ..... | |
| 53 | * 53 bits PI = 3.141592653589793 115997963 ..... , | |
| 54 | * 56 bits PI = 3.141592653589793 227020265 ..... , | |
| 55 | * | |
| 56 | * in hexadecimal: | |
| 57 | * pi = 3.243F6A8885A308D313198A2E.... | |
| 58 | * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps | |
| 59 | * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps | |
| 60 | * | |
| 61 | * In a test run with 356,000 random argument on [-1,1] * [-1,1] on a | |
| 62 | * VAX, the maximum observed error was 1.41 ulps (units of the last place) | |
| 63 | * compared with (PI/pi)*(the exact ARG(x+iy)). | |
| 64 | * | |
| 65 | * Note: | |
| 66 | * We use machine PI (the true pi rounded) in place of the actual | |
| 67 | * value of pi for all the trig and inverse trig functions. In general, | |
| 68 | * if trig is one of sin, cos, tan, then computed trig(y) returns the | |
| 69 | * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig | |
| 70 | * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the | |
| 71 | * trig functions have period PI, and trig(arctrig(x)) returns x for | |
| 72 | * all critical values x. | |
| 73 | * | |
| 74 | * Constants: | |
| 75 | * The hexadecimal values are the intended ones for the following constants. | |
| 76 | * The decimal values may be used, provided that the compiler will convert | |
| 77 | * from decimal to binary accurately enough to produce the hexadecimal values | |
| 78 | * shown. | |
| 79 | */ | |
| 80 | ||
| 81 | static double | |
| 82 | PIo4 = 7.8539816339744827900E-1 , /*Hex 2^ -1 * 1.921FB54442D18 */ | |
| 83 | PIo2 = 1.5707963267948965580E0 , /*Hex 2^ 0 * 1.921FB54442D18 */ | |
| 84 | PI = 3.1415926535897931160E0 ; /*Hex 2^ 1 * 1.921FB54442D18 */ | |
| 85 | ||
| 86 | double atan2(y,x) | |
| 87 | double y,x; | |
| 88 | { | |
| 89 | static double zero=0, one=1; | |
| 90 | double copysign(),atan2__A(),signy,signx; | |
| 91 | int finite(); | |
| 92 | ||
| 93 | /* if x or y is NAN */ | |
| 94 | if(x!=x) return(x); if(y!=y) return(y); | |
| 95 | ||
| 96 | /* copy down the sign of y and x */ | |
| 97 | signy = copysign(one,y); | |
| 98 | signx = copysign(one,x); | |
| 99 | ||
| 100 | /* when y = 0 */ | |
| 101 | if(y==zero) return((signx==one)?y:copysign(PI,signy)); | |
| 102 | ||
| 103 | /* when x = 0 */ | |
| 104 | if(x==zero) return(copysign(PIo2,signy)); | |
| 105 | ||
| 106 | /* when x is INF */ | |
| 107 | if(!finite(x)) | |
| 108 | if(!finite(y)) | |
| 109 | return(copysign((signx==one)?PIo4:3*PIo4,signy)); | |
| 110 | else | |
| 111 | return(copysign((signx==one)?zero:PI,signy)); | |
| 112 | ||
| 113 | /* when y is INF */ | |
| 114 | if(!finite(y)) return(copysign(PIo2,signy)); | |
| 115 | ||
| 116 | /* else let atan2__A do the work */ | |
| 117 | return(atan2__A(y,x)); | |
| 118 | } |