projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| f9fea09f 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 | |
| 27c51c7b | 15 | static char sccsid[] = |
| f77d20bd | 16 | "@(#)cabs.c 1.2 (Berkeley) 8/21/85; 1.4 (ucb.elefunt) %G%"; |
| f9fea09f ZAL |
17 | #endif not lint |
| 18 | ||
| 19 | /* CABS(Z) | |
| 20 | * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY | |
| 21 | * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) | |
| 22 | * CODED IN C BY K.C. NG, 11/28/84. | |
| 23 | * REVISED BY K.C. NG, 7/12/85. | |
| 24 | * | |
| 25 | * Required kernel function : | |
| 26 | * hypot(x,y) | |
| 27 | * | |
| 28 | * Method : | |
| 29 | * cabs(z) = hypot(x,y) . | |
| 30 | */ | |
| 31 | ||
| 32 | double cabs(z) | |
| 33 | struct { double x, y;} z; | |
| 34 | { | |
| 35 | double hypot(); | |
| 36 | return(hypot(z.x,z.y)); | |
| 37 | } | |
| 38 | ||
| 39 | ||
| 40 | /* HYPOT(X,Y) | |
| 41 | * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY | |
| 42 | * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) | |
| 43 | * CODED IN C BY K.C. NG, 11/28/84; | |
| 44 | * REVISED BY K.C. NG, 7/12/85. | |
| 45 | * | |
| 46 | * Required system supported functions : | |
| 47 | * copysign(x,y) | |
| 48 | * finite(x) | |
| 49 | * scalb(x,N) | |
| 50 | * sqrt(x) | |
| 51 | * | |
| 52 | * Method : | |
| 53 | * 1. replace x by |x| and y by |y|, and swap x and | |
| 54 | * y if y > x (hence x is never smaller than y). | |
| 55 | * 2. Hypot(x,y) is computed by: | |
| 56 | * Case I, x/y > 2 | |
| 57 | * | |
| 58 | * y | |
| 59 | * hypot = x + ----------------------------- | |
| 60 | * 2 | |
| 61 | * sqrt ( 1 + [x/y] ) + x/y | |
| 62 | * | |
| 63 | * Case II, x/y <= 2 | |
| 64 | * y | |
| 65 | * hypot = x + -------------------------------------------------- | |
| 66 | * 2 | |
| 67 | * [x/y] - 2 | |
| 68 | * (sqrt(2)+1) + (x-y)/y + ----------------------------- | |
| 69 | * 2 | |
| 70 | * sqrt ( 1 + [x/y] ) + sqrt(2) | |
| 71 | * | |
| 72 | * | |
| 73 | * | |
| 74 | * Special cases: | |
| 75 | * hypot(x,y) is INF if x or y is +INF or -INF; else | |
| 76 | * hypot(x,y) is NAN if x or y is NAN. | |
| 77 | * | |
| 78 | * Accuracy: | |
| 79 | * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units | |
| 80 | * in the last place). See Kahan's "Interval Arithmetic Options in the | |
| 81 | * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics | |
| 82 | * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate | |
| 83 | * code follows in comments.) In a test run with 500,000 random arguments | |
| 84 | * on a VAX, the maximum observed error was .959 ulps. | |
| 85 | * | |
| 86 | * Constants: | |
| 87 | * The hexadecimal values are the intended ones for the following constants. | |
| 88 | * The decimal values may be used, provided that the compiler will convert | |
| 89 | * from decimal to binary accurately enough to produce the hexadecimal values | |
| 90 | * shown. | |
| 91 | */ | |
| 92 | ||
| e0085737 | 93 | #if (defined(VAX)||defined(TAHOE)) /* VAX D format */ |
| f77d20bd ZAL |
94 | #ifdef VAX |
| 95 | #define _0x(A,B) 0x/**/A/**/B | |
| 96 | #else /* VAX */ | |
| 97 | #define _0x(A,B) 0x/**/B/**/A | |
| 98 | #endif /* VAX */ | |
| f9fea09f ZAL |
99 | /* static double */ |
| 100 | /* r2p1hi = 2.4142135623730950345E0 , Hex 2^ 2 * .9A827999FCEF32 */ | |
| 101 | /* r2p1lo = 1.4349369327986523769E-17 , Hex 2^-55 * .84597D89B3754B */ | |
| 102 | /* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */ | |
| f77d20bd ZAL |
103 | static long r2p1hix[] = { _0x(8279,411a), _0x(ef32,99fc)}; |
| 104 | static long r2p1lox[] = { _0x(597d,2484), _0x(754b,89b3)}; | |
| 105 | static long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)}; | |
| f9fea09f ZAL |
106 | #define r2p1hi (*(double*)r2p1hix) |
| 107 | #define r2p1lo (*(double*)r2p1lox) | |
| 108 | #define sqrt2 (*(double*)sqrt2x) | |
| 109 | #else /* IEEE double format */ | |
| 110 | static double | |
| 111 | r2p1hi = 2.4142135623730949234E0 , /*Hex 2^1 * 1.3504F333F9DE6 */ | |
| 112 | r2p1lo = 1.2537167179050217666E-16 , /*Hex 2^-53 * 1.21165F626CDD5 */ | |
| 113 | sqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */ | |
| 114 | #endif | |
| 115 | ||
| 116 | double hypot(x,y) | |
| 117 | double x, y; | |
| 118 | { | |
| 119 | static double zero=0, one=1, | |
| 120 | small=1.0E-18; /* fl(1+small)==1 */ | |
| 121 | static ibig=30; /* fl(1+2**(2*ibig))==1 */ | |
| 122 | double copysign(),scalb(),logb(),sqrt(),t,r; | |
| 123 | int finite(), exp; | |
| 124 | ||
| 125 | if(finite(x)) | |
| 126 | if(finite(y)) | |
| 127 | { | |
| 128 | x=copysign(x,one); | |
| 129 | y=copysign(y,one); | |
| 130 | if(y > x) | |
| 131 | { t=x; x=y; y=t; } | |
| 132 | if(x == zero) return(zero); | |
| 133 | if(y == zero) return(x); | |
| 134 | exp= logb(x); | |
| 135 | if(exp-(int)logb(y) > ibig ) | |
| 136 | /* raise inexact flag and return |x| */ | |
| 137 | { one+small; return(x); } | |
| 138 | ||
| 139 | /* start computing sqrt(x^2 + y^2) */ | |
| 140 | r=x-y; | |
| 141 | if(r>y) { /* x/y > 2 */ | |
| 142 | r=x/y; | |
| 143 | r=r+sqrt(one+r*r); } | |
| 144 | else { /* 1 <= x/y <= 2 */ | |
| 145 | r/=y; t=r*(r+2.0); | |
| 146 | r+=t/(sqrt2+sqrt(2.0+t)); | |
| 147 | r+=r2p1lo; r+=r2p1hi; } | |
| 148 | ||
| 149 | r=y/r; | |
| 150 | return(x+r); | |
| 151 | ||
| 152 | } | |
| 153 | ||
| 154 | else if(y==y) /* y is +-INF */ | |
| 155 | return(copysign(y,one)); | |
| 156 | else | |
| 157 | return(y); /* y is NaN and x is finite */ | |
| 158 | ||
| 159 | else if(x==x) /* x is +-INF */ | |
| 160 | return (copysign(x,one)); | |
| 161 | else if(finite(y)) | |
| 162 | return(x); /* x is NaN, y is finite */ | |
| 163 | else if(y!=y) return(y); /* x and y is NaN */ | |
| 164 | else return(copysign(y,one)); /* y is INF */ | |
| 165 | } | |
| 166 | ||
| 167 | /* A faster but less accurate version of cabs(x,y) */ | |
| 168 | #if 0 | |
| 169 | double hypot(x,y) | |
| 170 | double x, y; | |
| 171 | { | |
| 172 | static double zero=0, one=1; | |
| 173 | small=1.0E-18; /* fl(1+small)==1 */ | |
| 174 | static ibig=30; /* fl(1+2**(2*ibig))==1 */ | |
| 175 | double copysign(),scalb(),logb(),sqrt(),temp; | |
| 176 | int finite(), exp; | |
| 177 | ||
| 178 | if(finite(x)) | |
| 179 | if(finite(y)) | |
| 180 | { | |
| 181 | x=copysign(x,one); | |
| 182 | y=copysign(y,one); | |
| 183 | if(y > x) | |
| 184 | { temp=x; x=y; y=temp; } | |
| 185 | if(x == zero) return(zero); | |
| 186 | if(y == zero) return(x); | |
| 187 | exp= logb(x); | |
| 188 | x=scalb(x,-exp); | |
| 189 | if(exp-(int)logb(y) > ibig ) | |
| 190 | /* raise inexact flag and return |x| */ | |
| 191 | { one+small; return(scalb(x,exp)); } | |
| 192 | else y=scalb(y,-exp); | |
| 193 | return(scalb(sqrt(x*x+y*y),exp)); | |
| 194 | } | |
| 195 | ||
| 196 | else if(y==y) /* y is +-INF */ | |
| 197 | return(copysign(y,one)); | |
| 198 | else | |
| 199 | return(y); /* y is NaN and x is finite */ | |
| 200 | ||
| 201 | else if(x==x) /* x is +-INF */ | |
| 202 | return (copysign(x,one)); | |
| 203 | else if(finite(y)) | |
| 204 | return(x); /* x is NaN, y is finite */ | |
| 205 | else if(y!=y) return(y); /* x and y is NaN */ | |
| 206 | else return(copysign(y,one)); /* y is INF */ | |
| 207 | } | |
| 208 | #endif |