src/math/hypotl.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_hypotl.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /* long double version of hypot().  See comments in hypot.c. */
  13
  14 #include "libm.h"
  15
  16 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  17 long double hypotl(long double x, long double y)
  18 {
  19         return hypot(x, y);
  20 }
  21 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
  22
  23 #define GET_LDBL_EXPSIGN(i, v) do {     \
  24         union IEEEl2bits uv;            \
  25                                         \
  26         uv.e = v;                       \
  27         i = uv.xbits.expsign;           \
  28 } while (0)
  29
  30 #define GET_LDBL_MAN(h, l, v) do {      \
  31         union IEEEl2bits uv;            \
  32                                         \
  33         uv.e = v;                       \
  34         h = uv.bits.manh;               \
  35         l = uv.bits.manl;               \
  36 } while (0)
  37
  38 #define SET_LDBL_EXPSIGN(v, i) do {     \
  39         union IEEEl2bits uv;            \
  40                                         \
  41         uv.e = v;                       \
  42         uv.xbits.expsign = i;           \
  43         v = uv.e;                       \
  44 } while (0)
  45
  46 #undef GET_HIGH_WORD
  47 #define GET_HIGH_WORD(i, v)     GET_LDBL_EXPSIGN(i, v)
  48 #undef SET_HIGH_WORD
  49 #define SET_HIGH_WORD(v, i)     SET_LDBL_EXPSIGN(v, i)
  50
  51 #define DESW(exp)       (exp)           /* delta expsign word */
  52 #define ESW(exp)        (MAX_EXP - 1 + (exp))   /* expsign word */
  53 #define MANT_DIG        LDBL_MANT_DIG
  54 #define MAX_EXP         LDBL_MAX_EXP
  55
  56 #if LDBL_MANL_SIZE > 32
  57 typedef uint64_t man_t;
  58 #else
  59 typedef uint32_t man_t;
  60 #endif
  61
  62 long double hypotl(long double x, long double y)
  63 {
  64         long double a=x,b=y,t1,t2,y1,y2,w;
  65         int32_t j,k,ha,hb;
  66
  67         GET_HIGH_WORD(ha, x);
  68         ha &= 0x7fff;
  69         GET_HIGH_WORD(hb, y);
  70         hb &= 0x7fff;
  71         if (hb > ha) {
  72                 a = y;
  73                 b = x;
  74                 j=ha; ha=hb; hb=j;
  75         } else {
  76                 a = x;
  77                 b = y;
  78         }
  79         a = fabsl(a);
  80         b = fabsl(b);
  81         if (ha - hb > DESW(MANT_DIG+7))  /* x/y > 2**(MANT_DIG+7) */
  82                 return a+b;
  83         k = 0;
  84         if (ha > ESW(MAX_EXP/2-12)) {    /* a>2**(MAX_EXP/2-12) */
  85                 if (ha >= ESW(MAX_EXP)) {  /* Inf or NaN */
  86                         man_t manh, manl;
  87                         /* Use original arg order iff result is NaN; quieten sNaNs. */
  88                         w = fabsl(x+0.0)-fabsl(y+0.0);
  89                         GET_LDBL_MAN(manh,manl,a);
  90                         if (manh == LDBL_NBIT && manl == 0) w = a;
  91                         GET_LDBL_MAN(manh,manl,b);
  92                         if (hb >= ESW(MAX_EXP) && manh == LDBL_NBIT && manl == 0) w = b;
  93                         return w;
  94                 }
  95                 /* scale a and b by 2**-(MAX_EXP/2+88) */
  96                 ha -= DESW(MAX_EXP/2+88); hb -= DESW(MAX_EXP/2+88);
  97                 k += MAX_EXP/2+88;
  98                 SET_HIGH_WORD(a, ha);
  99                 SET_HIGH_WORD(b, hb);
 100         }
 101         if (hb < ESW(-(MAX_EXP/2-12))) {  /* b < 2**-(MAX_EXP/2-12) */
 102                 if (hb <= 0) {  /* subnormal b or 0 */
 103                         man_t manh, manl;
 104                         GET_LDBL_MAN(manh,manl,b);
 105                         if ((manh|manl) == 0)
 106                                 return a;
 107                         t1 = 0;
 108                         SET_HIGH_WORD(t1, ESW(MAX_EXP-2));  /* t1 = 2^(MAX_EXP-2) */
 109                         b *= t1;
 110                         a *= t1;
 111                         k -= MAX_EXP-2;
 112                 } else {            /* scale a and b by 2^(MAX_EXP/2+88) */
 113                         ha += DESW(MAX_EXP/2+88);
 114                         hb += DESW(MAX_EXP/2+88);
 115                         k -= MAX_EXP/2+88;
 116                         SET_HIGH_WORD(a, ha);
 117                         SET_HIGH_WORD(b, hb);
 118                 }
 119         }
 120         /* medium size a and b */
 121         w = a - b;
 122         if (w > b) {
 123                 t1 = a;
 124                 union IEEEl2bits uv;
 125                 uv.e = t1; uv.bits.manl = 0; t1 = uv.e;
 126                 t2 = a-t1;
 127                 w  = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
 128         } else {
 129                 a  = a+a;
 130                 y1 = b;
 131                 union IEEEl2bits uv;
 132                 uv.e = y1; uv.bits.manl = 0; y1 = uv.e;
 133                 y2 = b - y1;
 134                 t1 = a;
 135                 uv.e = t1; uv.bits.manl = 0; t1 = uv.e;
 136                 t2 = a - t1;
 137                 w  = sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 138         }
 139         if(k!=0) {
 140                 uint32_t high;
 141                 t1 = 1.0;
 142                 GET_HIGH_WORD(high, t1);
 143                 SET_HIGH_WORD(t1, high+DESW(k));
 144                 return t1*w;
 145         }
 146         return w;
 147 }
 148 #endif