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 10 logarithm of x.  See e_log.c and k_log.h for most
      19             :  * comments.
      20             :  *
      21             :  *    log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2)
      22             :  * in not-quite-routine extra precision.
      23             :  */
      24             : 
      25             : #include <float.h>
      26             : 
      27             : #include "math_private.h"
      28             : #include "k_log.h"
      29             :  
      30             : static const double
      31             : two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
      32             : ivln10hi   =  4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
      33             : ivln10lo   =  2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
      34             : log10_2hi  =  3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
      35             : log10_2lo  =  3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
      36             : 
      37             : static const double zero   =  0.0;
      38             : static volatile double vzero = 0.0;
      39             : 
      40             : double
      41           0 : __ieee754_log10(double x)
      42             : {
      43             :         double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
      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             :         /* See e_log2.c for most details. */
      71           0 :         hi = f - hfsq;
      72           0 :         SET_LOW_WORD(hi,0);
      73           0 :         lo = (f - hi) - hfsq + r;
      74           0 :         val_hi = hi*ivln10hi;
      75           0 :         y2 = y*log10_2hi;
      76           0 :         val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi;
      77             : 
      78             :         /*
      79             :          * Extra precision in for adding y*log10_2hi is not strictly needed
      80             :          * since there is no very large cancellation near x = sqrt(2) or
      81             :          * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
      82             :          * with some parallelism and it reduces the error for many args.
      83             :          */
      84           0 :         w = y2 + val_hi;
      85           0 :         val_lo += (y2 - w) + val_hi;
      86           0 :         val_hi = w;
      87             : 
      88           0 :         return val_lo + val_hi;
      89             : }