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