1 /* floatexp.c: exponential function and friends, based on floatnum. */
2 /*
3 Copyright (C) 2007, 2008 Wolf Lammen.
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License , or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; see the file COPYING. If not, write to:
17
18 The Free Software Foundation, Inc.
19 59 Temple Place, Suite 330
20 Boston, MA 02111-1307 USA.
21
22
23 You may contact the author by:
24 e-mail: ookami1 <at> gmx <dot> de
25 mail: Wolf Lammen
26 Oertzweg 45
27 22307 Hamburg
28 Germany
29
30 *************************************************************************/
31 #include "floatconst.h"
32 #include "floatcommon.h"
33 #include "floatseries.h"
34 #include "floatexp.h"
35
36 /* uses the addition theorem
37 cosh(2x)-1 == 2*(cosh x - 1)*(cosh x + 1)
38 to reduce the argument to the range |x| < 0.01.
39 Starting with x == 1, you need 7 reduction steps
40 to achieve the desired magnitude.
41 The relative error is < 8e-100 for a 100 digit result.
42 The return value is 0, if the result underflows.
43 |x| < 1, otherwise the final process, where
44 the reductions are unwinded, becomes too
45 unstable */
46 char
_coshminus1lt1(floatnum x,int digits)47 _coshminus1lt1(
48 floatnum x,
49 int digits)
50 {
51 floatstruct tmp;
52 int reductions;
53
54 if (float_iszero(x))
55 return 1;
56 float_abs(x);
57 reductions = 0;
58 while(float_getexponent(x) >= -2)
59 {
60 float_mul(x, x, &c1Div2, digits+1);
61 ++reductions;
62 }
63 if (!coshminus1near0(x, digits) && reductions == 0)
64 return !float_iszero(x);
65 float_create(&tmp);
66 for(; reductions-- > 0;)
67 {
68 float_mul(&tmp, x, x, digits);
69 float_add(x, x, x, digits+2);
70 float_add(x, x, &tmp, digits+2);
71 float_add(x, x, x, digits+2);
72 }
73 float_free(&tmp);
74 return 1;
75 }
76
77 /* sinh x == sqrt((cosh x - 1) * (cosh x + 1)) */
78 static void
_sinhfromcoshminus1(floatnum x,int digits)79 _sinhfromcoshminus1(
80 floatnum x,
81 int digits)
82 {
83 floatstruct tmp;
84
85 float_create(&tmp);
86 float_add(&tmp, x, &c2, digits);
87 float_mul(x, &tmp, x, digits+1);
88 float_sqrt(x, digits+1);
89 float_free(&tmp);
90 }
91
92 /* sinh x for |x| < 1. Derived from cosh x - 1.
93 The relative error is < 8e-100 for a 100 digit result */
94 void
_sinhlt1(floatnum x,int digits)95 _sinhlt1(
96 floatnum x,
97 int digits)
98 {
99 signed char sgn;
100 if (float_getexponent(x) < -digits)
101 /* for very small x: sinh(x) approx. == x. */
102 return;
103 sgn = float_getsign(x);
104 _coshminus1lt1(x, digits);
105 _sinhfromcoshminus1(x, digits);
106 float_setsign(x, sgn);
107 }
108
109 /* evaluates exp(x) - 1. This value can be obtained
110 by exp(x) - 1 == sinh(x) + cosh(x) - 1
111 relative error < 8e-100 for a 100 digit result */
112 void
_expminus1lt1(floatnum x,int digits)113 _expminus1lt1(
114 floatnum x,
115 int digits)
116 {
117 floatstruct tmp;
118 signed char sgn;
119
120 if (float_getexponent(x) < -digits || float_iszero(x))
121 /* for very small x: exp(x)-1 approx.== x */
122 return;
123 float_create(&tmp);
124 sgn = float_getsign(x);
125 _coshminus1lt1(x, digits);
126 float_copy(&tmp, x, EXACT);
127 _sinhfromcoshminus1(x, digits);
128 float_setsign(x, sgn);
129 float_add(x, x, &tmp, digits+1);
130 float_free(&tmp);
131 }
132
133 /* exp(x) for 0 <= x < ln 10
134 relative error < 5e-100 */
135 void
_expltln10(floatnum x,int digits)136 _expltln10(
137 floatnum x,
138 int digits)
139 {
140 int expx;
141 int factor;
142 char sgnf;
143
144 expx = float_getexponent(x);
145 factor = 1;
146 if (expx >= -1)
147 {
148 sgnf = leadingdigits(x, 2 + expx);
149 if (sgnf > 4)
150 {
151 if (sgnf < 9)
152 {
153 factor = 2;
154 float_sub(x, x, &cLn2, digits+1);
155 }
156 else if (sgnf < 14)
157 {
158 factor = 3;
159 float_sub(x, x, &cLn3, digits+1);
160 }
161 else if (sgnf < 21)
162 {
163 factor = 7;
164 float_sub(x, x, &cLn7, digits+1);
165 }
166 else
167 {
168 factor = 10;
169 float_sub(x, x, &cLn10, digits+1);
170 }
171 }
172 }
173 _expminus1lt1(x, digits);
174 float_add(x, x, &c1, digits+1);
175 if (factor != 1)
176 float_muli(x, x, factor, digits+1);
177 }
178
179 /* exp(x) for all x. Underflow or overflow is
180 indicated by the return value (0, if error)
181 relative error for 100 digit results is 5e-100 */
182 char
_exp(floatnum x,int digits)183 _exp(
184 floatnum x,
185 int digits)
186 {
187 floatstruct exp, tmp;
188 int expx, extra;
189 char ok;
190
191 if (float_iszero(x))
192 {
193 float_copy(x, &c1, EXACT);
194 return 1;
195 }
196 expx = float_getexponent(x);
197 if (expx >= (int)(BITS_IN_EXP >> 1))
198 /* obvious overflow or underflow */
199 return 0;
200
201 float_create(&exp);
202 float_create(&tmp);
203 float_setzero(&exp);
204
205 if (expx >= 0)
206 {
207 float_div(&exp, x, &cLn10, expx+1);
208 float_int(&exp);
209 extra = float_getexponent(&exp)+1;
210 float_mul(&tmp, &exp, &cLn10, digits+extra);
211 float_sub(x, x, &tmp, digits+extra);
212 if (float_cmp(x, &cLn10) >= 0)
213 {
214 float_add(&exp, &exp, &c1, EXACT);
215 float_sub(x, x, &cLn10, digits);
216 }
217 }
218 if (float_getsign(x) < 0)
219 {
220 float_sub(&exp, &exp, &c1, EXACT);
221 float_add(x, x, &cLn10, digits);
222 }
223 /* when we get here 0 <= x < ln 10 */
224 _expltln10(x, digits);
225 /* just in case rounding leads to a value >= 10 */
226 expx = float_getexponent(x);
227 if (expx != 0)
228 float_addi(&exp, &exp, expx, EXACT);
229
230 ok = 1;
231 if (!float_iszero(&exp))
232 {
233 expx = float_asinteger(&exp);
234 ok = expx != 0;
235 float_setexponent(x, expx);
236 }
237 float_free(&exp);
238 float_free(&tmp);
239 return ok && !float_isnan(x);
240 }
241
242 static char
_0_5exp(floatnum x,int digits)243 _0_5exp(
244 floatnum x,
245 int digits)
246 {
247 float_sub(x, x, &cLn2, digits + (3*logexp(x)/10)+1);
248 return _exp(x, digits);
249 }
250
251 /* exp(x)-1 for all x. Overflow is
252 indicated by the return value (0, if error)
253 relative error for 100 digit results is 8e-100 */
254 char
_expminus1(floatnum x,int digits)255 _expminus1(
256 floatnum x,
257 int digits)
258 {
259 int expr;
260
261 if (float_abscmp(x, &c1Div2) < 0)
262 {
263 _expminus1lt1(x, digits);
264 return 1;
265 }
266 if (float_getsign(x) < 0)
267 {
268 expr = (2*float_getexponent(x)/5);
269 if (expr >= digits)
270 float_setinteger(x, -1);
271 else
272 {
273 _exp(x, digits-expr);
274 float_sub(x, x, &c1, digits);
275 }
276 return 1;
277 }
278 if (!_exp(x, digits))
279 return 0;
280 float_sub(x, x, &c1, digits);
281 return 1;
282 }
283
284 static void
_addreciproc(floatnum x,int digits,signed char sgn)285 _addreciproc(
286 floatnum x,
287 int digits,
288 signed char sgn)
289 {
290 floatstruct tmp;
291 int expx;
292
293 expx = float_getexponent(x);
294 if (2*expx < digits)
295 {
296 float_create(&tmp);
297 float_muli(&tmp, x, 4, digits-2*expx);
298 float_reciprocal(&tmp, digits-2*expx);
299 float_setsign(&tmp, sgn);
300 float_add(x, x, &tmp, digits+1);
301 float_free(&tmp);
302 }
303 }
304
305 /* cosh(x)-1 for all x. Underflow or overflow is
306 indicated by the return value (0, if error)
307 relative error for 100 digit results is 6e-100 */
308
309 char
_coshminus1(floatnum x,int digits)310 _coshminus1(
311 floatnum x,
312 int digits)
313 {
314 if (float_getexponent(x) < 0 || float_iszero(x))
315 return _coshminus1lt1(x, digits);
316 if(!_0_5exp(x, digits))
317 return 0;
318 _addreciproc(x, digits, 1);
319 float_sub(x, x, &c1, digits);
320 return 1;
321 }
322
323 /* sinh(x) for all x. Overflow is
324 indicated by the return value (0, if error)
325 relative error for 100 digit results is < 8e-100 */
326 char
_sinh(floatnum x,int digits)327 _sinh(
328 floatnum x,
329 int digits)
330 {
331 if (float_getexponent(x) < 0 || float_iszero(x))
332 _sinhlt1(x, digits);
333 else
334 {
335 if(!_0_5exp(x, digits))
336 return 0;
337 _addreciproc(x, digits, -1);
338 }
339 return 1;
340 }
341
342 /* tanh(x) for |x| <= 0.5.
343 relative error for 100 digit results is < 7e-100 */
344 void
_tanhlt0_5(floatnum x,int digits)345 _tanhlt0_5(
346 floatnum x,
347 int digits)
348 {
349 floatstruct tmp;
350 signed char sgn;
351
352 float_create(&tmp);
353 sgn = float_getsign(x);
354 float_abs(x);
355 float_add(x, x, x, digits+1);
356 _expminus1lt1(x, digits);
357 float_add(&tmp, x, &c2, digits);
358 float_div(x, x, &tmp, digits);
359 float_setsign(x, sgn);
360 float_free(&tmp);
361 }
362
363 /* tanh(x)-1 for x > 0.
364 relative error for 100 digit results is < 9e-100 */
365 char
_tanhminus1gt0(floatnum x,int digits)366 _tanhminus1gt0(
367 floatnum x,
368 int digits)
369 {
370 if (float_add(x, x, x, digits+1) && _0_5exp(x, digits))
371 {
372 float_add(x, x, &c1Div2, digits+1);
373 float_reciprocal(x, digits+1);
374 float_setsign(x, -1);
375 return 1;
376 }
377 return 0;
378 }
379
380 void
_tanhgt0_5(floatnum x,int digits)381 _tanhgt0_5(
382 floatnum x,
383 int digits)
384 {
385 int expx;
386
387 expx = float_getexponent(x);
388 if (5*expx >= digits)
389 float_copy(x, &c1, EXACT);
390 else
391 {
392 _tanhminus1gt0(x, digits - 5*expx);
393 float_add(x, x, &c1, digits);
394 }
395 }
396
397 char
_power10(floatnum x,int digits)398 _power10(
399 floatnum x,
400 int digits)
401 {
402 int exp;
403
404 if (float_isinteger(x))
405 {
406 exp = float_asinteger(x);
407 if (exp == 0 && !float_iszero(x))
408 return 0;
409 float_copy(x, &c1, EXACT);
410 float_setexponent(x, exp);
411 return !float_isnan(x);
412 }
413 return float_mul(x, x, &cLn10, digits+2) && _exp(x, digits);
414 }
415