lib/log2.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 /* Best possible approximation of log(2) as a 'double'.  */
  23 #define LOG2 0.693147180559945309417232121458176568075
  24
  25 /* Best possible approximation of 1/log(2) as a 'double'.  */
  26 #define LOG2_INVERSE 1.44269504088896340735992468100189213743
  27
  28 /* sqrt(0.5).  */
  29 #define SQRT_HALF 0.707106781186547524400844362104849039284
  30
  31 double
  32 log2 (double x)
  33 {
  34   if (isnand (x))
  35     return x;
  36
  37   if (x <= 0.0)
  38     {
  39       if (x == 0.0)
  40         /* Return -Infinity.  */
  41         return - HUGE_VAL;
  42       else
  43         {
  44           /* Return NaN.  */
  45 #if defined _MSC_VER || (defined __sgi && !defined __GNUC__)
  46           static double zero;
  47           return zero / zero;
  48 #else
  49           return 0.0 / 0.0;
  50 #endif
  51         }
  52     }
  53
  54   /* Decompose x into
  55        x = 2^e * y
  56      where
  57        e is an integer,
  58        1/2 < y < 2.
  59      Then log2(x) = e + log2(y) = e + log(y)/log(2).  */
  60   {
  61     int e;
  62     double y;
  63
  64     y = frexp (x, &e);
  65     if (y < SQRT_HALF)
  66       {
  67         y = 2.0 * y;
  68         e = e - 1;
  69       }
  70
  71     return (double) e + log (y) * LOG2_INVERSE;
  72   }
  73 }