1 /* exp10.c
2 *
3 * Base 10 exponential function
4 * (Common antilogarithm)
5 *
6 *
7 *
8 * SYNOPSIS:
9 *
10 * double x, y, exp10();
11 *
12 * y = exp10( x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns 10 raised to the x power.
19 *
20 * Range reduction is accomplished by expressing the argument
21 * as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
22 * The Pade' form
23 *
24 * 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
25 *
26 * is used to approximate 10**f.
27 *
28 *
29 *
30 * ACCURACY:
31 *
32 * Relative error:
33 * arithmetic domain # trials peak rms
34 * IEEE -307,+307 30000 2.2e-16 5.5e-17
35 *
36 * ERROR MESSAGES:
37 *
38 * message condition value returned
39 * exp10 underflow x < -MAXL10 0.0
40 * exp10 overflow x > MAXL10 NPY_INFINITY
41 *
42 * IEEE arithmetic: MAXL10 = 308.2547155599167.
43 *
44 */
45
46 /*
47 * Cephes Math Library Release 2.2: January, 1991
48 * Copyright 1984, 1991 by Stephen L. Moshier
49 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
50 */
51
52
53 #include "mconf.h"
54
55 static double P[] = {
56 4.09962519798587023075E-2,
57 1.17452732554344059015E1,
58 4.06717289936872725516E2,
59 2.39423741207388267439E3,
60 };
61
62 static double Q[] = {
63 /* 1.00000000000000000000E0, */
64 8.50936160849306532625E1,
65 1.27209271178345121210E3,
66 2.07960819286001865907E3,
67 };
68
69 /* static double LOG102 = 3.01029995663981195214e-1; */
70 static double LOG210 = 3.32192809488736234787e0;
71 static double LG102A = 3.01025390625000000000E-1;
72 static double LG102B = 4.60503898119521373889E-6;
73
74 /* static double MAXL10 = 38.230809449325611792; */
75 static double MAXL10 = 308.2547155599167;
76
exp10(double x)77 double exp10(double x)
78 {
79 double px, xx;
80 short n;
81
82 if (cephes_isnan(x))
83 return (x);
84 if (x > MAXL10) {
85 return (NPY_INFINITY);
86 }
87
88 if (x < -MAXL10) { /* Would like to use MINLOG but can't */
89 sf_error("exp10", SF_ERROR_UNDERFLOW, NULL);
90 return (0.0);
91 }
92
93 /* Express 10**x = 10**g 2**n
94 * = 10**g 10**( n log10(2) )
95 * = 10**( g + n log10(2) )
96 */
97 px = floor(LOG210 * x + 0.5);
98 n = px;
99 x -= px * LG102A;
100 x -= px * LG102B;
101
102 /* rational approximation for exponential
103 * of the fractional part:
104 * 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
105 */
106 xx = x * x;
107 px = x * polevl(xx, P, 3);
108 x = px / (p1evl(xx, Q, 3) - px);
109 x = 1.0 + ldexp(x, 1);
110
111 /* multiply by power of 2 */
112 x = ldexp(x, n);
113
114 return (x);
115 }
116