1 /* 2 * hypot -- sqrt(p*p+q*q), but overflows only if the result does. 3 * See Cleve Moler and Donald Morrison, 4 * ``Replacing Square Roots by Pythagorean Sums,'' 5 * IBM Journal of Research and Development, 6 * Vol. 27, Number 6, pp. 577-581, Nov. 1983 7 */ 8 9 #include <u.h> 10 #include <libc.h> 11 12 double 13 hypot(double p, double q) 14 { 15 double r, s, pfac; 16 17 if(p < 0) 18 p = -p; 19 if(q < 0) 20 q = -q; 21 if(p < q) { 22 r = p; 23 p = q; 24 q = r; 25 } 26 if(p == 0) 27 return 0; 28 pfac = p; 29 r = q = q/p; 30 p = 1; 31 for(;;) { 32 r *= r; 33 s = r+4; 34 if(s == 4) 35 return p*pfac; 36 r /= s; 37 p += 2*r*p; 38 q *= r; 39 r = q/p; 40 } 41 return 0; 42 } 43