source: mainline/uspace/lib/math/generic/trig.c@ 3061bc1

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

Add single-precision variant for all functions. Allow generic implementations to call other functions while selecting the number of bits of precision, but not the implementation (generic or arch-specific).
  • Property mode set to 100644
File size: 8.4 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 <math.h>
37#include <trig.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/** Sine approximation by Taylor series (32-bit floating point)
51 *
52 * Compute the approximation of sine by a Taylor
53 * series (using the first TAYLOR_DEGREE terms).
54 * The approximation is reasonably accurate for
55 * arguments within the interval [-pi/4, pi/4].
56 *
57 * @param arg Sine argument.
58 *
59 * @return Sine value approximation.
60 *
61 */
62static float32_t taylor_sin_32(float32_t arg)
63{
64 float32_t ret = 0;
65 float32_t nom = 1;
66
67 for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
68 nom *= arg;
69
70 if ((i % 4) == 0)
71 ret += nom / factorials[i];
72 else if ((i % 4) == 2)
73 ret -= nom / factorials[i];
74 }
75
76 return ret;
77}
78
79/** Sine approximation by Taylor series (64-bit floating point)
80 *
81 * Compute the approximation of sine by a Taylor
82 * series (using the first TAYLOR_DEGREE terms).
83 * The approximation is reasonably accurate for
84 * arguments within the interval [-pi/4, pi/4].
85 *
86 * @param arg Sine argument.
87 *
88 * @return Sine value approximation.
89 *
90 */
91static float64_t taylor_sin_64(float64_t arg)
92{
93 float64_t ret = 0;
94 float64_t nom = 1;
95
96 for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
97 nom *= arg;
98
99 if ((i % 4) == 0)
100 ret += nom / factorials[i];
101 else if ((i % 4) == 2)
102 ret -= nom / factorials[i];
103 }
104
105 return ret;
106}
107
108/** Cosine approximation by Taylor series (32-bit floating point)
109 *
110 * Compute the approximation of cosine by a Taylor
111 * series (using the first TAYLOR_DEGREE terms).
112 * The approximation is reasonably accurate for
113 * arguments within the interval [-pi/4, pi/4].
114 *
115 * @param arg Cosine argument.
116 *
117 * @return Cosine value approximation.
118 *
119 */
120static float32_t taylor_cos_32(float32_t arg)
121{
122 float32_t ret = 1;
123 float32_t nom = 1;
124
125 for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
126 nom *= arg;
127
128 if ((i % 4) == 1)
129 ret -= nom / factorials[i];
130 else if ((i % 4) == 3)
131 ret += nom / factorials[i];
132 }
133
134 return ret;
135}
136
137/** Cosine approximation by Taylor series (64-bit floating point)
138 *
139 * Compute the approximation of cosine by a Taylor
140 * series (using the first TAYLOR_DEGREE terms).
141 * The approximation is reasonably accurate for
142 * arguments within the interval [-pi/4, pi/4].
143 *
144 * @param arg Cosine argument.
145 *
146 * @return Cosine value approximation.
147 *
148 */
149static float64_t taylor_cos_64(float64_t arg)
150{
151 float64_t ret = 1;
152 float64_t nom = 1;
153
154 for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
155 nom *= arg;
156
157 if ((i % 4) == 1)
158 ret -= nom / factorials[i];
159 else if ((i % 4) == 3)
160 ret += nom / factorials[i];
161 }
162
163 return ret;
164}
165
166/** Sine value for values within base period (32-bit floating point)
167 *
168 * Compute the value of sine for arguments within
169 * the base period [0, 2pi]. For arguments outside
170 * the base period the returned values can be
171 * very inaccurate or even completely wrong.
172 *
173 * @param arg Sine argument.
174 *
175 * @return Sine value.
176 *
177 */
178static float32_t base_sin_32(float32_t arg)
179{
180 unsigned int period = arg / (M_PI / 4);
181
182 switch (period) {
183 case 0:
184 return taylor_sin_32(arg);
185 case 1:
186 case 2:
187 return taylor_cos_32(arg - M_PI / 2);
188 case 3:
189 case 4:
190 return -taylor_sin_32(arg - M_PI);
191 case 5:
192 case 6:
193 return -taylor_cos_32(arg - 3 * M_PI / 2);
194 default:
195 return taylor_sin_32(arg - 2 * M_PI);
196 }
197}
198
199/** Sine value for values within base period (64-bit floating point)
200 *
201 * Compute the value of sine for arguments within
202 * the base period [0, 2pi]. For arguments outside
203 * the base period the returned values can be
204 * very inaccurate or even completely wrong.
205 *
206 * @param arg Sine argument.
207 *
208 * @return Sine value.
209 *
210 */
211static float64_t base_sin_64(float64_t arg)
212{
213 unsigned int period = arg / (M_PI / 4);
214
215 switch (period) {
216 case 0:
217 return taylor_sin_64(arg);
218 case 1:
219 case 2:
220 return taylor_cos_64(arg - M_PI / 2);
221 case 3:
222 case 4:
223 return -taylor_sin_64(arg - M_PI);
224 case 5:
225 case 6:
226 return -taylor_cos_64(arg - 3 * M_PI / 2);
227 default:
228 return taylor_sin_64(arg - 2 * M_PI);
229 }
230}
231
232/** Cosine value for values within base period (32-bit floating point)
233 *
234 * Compute the value of cosine for arguments within
235 * the base period [0, 2pi]. For arguments outside
236 * the base period the returned values can be
237 * very inaccurate or even completely wrong.
238 *
239 * @param arg Cosine argument.
240 *
241 * @return Cosine value.
242 *
243 */
244static float32_t base_cos_32(float32_t arg)
245{
246 unsigned int period = arg / (M_PI / 4);
247
248 switch (period) {
249 case 0:
250 return taylor_cos_32(arg);
251 case 1:
252 case 2:
253 return -taylor_sin_32(arg - M_PI / 2);
254 case 3:
255 case 4:
256 return -taylor_cos_32(arg - M_PI);
257 case 5:
258 case 6:
259 return taylor_sin_32(arg - 3 * M_PI / 2);
260 default:
261 return taylor_cos_32(arg - 2 * M_PI);
262 }
263}
264
265/** Cosine value for values within base period (64-bit floating point)
266 *
267 * Compute the value of cosine for arguments within
268 * the base period [0, 2pi]. For arguments outside
269 * the base period the returned values can be
270 * very inaccurate or even completely wrong.
271 *
272 * @param arg Cosine argument.
273 *
274 * @return Cosine value.
275 *
276 */
277static float64_t base_cos_64(float64_t arg)
278{
279 unsigned int period = arg / (M_PI / 4);
280
281 switch (period) {
282 case 0:
283 return taylor_cos_64(arg);
284 case 1:
285 case 2:
286 return -taylor_sin_64(arg - M_PI / 2);
287 case 3:
288 case 4:
289 return -taylor_cos_64(arg - M_PI);
290 case 5:
291 case 6:
292 return taylor_sin_64(arg - 3 * M_PI / 2);
293 default:
294 return taylor_cos_64(arg - 2 * M_PI);
295 }
296}
297
298/** Sine (32-bit floating point)
299 *
300 * Compute sine value.
301 *
302 * @param arg Sine argument.
303 *
304 * @return Sine value.
305 *
306 */
307float32_t float32_sin(float32_t arg)
308{
309 float32_t base_arg = fmod_f32(arg, 2 * M_PI);
310
311 if (base_arg < 0)
312 return -base_sin_32(-base_arg);
313
314 return base_sin_32(base_arg);
315}
316
317/** Sine (64-bit floating point)
318 *
319 * Compute sine value.
320 *
321 * @param arg Sine argument.
322 *
323 * @return Sine value.
324 *
325 */
326float64_t float64_sin(float64_t arg)
327{
328 float64_t base_arg = fmod_f64(arg, 2 * M_PI);
329
330 if (base_arg < 0)
331 return -base_sin_64(-base_arg);
332
333 return base_sin_64(base_arg);
334}
335
336/** Cosine (32-bit floating point)
337 *
338 * Compute cosine value.
339 *
340 * @param arg Cosine argument.
341 *
342 * @return Cosine value.
343 *
344 */
345float32_t float32_cos(float32_t arg)
346{
347 float32_t base_arg = fmod_f32(arg, 2 * M_PI);
348
349 if (base_arg < 0)
350 return base_cos_32(-base_arg);
351
352 return base_cos_32(base_arg);
353}
354
355/** Cosine (64-bit floating point)
356 *
357 * Compute cosine value.
358 *
359 * @param arg Cosine argument.
360 *
361 * @return Cosine value.
362 *
363 */
364float64_t float64_cos(float64_t arg)
365{
366 float64_t base_arg = fmod_f64(arg, 2 * M_PI);
367
368 if (base_arg < 0)
369 return base_cos_64(-base_arg);
370
371 return base_cos_64(base_arg);
372}
373
374/** @}
375 */
Note: See TracBrowser for help on using the repository browser.