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
hypot(double p,double q)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 }
42