summaryrefslogtreecommitdiff
path: root/lib/exp2l.c
blob: 14da3db049889d791f7538ad449f65f93d9c52de (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
/* Exponential base 2 function.
   Copyright (C) 2011-2021 Free Software Foundation, Inc.

   This file is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as
   published by the Free Software Foundation; either version 3 of the
   License, or (at your option) any later version.

   This file is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.  */

#include <config.h>

/* Specification.  */
#include <math.h>

#if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE

long double
exp2l (long double x)
{
  return exp2 (x);
}

#else

# include <float.h>

/* gl_expl_table[i] = exp((i - 128) * log(2)/256).  */
extern const long double gl_expl_table[257];

/* Best possible approximation of log(2) as a 'long double'.  */
#define LOG2 0.693147180559945309417232121458176568075L

/* Best possible approximation of 1/log(2) as a 'long double'.  */
#define LOG2_INVERSE 1.44269504088896340735992468100189213743L

/* Best possible approximation of log(2)/256 as a 'long double'.  */
#define LOG2_BY_256 0.00270760617406228636491106297444600221904L

/* Best possible approximation of 256/log(2) as a 'long double'.  */
#define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L

long double
exp2l (long double x)
{
  /* exp2(x) = exp(x*log(2)).
     If we would compute it like this, there would be rounding errors for
     integer or near-integer values of x.  To avoid these, we inline the
     algorithm for exp(), and the multiplication with log(2) cancels a
     division by log(2).  */

  if (isnanl (x))
    return x;

  if (x > (long double) LDBL_MAX_EXP)
    /* x > LDBL_MAX_EXP
       hence exp2(x) > 2^LDBL_MAX_EXP, overflows to Infinity.  */
    return HUGE_VALL;

  if (x < (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG))
    /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG)
       hence exp2(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG),
       underflows to zero.  */
    return 0.0L;

  /* Decompose x into
       x = n + m/256 + y/log(2)
     where
       n is an integer,
       m is an integer, -128 <= m <= 128,
       y is a number, |y| <= log(2)/512 + epsilon = 0.00135...
     Then
       exp2(x) = 2^n * exp(m * log(2)/256) * exp(y)
     The first factor is an ldexpl() call.
     The second factor is a table lookup.
     The third factor is computed
     - either as sinh(y) + cosh(y)
       where sinh(y) is computed through the power series:
         sinh(y) = y + y^3/3! + y^5/5! + ...
       and cosh(y) is computed as hypot(1, sinh(y)),
     - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z))
       where z = y/2
       and tanh(z) is computed through its power series:
         tanh(z) = z
                   - 1/3 * z^3
                   + 2/15 * z^5
                   - 17/315 * z^7
                   + 62/2835 * z^9
                   - 1382/155925 * z^11
                   + 21844/6081075 * z^13
                   - 929569/638512875 * z^15
                   + ...
       Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the
       z^13 term is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we
       can truncate the series after the z^11 term.  */

  {
    long double nm = roundl (x * 256.0L); /* = 256 * n + m */
    long double z = (x * 256.0L - nm) * (LOG2_BY_256 * 0.5L);

/* Coefficients of the power series for tanh(z).  */
#define TANH_COEFF_1   1.0L
#define TANH_COEFF_3  -0.333333333333333333333333333333333333334L
#define TANH_COEFF_5   0.133333333333333333333333333333333333334L
#define TANH_COEFF_7  -0.053968253968253968253968253968253968254L
#define TANH_COEFF_9   0.0218694885361552028218694885361552028218L
#define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L
#define TANH_COEFF_13  0.00359212803657248101692546136990581435026L
#define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L

    long double z2 = z * z;
    long double tanh_z =
      (((((TANH_COEFF_11
           * z2 + TANH_COEFF_9)
          * z2 + TANH_COEFF_7)
         * z2 + TANH_COEFF_5)
        * z2 + TANH_COEFF_3)
       * z2 + TANH_COEFF_1)
      * z;

    long double exp_y = (1.0L + tanh_z) / (1.0L - tanh_z);

    int n = (int) roundl (nm * (1.0L / 256.0L));
    int m = (int) nm - 256 * n;

    return ldexpl (gl_expl_table[128 + m] * exp_y, n);
  }
}

#endif