sys/src/libc/port/hypot.c

   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 }