math_ldbl.h

glibc 库, 不仅可以学习使用库函数,还可以学习函数的具体实现,是提高功力的好资料

#ifndef _MATH_PRIVATE_H_
#error "Never use <math_ldbl.h> directly; include <math_private.h> instead."
#endif

#include <sysdeps/ieee754/ldbl-128/math_ldbl.h>
#include <ieee754.h>
  
static inline void
ldbl_extract_mantissa (int64_t *hi64, u_int64_t *lo64, int *exp, long double x)
{
  /* We have 105 bits of mantissa plus one implicit digit.  Since
     106 bits are representable we use the first implicit digit for
     the number before the decimal point and the second implicit bit
     as bit 53 of the mantissa.  */
  unsigned long long hi, lo;
  int ediff;
  union ibm_extended_long_double eldbl;
  eldbl.d = x;
  *exp = eldbl.ieee.exponent - IBM_EXTENDED_LONG_DOUBLE_BIAS;

  lo = ((long long)eldbl.ieee.mantissa2 << 32) | eldbl.ieee.mantissa3;
  hi = ((long long)eldbl.ieee.mantissa0 << 32) | eldbl.ieee.mantissa1;
  /* If the lower double is not a denomal or zero then set the hidden
     53rd bit.  */
  if (eldbl.ieee.exponent2 > 0x001)
    {
      lo |= (1ULL << 52);
      lo = lo << 7; /* pre-shift lo to match ieee854.  */
      /* The lower double is normalized separately from the upper.  We
	 may need to adjust the lower manitissa to reflect this.  */
      ediff = eldbl.ieee.exponent - eldbl.ieee.exponent2;
      if (ediff > 53)
	lo = lo >> (ediff-53);
    }
  hi |= (1ULL << 52);
  
  if ((eldbl.ieee.negative != eldbl.ieee.negative2)
      && ((eldbl.ieee.exponent2 != 0) && (lo != 0LL)))
    {
      hi--;
      lo = (1ULL << 60) - lo;
      if (hi < (1ULL << 52))
	{
	  /* we have a borrow from the hidden bit, so shift left 1.  */
	  hi = (hi << 1) | (lo >> 59);
	  lo = 0xfffffffffffffffLL & (lo << 1);
	  *exp = *exp - 1;
	}
    }
  *lo64 = (hi << 60) | lo;
  *hi64 = hi >> 4;
}

static inline long double
ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
{
  union ibm_extended_long_double u;
  unsigned long hidden2, lzcount;
  unsigned long long hi, lo;

  u.ieee.negative = sign;
  u.ieee.negative2 = sign;
  u.ieee.exponent = exp + IBM_EXTENDED_LONG_DOUBLE_BIAS;
  u.ieee.exponent2 = exp-53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
  /* Expect 113 bits (112 bits + hidden) right justified in two longs.
     The low order 53 bits (52 + hidden) go into the lower double */ 
  lo = (lo64 >> 7)& ((1ULL << 53) - 1);
  hidden2 = (lo64 >> 59) &  1ULL;
  /* The high order 53 bits (52 + hidden) go into the upper double */
  hi = (lo64 >> 60) & ((1ULL << 11) - 1);
  hi |= (hi64 << 4);

  if (lo != 0LL)
    {
      /* hidden2 bit of low double controls rounding of the high double.
	 If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)
	 plus change the sign of the low double to compensate.  */
      if (hidden2)
	{
	  hi++;
	  u.ieee.negative2 = !sign;
	  lo = (1ULL << 53) - lo;
	}
      /* The hidden bit of the lo mantissa is zero so we need to
	 normalize the it for the low double.  Shift it left until the
	 hidden bit is '1' then adjust the 2nd exponent accordingly.  */ 

      if (sizeof (lo) == sizeof (long))
	lzcount = __builtin_clzl (lo);
      else if ((lo >> 32) != 0)
	lzcount = __builtin_clzl ((long) (lo >> 32));
      else
	lzcount = __builtin_clzl ((long) lo) + 32;
      lzcount = lzcount - 11;
      if (lzcount > 0)
	{
	  int expnt2 = u.ieee.exponent2 - lzcount;
	  if (expnt2 >= 1)
	    {
	      /* Not denormal.  Normalize and set low exponent.  */
	      lo = lo << lzcount;
	      u.ieee.exponent2 = expnt2;
	    }
	  else
	    {
	      /* Is denormal.  */
	      lo = lo << (lzcount + expnt2);
	      u.ieee.exponent2 = 0;
	    }
	}
    }
  else
    {
      u.ieee.negative2 = 0;
      u.ieee.exponent2 = 0;
    }

  u.ieee.mantissa3 = lo & ((1ULL << 32) - 1);
  u.ieee.mantissa2 = (lo >> 32) & ((1ULL << 20) - 1);
  u.ieee.mantissa1 = hi & ((1ULL << 32) - 1);
  u.ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1);
  return u.d;
}
  
/* gcc generates disgusting code to pack and unpack long doubles.
   This tells gcc that pack/unpack is really a nop.  We use fr1/fr2
   because those are the regs used to pass/return a single
   long double arg.  */
static inline long double
ldbl_pack (double a, double aa)
{
  register long double x __asm__ ("fr1");
  register double xh __asm__ ("fr1");
  register double xl __asm__ ("fr2");
  xh = a;
  xl = aa;
  __asm__ ("" : "=f" (x) : "f" (xh), "f" (xl));
  return x;
}

static inline void
ldbl_unpack (long double l, double *a, double *aa)
{
  register long double x __asm__ ("fr1");
  register double xh __asm__ ("fr1");
  register double xl __asm__ ("fr2");
  x = l;
  __asm__ ("" : "=f" (xh), "=f" (xl) : "f" (x));
  *a = xh;
  *aa = xl;
}


/* Convert a finite long double to canonical form.
   Does not handle +/-Inf properly.  */
static inline void
ldbl_canonicalize (double *a, double *aa)
{
  double xh, xl;

  xh = *a + *aa;
  xl = (*a - xh) + *aa;
  *a = xh;
  *aa = xl;
}

/* Simple inline nearbyint (double) function .
   Only works in the default rounding mode
   but is useful in long double rounding functions.  */
static inline double
ldbl_nearbyint (double a)
{
  double two52 = 0x10000000000000LL;

  if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
    {
      if (__builtin_expect ((a > 0.0), 1))
	{
	  a += two52;
	  a -= two52;
	}
      else if (__builtin_expect ((a < 0.0), 1))
	{
	  a = two52 - a;
	  a = -(a - two52);
	}
    }
  return a;
}