96     hx = __HI(x);
  97     k = ((hx&0x7ff00000)>>20) - 54;
  98     if (n< -50000) return tiny*x;       /*underflow*/
  99   }
 100   if (k==0x7ff) return x+x;             /* NaN or Inf */
 101   k = k+n;
 102   if (k >  0x7fe) return hugeX*copysign(hugeX,x); /* overflow  */
 103   if (k > 0)                            /* normal result */
 104     {__HI(x) = (hx&0x800fffff)|(k<<20); return x;}
 105   if (k <= -54) {
 106     if (n > 50000)      /* in case integer overflow in n+k */
 107       return hugeX*copysign(hugeX,x);   /*overflow*/
 108     else return tiny*copysign(tiny,x);  /*underflow*/
 109   }
 110   k += 54;                              /* subnormal result */
 111   __HI(x) = (hx&0x800fffff)|(k<<20);
 112   return x*twom54;
 113 }
 114 
 115 /* __ieee754_log(x)
 116  * Return the logrithm of x
 117  *
 118  * Method :
 119  *   1. Argument Reduction: find k and f such that
 120  *                    x = 2^k * (1+f),
 121  *       where  sqrt(2)/2 < 1+f < sqrt(2) .
 122  *
 123  *   2. Approximation of log(1+f).
 124  *    Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
 125  *             = 2s + 2/3 s**3 + 2/5 s**5 + .....,
 126  *             = 2s + s*R
 127  *      We use a special Reme algorithm on [0,0.1716] to generate
 128  *    a polynomial of degree 14 to approximate R The maximum error
 129  *    of this polynomial approximation is bounded by 2**-58.45. In
 130  *    other words,
 131  *                    2      4      6      8      10      12      14
 132  *        R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
 133  *    (the values of Lg1 to Lg7 are listed in the program)
 134  *    and
 135  *        |      2          14          |     -58.45
 136  *        | Lg1*s +...+Lg7*s    -  R(z) | <= 2

  96     hx = __HI(x);
  97     k = ((hx&0x7ff00000)>>20) - 54;
  98     if (n< -50000) return tiny*x;       /*underflow*/
  99   }
 100   if (k==0x7ff) return x+x;             /* NaN or Inf */
 101   k = k+n;
 102   if (k >  0x7fe) return hugeX*copysign(hugeX,x); /* overflow  */
 103   if (k > 0)                            /* normal result */
 104     {__HI(x) = (hx&0x800fffff)|(k<<20); return x;}
 105   if (k <= -54) {
 106     if (n > 50000)      /* in case integer overflow in n+k */
 107       return hugeX*copysign(hugeX,x);   /*overflow*/
 108     else return tiny*copysign(tiny,x);  /*underflow*/
 109   }
 110   k += 54;                              /* subnormal result */
 111   __HI(x) = (hx&0x800fffff)|(k<<20);
 112   return x*twom54;
 113 }
 114 
 115 /* __ieee754_log(x)
 116  * Return the logarithm of x
 117  *
 118  * Method :
 119  *   1. Argument Reduction: find k and f such that
 120  *                    x = 2^k * (1+f),
 121  *       where  sqrt(2)/2 < 1+f < sqrt(2) .
 122  *
 123  *   2. Approximation of log(1+f).
 124  *    Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
 125  *             = 2s + 2/3 s**3 + 2/5 s**5 + .....,
 126  *             = 2s + s*R
 127  *      We use a special Reme algorithm on [0,0.1716] to generate
 128  *    a polynomial of degree 14 to approximate R The maximum error
 129  *    of this polynomial approximation is bounded by 2**-58.45. In
 130  *    other words,
 131  *                    2      4      6      8      10      12      14
 132  *        R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
 133  *    (the values of Lg1 to Lg7 are listed in the program)
 134  *    and
 135  *        |      2          14          |     -58.45
 136  *        | Lg1*s +...+Lg7*s    -  R(z) | <= 2