nsz Git - musl/blob - src/math/log10f.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log10f.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /*
  13  * See comments in log10.c.
  14  */
  15
  16 #include "libm.h"
  17 #include "__log1pf.h"
  18
  19 static const float
  20 two25     =  3.3554432000e+07, /* 0x4c000000 */
  21 ivln10hi  =  4.3432617188e-01, /* 0x3ede6000 */
  22 ivln10lo  = -3.1689971365e-05, /* 0xb804ead9 */
  23 log10_2hi =  3.0102920532e-01, /* 0x3e9a2080 */
  24 log10_2lo =  7.9034151668e-07; /* 0x355427db */
  25
  26 static const float zero = 0.0;
  27
  28 float log10f(float x)
  29 {
  30         float f,hfsq,hi,lo,r,y;
  31         int32_t i,k,hx;
  32
  33         GET_FLOAT_WORD(hx, x);
  34
  35         k = 0;
  36         if (hx < 0x00800000) {  /* x < 2**-126  */
  37                 if ((hx&0x7fffffff) == 0)
  38                         return -two25/zero;  /* log(+-0)=-inf */
  39                 if (hx < 0)
  40                         return (x-x)/zero;   /* log(-#) = NaN */
  41                 /* subnormal number, scale up x */
  42                 k -= 25;
  43                 x *= two25;
  44                 GET_FLOAT_WORD(hx, x);
  45         }
  46         if (hx >= 0x7f800000)
  47                 return x+x;
  48         if (hx == 0x3f800000)
  49                 return zero;  /* log(1) = +0 */
  50         k += (hx>>23) - 127;
  51         hx &= 0x007fffff;
  52         i = (hx+(0x4afb0d))&0x800000;
  53         SET_FLOAT_WORD(x, hx|(i^0x3f800000));  /* normalize x or x/2 */
  54         k += i>>23;
  55         y = (float)k;
  56         f = x - (float)1.0;
  57         hfsq = (float)0.5*f*f;
  58         r = __log1pf(f);
  59
  60 // FIXME
  61 //      /* See log2f.c and log2.c for details. */
  62 //      if (sizeof(float_t) > sizeof(float))
  63 //              return (r - hfsq + f) * ((float_t)ivln10lo + ivln10hi) +
  64 //                  y * ((float_t)log10_2lo + log10_2hi);
  65         hi = f - hfsq;
  66         GET_FLOAT_WORD(hx, hi);
  67         SET_FLOAT_WORD(hi, hx&0xfffff000);
  68         lo = (f - hi) - hfsq + r;
  69         return y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi +
  70                 hi*ivln10hi + y*log10_2hi;
  71 }