source: mainline/uspace/lib/math/generic/exp.c@ c280d7e

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since c280d7e was b69786e, checked in by Jiri Svoboda <jiri@…>, 10 years ago

Add all functions required by C89 plus their single-precision variants from C99.

  • Property mode set to 100644
File size: 4.1 KB
Line 
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!) */
43static 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 */
62static 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 */
87static 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 */
109float32_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 */
139float64_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 */
Note: See TracBrowser for help on using the repository browser.