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