compiler/stdc/math/e_log2.c

   1
   2 /* @(#)e_log2.c 1.4 11/10/15 */
   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 #ifndef lint
  15 static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_log2.c,v 1.4 2011/10/15 05:23:28 das Exp $";
  16 #endif
  17
  18 #include "math.h"
  19 #include "math_private.h"
  20
  21 static const double
  22 Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
  23 Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
  24 Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
  25 Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
  26 Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
  27 Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
  28 Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
  29
  30 /*
  31  * We always inline k_log1p(), since doing so produces a
  32  * substantial performance improvement (~40% on amd64).
  33  */
  34 static inline double
  35 k_log1p(double f)
  36 {
  37         double hfsq,s,z,R,w,t1,t2;
  38
  39         s = f/(2.0+f);
  40         z = s*s;
  41         w = z*z;
  42         t1= w*(Lg2+w*(Lg4+w*Lg6));
  43         t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
  44         R = t2+t1;
  45         hfsq=0.5*f*f;
  46         return s*(hfsq+R);
  47 }
  48
  49 static const double
  50 two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  51 ivln2hi    =  1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
  52 ivln2lo    =  1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
  53
  54 static const double zero   =  0.0;
  55
  56 double
  57 __ieee754_log2(double x)
  58 {
  59         double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
  60         int32_t i,k,hx;
  61         uint32_t lx;
  62
  63         EXTRACT_WORDS(hx,lx,x);
  64
  65         k=0;
  66         if (hx < 0x00100000) {                  /* x < 2**-1022  */
  67             if (((hx&0x7fffffff)|lx)==0)
  68                 return -two54/zero;             /* log(+-0)=-inf */
  69             if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
  70             k -= 54; x *= two54; /* subnormal number, scale up x */
  71             GET_HIGH_WORD(hx,x);
  72         }
  73         if (hx >= 0x7ff00000) return x+x;
  74         if (hx == 0x3ff00000 && lx == 0)
  75             return zero;                        /* log(1) = +0 */
  76         k += (hx>>20)-1023;
  77         hx &= 0x000fffff;
  78         i = (hx+0x95f64)&0x100000;
  79         SET_HIGH_WORD(x,hx|(i^0x3ff00000));     /* normalize x or x/2 */
  80         k += (i>>20);
  81         y = (double)k;
  82         f = x - 1.0;
  83         hfsq = 0.5*f*f;
  84         r = k_log1p(f);
  85
  86         /*
  87          * f-hfsq must (for args near 1) be evaluated in extra precision
  88          * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
  89          * This is fairly efficient since f-hfsq only depends on f, so can
  90          * be evaluated in parallel with R.  Not combining hfsq with R also
  91          * keeps R small (though not as small as a true `lo' term would be),
  92          * so that extra precision is not needed for terms involving R.
  93          *
  94          * Compiler bugs involving extra precision used to break Dekker's
  95          * theorem for spitting f-hfsq as hi+lo, unless double_t was used
  96          * or the multi-precision calculations were avoided when double_t
  97          * has extra precision.  These problems are now automatically
  98          * avoided as a side effect of the optimization of combining the
  99          * Dekker splitting step with the clear-low-bits step.
 100          *
 101          * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
 102          * precision to avoid a very large cancellation when x is very near
 103          * these values.  Unlike the above cancellations, this problem is
 104          * specific to base 2.  It is strange that adding +-1 is so much
 105          * harder than adding +-ln2 or +-log10_2.
 106          *
 107          * This uses Dekker's theorem to normalize y+val_hi, so the
 108          * compiler bugs are back in some configurations, sigh.  And I
 109          * don't want to used double_t to avoid them, since that gives a
 110          * pessimization and the support for avoiding the pessimization
 111          * is not yet available.
 112          *
 113          * The multi-precision calculations for the multiplications are
 114          * routine.
 115          */
 116         hi = f - hfsq;
 117         SET_LOW_WORD(hi,0);
 118         lo = (f - hi) - hfsq + r;
 119         val_hi = hi*ivln2hi;
 120         val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;
 121
 122         /* spadd(val_hi, val_lo, y), except for not using double_t: */
 123         w = y + val_hi;
 124         val_lo += (y - w) + val_hi;
 125         val_hi = w;
 126
 127         return val_lo + val_hi;
 128 }