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