]>
Commit | Line | Data |
---|---|---|
9acda61d | 1 | /* Auxiliary routine for the Bessel functions (j0f, y0f, j1f, y1f). |
6d7e8eda | 2 | Copyright (C) 2021-2023 Free Software Foundation, Inc. |
9acda61d PZ |
3 | This file is part of the GNU C Library. |
4 | ||
5 | The GNU C Library is free software; you can redistribute it and/or | |
6 | modify it under the terms of the GNU Lesser General Public | |
7 | License as published by the Free Software Foundation; either | |
8 | version 2.1 of the License, or (at your option) any later version. | |
9 | ||
10 | The GNU C Library is distributed in the hope that it will be useful, | |
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 | Lesser General Public License for more details. | |
14 | ||
15 | You should have received a copy of the GNU Lesser General Public | |
16 | License along with the GNU C Library; if not, see | |
17 | <https://www.gnu.org/licenses/>. */ | |
18 | ||
19 | #ifndef _MATH_REDUCE_AUX_H | |
20 | #define _MATH_REDUCE_AUX_H | |
21 | ||
22 | #include <math.h> | |
23 | #include <math_private.h> | |
24 | #include <s_sincosf.h> | |
25 | ||
26 | /* Return h and update n such that: | |
27 | Now x - pi/4 - alpha = h + n*pi/2 mod (2*pi). */ | |
28 | static inline double | |
29 | reduce_aux (float x, int *n, double alpha) | |
30 | { | |
31 | double h; | |
32 | h = reduce_large (asuint (x), n); | |
33 | /* Now |x| = h+n*pi/2 mod 2*pi. */ | |
34 | /* Recover sign. */ | |
35 | if (x < 0) | |
36 | { | |
37 | h = -h; | |
38 | *n = -*n; | |
39 | } | |
40 | /* Subtract pi/4. */ | |
41 | double piover2 = 0xc.90fdaa22168cp-3; | |
42 | if (h >= 0) | |
43 | h -= piover2 / 2; | |
44 | else | |
45 | { | |
46 | h += piover2 / 2; | |
47 | (*n) --; | |
48 | } | |
49 | /* Subtract alpha and reduce if needed mod pi/2. */ | |
50 | h -= alpha; | |
51 | if (h > piover2) | |
52 | { | |
53 | h -= piover2; | |
54 | (*n) ++; | |
55 | } | |
56 | else if (h < -piover2) | |
57 | { | |
58 | h += piover2; | |
59 | (*n) --; | |
60 | } | |
61 | return h; | |
62 | } | |
63 | ||
64 | #endif |