1 /******************************** -*- C -*- ****************************
2 *
3 * Emulation for expl
4 *
5 *
6 ***********************************************************************/
7
8 /***********************************************************************
9 *
10 * Copyright 2002, 2006 Free Software Foundation, Inc.
11 * Written by Paolo Bonzini.
12 *
13 * This file is part of GNU Smalltalk.
14 *
15 * GNU Smalltalk is free software; you can redistribute it and/or modify it
16 * under the terms of the GNU General Public License as published by the Free
17 * Software Foundation; either version 2, or (at your option) any later
18 * version.
19 *
20 * Linking GNU Smalltalk statically or dynamically with other modules is
21 * making a combined work based on GNU Smalltalk. Thus, the terms and
22 * conditions of the GNU General Public License cover the whole
23 * combination.
24 *
25 * In addition, as a special exception, the Free Software Foundation
26 * give you permission to combine GNU Smalltalk with free software
27 * programs or libraries that are released under the GNU LGPL and with
28 * independent programs running under the GNU Smalltalk virtual machine.
29 *
30 * You may copy and distribute such a system following the terms of the
31 * GNU GPL for GNU Smalltalk and the licenses of the other code
32 * concerned, provided that you include the source code of that other
33 * code when and as the GNU GPL requires distribution of source code.
34 *
35 * Note that people who make modified versions of GNU Smalltalk are not
36 * obligated to grant this special exception for their modified
37 * versions; it is their choice whether to do so. The GNU General
38 * Public License gives permission to release a modified version without
39 * this exception; this exception also makes it possible to release a
40 * modified version which carries forward this exception.
41 *
42 * GNU Smalltalk is distributed in the hope that it will be useful, but WITHOUT
43 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
44 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
45 * more details.
46 *
47 * You should have received a copy of the GNU General Public License along with
48 * GNU Smalltalk; see the file COPYING. If not, write to the Free Software
49 * Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
50 *
51 ***********************************************************************/
52
53 #include <float.h>
54 #include <math.h>
55
56 #include "mathl.h"
57
58 static const long double C[] = {
59 /* Chebyshev polynom coeficients for (exp(x)-1)/x */
60 #define P1 C[0]
61 #define P2 C[1]
62 #define P3 C[2]
63 #define P4 C[3]
64 #define P5 C[4]
65 #define P6 C[5]
66 0.5L,
67 1.66666666666666666666666666666666683E-01L,
68 4.16666666666666666666654902320001674E-02L,
69 8.33333333333333333333314659767198461E-03L,
70 1.38888888889899438565058018857254025E-03L,
71 1.98412698413981650382436541785404286E-04L,
72
73 /* Smallest integer x for which e^x overflows. */
74 #define himark C[6]
75 11356.523406294143949491931077970765L,
76
77 /* Largest integer x for which e^x underflows. */
78 #define lomark C[7]
79 -11433.4627433362978788372438434526231L,
80
81 /* very small number */
82 #define TINY C[8]
83 1.0e-4900L,
84
85 /* 2^16383 */
86 #define TWO16383 C[9]
87 5.94865747678615882542879663314003565E+4931L};
88
89 long double
expl(long double x)90 expl (long double x)
91 {
92 /* Check for usual case. */
93 if (x < himark && x > lomark)
94 {
95 int exponent;
96 long double t, x22;
97 int k = 1;
98 long double result = 1.0;
99
100 /* Compute an integer power of e with a granularity of 0.125. */
101 exponent = (int) floorl (x * 8.0L);
102 x -= exponent / 8.0L;
103
104 if (x > 0.0625)
105 {
106 exponent++;
107 x -= 0.125L;
108 }
109
110 if (exponent < 0)
111 {
112 t = 0.8824969025845954028648921432290507362220L; /* e^-0.25 */
113 exponent = -exponent;
114 }
115 else
116 t = 1.1331484530668263168290072278117938725655L; /* e^0.25 */
117
118 while (exponent)
119 {
120 if (exponent & k)
121 {
122 result *= t;
123 exponent ^= k;
124 }
125 t *= t;
126 k <<= 1;
127 }
128
129 /* Approximate (e^x - 1)/x, using a seventh-degree polynomial,
130 with maximum error in [-2^-16-2^-53,2^-16+2^-53]
131 less than 4.8e-39. */
132 x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
133
134 return result + result * x22;
135 }
136 /* Exceptional cases: */
137 else if (x < himark)
138 {
139 if (x + x == x)
140 /* e^-inf == 0, with no error. */
141 return 0;
142 else
143 /* Underflow */
144 return TINY * TINY;
145 }
146 else
147 /* Return x, if x is a NaN or Inf; or overflow, otherwise. */
148 return TWO16383*x;
149 }
150
151 #if 0
152 int
153 main ()
154 {
155 printf ("%.16Lg\n", expl(1.0L));
156 printf ("%.16Lg\n", expl(-1.0L));
157 printf ("%.16Lg\n", expl(2.0L));
158 printf ("%.16Lg\n", expl(4.0L));
159 printf ("%.16Lg\n", expl(-2.0L));
160 printf ("%.16Lg\n", expl(-4.0L));
161 printf ("%.16Lg\n", expl(0.0625L));
162 printf ("%.16Lg\n", expl(0.3L));
163 printf ("%.16Lg\n", expl(0.6L));
164 }
165 #endif
166