lib/log2f.c

   1 /* Base 2 logarithm.
   2    Copyright (C) 2012-2014 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
  16
  17 #include <config.h>
  18
  19 /* Specification.  */
  20 #include <math.h>
  21
  22 #if HAVE_LOG2
  23
  24 float
  25 log2f (float x)
  26 {
  27   return (float) log2 ((double) x);
  28 }
  29
  30 #else
  31
  32 /* Best possible approximation of log(2) as a 'float'.  */
  33 #define LOG2 0.693147180559945309417232121458176568075f
  34
  35 /* Best possible approximation of 1/log(2) as a 'float'.  */
  36 #define LOG2_INVERSE 1.44269504088896340735992468100189213743f
  37
  38 /* sqrt(0.5).  */
  39 #define SQRT_HALF 0.707106781186547524400844362104849039284f
  40
  41 float
  42 log2f (float x)
  43 {
  44   if (isnanf (x))
  45     return x;
  46
  47   if (x <= 0.0f)
  48     {
  49       if (x == 0.0f)
  50         /* Return -Infinity.  */
  51         return - HUGE_VALF;
  52       else
  53         {
  54           /* Return NaN.  */
  55 #if defined _MSC_VER
  56           static float zero;
  57           return zero / zero;
  58 #else
  59           return 0.0f / 0.0f;
  60 #endif
  61         }
  62     }
  63
  64   /* Decompose x into
  65        x = 2^e * y
  66      where
  67        e is an integer,
  68        1/2 < y < 2.
  69      Then log2(x) = e + log2(y) = e + log(y)/log(2).  */
  70   {
  71     int e;
  72     float y;
  73
  74     y = frexpf (x, &e);
  75     if (y < SQRT_HALF)
  76       {
  77         y = 2.0f * y;
  78         e = e - 1;
  79       }
  80
  81     return (float) e + logf (y) * LOG2_INVERSE;
  82   }
  83 }
  84
  85 #endif