Mercurial > gnulib

annotate lib/expl.c @ 40186:8964917f9574

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
autoupdate
author Karl Berry <karl@freefriends.org>
date Mon, 18 Feb 2019 08:02:49 -0800
parents b06060465f09
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
1 /* Exponential function.
40057
b06060465f09 maint: Run 'make update-copyright'
Paul Eggert <eggert@cs.ucla.edu>
parents: 19484
diff changeset
2 Copyright (C) 2011-2019 Free Software Foundation, Inc.
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
3
9309
bbbbbf4cd1c5 Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents: 8186
diff changeset
4 This program is free software: you can redistribute it and/or modify
4435
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
5 it under the terms of the GNU General Public License as published by
9309
bbbbbf4cd1c5 Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents: 8186
diff changeset
6 the Free Software Foundation; either version 3 of the License, or
bbbbbf4cd1c5 Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents: 8186
diff changeset
7 (at your option) any later version.
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
8
4435
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
9 This program is distributed in the hope that it will be useful,
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
d4efdb886b2a Switch from LGPL to GPL.
Paul Eggert <eggert@cs.ucla.edu>
parents: 4239
diff changeset
12 GNU General Public License for more details.
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
13
9309
bbbbbf4cd1c5 Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents: 8186
diff changeset
14 You should have received a copy of the GNU General Public License
19190
9759915b2aca all: prefer https: URLs
Paul Eggert <eggert@cs.ucla.edu>
parents: 18626
diff changeset
15 along with this program. If not, see <https://www.gnu.org/licenses/>. */
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
16
8186
c115624bf99a New module 'math'. <math.h> replaces mathl.h.
Bruno Haible <bruno@clisp.org>
parents: 5848
diff changeset
17 #include <config.h>
c115624bf99a New module 'math'. <math.h> replaces mathl.h.
Bruno Haible <bruno@clisp.org>
parents: 5848
diff changeset
18
c115624bf99a New module 'math'. <math.h> replaces mathl.h.
Bruno Haible <bruno@clisp.org>
parents: 5848
diff changeset
19 /* Specification. */
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
20 #include <math.h>
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
21
15937
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
23
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
24 long double
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
25 expl (long double x)
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
26 {
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
27 return exp (x);
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
28 }
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
29
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
30 #else
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
31
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
32 # include <float.h>
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
33
16620
4a578c3945b8 New module 'exp2l'.
Bruno Haible <bruno@clisp.org>
parents: 16578
diff changeset
34 /* gl_expl_table[i] = exp((i - 128) * log(2)/256). */
4a578c3945b8 New module 'exp2l'.
Bruno Haible <bruno@clisp.org>
parents: 16578
diff changeset
35 extern const long double gl_expl_table[257];
4a578c3945b8 New module 'exp2l'.
Bruno Haible <bruno@clisp.org>
parents: 16578
diff changeset
36
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
37 /* A value slightly larger than log(2). */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
38 #define LOG2_PLUS_EPSILON 0.6931471805599454L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
39
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
40 /* Best possible approximation of log(2) as a 'long double'. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
41 #define LOG2 0.693147180559945309417232121458176568075L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
42
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
43 /* Best possible approximation of 1/log(2) as a 'long double'. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
44 #define LOG2_INVERSE 1.44269504088896340735992468100189213743L
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
45
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
46 /* Best possible approximation of log(2)/256 as a 'long double'. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
47 #define LOG2_BY_256 0.00270760617406228636491106297444600221904L
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
48
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
49 /* Best possible approximation of 256/log(2) as a 'long double'. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
50 #define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
51
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
52 /* The upper 32 bits of log(2)/256. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
53 #define LOG2_BY_256_HI_PART 0.0027076061733168899081647396087646484375L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
54 /* log(2)/256 - LOG2_HI_PART. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
55 #define LOG2_BY_256_LO_PART \
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
56 0.000000000000745396456746323365681353781544922399845L
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
57
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
58 long double
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
59 expl (long double x)
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
60 {
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
61 if (isnanl (x))
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
62 return x;
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
63
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
64 if (x >= (long double) LDBL_MAX_EXP * LOG2_PLUS_EPSILON)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
65 /* x > LDBL_MAX_EXP * log(2)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
66 hence exp(x) > 2^LDBL_MAX_EXP, overflows to Infinity. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
67 return HUGE_VALL;
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
68
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
69 if (x <= (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * LOG2_PLUS_EPSILON)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
70 /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * log(2)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
71 hence exp(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG),
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
72 underflows to zero. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
73 return 0.0L;
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
74
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
75 /* Decompose x into
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
76 x = n * log(2) + m * log(2)/256 + y
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
77 where
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
78 n is an integer,
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
79 m is an integer, -128 <= m <= 128,
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
80 y is a number, |y| <= log(2)/512 + epsilon = 0.00135...
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
81 Then
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
82 exp(x) = 2^n * exp(m * log(2)/256) * exp(y)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
83 The first factor is an ldexpl() call.
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
84 The second factor is a table lookup.
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
85 The third factor is computed
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
86 - either as sinh(y) + cosh(y)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
87 where sinh(y) is computed through the power series:
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
88 sinh(y) = y + y^3/3! + y^5/5! + ...
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
89 and cosh(y) is computed as hypot(1, sinh(y)),
16578
5626e55f7f0b expl: Simplify computation.
Bruno Haible <bruno@clisp.org>
parents: 16576
diff changeset
90 - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z))
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
91 where z = y/2
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
92 and tanh(z) is computed through its power series:
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
93 tanh(z) = z
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
94 - 1/3 * z^3
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
95 + 2/15 * z^5
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
96 - 17/315 * z^7
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
97 + 62/2835 * z^9
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
98 - 1382/155925 * z^11
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
99 + 21844/6081075 * z^13
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
100 - 929569/638512875 * z^15
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
101 + ...
16651
6c0da1a4068d Fix some comments.
Bruno Haible <bruno@clisp.org>
parents: 16620
diff changeset
102 Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the
6c0da1a4068d Fix some comments.
Bruno Haible <bruno@clisp.org>
parents: 16620
diff changeset
103 z^13 term is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we
6c0da1a4068d Fix some comments.
Bruno Haible <bruno@clisp.org>
parents: 16620
diff changeset
104 can truncate the series after the z^11 term.
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
105
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
106 Given the usual bounds LDBL_MAX_EXP <= 16384, LDBL_MIN_EXP >= -16381,
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
107 LDBL_MANT_DIG <= 120, we can estimate x: -11440 <= x <= 11357.
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
108 This means, when dividing x by log(2), where we want x mod log(2)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
109 to be precise to LDBL_MANT_DIG bits, we have to use an approximation
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
110 to log(2) that has 14+LDBL_MANT_DIG bits. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
111
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
112 {
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
113 long double nm = roundl (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
114 /* n has at most 15 bits, nm therefore has at most 23 bits, therefore
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
115 n * LOG2_HI_PART is computed exactly, and n * LOG2_LO_PART is computed
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
116 with an absolute error < 2^15 * 2e-10 * 2^-LDBL_MANT_DIG. */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
117 long double y_tmp = x - nm * LOG2_BY_256_HI_PART;
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
118 long double y = y_tmp - nm * LOG2_BY_256_LO_PART;
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
119 long double z = 0.5L * y;
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
120
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
121 /* Coefficients of the power series for tanh(z). */
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
122 #define TANH_COEFF_1 1.0L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
123 #define TANH_COEFF_3 -0.333333333333333333333333333333333333334L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
124 #define TANH_COEFF_5 0.133333333333333333333333333333333333334L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
125 #define TANH_COEFF_7 -0.053968253968253968253968253968253968254L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
126 #define TANH_COEFF_9 0.0218694885361552028218694885361552028218L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
127 #define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
128 #define TANH_COEFF_13 0.00359212803657248101692546136990581435026L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
129 #define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
130
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
131 long double z2 = z * z;
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
132 long double tanh_z =
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
133 (((((TANH_COEFF_11
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
134 * z2 + TANH_COEFF_9)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
135 * z2 + TANH_COEFF_7)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
136 * z2 + TANH_COEFF_5)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
137 * z2 + TANH_COEFF_3)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
138 * z2 + TANH_COEFF_1)
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
139 * z;
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
140
16578
5626e55f7f0b expl: Simplify computation.
Bruno Haible <bruno@clisp.org>
parents: 16576
diff changeset
141 long double exp_y = (1.0L + tanh_z) / (1.0L - tanh_z);
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
142
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
143 int n = (int) roundl (nm * (1.0L / 256.0L));
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
144 int m = (int) nm - 256 * n;
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
145
16620
4a578c3945b8 New module 'exp2l'.
Bruno Haible <bruno@clisp.org>
parents: 16578
diff changeset
146 return ldexpl (gl_expl_table[128 + m] * exp_y, n);
16576
6c7c671909d3 expl: Fix precision of computed result.
Bruno Haible <bruno@clisp.org>
parents: 16511
diff changeset
147 }
4239
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
148 }
b33b8603715c Transcendental functions for 'long double', from Paolo Bonzini.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
149
15937
16cba28a5e86 expl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents: 14079
diff changeset
150 #endif