compiler/stdc/math/e_log2f.c

   1 /* e_log2f.c -- float version of e_log2.c.
   2  */
   3
   4 /*
   5  * ====================================================
   6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   7  *
   8  * Developed at SunPro, a Sun Microsystems, Inc. business.
   9  * Permission to use, copy, modify, and distribute this
  10  * software is freely granted, provided that this notice
  11  * is preserved.
  12  * ====================================================
  13  */
  14
  15 #ifndef lint
  16 static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_log2f.c,v 1.5 2011/10/15 05:23:28 das Exp $";
  17 #endif
  18
  19 #include "math.h"
  20 #include "math_private.h"
  21
  22 static const float
  23 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
  24 Lg1 =      0xaaaaaa.0p-24,      /* 0.66666662693 */
  25 Lg2 =      0xccce13.0p-25,      /* 0.40000972152 */
  26 Lg3 =      0x91e9ee.0p-25,      /* 0.28498786688 */
  27 Lg4 =      0xf89e26.0p-26;      /* 0.24279078841 */
  28
  29 static inline float
  30 k_log1pf(float f)
  31 {
  32         float hfsq,s,z,R,w,t1,t2;
  33
  34         s = f/((float)2.0+f);
  35         z = s*s;
  36         w = z*z;
  37         t1= w*(Lg2+w*Lg4);
  38         t2= z*(Lg1+w*Lg3);
  39         R = t2+t1;
  40         hfsq=(float)0.5*f*f;
  41         return s*(hfsq+R);
  42 }
  43
  44 static const float
  45 two25      =  3.3554432000e+07, /* 0x4c000000 */
  46 ivln2hi    =  1.4428710938e+00, /* 0x3fb8b000 */
  47 ivln2lo    = -1.7605285393e-04; /* 0xb9389ad4 */
  48
  49 static const float zero   =  0.0;
  50
  51 float
  52 __ieee754_log2f(float x)
  53 {
  54         float f,hfsq,hi,lo,r,y;
  55         int32_t i,k,hx;
  56
  57         GET_FLOAT_WORD(hx,x);
  58
  59         k=0;
  60         if (hx < 0x00800000) {                  /* x < 2**-126  */
  61             if ((hx&0x7fffffff)==0)
  62                 return -two25/zero;             /* log(+-0)=-inf */
  63             if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
  64             k -= 25; x *= two25; /* subnormal number, scale up x */
  65             GET_FLOAT_WORD(hx,x);
  66         }
  67         if (hx >= 0x7f800000) return x+x;
  68         if (hx == 0x3f800000)
  69             return zero;                        /* log(1) = +0 */
  70         k += (hx>>23)-127;
  71         hx &= 0x007fffff;
  72         i = (hx+(0x4afb0d))&0x800000;
  73         SET_FLOAT_WORD(x,hx|(i^0x3f800000));    /* normalize x or x/2 */
  74         k += (i>>23);
  75         y = (float)k;
  76         f = x - (float)1.0;
  77         hfsq = (float)0.5*f*f;
  78         r = k_log1pf(f);
  79
  80         /*
  81          * We no longer need to avoid falling into the multi-precision
  82          * calculations due to compiler bugs breaking Dekker's theorem.
  83          * Keep avoiding this as an optimization.  See e_log2.c for more
  84          * details (some details are here only because the optimization
  85          * is not yet available in double precision).
  86          *
  87          * Another compiler bug turned up.  With gcc on i386,
  88          * (ivln2lo + ivln2hi) would be evaluated in float precision
  89          * despite runtime evaluations using double precision.  So we
  90          * must cast one of its terms to float_t.  This makes the whole
  91          * expression have type float_t, so return is forced to waste
  92          * time clobbering its extra precision.
  93          */
  94         if (sizeof(float_t) > sizeof(float))
  95                 return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;
  96
  97         hi = f - hfsq;
  98         GET_FLOAT_WORD(hx,hi);
  99         SET_FLOAT_WORD(hi,hx&0xfffff000);
 100         lo = (f - hi) - hfsq + r;
 101         return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y;
 102 }