|
19 | 19 | * |
20 | 20 | * The original C code, long comment, copyright, license, and constants are from [Cephes]{@link http://www.netlib.org/cephes}. |
21 | 21 | * The implementation follows the original, but has been modified according to project conventions. |
22 | | -* |
23 | | -* Copyright 1984, 1991, 2000 by Stephen L. Moshier |
24 | | -* |
25 | | -* Some software in this archive may be from the book _Methods and Programs for Mathematical Functions_ |
26 | | -* (Prentice-Hall or Simon & Schuster International, 1989) or from the Cephes Mathematical Library, |
27 | | -* a commercial product. In either event, it is copyrighted by the author. |
28 | | -* What you see here may be used freely but it comes with no support or guarantee. |
29 | | -* |
30 | | -* Stephen L. Moshier |
31 | | - |
32 | 22 | */ |
33 | 23 |
|
34 | 24 | #include "stdlib/math/base/special/exp10f.h" |
|
40 | 30 | #include "stdlib/math/base/special/ldexpf.h" |
41 | 31 | #include <stdint.h> |
42 | 32 |
|
43 | | -// Constants: |
44 | | -static const float LOG210F = 3.3219281f; |
45 | | -static const float LG102AF = 0.30102539f; |
46 | | -static const float LG102BF = 4.6050378e-06f; |
| 33 | +// Cephes constants: |
47 | 34 |
|
48 | | -// Evaluate polynomial P(x): |
49 | | -static float polyval_p( const float x ) { |
50 | | - return 2394.2374f + (x * (406.7173f + (x * (11.745273f + (x * 0.04099625f))))); |
| 35 | +static float LOG210 = 3.32192809488736234787e0; |
| 36 | +static float LG102A = 3.00781250000000000000E-1; |
| 37 | +static float LG102B = 2.48745663981195213739E-4; |
| 38 | +static float MAXL10 = 38.230809449325611792; |
| 39 | +/** |
| 40 | +* Evaluates a polynomial using Horner's rule. |
| 41 | +* P(x)=a5+x(a4+x(a3+x(a2+x(a1+xa0)))) |
| 42 | +* |
| 43 | +* @param x input value |
| 44 | +* @return result of polynomial evaluation |
| 45 | +*/ |
| 46 | +static float polyval( const float x ) { |
| 47 | + return 2.302585167056758e+0f + x * (2.650948748208892e+0f + x * (2.034649854009453e+0f + x * (1.171292686296281e+0f + x * (5.420251702225484e-1f + x * 2.063216740311022e-1f)))); |
51 | 48 | } |
52 | 49 |
|
53 | | -// Evaluate polynomial Q(x): |
54 | | -static float polyval_q( const float x ) { |
55 | | - return 2079.6082f + (x * (1272.0927f + (x * (85.09362f + (x * 1.0f))))); |
56 | | -} |
57 | 50 |
|
58 | 51 | /** |
59 | | -* Returns 10 raised to the x power (single-precision). |
| 52 | +* Computes 10^x for a single-precision floating-point number. |
60 | 53 | * |
61 | 54 | * ## Method |
62 | 55 | * |
63 | | -* - Range reduction is done using 10^x = 2^n * 10^f, with |f| < 0.5 * log10(2) |
64 | | -* - The reduced part is evaluated with a rational approximation: |
65 | | -* |
66 | | -* 1 + 2x P(x^2) / (Q(x^2) - P(x^2)) |
| 56 | +* - Performs range reduction: 10^x = 2^n * 10^f, where |f| < 0.5*log10(2) |
| 57 | +* - Approximates 10^f using a polynomial expansion from Cephes |
67 | 58 | * |
68 | 59 | * @param x input value |
69 | | -* @return 10^x |
| 60 | +* @return single-precision result of 10^x |
70 | 61 | */ |
71 | 62 | float stdlib_base_exp10f( const float x ) { |
72 | 63 | float xc; |
73 | 64 | float px; |
74 | | - float xx; |
75 | 65 | int32_t n; |
76 | 66 |
|
77 | 67 | if ( stdlib_base_is_nanf( x ) ) { |
78 | 68 | return 0.0f / 0.0f; |
79 | 69 | } |
80 | | - if ( x > STDLIB_CONSTANT_FLOAT32_MAX_BASE10_EXPONENT ) { |
| 70 | + if ( x > MAXL10 ) { |
81 | 71 | return STDLIB_CONSTANT_FLOAT32_PINF; |
82 | 72 | } |
| 73 | + if ( x < -MAXL10 ) { |
| 74 | + return 0.0f; |
| 75 | + } |
| 76 | + if ( x == 0.0f ) { |
| 77 | + return 1.0f; |
| 78 | + } |
83 | 79 |
|
84 | | - // Compute n = round(x * log2(10)): |
85 | | - px = stdlib_base_floorf( (LOG210F * x) + 0.5f ); |
| 80 | + // Range reduction: |
| 81 | + px = stdlib_base_floorf( (LOG210 * x) + 0.5f ); |
86 | 82 | n = (int32_t)px; |
87 | | - |
88 | | - // Compute f = x - n * log10(2): |
89 | 83 | xc = x; |
90 | | - xc -= (px * LG102AF); |
91 | | - xc -= (px * LG102BF); |
| 84 | + xc -= px * LG102A; |
| 85 | + xc -= px * LG102B; |
92 | 86 |
|
93 | | - // Evaluate rational approximation for 10^f: |
94 | | - xx = xc * xc; |
95 | | - px = xc * polyval_p( xx ); |
96 | | - xc = px / (polyval_q( xx ) - px); |
97 | | - xc = 1.0f + stdlib_base_ldexpf( xc, 1 ); |
| 87 | + // Polynomial approximation: 10^f ≈ 1 + f * P(f) |
| 88 | + px = xc * polyval( xc ); |
| 89 | + xc = 1.0f + px; |
98 | 90 |
|
99 | | - // Return result: |
| 91 | + // Multiply by power of 2: 10^x = 10^f * 2^n |
100 | 92 | return stdlib_base_ldexpf( xc, n ); |
101 | 93 | } |
0 commit comments