Branch data     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 "cdefs-compat.h"
      15                 :             : //__FBSDID("$FreeBSD: src/lib/msun/src/e_log10.c,v 1.15 2011/10/15 05:23:28 das Exp $");
      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                 :             : #include <openlibm_math.h>
      27                 :             : 
      28                 :             : #include "math_private.h"
      29                 :             : #include "k_log.h"
      30                 :             : 
      31                 :             : static const double
      32                 :             : two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
      33                 :             : ivln10hi   =  4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
      34                 :             : ivln10lo   =  2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
      35                 :             : log10_2hi  =  3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
      36                 :             : log10_2lo  =  3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
      37                 :             : 
      38                 :             : static const double zero   =  0.0;
      39                 :             : 
      40                 :             : OLM_DLLEXPORT double
      41                 :          13 : __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                 :          13 :         EXTRACT_WORDS(hx,lx,x);
      48                 :             : 
      49                 :          13 :         k=0;
      50         [ +  + ]:          13 :         if (hx < 0x00100000) {                       /* x < 2**-1022  */
      51         [ +  + ]:           4 :             if (((hx&0x7fffffff)|lx)==0)
      52                 :           2 :                 return -two54/zero;             /* log(+-0)=-inf */
      53         [ +  - ]:           2 :             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         [ +  + ]:           9 :         if (hx >= 0x7ff00000) return x+x;
      58   [ +  +  +  - ]:           8 :         if (hx == 0x3ff00000 && lx == 0)
      59                 :           2 :             return zero;                        /* log(1) = +0 */
      60                 :           6 :         k += (hx>>20)-1023;
      61                 :           6 :         hx &= 0x000fffff;
      62                 :           6 :         i = (hx+0x95f64)&0x100000;
      63                 :           6 :         SET_HIGH_WORD(x,hx|(i^0x3ff00000));     /* normalize x or x/2 */
      64                 :           6 :         k += (i>>20);
      65                 :           6 :         y = (double)k;
      66                 :           6 :         f = x - 1.0;
      67                 :           6 :         hfsq = 0.5*f*f;
      68                 :           6 :         r = k_log1p(f);
      69                 :             : 
      70                 :             :         /* See e_log2.c for most details. */
      71                 :           6 :         hi = f - hfsq;
      72                 :           6 :         SET_LOW_WORD(hi,0);
      73                 :           6 :         lo = (f - hi) - hfsq + r;
      74                 :           6 :         val_hi = hi*ivln10hi;
      75                 :           6 :         y2 = y*log10_2hi;
      76                 :           6 :         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                 :           6 :         w = y2 + val_hi;
      85                 :           6 :         val_lo += (y2 - w) + val_hi;
      86                 :           6 :         val_hi = w;
      87                 :             : 
      88                 :           6 :         return val_lo + val_hi;
      89                 :             : }
      90                 :             : 
      91                 :             : #if (LDBL_MANT_DIG == 53)
      92                 :             : openlibm_weak_reference(log10, log10l);
      93                 :             : #endif