Bump to 1.14.1
[platform/upstream/augeas.git] / lib / log2l.c
1 /* Base 2 logarithm.
2    Copyright (C) 2011-2016 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_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