Back to index

glibc  2.9
math_ldbl.h
Go to the documentation of this file.
00001 #ifndef _MATH_PRIVATE_H_
00002 #error "Never use <math_ldbl.h> directly; include <math_private.h> instead."
00003 #endif
00004 
00005 #include <sysdeps/ieee754/ldbl-128/math_ldbl.h>
00006 #include <ieee754.h>
00007   
00008 static inline void
00009 ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x)
00010 {
00011   /* We have 105 bits of mantissa plus one implicit digit.  Since
00012      106 bits are representable we use the first implicit digit for
00013      the number before the decimal point and the second implicit bit
00014      as bit 53 of the mantissa.  */
00015   unsigned long long hi, lo;
00016   int ediff;
00017   union ibm_extended_long_double eldbl;
00018   eldbl.d = x;
00019   *exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;
00020 
00021   lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3;
00022   hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1;
00023   /* If the lower double is not a denomal or zero then set the hidden
00024      53rd bit.  */
00025   if (eldbl.ieee.exponent2 > 0x001)
00026     {
00027       lo |= (1ULL << 52);
00028       lo = lo << 7; /* pre-shift lo to match ieee854.  */
00029       /* The lower double is normalized separately from the upper.  We
00030         may need to adjust the lower manitissa to reflect this.  */
00031       ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2;
00032       if (ediff > 53)
00033        lo = lo >> (ediff-53);
00034     }
00035   hi |= (1ULL << 52);
00036   
00037   if ((eldbl.ieee.negative != eldbl.ieee.negative2)
00038       && ((eldbl.ieee.exponent2 != 0) && (lo != 0LL)))
00039     {
00040       hi--;
00041       lo = (1ULL << 60) - lo;
00042       if (hi < (1ULL << 52))
00043        {
00044          /* we have a borrow from the hidden bit, so shift left 1.  */
00045          hi = (hi << 1) | (lo >> 59);
00046          lo = 0xfffffffffffffffLL & (lo << 1);
00047          *exp = *exp - 1;
00048        }
00049     }
00050   *lo64 = (hi << 60) | lo;
00051   *hi64 = hi >> 4;
00052 }
00053 
00054 static inline long double
00055 ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
00056 {
00057   union ibm_extended_long_double u;
00058   unsigned long hidden2, lzcount;
00059   unsigned long long hi, lo;
00060 
00061   u.ieee.negative = sign;
00062   u.ieee.negative2 = sign;
00063   u.ieee.exponent = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS;
00064   u.ieee.exponent2 = exp-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
00065   /* Expect 113 bits (112 bits + hidden) right justified in two longs.
00066      The low order 53 bits (52 + hidden) go into the lower double */ 
00067   lo = (lo64 >> 7)& ((1ULL << 53) - 1);
00068   hidden2 = (lo64 >> 59) &  1ULL;
00069   /* The high order 53 bits (52 + hidden) go into the upper double */
00070   hi = (lo64 >> 60) & ((1ULL << 11) - 1);
00071   hi |= (hi64 << 4);
00072 
00073   if (lo != 0LL)
00074     {
00075       /* hidden2 bit of low double controls rounding of the high double.
00076         If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)
00077         plus change the sign of the low double to compensate.  */
00078       if (hidden2)
00079        {
00080          hi++;
00081          u.ieee.negative2 = !sign;
00082          lo = (1ULL << 53) - lo;
00083        }
00084       /* The hidden bit of the lo mantissa is zero so we need to
00085         normalize the it for the low double.  Shift it left until the
00086         hidden bit is '1' then adjust the 2nd exponent accordingly.  */ 
00087 
00088       if (sizeof (lo) == sizeof (long))
00089        lzcount = __builtin_clzl (lo);
00090       else if ((lo >> 32) != 0)
00091        lzcount = __builtin_clzl ((long) (lo >> 32));
00092       else
00093        lzcount = __builtin_clzl ((long) lo) + 32;
00094       lzcount = lzcount - 11;
00095       if (lzcount > 0)
00096        {
00097          int expnt2 = u.ieee.exponent2 - lzcount;
00098          if (expnt2 >= 1)
00099            {
00100              /* Not denormal.  Normalize and set low exponent.  */
00101              lo = lo << lzcount;
00102              u.ieee.exponent2 = expnt2;
00103            }
00104          else
00105            {
00106              /* Is denormal.  */
00107              lo = lo << (lzcount + expnt2);
00108              u.ieee.exponent2 = 0;
00109            }
00110        }
00111     }
00112   else
00113     {
00114       u.ieee.negative2 = 0;
00115       u.ieee.exponent2 = 0;
00116     }
00117 
00118   u.ieee.mantissa3 = lo & ((1ULL << 32) - 1);
00119   u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1);
00120   u.ieee.mantissa1 = hi & ((1ULL << 32) - 1);
00121   u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1);
00122   return u.d;
00123 }
00124   
00125 /* gcc generates disgusting code to pack and unpack long doubles.
00126    This tells gcc that pack/unpack is really a nop.  We use fr1/fr2
00127    because those are the regs used to pass/return a single
00128    long double arg.  */
00129 static inline long double
00130 ldbl_pack (double a, double aa)
00131 {
00132   register long double x __asm__ ("fr1");
00133   register double xh __asm__ ("fr1");
00134   register double xl __asm__ ("fr2");
00135   xh = a;
00136   xl = aa;
00137   __asm__ ("" : "=f" (x) : "f" (xh), "f" (xl));
00138   return x;
00139 }
00140 
00141 static inline void
00142 ldbl_unpack (long double l, double *a, double *aa)
00143 {
00144   register long double x __asm__ ("fr1");
00145   register double xh __asm__ ("fr1");
00146   register double xl __asm__ ("fr2");
00147   x = l;
00148   __asm__ ("" : "=f" (xh), "=f" (xl) : "f" (x));
00149   *a = xh;
00150   *aa = xl;
00151 }
00152 
00153 
00154 /* Convert a finite long double to canonical form.
00155    Does not handle +/-Inf properly.  */
00156 static inline void
00157 ldbl_canonicalize (double *a, double *aa)
00158 {
00159   double xh, xl;
00160 
00161   xh = *a + *aa;
00162   xl = (*a - xh) + *aa;
00163   *a = xh;
00164   *aa = xl;
00165 }
00166 
00167 /* Simple inline nearbyint (double) function .
00168    Only works in the default rounding mode
00169    but is useful in long double rounding functions.  */
00170 static inline double
00171 ldbl_nearbyint (double a)
00172 {
00173   double two52 = 0x10000000000000LL;
00174 
00175   if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
00176     {
00177       if (__builtin_expect ((a > 0.0), 1))
00178        {
00179          a += two52;
00180          a -= two52;
00181        }
00182       else if (__builtin_expect ((a < 0.0), 1))
00183        {
00184          a = two52 - a;
00185          a = -(a - two52);
00186        }
00187     }
00188   return a;
00189 }