sysdeps/ieee754/ldbl-128ibm/math_ldbl.h

   1 #ifndef _MATH_PRIVATE_H_
   2 #error "Never use <math_ldbl.h> directly; include <math_private.h> instead."
   3 #endif
   4
   5 #include <ieee754.h>
   6 #include <stdint.h>
   7
   8 /* To suit our callers we return *hi64 and *lo64 as if they came from
   9    an ieee854 112 bit mantissa, that is, 48 bits in *hi64 (plus one
  10    implicit bit) and 64 bits in *lo64.  */
  11
  12 static inline void
  13 ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x)
  14 {
  15   /* We have 105 bits of mantissa plus one implicit digit.  Since
  16      106 bits are representable we use the first implicit digit for
  17      the number before the decimal point and the second implicit bit
  18      as bit 53 of the mantissa.  */
  19   uint64_t hi, lo;
  20   union ibm_extended_long_double u;
  21
  22   u.ld = x;
  23   *exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
  24
  25   lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
  26   hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
  27
  28   if (u.d[0].ieee.exponent != 0)
  29     {
  30       int ediff;
  31
  32       /* If not a denormal or zero then we have an implicit 53rd bit.  */
  33       hi |= (uint64_t) 1 << 52;
  34
  35       if (u.d[1].ieee.exponent != 0)
  36         lo |= (uint64_t) 1 << 52;
  37       else
  38         /* A denormal is to be interpreted as having a biased exponent
  39            of 1.  */
  40         lo = lo << 1;
  41
  42       /* We are going to shift 4 bits out of hi later, because we only
  43          want 48 bits in *hi64.  That means we want 60 bits in lo, but
  44          we currently only have 53.  Shift the value up.  */
  45       lo = lo << 7;
  46
  47       /* The lower double is normalized separately from the upper.
  48          We may need to adjust the lower mantissa to reflect this.
  49          The difference between the exponents can be larger than 53
  50          when the low double is much less than 1ULP of the upper
  51          (in which case there are significant bits, all 0's or all
  52          1's, between the two significands).  The difference between
  53          the exponents can be less than 53 when the upper double
  54          exponent is nearing its minimum value (in which case the low
  55          double is denormal ie. has an exponent of zero).  */
  56       ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
  57       if (ediff > 0)
  58         {
  59           if (ediff < 64)
  60             lo = lo >> ediff;
  61           else
  62             lo = 0;
  63         }
  64       else if (ediff < 0)
  65         lo = lo << -ediff;
  66
  67       if (u.d[0].ieee.negative != u.d[1].ieee.negative
  68           && lo != 0)
  69         {
  70           hi--;
  71           lo = ((uint64_t) 1 << 60) - lo;
  72           if (hi < (uint64_t) 1 << 52)
  73             {
  74               /* We have a borrow from the hidden bit, so shift left 1.  */
  75               hi = (hi << 1) | (lo >> 59);
  76               lo = (((uint64_t) 1 << 60) - 1) & (lo << 1);
  77               *exp = *exp - 1;
  78             }
  79         }
  80     }
  81   else
  82     /* If the larger magnitude double is denormal then the smaller
  83        one must be zero.  */
  84     hi = hi << 1;
  85
  86   *lo64 = (hi << 60) | lo;
  87   *hi64 = hi >> 4;
  88 }
  89
  90 static inline long double
  91 ldbl_insert_mantissa (int sign, int exp, int64_t hi64, uint64_t lo64)
  92 {
  93   union ibm_extended_long_double u;
  94   int expnt2;
  95   uint64_t hi, lo;
  96
  97   u.d[0].ieee.negative = sign;
  98   u.d[1].ieee.negative = sign;
  99   u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS;
 100   u.d[1].ieee.exponent = 0;
 101   expnt2 = exp - 53 + IEEE754_DOUBLE_BIAS;
 102
 103   /* Expect 113 bits (112 bits + hidden) right justified in two longs.
 104      The low order 53 bits (52 + hidden) go into the lower double */
 105   lo = (lo64 >> 7) & (((uint64_t) 1 << 53) - 1);
 106   /* The high order 53 bits (52 + hidden) go into the upper double */
 107   hi = lo64 >> 60;
 108   hi |= hi64 << 4;
 109
 110   if (lo != 0)
 111     {
 112       int lzcount;
 113
 114       /* hidden bit of low double controls rounding of the high double.
 115          If hidden is '1' and either the explicit mantissa is non-zero
 116          or hi is odd, then round up hi and adjust lo (2nd mantissa)
 117          plus change the sign of the low double to compensate.  */
 118       if ((lo & ((uint64_t) 1 << 52)) != 0
 119           && ((hi & 1) != 0 || (lo & (((uint64_t) 1 << 52) - 1)) != 0))
 120         {
 121           hi++;
 122           if ((hi & ((uint64_t) 1 << 53)) != 0)
 123             {
 124               hi = hi >> 1;
 125               u.d[0].ieee.exponent++;
 126             }
 127           u.d[1].ieee.negative = !sign;
 128           lo = ((uint64_t) 1 << 53) - lo;
 129         }
 130
 131       /* Normalize the low double.  Shift the mantissa left until
 132          the hidden bit is '1' and adjust the exponent accordingly.  */
 133
 134       if (sizeof (lo) == sizeof (long))
 135         lzcount = __builtin_clzl (lo);
 136       else if ((lo >> 32) != 0)
 137         lzcount = __builtin_clzl ((long) (lo >> 32));
 138       else
 139         lzcount = __builtin_clzl ((long) lo) + 32;
 140       lzcount = lzcount - (64 - 53);
 141       lo <<= lzcount;
 142       expnt2 -= lzcount;
 143
 144       if (expnt2 >= 1)
 145         /* Not denormal.  */
 146         u.d[1].ieee.exponent = expnt2;
 147       else
 148         {
 149           /* Is denormal.  Note that biased exponent of 0 is treated
 150              as if it was 1, hence the extra shift.  */
 151           if (expnt2 > -53)
 152             lo >>= 1 - expnt2;
 153           else
 154             lo = 0;
 155         }
 156     }
 157   else
 158     u.d[1].ieee.negative = 0;
 159
 160   u.d[1].ieee.mantissa1 = lo;
 161   u.d[1].ieee.mantissa0 = lo >> 32;
 162   u.d[0].ieee.mantissa1 = hi;
 163   u.d[0].ieee.mantissa0 = hi >> 32;
 164   return u.ld;
 165 }
 166
 167 /* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
 168    of long double implemented as double double.  */
 169 static inline long double
 170 default_ldbl_pack (double a, double aa)
 171 {
 172   union ibm_extended_long_double u;
 173   u.d[0].d = a;
 174   u.d[1].d = aa;
 175   return u.ld;
 176 }
 177
 178 static inline void
 179 default_ldbl_unpack (long double l, double *a, double *aa)
 180 {
 181   union ibm_extended_long_double u;
 182   u.ld = l;
 183   *a = u.d[0].d;
 184   *aa = u.d[1].d;
 185 }
 186
 187 #ifndef ldbl_pack
 188 # define ldbl_pack   default_ldbl_pack
 189 #endif
 190 #ifndef ldbl_unpack
 191 # define ldbl_unpack default_ldbl_unpack
 192 #endif
 193
 194 /* Extract high double.  */
 195 #define ldbl_high(x) ((double) x)
 196
 197 /* Convert a finite long double to canonical form.
 198    Does not handle +/-Inf properly.  */
 199 static inline void
 200 ldbl_canonicalize (double *a, double *aa)
 201 {
 202   double xh, xl;
 203
 204   xh = *a + *aa;
 205   xl = (*a - xh) + *aa;
 206   *a = xh;
 207   *aa = xl;
 208 }
 209
 210 /* Simple inline nearbyint (double) function.
 211    Only works in the default rounding mode
 212    but is useful in long double rounding functions.  */
 213 static inline double
 214 ldbl_nearbyint (double a)
 215 {
 216   double two52 = 0x1p52;
 217
 218   if (__glibc_likely ((__builtin_fabs (a) < two52)))
 219     {
 220       if (__glibc_likely ((a > 0.0)))
 221         {
 222           a += two52;
 223           a -= two52;
 224         }
 225       else if (__glibc_likely ((a < 0.0)))
 226         {
 227           a = two52 - a;
 228           a = -(a - two52);
 229         }
 230     }
 231   return a;
 232 }