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