1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
5 | * |
6 | * This program is free software; you can redistribute it and/or modify |
7 | * it under the terms of the GNU Lesser General Public License as published by |
8 | * the Free Software Foundation; either version 2.1 of the License, or |
9 | * (at your option) any later version. |
10 | * |
11 | * This program is distributed in the hope that it will be useful, |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 | * GNU Lesser General Public License for more details. |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License |
17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
18 | */ |
19 | /**************************************************************************/ |
20 | /* MODULE_NAME:slowexp.c */ |
21 | /* */ |
22 | /* FUNCTION:slowexp */ |
23 | /* */ |
24 | /* FILES NEEDED:mpa.h */ |
25 | /* mpa.c mpexp.c */ |
26 | /* */ |
27 | /*Converting from double precision to Multi-precision and calculating */ |
28 | /* e^x */ |
29 | /**************************************************************************/ |
30 | #include <math_private.h> |
31 | |
32 | #include <stap-probe.h> |
33 | |
34 | #ifndef USE_LONG_DOUBLE_FOR_MP |
35 | # include "mpa.h" |
36 | void __mpexp (mp_no *x, mp_no *y, int p); |
37 | #endif |
38 | |
39 | #ifndef SECTION |
40 | # define SECTION |
41 | #endif |
42 | |
43 | /*Converting from double precision to Multi-precision and calculating e^x */ |
44 | double |
45 | SECTION |
46 | __slowexp (double x) |
47 | { |
48 | #ifndef USE_LONG_DOUBLE_FOR_MP |
49 | double w, z, res, eps = 3.0e-26; |
50 | int p; |
51 | mp_no mpx, mpy, mpz, mpw, mpeps, mpcor; |
52 | |
53 | /* Use the multiple precision __MPEXP function to compute the exponential |
54 | First at 144 bits and if it is not accurate enough, at 768 bits. */ |
55 | p = 6; |
56 | __dbl_mp (x, &mpx, p); |
57 | __mpexp (&mpx, &mpy, p); |
58 | __dbl_mp (eps, &mpeps, p); |
59 | __mul (&mpeps, &mpy, &mpcor, p); |
60 | __add (&mpy, &mpcor, &mpw, p); |
61 | __sub (&mpy, &mpcor, &mpz, p); |
62 | __mp_dbl (&mpw, &w, p); |
63 | __mp_dbl (&mpz, &z, p); |
64 | if (w == z) |
65 | { |
66 | /* Track how often we get to the slow exp code plus |
67 | its input/output values. */ |
68 | LIBC_PROBE (slowexp_p6, 2, &x, &w); |
69 | return w; |
70 | } |
71 | else |
72 | { |
73 | p = 32; |
74 | __dbl_mp (x, &mpx, p); |
75 | __mpexp (&mpx, &mpy, p); |
76 | __mp_dbl (&mpy, &res, p); |
77 | |
78 | /* Track how often we get to the uber-slow exp code plus |
79 | its input/output values. */ |
80 | LIBC_PROBE (slowexp_p32, 2, &x, &res); |
81 | return res; |
82 | } |
83 | #else |
84 | return (double) __ieee754_expl((long double)x); |
85 | #endif |
86 | } |
87 | |