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

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_log2f.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 log2.c.
  14  */
  15
  16 #include "libm.h"
  17 #include "__log1pf.h"
  18
  19 static const float
  20 two25   =  3.3554432000e+07, /* 0x4c000000 */
  21 ivln2hi =  1.4428710938e+00, /* 0x3fb8b000 */
  22 ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
  23
  24 float log2f(float x)
  25 {
  26         float f,hfsq,hi,lo,r,y;
  27         int32_t i,k,hx;
  28
  29         GET_FLOAT_WORD(hx, x);
  30
  31         k = 0;
  32         if (hx < 0x00800000) {  /* x < 2**-126  */
  33                 if ((hx&0x7fffffff) == 0)
  34                         return -two25/0.0f;  /* log(+-0)=-inf */
  35                 if (hx < 0)
  36                         return (x-x)/0.0f;   /* log(-#) = NaN */
  37                 /* subnormal number, scale up x */
  38                 k -= 25;
  39                 x *= two25;
  40                 GET_FLOAT_WORD(hx, x);
  41         }
  42         if (hx >= 0x7f800000)
  43                 return x+x;
  44         if (hx == 0x3f800000)
  45                 return 0.0f;  /* log(1) = +0 */
  46         k += (hx>>23) - 127;
  47         hx &= 0x007fffff;
  48         i = (hx+(0x4afb0d))&0x800000;
  49         SET_FLOAT_WORD(x, hx|(i^0x3f800000));  /* normalize x or x/2 */
  50         k += i>>23;
  51         y = (float)k;
  52         f = x - 1.0f;
  53         hfsq = 0.5f * f * f;
  54         r = __log1pf(f);
  55
  56         /*
  57          * We no longer need to avoid falling into the multi-precision
  58          * calculations due to compiler bugs breaking Dekker's theorem.
  59          * Keep avoiding this as an optimization.  See log2.c for more
  60          * details (some details are here only because the optimization
  61          * is not yet available in double precision).
  62          *
  63          * Another compiler bug turned up.  With gcc on i386,
  64          * (ivln2lo + ivln2hi) would be evaluated in float precision
  65          * despite runtime evaluations using double precision.  So we
  66          * must cast one of its terms to float_t.  This makes the whole
  67          * expression have type float_t, so return is forced to waste
  68          * time clobbering its extra precision.
  69          */
  70 // FIXME
  71 //      if (sizeof(float_t) > sizeof(float))
  72 //              return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;
  73
  74         hi = f - hfsq;
  75         GET_FLOAT_WORD(hx,hi);
  76         SET_FLOAT_WORD(hi,hx&0xfffff000);
  77         lo = (f - hi) - hfsq + r;
  78         return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y;
  79 }