1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2016 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"
36void __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 */
44double
45SECTION
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