lib/log2l.c

   1 /* Base 2 logarithm.
   2    Copyright (C) 2011-2017 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 <https://www.gnu.org/licenses/>.  */
  16
  17 #include <config.h>
  18
  19 /* Specification.  */
  20 #include <math.h>
  21
  22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
  23
  24 long double
  25 log2l (long double x)
  26 {
  27   return log2 (x);
  28 }
  29
  30 #else
  31
  32 /* Best possible approximation of log(2) as a 'long double'.  */
  33 #define LOG2 0.693147180559945309417232121458176568075L
  34
  35 /* Best possible approximation of 1/log(2) as a 'long double'.  */
  36 #define LOG2_INVERSE 1.44269504088896340735992468100189213743L
  37
  38 /* sqrt(0.5).  */
  39 #define SQRT_HALF 0.707106781186547524400844362104849039284L
  40
  41 long double
  42 log2l (long double x)
  43 {
  44   if (isnanl (x))
  45     return x;
  46
  47   if (x <= 0.0L)
  48     {
  49       if (x == 0.0L)
  50         /* Return -Infinity.  */
  51         return - HUGE_VALL;
  52       else
  53         {
  54           /* Return NaN.  */
  55 #if defined _MSC_VER || (defined __sgi && !defined __GNUC__)
  56           static long double zero;
  57           return zero / zero;
  58 #else
  59           return 0.0L / 0.0L;
  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     long double y;
  73
  74     y = frexpl (x, &e);
  75     if (y < SQRT_HALF)
  76       {
  77         y = 2.0L * y;
  78         e = e - 1;
  79       }
  80
  81     return (long double) e + logl (y) * LOG2_INVERSE;
  82   }
  83 }
  84
  85 #endif