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