Line data    Source code

       1             : 
       2             : /* @(#)e_log10.c 1.3 95/01/18 */
       3             : /*
       4             :  * ====================================================
       5             :  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       6             :  *
       7             :  * Developed at SunSoft, a Sun Microsystems, Inc. business.
       8             :  * Permission to use, copy, modify, and distribute this
       9             :  * software is freely granted, provided that this notice 
      10             :  * is preserved.
      11             :  * ====================================================
      12             :  */
      13             : 
      14             : //#include <sys/cdefs.h>
      15             : //__FBSDID("$FreeBSD$");
      16             : 
      17             : /*
      18             :  * Return the base 2 logarithm of x.  See e_log.c and k_log.h for most
      19             :  * comments.
      20             :  *
      21             :  * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
      22             :  * then does the combining and scaling steps
      23             :  *    log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
      24             :  * in not-quite-routine extra precision.
      25             :  */
      26             : 
      27             : #include <float.h>
      28             : 
      29             : #include "math_private.h"
      30             : #include "k_log.h"
      31             : 
      32             : static const double
      33             : two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
      34             : ivln2hi    =  1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
      35             : ivln2lo    =  1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
      36             : 
      37             : static const double zero   =  0.0;
      38             : static volatile double vzero = 0.0;
      39             : 
      40             : double
      41           0 : __ieee754_log2(double x)
      42             : {
      43             :         double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
      44             :         int32_t i,k,hx;
      45             :         u_int32_t lx;
      46             : 
      47           0 :         EXTRACT_WORDS(hx,lx,x);
      48             : 
      49           0 :         k=0;
      50           0 :         if (hx < 0x00100000) {                       /* x < 2**-1022  */
      51           0 :             if (((hx&0x7fffffff)|lx)==0)
      52           0 :                 return -two54/vzero;            /* log(+-0)=-inf */
      53           0 :             if (hx<0) return (x-x)/zero;     /* log(-#) = NaN */
      54           0 :             k -= 54; x *= two54; /* subnormal number, scale up x */
      55           0 :             GET_HIGH_WORD(hx,x);
      56             :         }
      57           0 :         if (hx >= 0x7ff00000) return x+x;
      58           0 :         if (hx == 0x3ff00000 && lx == 0)
      59           0 :             return zero;                        /* log(1) = +0 */
      60           0 :         k += (hx>>20)-1023;
      61           0 :         hx &= 0x000fffff;
      62           0 :         i = (hx+0x95f64)&0x100000;
      63           0 :         SET_HIGH_WORD(x,hx|(i^0x3ff00000));     /* normalize x or x/2 */
      64           0 :         k += (i>>20);
      65           0 :         y = (double)k;
      66           0 :         f = x - 1.0;
      67           0 :         hfsq = 0.5*f*f;
      68           0 :         r = k_log1p(f);
      69             : 
      70             :         /*
      71             :          * f-hfsq must (for args near 1) be evaluated in extra precision
      72             :          * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
      73             :          * This is fairly efficient since f-hfsq only depends on f, so can
      74             :          * be evaluated in parallel with R.  Not combining hfsq with R also
      75             :          * keeps R small (though not as small as a true `lo' term would be),
      76             :          * so that extra precision is not needed for terms involving R.
      77             :          *
      78             :          * Compiler bugs involving extra precision used to break Dekker's
      79             :          * theorem for spitting f-hfsq as hi+lo, unless double_t was used
      80             :          * or the multi-precision calculations were avoided when double_t
      81             :          * has extra precision.  These problems are now automatically
      82             :          * avoided as a side effect of the optimization of combining the
      83             :          * Dekker splitting step with the clear-low-bits step.
      84             :          *
      85             :          * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
      86             :          * precision to avoid a very large cancellation when x is very near
      87             :          * these values.  Unlike the above cancellations, this problem is
      88             :          * specific to base 2.  It is strange that adding +-1 is so much
      89             :          * harder than adding +-ln2 or +-log10_2.
      90             :          *
      91             :          * This uses Dekker's theorem to normalize y+val_hi, so the
      92             :          * compiler bugs are back in some configurations, sigh.  And I
      93             :          * don't want to used double_t to avoid them, since that gives a
      94             :          * pessimization and the support for avoiding the pessimization
      95             :          * is not yet available.
      96             :          *
      97             :          * The multi-precision calculations for the multiplications are
      98             :          * routine.
      99             :          */
     100           0 :         hi = f - hfsq;
     101           0 :         SET_LOW_WORD(hi,0);
     102           0 :         lo = (f - hi) - hfsq + r;
     103           0 :         val_hi = hi*ivln2hi;
     104           0 :         val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;
     105             : 
     106             :         /* spadd(val_hi, val_lo, y), except for not using double_t: */
     107           0 :         w = y + val_hi;
     108           0 :         val_lo += (y - w) + val_hi;
     109           0 :         val_hi = w;
     110             : 
     111           0 :         return val_lo + val_hi;
     112             : }