1 /*
2 * Copyright (c) 2003, 2007-14 Matteo Frigo
3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program 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
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 *
19 */
20
21
22 /* trigonometric functions */
23 #include "kernel/ifftw.h"
24 #include <math.h>
25
26 #if defined(TRIGREAL_IS_LONG_DOUBLE)
27 # define COS cosl
28 # define SIN sinl
29 # define KTRIG(x) (x##L)
30 # if defined(HAVE_DECL_SINL) && !HAVE_DECL_SINL
31 extern long double sinl(long double x);
32 # endif
33 # if defined(HAVE_DECL_COSL) && !HAVE_DECL_COSL
34 extern long double cosl(long double x);
35 # endif
36 #elif defined(TRIGREAL_IS_QUAD)
37 # define COS cosq
38 # define SIN sinq
39 # define KTRIG(x) (x##Q)
40 extern __float128 sinq(__float128 x);
41 extern __float128 cosq(__float128 x);
42 #else
43 # define COS cos
44 # define SIN sin
45 # define KTRIG(x) (x)
46 #endif
47
48 static const trigreal K2PI =
49 KTRIG(6.2831853071795864769252867665590057683943388);
50 #define by2pi(m, n) ((K2PI * (m)) / (n))
51
52 /*
53 * Improve accuracy by reducing x to range [0..1/8]
54 * before multiplication by 2 * PI.
55 */
56
real_cexp(INT m,INT n,trigreal * out)57 static void real_cexp(INT m, INT n, trigreal *out)
58 {
59 trigreal theta, c, s, t;
60 unsigned octant = 0;
61 INT quarter_n = n;
62
63 n += n; n += n;
64 m += m; m += m;
65
66 if (m < 0) m += n;
67 if (m > n - m) { m = n - m; octant |= 4; }
68 if (m - quarter_n > 0) { m = m - quarter_n; octant |= 2; }
69 if (m > quarter_n - m) { m = quarter_n - m; octant |= 1; }
70
71 theta = by2pi(m, n);
72 c = COS(theta); s = SIN(theta);
73
74 if (octant & 1) { t = c; c = s; s = t; }
75 if (octant & 2) { t = c; c = -s; s = t; }
76 if (octant & 4) { s = -s; }
77
78 out[0] = c;
79 out[1] = s;
80 }
81
choose_twshft(INT n)82 static INT choose_twshft(INT n)
83 {
84 INT log2r = 0;
85 while (n > 0) {
86 ++log2r;
87 n /= 4;
88 }
89 return log2r;
90 }
91
cexpl_sqrtn_table(triggen * p,INT m,trigreal * res)92 static void cexpl_sqrtn_table(triggen *p, INT m, trigreal *res)
93 {
94 m += p->n * (m < 0);
95
96 {
97 INT m0 = m & p->twmsk;
98 INT m1 = m >> p->twshft;
99 trigreal wr0 = p->W0[2 * m0];
100 trigreal wi0 = p->W0[2 * m0 + 1];
101 trigreal wr1 = p->W1[2 * m1];
102 trigreal wi1 = p->W1[2 * m1 + 1];
103
104 res[0] = wr1 * wr0 - wi1 * wi0;
105 res[1] = wi1 * wr0 + wr1 * wi0;
106 }
107 }
108
109 /* multiply (xr, xi) by exp(FFT_SIGN * 2*pi*i*m/n) */
rotate_sqrtn_table(triggen * p,INT m,R xr,R xi,R * res)110 static void rotate_sqrtn_table(triggen *p, INT m, R xr, R xi, R *res)
111 {
112 m += p->n * (m < 0);
113
114 {
115 INT m0 = m & p->twmsk;
116 INT m1 = m >> p->twshft;
117 trigreal wr0 = p->W0[2 * m0];
118 trigreal wi0 = p->W0[2 * m0 + 1];
119 trigreal wr1 = p->W1[2 * m1];
120 trigreal wi1 = p->W1[2 * m1 + 1];
121 trigreal wr = wr1 * wr0 - wi1 * wi0;
122 trigreal wi = wi1 * wr0 + wr1 * wi0;
123
124 #if FFT_SIGN == -1
125 res[0] = xr * wr + xi * wi;
126 res[1] = xi * wr - xr * wi;
127 #else
128 res[0] = xr * wr - xi * wi;
129 res[1] = xi * wr + xr * wi;
130 #endif
131 }
132 }
133
cexpl_sincos(triggen * p,INT m,trigreal * res)134 static void cexpl_sincos(triggen *p, INT m, trigreal *res)
135 {
136 real_cexp(m, p->n, res);
137 }
138
cexp_zero(triggen * p,INT m,R * res)139 static void cexp_zero(triggen *p, INT m, R *res)
140 {
141 UNUSED(p); UNUSED(m);
142 res[0] = 0;
143 res[1] = 0;
144 }
145
cexpl_zero(triggen * p,INT m,trigreal * res)146 static void cexpl_zero(triggen *p, INT m, trigreal *res)
147 {
148 UNUSED(p); UNUSED(m);
149 res[0] = 0;
150 res[1] = 0;
151 }
152
cexp_generic(triggen * p,INT m,R * res)153 static void cexp_generic(triggen *p, INT m, R *res)
154 {
155 trigreal resl[2];
156 p->cexpl(p, m, resl);
157 res[0] = (R)resl[0];
158 res[1] = (R)resl[1];
159 }
160
rotate_generic(triggen * p,INT m,R xr,R xi,R * res)161 static void rotate_generic(triggen *p, INT m, R xr, R xi, R *res)
162 {
163 trigreal w[2];
164 p->cexpl(p, m, w);
165 res[0] = xr * w[0] - xi * (FFT_SIGN * w[1]);
166 res[1] = xi * w[0] + xr * (FFT_SIGN * w[1]);
167 }
168
X(mktriggen)169 triggen *X(mktriggen)(enum wakefulness wakefulness, INT n)
170 {
171 INT i, n0, n1;
172 triggen *p = (triggen *)MALLOC(sizeof(*p), TWIDDLES);
173
174 p->n = n;
175 p->W0 = p->W1 = 0;
176 p->cexp = 0;
177 p->rotate = 0;
178
179 switch (wakefulness) {
180 case SLEEPY:
181 A(0 /* can't happen */);
182 break;
183
184 case AWAKE_SQRTN_TABLE: {
185 INT twshft = choose_twshft(n);
186
187 p->twshft = twshft;
188 p->twradix = ((INT)1) << twshft;
189 p->twmsk = p->twradix - 1;
190
191 n0 = p->twradix;
192 n1 = (n + n0 - 1) / n0;
193
194 p->W0 = (trigreal *)MALLOC(n0 * 2 * sizeof(trigreal), TWIDDLES);
195 p->W1 = (trigreal *)MALLOC(n1 * 2 * sizeof(trigreal), TWIDDLES);
196
197 for (i = 0; i < n0; ++i)
198 real_cexp(i, n, p->W0 + 2 * i);
199
200 for (i = 0; i < n1; ++i)
201 real_cexp(i * p->twradix, n, p->W1 + 2 * i);
202
203 p->cexpl = cexpl_sqrtn_table;
204 p->rotate = rotate_sqrtn_table;
205 break;
206 }
207
208 case AWAKE_SINCOS:
209 p->cexpl = cexpl_sincos;
210 break;
211
212 case AWAKE_ZERO:
213 p->cexp = cexp_zero;
214 p->cexpl = cexpl_zero;
215 break;
216 }
217
218 if (!p->cexp) {
219 if (sizeof(trigreal) == sizeof(R))
220 p->cexp = (void (*)(triggen *, INT, R *))p->cexpl;
221 else
222 p->cexp = cexp_generic;
223 }
224 if (!p->rotate)
225 p->rotate = rotate_generic;
226 return p;
227 }
228
X(triggen_destroy)229 void X(triggen_destroy)(triggen *p)
230 {
231 X(ifree0)(p->W0);
232 X(ifree0)(p->W1);
233 X(ifree)(p);
234 }
235