sysdeps/powerpc/fpu/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
   8 static inline void
   9 ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x)
  10 {
  11   /* We have 105 bits of mantissa plus one implicit digit.  Since
  12      106 bits are representable we use the first implicit digit for
  13      the number before the decimal point and the second implicit bit
  14      as bit 53 of the mantissa.  */
  15   unsigned long long hi, lo;
  16   int ediff;
  17   union ibm_extended_long_double eldbl;
  18   eldbl.d = x;
  19   *exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;
  20
  21   lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3;
  22   hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1;
  23   /* If the lower double is not a denomal or zero then set the hidden
  24      53rd bit.  */
  25   if (eldbl.ieee.exponent2 > 0x001)
  26     {
  27       lo |= (1ULL << 52);
  28       lo = lo << 7; /* pre-shift lo to match ieee854.  */
  29       /* The lower double is normalized separately from the upper.  We
  30          may need to adjust the lower manitissa to reflect this.  */
  31       ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2;
  32       if (ediff > 53)
  33         lo = lo >> (ediff-53);
  34     }
  35   hi |= (1ULL << 52);
  36
  37   if ((eldbl.ieee.negative != eldbl.ieee.negative2)
  38       && ((eldbl.ieee.exponent2 != 0) && (lo != 0LL)))
  39     {
  40       hi--;
  41       lo = (1ULL << 60) - lo;
  42       if (hi < (1ULL << 52))
  43         {
  44           /* we have a borrow from the hidden bit, so shift left 1.  */
  45           hi = (hi << 1) | (lo >> 59);
  46           lo = 0xfffffffffffffffLL & (lo << 1);
  47           *exp = *exp - 1;
  48         }
  49     }
  50   *lo64 = (hi << 60) | lo;
  51   *hi64 = hi >> 4;
  52 }
  53
  54 static inline long double
  55 ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
  56 {
  57   union ibm_extended_long_double u;
  58   unsigned long hidden2, lzcount;
  59   unsigned long long hi, lo;
  60
  61   u.ieee.negative = sign;
  62   u.ieee.negative2 = sign;
  63   u.ieee.exponent = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS;
  64   u.ieee.exponent2 = exp-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
  65   /* Expect 113 bits (112 bits + hidden) right justified in two longs.
  66      The low order 53 bits (52 + hidden) go into the lower double */
  67   lo = (lo64 >> 7)& ((1ULL << 53) - 1);
  68   hidden2 = (lo64 >> 59) &  1ULL;
  69   /* The high order 53 bits (52 + hidden) go into the upper double */
  70   hi = (lo64 >> 60) & ((1ULL << 11) - 1);
  71   hi |= (hi64 << 4);
  72
  73   if (lo != 0LL)
  74     {
  75       /* hidden2 bit of low double controls rounding of the high double.
  76          If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)
  77          plus change the sign of the low double to compensate.  */
  78       if (hidden2)
  79         {
  80           hi++;
  81           u.ieee.negative2 = !sign;
  82           lo = (1ULL << 53) - lo;
  83         }
  84       /* The hidden bit of the lo mantissa is zero so we need to
  85          normalize the it for the low double.  Shift it left until the
  86          hidden bit is '1' then adjust the 2nd exponent accordingly.  */
  87
  88       if (sizeof (lo) == sizeof (long))
  89         lzcount = __builtin_clzl (lo);
  90       else if ((lo >> 32) != 0)
  91         lzcount = __builtin_clzl ((long) (lo >> 32));
  92       else
  93         lzcount = __builtin_clzl ((long) lo) + 32;
  94       lzcount = lzcount - 11;
  95       if (lzcount > 0)
  96         {
  97           int expnt2 = u.ieee.exponent2 - lzcount;
  98           if (expnt2 >= 1)
  99             {
 100               /* Not denormal.  Normalize and set low exponent.  */
 101               lo = lo << lzcount;
 102               u.ieee.exponent2 = expnt2;
 103             }
 104           else
 105             {
 106               /* Is denormal.  */
 107               lo = lo << (lzcount + expnt2);
 108               u.ieee.exponent2 = 0;
 109             }
 110         }
 111     }
 112   else
 113     {
 114       u.ieee.negative2 = 0;
 115       u.ieee.exponent2 = 0;
 116     }
 117
 118   u.ieee.mantissa3 = lo & ((1ULL << 32) - 1);
 119   u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1);
 120   u.ieee.mantissa1 = hi & ((1ULL << 32) - 1);
 121   u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1);
 122   return u.d;
 123 }
 124
 125 /* gcc generates disgusting code to pack and unpack long doubles.
 126    This tells gcc that pack/unpack is really a nop.  We use fr1/fr2
 127    because those are the regs used to pass/return a single
 128    long double arg.  */
 129 static inline long double
 130 ldbl_pack (double a, double aa)
 131 {
 132   register long double x __asm__ ("fr1");
 133   register double xh __asm__ ("fr1");
 134   register double xl __asm__ ("fr2");
 135   xh = a;
 136   xl = aa;
 137   __asm__ ("" : "=f" (x) : "f" (xh), "f" (xl));
 138   return x;
 139 }
 140
 141 static inline void
 142 ldbl_unpack (long double l, double *a, double *aa)
 143 {
 144   register long double x __asm__ ("fr1");
 145   register double xh __asm__ ("fr1");
 146   register double xl __asm__ ("fr2");
 147   x = l;
 148   __asm__ ("" : "=f" (xh), "=f" (xl) : "f" (x));
 149   *a = xh;
 150   *aa = xl;
 151 }
 152
 153
 154 /* Convert a finite long double to canonical form.
 155    Does not handle +/-Inf properly.  */
 156 static inline void
 157 ldbl_canonicalize (double *a, double *aa)
 158 {
 159   double xh, xl;
 160
 161   xh = *a + *aa;
 162   xl = (*a - xh) + *aa;
 163   *a = xh;
 164   *aa = xl;
 165 }
 166
 167 /* Simple inline nearbyint (double) function .
 168    Only works in the default rounding mode
 169    but is useful in long double rounding functions.  */
 170 static inline double
 171 ldbl_nearbyint (double a)
 172 {
 173   double two52 = 0x10000000000000LL;
 174
 175   if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
 176     {
 177       if (__builtin_expect ((a > 0.0), 1))
 178         {
 179           a += two52;
 180           a -= two52;
 181         }
 182       else if (__builtin_expect ((a < 0.0), 1))
 183         {
 184           a = two52 - a;
 185           a = -(a - two52);
 186         }
 187     }
 188   return a;
 189 }