lib/msun/src/e_log10.c

   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
  22 #include "math.h"
  23 #include "math_private.h"
  24 #include "k_log.h"
  25
  26 static const double
  27 two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  28 ivln10hi   =  4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
  29 ivln10lo   =  2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
  30 log10_2hi  =  3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
  31 log10_2lo  =  3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
  32
  33 static const double zero   =  0.0;
  34
  35 double
  36 __ieee754_log10(double x)
  37 {
  38         double f,hi,lo,y,z;
  39         int32_t i,k,hx;
  40         u_int32_t lx;
  41
  42         EXTRACT_WORDS(hx,lx,x);
  43
  44         k=0;
  45         if (hx < 0x00100000) {                  /* x < 2**-1022  */
  46             if (((hx&0x7fffffff)|lx)==0)
  47                 return -two54/zero;             /* log(+-0)=-inf */
  48             if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
  49             k -= 54; x *= two54; /* subnormal number, scale up x */
  50             GET_HIGH_WORD(hx,x);
  51         }
  52         if (hx >= 0x7ff00000) return x+x;
  53         k += (hx>>20)-1023;
  54         hx &= 0x000fffff;
  55         i = (hx+0x95f64)&0x100000;
  56         SET_HIGH_WORD(x,hx|(i^0x3ff00000));     /* normalize x or x/2 */
  57         k += (i>>20);
  58         y = (double)k;
  59         f = __kernel_log(x);
  60         hi = x = x - 1;
  61         SET_LOW_WORD(hi,0);
  62         lo = x - hi;
  63         z = y*log10_2lo + (x+f)*ivln10lo + (lo+f)*ivln10hi + hi*ivln10hi;
  64         return  z+y*log10_2hi;
  65 }