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