1 | /*
|
---|
2 | * Copyright (c) 2015 Jiri Svoboda
|
---|
3 | * Copyright (c) 2014 Martin Decky
|
---|
4 | * All rights reserved.
|
---|
5 | *
|
---|
6 | * Redistribution and use in source and binary forms, with or without
|
---|
7 | * modification, are permitted provided that the following conditions
|
---|
8 | * are met:
|
---|
9 | *
|
---|
10 | * - Redistributions of source code must retain the above copyright
|
---|
11 | * notice, this list of conditions and the following disclaimer.
|
---|
12 | * - Redistributions in binary form must reproduce the above copyright
|
---|
13 | * notice, this list of conditions and the following disclaimer in the
|
---|
14 | * documentation and/or other materials provided with the distribution.
|
---|
15 | * - The name of the author may not be used to endorse or promote products
|
---|
16 | * derived from this software without specific prior written permission.
|
---|
17 | *
|
---|
18 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
|
---|
19 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
|
---|
20 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
|
---|
21 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
|
---|
22 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
|
---|
23 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
---|
24 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
---|
25 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
---|
26 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
|
---|
27 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
---|
28 | */
|
---|
29 |
|
---|
30 | /** @addtogroup libmath
|
---|
31 | * @{
|
---|
32 | */
|
---|
33 | /** @file
|
---|
34 | */
|
---|
35 |
|
---|
36 | #include <exp.h>
|
---|
37 | #include <math.h>
|
---|
38 |
|
---|
39 | #define TAYLOR_DEGREE_32 13
|
---|
40 | #define TAYLOR_DEGREE_64 21
|
---|
41 |
|
---|
42 | /** Precomputed values for factorial (starting from 1!) */
|
---|
43 | static float64_t factorials[TAYLOR_DEGREE_64] = {
|
---|
44 | 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
|
---|
45 | 479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L,
|
---|
46 | 20922789888000.0L, 355687428096000.0L, 6402373705728000.0L,
|
---|
47 | 121645100408832000.0L, 2432902008176640000.0L, 51090942171709440000.0L
|
---|
48 | };
|
---|
49 |
|
---|
50 | /** Exponential approximation by Taylor series (32-bit floating point)
|
---|
51 | *
|
---|
52 | * Compute the approximation of exponential by a Taylor
|
---|
53 | * series (using the first TAYLOR_DEGREE terms).
|
---|
54 | * The approximation is reasonably accurate for
|
---|
55 | * arguments within the interval XXXX.
|
---|
56 | *
|
---|
57 | * @param arg Argument.
|
---|
58 | *
|
---|
59 | * @return Exponential value approximation.
|
---|
60 | *
|
---|
61 | */
|
---|
62 | static float32_t taylor_exp_32(float32_t arg)
|
---|
63 | {
|
---|
64 | float32_t ret = 1;
|
---|
65 | float32_t nom = 1;
|
---|
66 |
|
---|
67 | for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
|
---|
68 | nom *= arg;
|
---|
69 | ret += nom / factorials[i];
|
---|
70 | }
|
---|
71 |
|
---|
72 | return ret;
|
---|
73 | }
|
---|
74 |
|
---|
75 | /** Exponential approximation by Taylor series (64-bit floating point)
|
---|
76 | *
|
---|
77 | * Compute the approximation of exponential by a Taylor
|
---|
78 | * series (using the first TAYLOR_DEGREE terms).
|
---|
79 | * The approximation is reasonably accurate for
|
---|
80 | * arguments within the interval XXXX.
|
---|
81 | *
|
---|
82 | * @param arg Argument.
|
---|
83 | *
|
---|
84 | * @return Exponential value approximation.
|
---|
85 | *
|
---|
86 | */
|
---|
87 | static float64_t taylor_exp_64(float64_t arg)
|
---|
88 | {
|
---|
89 | float64_t ret = 1;
|
---|
90 | float64_t nom = 1;
|
---|
91 |
|
---|
92 | for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
|
---|
93 | nom *= arg;
|
---|
94 | ret += nom / factorials[i];
|
---|
95 | }
|
---|
96 |
|
---|
97 | return ret;
|
---|
98 | }
|
---|
99 |
|
---|
100 | /** Exponential (32-bit floating point)
|
---|
101 | *
|
---|
102 | * Compute exponential value.
|
---|
103 | *
|
---|
104 | * @param arg Exponential argument.
|
---|
105 | *
|
---|
106 | * @return Exponential value.
|
---|
107 | *
|
---|
108 | */
|
---|
109 | float32_t float32_exp(float32_t arg)
|
---|
110 | {
|
---|
111 | float32_t f;
|
---|
112 | float32_t i;
|
---|
113 | float32_u r;
|
---|
114 |
|
---|
115 | /*
|
---|
116 | * e^a = (2 ^ log2(e))^a = 2 ^ (log2(e) * a)
|
---|
117 | * log2(e) * a = i + f | f in [0, 1]
|
---|
118 | * e ^ a = 2 ^ (i + f) = 2^f * 2^i = (e ^ log(2))^f * 2^i =
|
---|
119 | * e^(log(2)*f) * 2^i
|
---|
120 | */
|
---|
121 |
|
---|
122 | i = trunc_f32(arg * M_LOG2E);
|
---|
123 | f = arg * M_LOG2E - i;
|
---|
124 |
|
---|
125 | r.val = taylor_exp_32(M_LN2 * f);
|
---|
126 | r.data.parts.exp += i;
|
---|
127 | return r.val;
|
---|
128 | }
|
---|
129 |
|
---|
130 | /** Exponential (64-bit floating point)
|
---|
131 | *
|
---|
132 | * Compute exponential value.
|
---|
133 | *
|
---|
134 | * @param arg Exponential argument.
|
---|
135 | *
|
---|
136 | * @return Exponential value.
|
---|
137 | *
|
---|
138 | */
|
---|
139 | float64_t float64_exp(float64_t arg)
|
---|
140 | {
|
---|
141 | float64_t f;
|
---|
142 | float64_t i;
|
---|
143 | float64_u r;
|
---|
144 |
|
---|
145 | /*
|
---|
146 | * e^a = (2 ^ log2(e))^a = 2 ^ (log2(e) * a)
|
---|
147 | * log2(e) * a = i + f | f in [0, 1]
|
---|
148 | * e ^ a = 2 ^ (i + f) = 2^f * 2^i = (e ^ log(2))^f * 2^i =
|
---|
149 | * e^(log(2)*f) * 2^i
|
---|
150 | */
|
---|
151 |
|
---|
152 | i = trunc_f64(arg * M_LOG2E);
|
---|
153 | f = arg * M_LOG2E - i;
|
---|
154 |
|
---|
155 | r.val = taylor_exp_64(M_LN2 * f);
|
---|
156 | r.data.parts.exp += i;
|
---|
157 | return r.val;
|
---|
158 | }
|
---|
159 |
|
---|
160 | /** @}
|
---|
161 | */
|
---|