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