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 /* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
00126    of long double implemented as double double.  */
00127 static inline long double
00128 ldbl_pack (double a, double aa)
00129 {
00130   union ibm_extended_long_double u;
00131   u.dd[0] = a;
00132   u.dd[1] = aa;
00133   return u.d;
00134 }
00135 
00136 static inline void
00137 ldbl_unpack (long double l, double *a, double *aa)
00138 {
00139   union ibm_extended_long_double u;
00140   u.d = l;
00141   *a = u.dd[0];
00142   *aa = u.dd[1];
00143 }
00144 
00145 
00146 /* Convert a finite long double to canonical form.
00147    Does not handle +/-Inf properly.  */
00148 static inline void
00149 ldbl_canonicalize (double *a, double *aa)
00150 {
00151   double xh, xl;
00152 
00153   xh = *a + *aa;
00154   xl = (*a - xh) + *aa;
00155   *a = xh;
00156   *aa = xl;
00157 }
00158 
00159 /* Simple inline nearbyint (double) function .
00160    Only works in the default rounding mode
00161    but is useful in long double rounding functions.  */
00162 static inline double
00163 ldbl_nearbyint (double a)
00164 {
00165   double two52 = 0x10000000000000LL;
00166 
00167   if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
00168     {
00169       if (__builtin_expect ((a > 0.0), 1))
00170        {
00171          a += two52;
00172          a -= two52;
00173        }
00174       else if (__builtin_expect ((a < 0.0), 1))
00175        {
00176          a = two52 - a;
00177          a = -(a - two52);
00178        }
00179     }
00180   return a;
00181 }