src/math/hypotl.c

   1 #include "libm.h"
   2
   3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
   4 long double hypotl(long double x, long double y)
   5 {
   6         return hypot(x, y);
   7 }
   8 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
   9 #if LDBL_MANT_DIG == 64
  10 #define SPLIT (0x1p32L+1)
  11 #elif LDBL_MANT_DIG == 113
  12 #define SPLIT (0x1p57L+1)
  13 #endif
  14
  15 static void sq(long double *hi, long double *lo, long double x)
  16 {
  17         long double xh, xl, xc;
  18         xc = x*SPLIT;
  19         xh = x - xc + xc;
  20         xl = x - xh;
  21         *hi = x*x;
  22         *lo = xh*xh - *hi + 2*xh*xl + xl*xl;
  23 }
  24
  25 long double hypotl(long double x, long double y)
  26 {
  27         union ldshape ux = {x}, uy = {y};
  28         int ex, ey;
  29         long double hx, lx, hy, ly, z;
  30
  31         ux.i.se &= 0x7fff;
  32         uy.i.se &= 0x7fff;
  33         if (ux.i.se < uy.i.se) {
  34                 ex = uy.i.se;
  35                 ey = ux.i.se;
  36                 x = uy.f;
  37                 y = ux.f;
  38         } else {
  39                 ex = ux.i.se;
  40                 ey = uy.i.se;
  41                 x = ux.f;
  42                 y = uy.f;
  43         }
  44
  45         if (ex == 0x7fff && isinf(y))
  46                 return y;
  47         if (ex == 0x7fff || y == 0)
  48                 return x;
  49         if (ex - ey > LDBL_MANT_DIG)
  50                 return x + y;
  51
  52         z = 1;
  53         if (ex > 0x3fff+8000) {
  54                 z = 0x1p10000L;
  55                 x *= 0x1p-10000L;
  56                 y *= 0x1p-10000L;
  57         } else if (ey < 0x3fff-8000) {
  58                 z = 0x1p-10000L;
  59                 x *= 0x1p10000L;
  60                 y *= 0x1p10000L;
  61         }
  62         sq(&hx, &lx, x);
  63         sq(&hy, &ly, y);
  64         return z*sqrtl(ly+lx+hy+hx);
  65 }
  66 #endif