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

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since a3c6a85 was a3c6a85, checked in by Jiří Zárevúcky <jiri.zarevucky@…>, 7 years ago

Document sincos() and sincosf().
  • Property mode set to 100644
File size: 9.2 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
38#define TAYLOR_DEGREE_32 13
39#define TAYLOR_DEGREE_64 21
40
41/** Precomputed values for factorial (starting from 1!) */
42static double factorials[TAYLOR_DEGREE_64] = {
43 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
44 479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L,
45 20922789888000.0L, 355687428096000.0L, 6402373705728000.0L,
46 121645100408832000.0L, 2432902008176640000.0L, 51090942171709440000.0L
47};
48
49/** Sine approximation by Taylor series (32-bit floating point)
50 *
51 * Compute the approximation of sine by a Taylor
52 * series (using the first TAYLOR_DEGREE terms).
53 * The approximation is reasonably accurate for
54 * arguments within the interval [-pi/4, pi/4].
55 *
56 * @param arg Sine argument.
57 *
58 * @return Sine value approximation.
59 *
60 */
61static float taylor_sin_32(float arg)
62{
63 float ret = 0;
64 float nom = 1;
65
66 for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
67 nom *= arg;
68
69 if ((i % 4) == 0)
70 ret += nom / factorials[i];
71 else if ((i % 4) == 2)
72 ret -= nom / factorials[i];
73 }
74
75 return ret;
76}
77
78/** Sine approximation by Taylor series (64-bit floating point)
79 *
80 * Compute the approximation of sine by a Taylor
81 * series (using the first TAYLOR_DEGREE terms).
82 * The approximation is reasonably accurate for
83 * arguments within the interval [-pi/4, pi/4].
84 *
85 * @param arg Sine argument.
86 *
87 * @return Sine value approximation.
88 *
89 */
90static double taylor_sin_64(double arg)
91{
92 double ret = 0;
93 double nom = 1;
94
95 for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
96 nom *= arg;
97
98 if ((i % 4) == 0)
99 ret += nom / factorials[i];
100 else if ((i % 4) == 2)
101 ret -= nom / factorials[i];
102 }
103
104 return ret;
105}
106
107/** Cosine approximation by Taylor series (32-bit floating point)
108 *
109 * Compute the approximation of cosine by a Taylor
110 * series (using the first TAYLOR_DEGREE terms).
111 * The approximation is reasonably accurate for
112 * arguments within the interval [-pi/4, pi/4].
113 *
114 * @param arg Cosine argument.
115 *
116 * @return Cosine value approximation.
117 *
118 */
119static float taylor_cos_32(float arg)
120{
121 float ret = 1;
122 float nom = 1;
123
124 for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
125 nom *= arg;
126
127 if ((i % 4) == 1)
128 ret -= nom / factorials[i];
129 else if ((i % 4) == 3)
130 ret += nom / factorials[i];
131 }
132
133 return ret;
134}
135
136/** Cosine approximation by Taylor series (64-bit floating point)
137 *
138 * Compute the approximation of cosine by a Taylor
139 * series (using the first TAYLOR_DEGREE terms).
140 * The approximation is reasonably accurate for
141 * arguments within the interval [-pi/4, pi/4].
142 *
143 * @param arg Cosine argument.
144 *
145 * @return Cosine value approximation.
146 *
147 */
148static double taylor_cos_64(double arg)
149{
150 double ret = 1;
151 double nom = 1;
152
153 for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
154 nom *= arg;
155
156 if ((i % 4) == 1)
157 ret -= nom / factorials[i];
158 else if ((i % 4) == 3)
159 ret += nom / factorials[i];
160 }
161
162 return ret;
163}
164
165/** Sine value for values within base period (32-bit floating point)
166 *
167 * Compute the value of sine for arguments within
168 * the base period [0, 2pi]. For arguments outside
169 * the base period the returned values can be
170 * very inaccurate or even completely wrong.
171 *
172 * @param arg Sine argument.
173 *
174 * @return Sine value.
175 *
176 */
177static float base_sin_32(float arg)
178{
179 unsigned int period = arg / (M_PI / 4);
180
181 switch (period) {
182 case 0:
183 return taylor_sin_32(arg);
184 case 1:
185 case 2:
186 return taylor_cos_32(arg - M_PI / 2);
187 case 3:
188 case 4:
189 return -taylor_sin_32(arg - M_PI);
190 case 5:
191 case 6:
192 return -taylor_cos_32(arg - 3 * M_PI / 2);
193 default:
194 return taylor_sin_32(arg - 2 * M_PI);
195 }
196}
197
198/** Sine value for values within base period (64-bit floating point)
199 *
200 * Compute the value of sine for arguments within
201 * the base period [0, 2pi]. For arguments outside
202 * the base period the returned values can be
203 * very inaccurate or even completely wrong.
204 *
205 * @param arg Sine argument.
206 *
207 * @return Sine value.
208 *
209 */
210static double base_sin_64(double arg)
211{
212 unsigned int period = arg / (M_PI / 4);
213
214 switch (period) {
215 case 0:
216 return taylor_sin_64(arg);
217 case 1:
218 case 2:
219 return taylor_cos_64(arg - M_PI / 2);
220 case 3:
221 case 4:
222 return -taylor_sin_64(arg - M_PI);
223 case 5:
224 case 6:
225 return -taylor_cos_64(arg - 3 * M_PI / 2);
226 default:
227 return taylor_sin_64(arg - 2 * M_PI);
228 }
229}
230
231/** Cosine value for values within base period (32-bit floating point)
232 *
233 * Compute the value of cosine for arguments within
234 * the base period [0, 2pi]. For arguments outside
235 * the base period the returned values can be
236 * very inaccurate or even completely wrong.
237 *
238 * @param arg Cosine argument.
239 *
240 * @return Cosine value.
241 *
242 */
243static float base_cos_32(float arg)
244{
245 unsigned int period = arg / (M_PI / 4);
246
247 switch (period) {
248 case 0:
249 return taylor_cos_32(arg);
250 case 1:
251 case 2:
252 return -taylor_sin_32(arg - M_PI / 2);
253 case 3:
254 case 4:
255 return -taylor_cos_32(arg - M_PI);
256 case 5:
257 case 6:
258 return taylor_sin_32(arg - 3 * M_PI / 2);
259 default:
260 return taylor_cos_32(arg - 2 * M_PI);
261 }
262}
263
264/** Cosine value for values within base period (64-bit floating point)
265 *
266 * Compute the value of cosine for arguments within
267 * the base period [0, 2pi]. For arguments outside
268 * the base period the returned values can be
269 * very inaccurate or even completely wrong.
270 *
271 * @param arg Cosine argument.
272 *
273 * @return Cosine value.
274 *
275 */
276static double base_cos_64(double arg)
277{
278 unsigned int period = arg / (M_PI / 4);
279
280 switch (period) {
281 case 0:
282 return taylor_cos_64(arg);
283 case 1:
284 case 2:
285 return -taylor_sin_64(arg - M_PI / 2);
286 case 3:
287 case 4:
288 return -taylor_cos_64(arg - M_PI);
289 case 5:
290 case 6:
291 return taylor_sin_64(arg - 3 * M_PI / 2);
292 default:
293 return taylor_cos_64(arg - 2 * M_PI);
294 }
295}
296
297/** Sine (32-bit floating point)
298 *
299 * Compute sine value.
300 *
301 * @param arg Sine argument.
302 *
303 * @return Sine value.
304 *
305 */
306float sinf(float arg)
307{
308 float base_arg = fmodf(arg, 2 * M_PI);
309
310 if (base_arg < 0)
311 return -base_sin_32(-base_arg);
312
313 return base_sin_32(base_arg);
314}
315
316/** Sine (64-bit floating point)
317 *
318 * Compute sine value.
319 *
320 * @param arg Sine argument.
321 *
322 * @return Sine value.
323 *
324 */
325double sin(double arg)
326{
327 double base_arg = fmod(arg, 2 * M_PI);
328
329 if (base_arg < 0)
330 return -base_sin_64(-base_arg);
331
332 return base_sin_64(base_arg);
333}
334
335/** Cosine (32-bit floating point)
336 *
337 * Compute cosine value.
338 *
339 * @param arg Cosine argument.
340 *
341 * @return Cosine value.
342 *
343 */
344float cosf(float arg)
345{
346 float base_arg = fmodf(arg, 2 * M_PI);
347
348 if (base_arg < 0)
349 return base_cos_32(-base_arg);
350
351 return base_cos_32(base_arg);
352}
353
354/** Cosine (64-bit floating point)
355 *
356 * Compute cosine value.
357 *
358 * @param arg Cosine argument.
359 *
360 * @return Cosine value.
361 *
362 */
363double cos(double arg)
364{
365 double base_arg = fmod(arg, 2 * M_PI);
366
367 if (base_arg < 0)
368 return base_cos_64(-base_arg);
369
370 return base_cos_64(base_arg);
371}
372
373/**
374 * Computes sine and cosine at the same time, which might be more efficient than
375 * computing each separately.
376 *
377 * @param x Input value.
378 * @param s Output sine value, *s = sinf(x).
379 * @param c Output cosine value, *c = cosf(x).
380 */
381void sincosf(float x, float *s, float *c)
382{
383 float base_arg = fmodf(x, 2 * M_PI);
384
385 if (base_arg < 0) {
386 *s = -base_sin_32(-base_arg);
387 *c = base_cos_32(-base_arg);
388 } else {
389 *s = base_sin_32(base_arg);
390 *c = base_cos_32(base_arg);
391 }
392}
393
394/**
395 * Computes sine and cosine at the same time, which might be more efficient than
396 * computing each separately.
397 *
398 * @param x Input value.
399 * @param s Output sine value, *s = sin(x).
400 * @param c Output cosine value, *c = cos(x).
401 */
402void sincos(double x, double *s, double *c)
403{
404 double base_arg = fmod(x, 2 * M_PI);
405
406 if (base_arg < 0) {
407 *s = -base_sin_64(-base_arg);
408 *c = base_cos_64(-base_arg);
409 } else {
410 *s = base_sin_64(base_arg);
411 *c = base_cos_64(base_arg);
412 }
413}
414
415/** @}
416 */
Note: See TracBrowser for help on using the repository browser.