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