1 /*
2 compile with cc -DTESTING_FFTPACK fftpack.c in order to build the
3 test application.
4
5 This is an f2c translation of the full fftpack sources as found on
6 http://www.netlib.org/fftpack/ The translated code has been
7 slightlty edited to remove the ugliest artefacts of the translation
8 (a hundred of wild GOTOs were wiped during that operation).
9
10 The original fftpack file was written by Paul N. Swarztrauber
11 (Version 4, 1985), in fortran 77.
12
13 FFTPACK license:
14
15 http://www.cisl.ucar.edu/css/software/fftpack5/ftpk.html
16
17 Copyright (c) 2004 the University Corporation for Atmospheric
18 Research ("UCAR"). All rights reserved. Developed by NCAR's
19 Computational and Information Systems Laboratory, UCAR,
20 www.cisl.ucar.edu.
21
22 Redistribution and use of the Software in source and binary forms,
23 with or without modification, is permitted provided that the
24 following conditions are met:
25
26 - Neither the names of NCAR's Computational and Information Systems
27 Laboratory, the University Corporation for Atmospheric Research,
28 nor the names of its sponsors or contributors may be used to
29 endorse or promote products derived from this Software without
30 specific prior written permission.
31
32 - Redistributions of source code must retain the above copyright
33 notices, this list of conditions, and the disclaimer below.
34
35 - Redistributions in binary form must reproduce the above copyright
36 notice, this list of conditions, and the disclaimer below in the
37 documentation and/or other materials provided with the
38 distribution.
39
40 THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
41 EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF
42 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
43 NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT
44 HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL,
45 EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN
46 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
47 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE
48 SOFTWARE.
49
50 ChangeLog:
51 2011/10/02: this is my first release of this file.
52 */
53
54 #include "fftpack.h"
55 #include <math.h>
56
57 typedef fftpack_real real;
58 typedef fftpack_int integer;
59
60 #ifndef FFTPACK_DOUBLE_PRECISION
61 #define FFTPACK_COS cosf
62 #define FFTPACK_SIN sinf
63 #else
64 #define FFTPACK_COS cos
65 #define FFTPACK_SIN sin
66 #endif
67
68
69 typedef struct f77complex {
70 real r, i;
71 } f77complex;
72
73 #ifdef TESTING_FFTPACK
c_abs(f77complex * c)74 static real c_abs(f77complex *c) { return sqrt(c->r*c->r + c->i*c->i); }
dmax(double a,double b)75 static double dmax(double a, double b) { return a < b ? b : a; }
76 #endif
77
78 /* define own constants required to turn off g++ extensions .. */
79 #ifndef M_PI
80 #define M_PI 3.14159265358979323846 /* pi */
81 #endif
82
83 #ifndef M_SQRT2
84 #define M_SQRT2 1.41421356237309504880 /* sqrt(2) */
85 #endif
86
87
88 /* translated by f2c (version 20061008), and slightly edited */
89
passfb(integer * nac,integer ido,integer ip,integer l1,integer idl1,real * cc,real * c1,real * c2,real * ch,real * ch2,const real * wa,real fsign)90 static void passfb(integer *nac, integer ido, integer ip, integer l1, integer idl1,
91 real *cc, real *c1, real *c2, real *ch, real *ch2, const real *wa, real fsign)
92 {
93 /* System generated locals */
94 integer ch_offset, cc_offset,
95 c1_offset, c2_offset, ch2_offset;
96
97 /* Local variables */
98 integer i, j, k, l, jc, lc, ik, idj, idl, inc, idp;
99 real wai, war;
100 integer ipp2, idij, idlj, idot, ipph;
101
102
103 #define c1_ref(a_1,a_2,a_3) c1[((a_3)*l1 + (a_2))*ido + a_1]
104 #define c2_ref(a_1,a_2) c2[(a_2)*idl1 + a_1]
105 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*ip + (a_2))*ido + a_1]
106 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
107 #define ch2_ref(a_1,a_2) ch2[(a_2)*idl1 + a_1]
108
109 /* Parameter adjustments */
110 ch_offset = 1 + ido * (1 + l1);
111 ch -= ch_offset;
112 c1_offset = 1 + ido * (1 + l1);
113 c1 -= c1_offset;
114 cc_offset = 1 + ido * (1 + ip);
115 cc -= cc_offset;
116 ch2_offset = 1 + idl1;
117 ch2 -= ch2_offset;
118 c2_offset = 1 + idl1;
119 c2 -= c2_offset;
120 --wa;
121
122 /* Function Body */
123 idot = ido / 2;
124 ipp2 = ip + 2;
125 ipph = (ip + 1) / 2;
126 idp = ip * ido;
127
128 if (ido >= l1) {
129 for (j = 2; j <= ipph; ++j) {
130 jc = ipp2 - j;
131 for (k = 1; k <= l1; ++k) {
132 for (i = 1; i <= ido; ++i) {
133 ch_ref(i, k, j) = cc_ref(i, j, k) + cc_ref(i, jc, k);
134 ch_ref(i, k, jc) = cc_ref(i, j, k) - cc_ref(i, jc, k);
135 }
136 }
137 }
138 for (k = 1; k <= l1; ++k) {
139 for (i = 1; i <= ido; ++i) {
140 ch_ref(i, k, 1) = cc_ref(i, 1, k);
141 }
142 }
143 } else {
144 for (j = 2; j <= ipph; ++j) {
145 jc = ipp2 - j;
146 for (i = 1; i <= ido; ++i) {
147 for (k = 1; k <= l1; ++k) {
148 ch_ref(i, k, j) = cc_ref(i, j, k) + cc_ref(i, jc, k);
149 ch_ref(i, k, jc) = cc_ref(i, j, k) - cc_ref(i, jc, k);
150 }
151 }
152 }
153 for (i = 1; i <= ido; ++i) {
154 for (k = 1; k <= l1; ++k) {
155 ch_ref(i, k, 1) = cc_ref(i, 1, k);
156 }
157 }
158 }
159 idl = 2 - ido;
160 inc = 0;
161 for (l = 2; l <= ipph; ++l) {
162 lc = ipp2 - l;
163 idl += ido;
164 for (ik = 1; ik <= idl1; ++ik) {
165 c2_ref(ik, l) = ch2_ref(ik, 1) + wa[idl - 1] * ch2_ref(ik, 2);
166 c2_ref(ik, lc) = fsign*wa[idl] * ch2_ref(ik, ip);
167 }
168 idlj = idl;
169 inc += ido;
170 for (j = 3; j <= ipph; ++j) {
171 jc = ipp2 - j;
172 idlj += inc;
173 if (idlj > idp) {
174 idlj -= idp;
175 }
176 war = wa[idlj - 1];
177 wai = wa[idlj];
178 for (ik = 1; ik <= idl1; ++ik) {
179 c2_ref(ik, l) = c2_ref(ik, l) + war * ch2_ref(ik, j);
180 c2_ref(ik, lc) = c2_ref(ik, lc) + fsign*wai * ch2_ref(ik, jc);
181 }
182 }
183 }
184 for (j = 2; j <= ipph; ++j) {
185 for (ik = 1; ik <= idl1; ++ik) {
186 ch2_ref(ik, 1) = ch2_ref(ik, 1) + ch2_ref(ik, j);
187 }
188 }
189 for (j = 2; j <= ipph; ++j) {
190 jc = ipp2 - j;
191 for (ik = 2; ik <= idl1; ik += 2) {
192 ch2_ref(ik - 1, j) = c2_ref(ik - 1, j) - c2_ref(ik, jc);
193 ch2_ref(ik - 1, jc) = c2_ref(ik - 1, j) + c2_ref(ik, jc);
194 ch2_ref(ik, j) = c2_ref(ik, j) + c2_ref(ik - 1, jc);
195 ch2_ref(ik, jc) = c2_ref(ik, j) - c2_ref(ik - 1, jc);
196 }
197 }
198 *nac = 1;
199 if (ido == 2) {
200 return;
201 }
202 *nac = 0;
203 for (ik = 1; ik <= idl1; ++ik) {
204 c2_ref(ik, 1) = ch2_ref(ik, 1);
205 }
206 for (j = 2; j <= ip; ++j) {
207 for (k = 1; k <= l1; ++k) {
208 c1_ref(1, k, j) = ch_ref(1, k, j);
209 c1_ref(2, k, j) = ch_ref(2, k, j);
210 }
211 }
212 if (idot <= l1) {
213 idij = 0;
214 for (j = 2; j <= ip; ++j) {
215 idij += 2;
216 for (i = 4; i <= ido; i += 2) {
217 idij += 2;
218 for (k = 1; k <= l1; ++k) {
219 c1_ref(i - 1, k, j) = wa[idij - 1] * ch_ref(i - 1, k, j) - fsign*wa[idij] * ch_ref(i, k, j);
220 c1_ref(i, k, j) = wa[idij - 1] * ch_ref(i, k, j) + fsign*wa[idij] * ch_ref(i - 1, k, j);
221 }
222 }
223 }
224 return;
225 }
226 idj = 2 - ido;
227 for (j = 2; j <= ip; ++j) {
228 idj += ido;
229 for (k = 1; k <= l1; ++k) {
230 idij = idj;
231 for (i = 4; i <= ido; i += 2) {
232 idij += 2;
233 c1_ref(i - 1, k, j) = wa[idij - 1] * ch_ref(i - 1, k, j) - fsign*wa[idij] * ch_ref(i, k, j);
234 c1_ref(i, k, j) = wa[idij - 1] * ch_ref(i, k, j) + fsign*wa[idij] * ch_ref(i - 1, k, j);
235 }
236 }
237 }
238 } /* passb */
239
240 #undef ch2_ref
241 #undef ch_ref
242 #undef cc_ref
243 #undef c2_ref
244 #undef c1_ref
245
246
passb2(integer ido,integer l1,const real * cc,real * ch,const real * wa1)247 static void passb2(integer ido, integer l1, const real *cc, real *ch, const real *wa1)
248 {
249 /* System generated locals */
250 integer cc_offset, ch_offset;
251
252 /* Local variables */
253 integer i, k;
254 real ti2, tr2;
255
256
257 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*2 + (a_2))*ido + a_1]
258 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
259
260 /* Parameter adjustments */
261 ch_offset = 1 + ido * (1 + l1);
262 ch -= ch_offset;
263 cc_offset = 1 + ido * 3;
264 cc -= cc_offset;
265 --wa1;
266
267 /* Function Body */
268 if (ido <= 2) {
269 for (k = 1; k <= l1; ++k) {
270 ch_ref(1, k, 1) = cc_ref(1, 1, k) + cc_ref(1, 2, k);
271 ch_ref(1, k, 2) = cc_ref(1, 1, k) - cc_ref(1, 2, k);
272 ch_ref(2, k, 1) = cc_ref(2, 1, k) + cc_ref(2, 2, k);
273 ch_ref(2, k, 2) = cc_ref(2, 1, k) - cc_ref(2, 2, k);
274 }
275 return;
276 }
277 for (k = 1; k <= l1; ++k) {
278 for (i = 2; i <= ido; i += 2) {
279 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + cc_ref(i - 1, 2, k);
280 tr2 = cc_ref(i - 1, 1, k) - cc_ref(i - 1, 2, k);
281 ch_ref(i, k, 1) = cc_ref(i, 1, k) + cc_ref(i, 2, k);
282 ti2 = cc_ref(i, 1, k) - cc_ref(i, 2, k);
283 ch_ref(i, k, 2) = wa1[i - 1] * ti2 + wa1[i] * tr2;
284 ch_ref(i - 1, k, 2) = wa1[i - 1] * tr2 - wa1[i] * ti2;
285 }
286 }
287 } /* passb2 */
288
289 #undef ch_ref
290 #undef cc_ref
291
292
passb3(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2)293 static void passb3(integer ido, integer l1, const real *cc, real *ch, const real *wa1, const real *wa2)
294 {
295 static const real taur = -.5f;
296 static const real taui = .866025403784439f;
297
298 /* System generated locals */
299 integer cc_offset, ch_offset;
300
301 /* Local variables */
302 integer i, k;
303 real ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
304
305
306 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*3 + (a_2))*ido + a_1]
307 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
308
309 /* Parameter adjustments */
310 ch_offset = 1 + ido * (1 + l1);
311 ch -= ch_offset;
312 cc_offset = 1 + (ido << 2);
313 cc -= cc_offset;
314 --wa1;
315 --wa2;
316
317 /* Function Body */
318 if (ido == 2) {
319 for (k = 1; k <= l1; ++k) {
320 tr2 = cc_ref(1, 2, k) + cc_ref(1, 3, k);
321 cr2 = cc_ref(1, 1, k) + taur * tr2;
322 ch_ref(1, k, 1) = cc_ref(1, 1, k) + tr2;
323 ti2 = cc_ref(2, 2, k) + cc_ref(2, 3, k);
324 ci2 = cc_ref(2, 1, k) + taur * ti2;
325 ch_ref(2, k, 1) = cc_ref(2, 1, k) + ti2;
326 cr3 = taui * (cc_ref(1, 2, k) - cc_ref(1, 3, k));
327 ci3 = taui * (cc_ref(2, 2, k) - cc_ref(2, 3, k));
328 ch_ref(1, k, 2) = cr2 - ci3;
329 ch_ref(1, k, 3) = cr2 + ci3;
330 ch_ref(2, k, 2) = ci2 + cr3;
331 ch_ref(2, k, 3) = ci2 - cr3;
332 }
333 } else {
334 for (k = 1; k <= l1; ++k) {
335 for (i = 2; i <= ido; i += 2) {
336 tr2 = cc_ref(i - 1, 2, k) + cc_ref(i - 1, 3, k);
337 cr2 = cc_ref(i - 1, 1, k) + taur * tr2;
338 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + tr2;
339 ti2 = cc_ref(i, 2, k) + cc_ref(i, 3, k);
340 ci2 = cc_ref(i, 1, k) + taur * ti2;
341 ch_ref(i, k, 1) = cc_ref(i, 1, k) + ti2;
342 cr3 = taui * (cc_ref(i - 1, 2, k) - cc_ref(i - 1, 3, k));
343 ci3 = taui * (cc_ref(i, 2, k) - cc_ref(i, 3, k));
344 dr2 = cr2 - ci3;
345 dr3 = cr2 + ci3;
346 di2 = ci2 + cr3;
347 di3 = ci2 - cr3;
348 ch_ref(i, k, 2) = wa1[i - 1] * di2 + wa1[i] * dr2;
349 ch_ref(i - 1, k, 2) = wa1[i - 1] * dr2 - wa1[i] * di2;
350 ch_ref(i, k, 3) = wa2[i - 1] * di3 + wa2[i] * dr3;
351 ch_ref(i - 1, k, 3) = wa2[i - 1] * dr3 - wa2[i] * di3;
352 }
353 }
354 }
355 } /* passb3 */
356
357 #undef ch_ref
358 #undef cc_ref
359
360
passb4(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2,const real * wa3)361 static void passb4(integer ido, integer l1, const real *cc, real *ch,
362 const real *wa1, const real *wa2, const real *wa3)
363 {
364 /* System generated locals */
365 integer cc_offset, ch_offset;
366
367 /* Local variables */
368 integer i, k;
369 real ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
370
371
372 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*4 + (a_2))*ido + a_1]
373 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
374
375 /* Parameter adjustments */
376 ch_offset = 1 + ido * (1 + l1);
377 ch -= ch_offset;
378 cc_offset = 1 + ido * 5;
379 cc -= cc_offset;
380 --wa1;
381 --wa2;
382 --wa3;
383
384 /* Function Body */
385 if (ido == 2) {
386 for (k = 1; k <= l1; ++k) {
387 ti1 = cc_ref(2, 1, k) - cc_ref(2, 3, k);
388 ti2 = cc_ref(2, 1, k) + cc_ref(2, 3, k);
389 tr4 = cc_ref(2, 4, k) - cc_ref(2, 2, k);
390 ti3 = cc_ref(2, 2, k) + cc_ref(2, 4, k);
391 tr1 = cc_ref(1, 1, k) - cc_ref(1, 3, k);
392 tr2 = cc_ref(1, 1, k) + cc_ref(1, 3, k);
393 ti4 = cc_ref(1, 2, k) - cc_ref(1, 4, k);
394 tr3 = cc_ref(1, 2, k) + cc_ref(1, 4, k);
395 ch_ref(1, k, 1) = tr2 + tr3;
396 ch_ref(1, k, 3) = tr2 - tr3;
397 ch_ref(2, k, 1) = ti2 + ti3;
398 ch_ref(2, k, 3) = ti2 - ti3;
399 ch_ref(1, k, 2) = tr1 + tr4;
400 ch_ref(1, k, 4) = tr1 - tr4;
401 ch_ref(2, k, 2) = ti1 + ti4;
402 ch_ref(2, k, 4) = ti1 - ti4;
403 }
404 } else {
405 for (k = 1; k <= l1; ++k) {
406 for (i = 2; i <= ido; i += 2) {
407 ti1 = cc_ref(i, 1, k) - cc_ref(i, 3, k);
408 ti2 = cc_ref(i, 1, k) + cc_ref(i, 3, k);
409 ti3 = cc_ref(i, 2, k) + cc_ref(i, 4, k);
410 tr4 = cc_ref(i, 4, k) - cc_ref(i, 2, k);
411 tr1 = cc_ref(i - 1, 1, k) - cc_ref(i - 1, 3, k);
412 tr2 = cc_ref(i - 1, 1, k) + cc_ref(i - 1, 3, k);
413 ti4 = cc_ref(i - 1, 2, k) - cc_ref(i - 1, 4, k);
414 tr3 = cc_ref(i - 1, 2, k) + cc_ref(i - 1, 4, k);
415 ch_ref(i - 1, k, 1) = tr2 + tr3;
416 cr3 = tr2 - tr3;
417 ch_ref(i, k, 1) = ti2 + ti3;
418 ci3 = ti2 - ti3;
419 cr2 = tr1 + tr4;
420 cr4 = tr1 - tr4;
421 ci2 = ti1 + ti4;
422 ci4 = ti1 - ti4;
423 ch_ref(i - 1, k, 2) = wa1[i - 1] * cr2 - wa1[i] * ci2;
424 ch_ref(i, k, 2) = wa1[i - 1] * ci2 + wa1[i] * cr2;
425 ch_ref(i - 1, k, 3) = wa2[i - 1] * cr3 - wa2[i] * ci3;
426 ch_ref(i, k, 3) = wa2[i - 1] * ci3 + wa2[i] * cr3;
427 ch_ref(i - 1, k, 4) = wa3[i - 1] * cr4 - wa3[i] * ci4;
428 ch_ref(i, k, 4) = wa3[i - 1] * ci4 + wa3[i] * cr4;
429 }
430 }
431 }
432 } /* passb4 */
433
434 #undef ch_ref
435 #undef cc_ref
436
437 /* passf5 and passb5 merged */
passfb5(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2,const real * wa3,const real * wa4,real fsign)438 static void passfb5(integer ido, integer l1, const real *cc, real *ch,
439 const real *wa1, const real *wa2, const real *wa3, const real *wa4, real fsign)
440 {
441 const real tr11 = .309016994374947f;
442 const real ti11 = .951056516295154f*fsign;
443 const real tr12 = -.809016994374947f;
444 const real ti12 = .587785252292473f*fsign;
445
446 /* System generated locals */
447 integer cc_offset, ch_offset;
448
449 /* Local variables */
450 integer i, k;
451 real ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
452 ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
453
454
455 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*5 + (a_2))*ido + a_1]
456 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
457
458 /* Parameter adjustments */
459 ch_offset = 1 + ido * (1 + l1);
460 ch -= ch_offset;
461 cc_offset = 1 + ido * 6;
462 cc -= cc_offset;
463 --wa1;
464 --wa2;
465 --wa3;
466 --wa4;
467
468 /* Function Body */
469 if (ido == 2) {
470 for (k = 1; k <= l1; ++k) {
471 ti5 = cc_ref(2, 2, k) - cc_ref(2, 5, k);
472 ti2 = cc_ref(2, 2, k) + cc_ref(2, 5, k);
473 ti4 = cc_ref(2, 3, k) - cc_ref(2, 4, k);
474 ti3 = cc_ref(2, 3, k) + cc_ref(2, 4, k);
475 tr5 = cc_ref(1, 2, k) - cc_ref(1, 5, k);
476 tr2 = cc_ref(1, 2, k) + cc_ref(1, 5, k);
477 tr4 = cc_ref(1, 3, k) - cc_ref(1, 4, k);
478 tr3 = cc_ref(1, 3, k) + cc_ref(1, 4, k);
479 ch_ref(1, k, 1) = cc_ref(1, 1, k) + tr2 + tr3;
480 ch_ref(2, k, 1) = cc_ref(2, 1, k) + ti2 + ti3;
481 cr2 = cc_ref(1, 1, k) + tr11 * tr2 + tr12 * tr3;
482 ci2 = cc_ref(2, 1, k) + tr11 * ti2 + tr12 * ti3;
483 cr3 = cc_ref(1, 1, k) + tr12 * tr2 + tr11 * tr3;
484 ci3 = cc_ref(2, 1, k) + tr12 * ti2 + tr11 * ti3;
485 cr5 = ti11 * tr5 + ti12 * tr4;
486 ci5 = ti11 * ti5 + ti12 * ti4;
487 cr4 = ti12 * tr5 - ti11 * tr4;
488 ci4 = ti12 * ti5 - ti11 * ti4;
489 ch_ref(1, k, 2) = cr2 - ci5;
490 ch_ref(1, k, 5) = cr2 + ci5;
491 ch_ref(2, k, 2) = ci2 + cr5;
492 ch_ref(2, k, 3) = ci3 + cr4;
493 ch_ref(1, k, 3) = cr3 - ci4;
494 ch_ref(1, k, 4) = cr3 + ci4;
495 ch_ref(2, k, 4) = ci3 - cr4;
496 ch_ref(2, k, 5) = ci2 - cr5;
497 }
498 } else {
499 for (k = 1; k <= l1; ++k) {
500 for (i = 2; i <= ido; i += 2) {
501 ti5 = cc_ref(i, 2, k) - cc_ref(i, 5, k);
502 ti2 = cc_ref(i, 2, k) + cc_ref(i, 5, k);
503 ti4 = cc_ref(i, 3, k) - cc_ref(i, 4, k);
504 ti3 = cc_ref(i, 3, k) + cc_ref(i, 4, k);
505 tr5 = cc_ref(i - 1, 2, k) - cc_ref(i - 1, 5, k);
506 tr2 = cc_ref(i - 1, 2, k) + cc_ref(i - 1, 5, k);
507 tr4 = cc_ref(i - 1, 3, k) - cc_ref(i - 1, 4, k);
508 tr3 = cc_ref(i - 1, 3, k) + cc_ref(i - 1, 4, k);
509 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + tr2 + tr3;
510 ch_ref(i, k, 1) = cc_ref(i, 1, k) + ti2 + ti3;
511 cr2 = cc_ref(i - 1, 1, k) + tr11 * tr2 + tr12 * tr3;
512 ci2 = cc_ref(i, 1, k) + tr11 * ti2 + tr12 * ti3;
513 cr3 = cc_ref(i - 1, 1, k) + tr12 * tr2 + tr11 * tr3;
514 ci3 = cc_ref(i, 1, k) + tr12 * ti2 + tr11 * ti3;
515 cr5 = ti11 * tr5 + ti12 * tr4;
516 ci5 = ti11 * ti5 + ti12 * ti4;
517 cr4 = ti12 * tr5 - ti11 * tr4;
518 ci4 = ti12 * ti5 - ti11 * ti4;
519 dr3 = cr3 - ci4;
520 dr4 = cr3 + ci4;
521 di3 = ci3 + cr4;
522 di4 = ci3 - cr4;
523 dr5 = cr2 + ci5;
524 dr2 = cr2 - ci5;
525 di5 = ci2 - cr5;
526 di2 = ci2 + cr5;
527 ch_ref(i - 1, k, 2) = wa1[i - 1] * dr2 - fsign*wa1[i] * di2;
528 ch_ref(i, k, 2) = wa1[i - 1] * di2 + fsign*wa1[i] * dr2;
529 ch_ref(i - 1, k, 3) = wa2[i - 1] * dr3 - fsign*wa2[i] * di3;
530 ch_ref(i, k, 3) = wa2[i - 1] * di3 + fsign*wa2[i] * dr3;
531 ch_ref(i - 1, k, 4) = wa3[i - 1] * dr4 - fsign*wa3[i] * di4;
532 ch_ref(i, k, 4) = wa3[i - 1] * di4 + fsign*wa3[i] * dr4;
533 ch_ref(i - 1, k, 5) = wa4[i - 1] * dr5 - fsign*wa4[i] * di5;
534 ch_ref(i, k, 5) = wa4[i - 1] * di5 + fsign*wa4[i] * dr5;
535 }
536 }
537 }
538 } /* passb5 */
539
540 #undef ch_ref
541 #undef cc_ref
542
passf2(integer ido,integer l1,const real * cc,real * ch,const real * wa1)543 static void passf2(integer ido, integer l1, const real *cc, real *ch, const real *wa1)
544 {
545 /* System generated locals */
546 integer cc_offset, ch_offset;
547
548 /* Local variables */
549 integer i, k;
550 real ti2, tr2;
551
552
553 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*2 + (a_2))*ido + a_1]
554 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
555
556 /* Parameter adjustments */
557 ch_offset = 1 + ido * (1 + l1);
558 ch -= ch_offset;
559 cc_offset = 1 + ido * 3;
560 cc -= cc_offset;
561 --wa1;
562
563 /* Function Body */
564 if (ido == 2) {
565 for (k = 1; k <= l1; ++k) {
566 ch_ref(1, k, 1) = cc_ref(1, 1, k) + cc_ref(1, 2, k);
567 ch_ref(1, k, 2) = cc_ref(1, 1, k) - cc_ref(1, 2, k);
568 ch_ref(2, k, 1) = cc_ref(2, 1, k) + cc_ref(2, 2, k);
569 ch_ref(2, k, 2) = cc_ref(2, 1, k) - cc_ref(2, 2, k);
570 }
571 } else {
572 for (k = 1; k <= l1; ++k) {
573 for (i = 2; i <= ido; i += 2) {
574 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + cc_ref(i - 1, 2,
575 k);
576 tr2 = cc_ref(i - 1, 1, k) - cc_ref(i - 1, 2, k);
577 ch_ref(i, k, 1) = cc_ref(i, 1, k) + cc_ref(i, 2, k);
578 ti2 = cc_ref(i, 1, k) - cc_ref(i, 2, k);
579 ch_ref(i, k, 2) = wa1[i - 1] * ti2 - wa1[i] * tr2;
580 ch_ref(i - 1, k, 2) = wa1[i - 1] * tr2 + wa1[i] * ti2;
581 }
582 }
583 }
584 } /* passf2 */
585
586 #undef ch_ref
587 #undef cc_ref
588
589
passf3(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2)590 static void passf3(integer ido, integer l1, const real *cc, real *ch,
591 const real *wa1, const real *wa2)
592 {
593 static const real taur = -.5f;
594 static const real taui = -.866025403784439f;
595
596 /* System generated locals */
597 integer cc_offset, ch_offset;
598
599 /* Local variables */
600 integer i, k;
601 real ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
602
603
604 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*3 + (a_2))*ido + a_1]
605 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
606
607 /* Parameter adjustments */
608 ch_offset = 1 + ido * (1 + l1);
609 ch -= ch_offset;
610 cc_offset = 1 + (ido << 2);
611 cc -= cc_offset;
612 --wa1;
613 --wa2;
614
615 /* Function Body */
616 if (ido == 2) {
617 for (k = 1; k <= l1; ++k) {
618 tr2 = cc_ref(1, 2, k) + cc_ref(1, 3, k);
619 cr2 = cc_ref(1, 1, k) + taur * tr2;
620 ch_ref(1, k, 1) = cc_ref(1, 1, k) + tr2;
621 ti2 = cc_ref(2, 2, k) + cc_ref(2, 3, k);
622 ci2 = cc_ref(2, 1, k) + taur * ti2;
623 ch_ref(2, k, 1) = cc_ref(2, 1, k) + ti2;
624 cr3 = taui * (cc_ref(1, 2, k) - cc_ref(1, 3, k));
625 ci3 = taui * (cc_ref(2, 2, k) - cc_ref(2, 3, k));
626 ch_ref(1, k, 2) = cr2 - ci3;
627 ch_ref(1, k, 3) = cr2 + ci3;
628 ch_ref(2, k, 2) = ci2 + cr3;
629 ch_ref(2, k, 3) = ci2 - cr3;
630 }
631 } else {
632 for (k = 1; k <= l1; ++k) {
633 for (i = 2; i <= ido; i += 2) {
634 tr2 = cc_ref(i - 1, 2, k) + cc_ref(i - 1, 3, k);
635 cr2 = cc_ref(i - 1, 1, k) + taur * tr2;
636 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + tr2;
637 ti2 = cc_ref(i, 2, k) + cc_ref(i, 3, k);
638 ci2 = cc_ref(i, 1, k) + taur * ti2;
639 ch_ref(i, k, 1) = cc_ref(i, 1, k) + ti2;
640 cr3 = taui * (cc_ref(i - 1, 2, k) - cc_ref(i - 1, 3, k));
641 ci3 = taui * (cc_ref(i, 2, k) - cc_ref(i, 3, k));
642 dr2 = cr2 - ci3;
643 dr3 = cr2 + ci3;
644 di2 = ci2 + cr3;
645 di3 = ci2 - cr3;
646 ch_ref(i, k, 2) = wa1[i - 1] * di2 - wa1[i] * dr2;
647 ch_ref(i - 1, k, 2) = wa1[i - 1] * dr2 + wa1[i] * di2;
648 ch_ref(i, k, 3) = wa2[i - 1] * di3 - wa2[i] * dr3;
649 ch_ref(i - 1, k, 3) = wa2[i - 1] * dr3 + wa2[i] * di3;
650 }
651 }
652 }
653 } /* passf3 */
654
655 #undef ch_ref
656 #undef cc_ref
657
658
passf4(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2,const real * wa3)659 static void passf4(integer ido, integer l1, const real *cc, real *ch,
660 const real *wa1, const real *wa2, const real *wa3)
661 {
662 /* System generated locals */
663 integer cc_offset, ch_offset;
664
665 /* Local variables */
666 integer i, k;
667 real ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
668
669
670 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*4 + (a_2))*ido + a_1]
671 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
672
673 /* Parameter adjustments */
674 ch_offset = 1 + ido * (1 + l1);
675 ch -= ch_offset;
676 cc_offset = 1 + ido * 5;
677 cc -= cc_offset;
678 --wa1;
679 --wa2;
680 --wa3;
681
682 /* Function Body */
683 if (ido == 2) {
684 for (k = 1; k <= l1; ++k) {
685 ti1 = cc_ref(2, 1, k) - cc_ref(2, 3, k);
686 ti2 = cc_ref(2, 1, k) + cc_ref(2, 3, k);
687 tr4 = cc_ref(2, 2, k) - cc_ref(2, 4, k);
688 ti3 = cc_ref(2, 2, k) + cc_ref(2, 4, k);
689 tr1 = cc_ref(1, 1, k) - cc_ref(1, 3, k);
690 tr2 = cc_ref(1, 1, k) + cc_ref(1, 3, k);
691 ti4 = cc_ref(1, 4, k) - cc_ref(1, 2, k);
692 tr3 = cc_ref(1, 2, k) + cc_ref(1, 4, k);
693 ch_ref(1, k, 1) = tr2 + tr3;
694 ch_ref(1, k, 3) = tr2 - tr3;
695 ch_ref(2, k, 1) = ti2 + ti3;
696 ch_ref(2, k, 3) = ti2 - ti3;
697 ch_ref(1, k, 2) = tr1 + tr4;
698 ch_ref(1, k, 4) = tr1 - tr4;
699 ch_ref(2, k, 2) = ti1 + ti4;
700 ch_ref(2, k, 4) = ti1 - ti4;
701 }
702 } else {
703 for (k = 1; k <= l1; ++k) {
704 for (i = 2; i <= ido; i += 2) {
705 ti1 = cc_ref(i, 1, k) - cc_ref(i, 3, k);
706 ti2 = cc_ref(i, 1, k) + cc_ref(i, 3, k);
707 ti3 = cc_ref(i, 2, k) + cc_ref(i, 4, k);
708 tr4 = cc_ref(i, 2, k) - cc_ref(i, 4, k);
709 tr1 = cc_ref(i - 1, 1, k) - cc_ref(i - 1, 3, k);
710 tr2 = cc_ref(i - 1, 1, k) + cc_ref(i - 1, 3, k);
711 ti4 = cc_ref(i - 1, 4, k) - cc_ref(i - 1, 2, k);
712 tr3 = cc_ref(i - 1, 2, k) + cc_ref(i - 1, 4, k);
713 ch_ref(i - 1, k, 1) = tr2 + tr3;
714 cr3 = tr2 - tr3;
715 ch_ref(i, k, 1) = ti2 + ti3;
716 ci3 = ti2 - ti3;
717 cr2 = tr1 + tr4;
718 cr4 = tr1 - tr4;
719 ci2 = ti1 + ti4;
720 ci4 = ti1 - ti4;
721 ch_ref(i - 1, k, 2) = wa1[i - 1] * cr2 + wa1[i] * ci2;
722 ch_ref(i, k, 2) = wa1[i - 1] * ci2 - wa1[i] * cr2;
723 ch_ref(i - 1, k, 3) = wa2[i - 1] * cr3 + wa2[i] * ci3;
724 ch_ref(i, k, 3) = wa2[i - 1] * ci3 - wa2[i] * cr3;
725 ch_ref(i - 1, k, 4) = wa3[i - 1] * cr4 + wa3[i] * ci4;
726 ch_ref(i, k, 4) = wa3[i - 1] * ci4 - wa3[i] * cr4;
727 }
728 }
729 }
730 } /* passf4 */
731
732 #undef ch_ref
733 #undef cc_ref
734
radb2(integer ido,integer l1,const real * cc,real * ch,const real * wa1)735 static void radb2(integer ido, integer l1, const real *cc, real *ch, const real *wa1)
736 {
737 /* System generated locals */
738 integer cc_offset, ch_offset;
739
740 /* Local variables */
741 integer i, k, ic;
742 real ti2, tr2;
743 integer idp2;
744
745
746 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*2 + (a_2))*ido + a_1]
747 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
748
749 /* Parameter adjustments */
750 ch_offset = 1 + ido * (1 + l1);
751 ch -= ch_offset;
752 cc_offset = 1 + ido * 3;
753 cc -= cc_offset;
754 --wa1;
755
756 /* Function Body */
757 for (k = 1; k <= l1; ++k) {
758 ch_ref(1, k, 1) = cc_ref(1, 1, k) + cc_ref(ido, 2, k);
759 ch_ref(1, k, 2) = cc_ref(1, 1, k) - cc_ref(ido, 2, k);
760 }
761 if (ido < 2) return;
762 else if (ido != 2) {
763 idp2 = ido + 2;
764 for (k = 1; k <= l1; ++k) {
765 for (i = 3; i <= ido; i += 2) {
766 ic = idp2 - i;
767 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + cc_ref(ic - 1, 2,
768 k);
769 tr2 = cc_ref(i - 1, 1, k) - cc_ref(ic - 1, 2, k);
770 ch_ref(i, k, 1) = cc_ref(i, 1, k) - cc_ref(ic, 2, k);
771 ti2 = cc_ref(i, 1, k) + cc_ref(ic, 2, k);
772 ch_ref(i - 1, k, 2) = wa1[i - 2] * tr2 - wa1[i - 1] * ti2;
773 ch_ref(i, k, 2) = wa1[i - 2] * ti2 + wa1[i - 1] * tr2;
774 }
775 }
776 if (ido % 2 == 1) return;
777 }
778 for (k = 1; k <= l1; ++k) {
779 ch_ref(ido, k, 1) = cc_ref(ido, 1, k) + cc_ref(ido, 1, k);
780 ch_ref(ido, k, 2) = -(cc_ref(1, 2, k) + cc_ref(1, 2, k));
781 }
782 } /* radb2 */
783
784 #undef ch_ref
785 #undef cc_ref
786
787
radb3(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2)788 static void radb3(integer ido, integer l1, const real *cc, real *ch,
789 const real *wa1, const real *wa2)
790 {
791 /* Initialized data */
792
793 static const real taur = -.5f;
794 static const real taui = .866025403784439f;
795
796 /* System generated locals */
797 integer cc_offset, ch_offset;
798
799 /* Local variables */
800 integer i, k, ic;
801 real ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
802 integer idp2;
803
804
805 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*3 + (a_2))*ido + a_1]
806 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
807
808 /* Parameter adjustments */
809 ch_offset = 1 + ido * (1 + l1);
810 ch -= ch_offset;
811 cc_offset = 1 + (ido << 2);
812 cc -= cc_offset;
813 --wa1;
814 --wa2;
815
816 /* Function Body */
817 for (k = 1; k <= l1; ++k) {
818 tr2 = cc_ref(ido, 2, k) + cc_ref(ido, 2, k);
819 cr2 = cc_ref(1, 1, k) + taur * tr2;
820 ch_ref(1, k, 1) = cc_ref(1, 1, k) + tr2;
821 ci3 = taui * (cc_ref(1, 3, k) + cc_ref(1, 3, k));
822 ch_ref(1, k, 2) = cr2 - ci3;
823 ch_ref(1, k, 3) = cr2 + ci3;
824 }
825 if (ido == 1) {
826 return;
827 }
828 idp2 = ido + 2;
829 for (k = 1; k <= l1; ++k) {
830 for (i = 3; i <= ido; i += 2) {
831 ic = idp2 - i;
832 tr2 = cc_ref(i - 1, 3, k) + cc_ref(ic - 1, 2, k);
833 cr2 = cc_ref(i - 1, 1, k) + taur * tr2;
834 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + tr2;
835 ti2 = cc_ref(i, 3, k) - cc_ref(ic, 2, k);
836 ci2 = cc_ref(i, 1, k) + taur * ti2;
837 ch_ref(i, k, 1) = cc_ref(i, 1, k) + ti2;
838 cr3 = taui * (cc_ref(i - 1, 3, k) - cc_ref(ic - 1, 2, k));
839 ci3 = taui * (cc_ref(i, 3, k) + cc_ref(ic, 2, k));
840 dr2 = cr2 - ci3;
841 dr3 = cr2 + ci3;
842 di2 = ci2 + cr3;
843 di3 = ci2 - cr3;
844 ch_ref(i - 1, k, 2) = wa1[i - 2] * dr2 - wa1[i - 1] * di2;
845 ch_ref(i, k, 2) = wa1[i - 2] * di2 + wa1[i - 1] * dr2;
846 ch_ref(i - 1, k, 3) = wa2[i - 2] * dr3 - wa2[i - 1] * di3;
847 ch_ref(i, k, 3) = wa2[i - 2] * di3 + wa2[i - 1] * dr3;
848 }
849 }
850 } /* radb3 */
851
852 #undef ch_ref
853 #undef cc_ref
854
855
radb4(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2,const real * wa3)856 static void radb4(integer ido, integer l1, const real *cc, real *ch,
857 const real *wa1, const real *wa2, const real *wa3)
858 {
859 /* Initialized data */
860
861 static const real sqrt2 = 1.414213562373095f;
862
863 /* System generated locals */
864 integer cc_offset, ch_offset;
865
866 /* Local variables */
867 integer i, k, ic;
868 real ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
869 integer idp2;
870
871
872 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*4 + (a_2))*ido + a_1]
873 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
874
875 /* Parameter adjustments */
876 ch_offset = 1 + ido * (1 + l1);
877 ch -= ch_offset;
878 cc_offset = 1 + ido * 5;
879 cc -= cc_offset;
880 --wa1;
881 --wa2;
882 --wa3;
883
884 /* Function Body */
885 for (k = 1; k <= l1; ++k) {
886 tr1 = cc_ref(1, 1, k) - cc_ref(ido, 4, k);
887 tr2 = cc_ref(1, 1, k) + cc_ref(ido, 4, k);
888 tr3 = cc_ref(ido, 2, k) + cc_ref(ido, 2, k);
889 tr4 = cc_ref(1, 3, k) + cc_ref(1, 3, k);
890 ch_ref(1, k, 1) = tr2 + tr3;
891 ch_ref(1, k, 2) = tr1 - tr4;
892 ch_ref(1, k, 3) = tr2 - tr3;
893 ch_ref(1, k, 4) = tr1 + tr4;
894 }
895 if (ido < 2) return;
896 if (ido != 2) {
897 idp2 = ido + 2;
898 for (k = 1; k <= l1; ++k) {
899 for (i = 3; i <= ido; i += 2) {
900 ic = idp2 - i;
901 ti1 = cc_ref(i, 1, k) + cc_ref(ic, 4, k);
902 ti2 = cc_ref(i, 1, k) - cc_ref(ic, 4, k);
903 ti3 = cc_ref(i, 3, k) - cc_ref(ic, 2, k);
904 tr4 = cc_ref(i, 3, k) + cc_ref(ic, 2, k);
905 tr1 = cc_ref(i - 1, 1, k) - cc_ref(ic - 1, 4, k);
906 tr2 = cc_ref(i - 1, 1, k) + cc_ref(ic - 1, 4, k);
907 ti4 = cc_ref(i - 1, 3, k) - cc_ref(ic - 1, 2, k);
908 tr3 = cc_ref(i - 1, 3, k) + cc_ref(ic - 1, 2, k);
909 ch_ref(i - 1, k, 1) = tr2 + tr3;
910 cr3 = tr2 - tr3;
911 ch_ref(i, k, 1) = ti2 + ti3;
912 ci3 = ti2 - ti3;
913 cr2 = tr1 - tr4;
914 cr4 = tr1 + tr4;
915 ci2 = ti1 + ti4;
916 ci4 = ti1 - ti4;
917 ch_ref(i - 1, k, 2) = wa1[i - 2] * cr2 - wa1[i - 1] * ci2;
918 ch_ref(i, k, 2) = wa1[i - 2] * ci2 + wa1[i - 1] * cr2;
919 ch_ref(i - 1, k, 3) = wa2[i - 2] * cr3 - wa2[i - 1] * ci3;
920 ch_ref(i, k, 3) = wa2[i - 2] * ci3 + wa2[i - 1] * cr3;
921 ch_ref(i - 1, k, 4) = wa3[i - 2] * cr4 - wa3[i - 1] * ci4;
922 ch_ref(i, k, 4) = wa3[i - 2] * ci4 + wa3[i - 1] * cr4;
923 }
924 }
925 if (ido % 2 == 1) return;
926 }
927 for (k = 1; k <= l1; ++k) {
928 ti1 = cc_ref(1, 2, k) + cc_ref(1, 4, k);
929 ti2 = cc_ref(1, 4, k) - cc_ref(1, 2, k);
930 tr1 = cc_ref(ido, 1, k) - cc_ref(ido, 3, k);
931 tr2 = cc_ref(ido, 1, k) + cc_ref(ido, 3, k);
932 ch_ref(ido, k, 1) = tr2 + tr2;
933 ch_ref(ido, k, 2) = sqrt2 * (tr1 - ti1);
934 ch_ref(ido, k, 3) = ti2 + ti2;
935 ch_ref(ido, k, 4) = -sqrt2 * (tr1 + ti1);
936 }
937 } /* radb4 */
938
939 #undef ch_ref
940 #undef cc_ref
941
942
radb5(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2,const real * wa3,const real * wa4)943 static void radb5(integer ido, integer l1, const real *cc, real *ch,
944 const real *wa1, const real *wa2, const real *wa3, const real *wa4)
945 {
946 /* Initialized data */
947
948 static const real tr11 = .309016994374947f;
949 static const real ti11 = .951056516295154f;
950 static const real tr12 = -.809016994374947f;
951 static const real ti12 = .587785252292473f;
952
953 /* System generated locals */
954 integer cc_offset, ch_offset;
955
956 /* Local variables */
957 integer i, k, ic;
958 real ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
959 ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
960 integer idp2;
961
962
963 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*5 + (a_2))*ido + a_1]
964 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
965
966 /* Parameter adjustments */
967 ch_offset = 1 + ido * (1 + l1);
968 ch -= ch_offset;
969 cc_offset = 1 + ido * 6;
970 cc -= cc_offset;
971 --wa1;
972 --wa2;
973 --wa3;
974 --wa4;
975
976 /* Function Body */
977 for (k = 1; k <= l1; ++k) {
978 ti5 = cc_ref(1, 3, k) + cc_ref(1, 3, k);
979 ti4 = cc_ref(1, 5, k) + cc_ref(1, 5, k);
980 tr2 = cc_ref(ido, 2, k) + cc_ref(ido, 2, k);
981 tr3 = cc_ref(ido, 4, k) + cc_ref(ido, 4, k);
982 ch_ref(1, k, 1) = cc_ref(1, 1, k) + tr2 + tr3;
983 cr2 = cc_ref(1, 1, k) + tr11 * tr2 + tr12 * tr3;
984 cr3 = cc_ref(1, 1, k) + tr12 * tr2 + tr11 * tr3;
985 ci5 = ti11 * ti5 + ti12 * ti4;
986 ci4 = ti12 * ti5 - ti11 * ti4;
987 ch_ref(1, k, 2) = cr2 - ci5;
988 ch_ref(1, k, 3) = cr3 - ci4;
989 ch_ref(1, k, 4) = cr3 + ci4;
990 ch_ref(1, k, 5) = cr2 + ci5;
991 }
992 if (ido == 1) {
993 return;
994 }
995 idp2 = ido + 2;
996 for (k = 1; k <= l1; ++k) {
997 for (i = 3; i <= ido; i += 2) {
998 ic = idp2 - i;
999 ti5 = cc_ref(i, 3, k) + cc_ref(ic, 2, k);
1000 ti2 = cc_ref(i, 3, k) - cc_ref(ic, 2, k);
1001 ti4 = cc_ref(i, 5, k) + cc_ref(ic, 4, k);
1002 ti3 = cc_ref(i, 5, k) - cc_ref(ic, 4, k);
1003 tr5 = cc_ref(i - 1, 3, k) - cc_ref(ic - 1, 2, k);
1004 tr2 = cc_ref(i - 1, 3, k) + cc_ref(ic - 1, 2, k);
1005 tr4 = cc_ref(i - 1, 5, k) - cc_ref(ic - 1, 4, k);
1006 tr3 = cc_ref(i - 1, 5, k) + cc_ref(ic - 1, 4, k);
1007 ch_ref(i - 1, k, 1) = cc_ref(i - 1, 1, k) + tr2 + tr3;
1008 ch_ref(i, k, 1) = cc_ref(i, 1, k) + ti2 + ti3;
1009 cr2 = cc_ref(i - 1, 1, k) + tr11 * tr2 + tr12 * tr3;
1010 ci2 = cc_ref(i, 1, k) + tr11 * ti2 + tr12 * ti3;
1011 cr3 = cc_ref(i - 1, 1, k) + tr12 * tr2 + tr11 * tr3;
1012 ci3 = cc_ref(i, 1, k) + tr12 * ti2 + tr11 * ti3;
1013 cr5 = ti11 * tr5 + ti12 * tr4;
1014 ci5 = ti11 * ti5 + ti12 * ti4;
1015 cr4 = ti12 * tr5 - ti11 * tr4;
1016 ci4 = ti12 * ti5 - ti11 * ti4;
1017 dr3 = cr3 - ci4;
1018 dr4 = cr3 + ci4;
1019 di3 = ci3 + cr4;
1020 di4 = ci3 - cr4;
1021 dr5 = cr2 + ci5;
1022 dr2 = cr2 - ci5;
1023 di5 = ci2 - cr5;
1024 di2 = ci2 + cr5;
1025 ch_ref(i - 1, k, 2) = wa1[i - 2] * dr2 - wa1[i - 1] * di2;
1026 ch_ref(i, k, 2) = wa1[i - 2] * di2 + wa1[i - 1] * dr2;
1027 ch_ref(i - 1, k, 3) = wa2[i - 2] * dr3 - wa2[i - 1] * di3;
1028 ch_ref(i, k, 3) = wa2[i - 2] * di3 + wa2[i - 1] * dr3;
1029 ch_ref(i - 1, k, 4) = wa3[i - 2] * dr4 - wa3[i - 1] * di4;
1030 ch_ref(i, k, 4) = wa3[i - 2] * di4 + wa3[i - 1] * dr4;
1031 ch_ref(i - 1, k, 5) = wa4[i - 2] * dr5 - wa4[i - 1] * di5;
1032 ch_ref(i, k, 5) = wa4[i - 2] * di5 + wa4[i - 1] * dr5;
1033 }
1034 }
1035 } /* radb5 */
1036
1037 #undef ch_ref
1038 #undef cc_ref
1039
1040
radbg(integer ido,integer ip,integer l1,integer idl1,const real * cc,real * c1,real * c2,real * ch,real * ch2,const real * wa)1041 static void radbg(integer ido, integer ip, integer l1, integer idl1,
1042 const real *cc, real *c1, real *c2, real *ch, real *ch2, const real *wa)
1043 {
1044 /* System generated locals */
1045 integer ch_offset, cc_offset,
1046 c1_offset, c2_offset, ch2_offset;
1047
1048 /* Local variables */
1049 integer i, j, k, l, j2, ic, jc, lc, ik, is;
1050 real dc2, ai1, ai2, ar1, ar2, ds2;
1051 integer nbd;
1052 real dcp, arg, dsp, ar1h, ar2h;
1053 integer idp2, ipp2, idij, ipph;
1054
1055
1056 #define c1_ref(a_1,a_2,a_3) c1[((a_3)*l1 + (a_2))*ido + a_1]
1057 #define c2_ref(a_1,a_2) c2[(a_2)*idl1 + a_1]
1058 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*ip + (a_2))*ido + a_1]
1059 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
1060 #define ch2_ref(a_1,a_2) ch2[(a_2)*idl1 + a_1]
1061
1062 /* Parameter adjustments */
1063 ch_offset = 1 + ido * (1 + l1);
1064 ch -= ch_offset;
1065 c1_offset = 1 + ido * (1 + l1);
1066 c1 -= c1_offset;
1067 cc_offset = 1 + ido * (1 + ip);
1068 cc -= cc_offset;
1069 ch2_offset = 1 + idl1;
1070 ch2 -= ch2_offset;
1071 c2_offset = 1 + idl1;
1072 c2 -= c2_offset;
1073 --wa;
1074
1075 /* Function Body */
1076 arg = (2*M_PI) / (real) (ip);
1077 dcp = FFTPACK_COS(arg);
1078 dsp = FFTPACK_SIN(arg);
1079 idp2 = ido + 2;
1080 nbd = (ido - 1) / 2;
1081 ipp2 = ip + 2;
1082 ipph = (ip + 1) / 2;
1083 if (ido >= l1) {
1084 for (k = 1; k <= l1; ++k) {
1085 for (i = 1; i <= ido; ++i) {
1086 ch_ref(i, k, 1) = cc_ref(i, 1, k);
1087 }
1088 }
1089 } else {
1090 for (i = 1; i <= ido; ++i) {
1091 for (k = 1; k <= l1; ++k) {
1092 ch_ref(i, k, 1) = cc_ref(i, 1, k);
1093 }
1094 }
1095 }
1096 for (j = 2; j <= ipph; ++j) {
1097 jc = ipp2 - j;
1098 j2 = j + j;
1099 for (k = 1; k <= l1; ++k) {
1100 ch_ref(1, k, j) = cc_ref(ido, j2 - 2, k) + cc_ref(ido, j2 - 2, k);
1101 ch_ref(1, k, jc) = cc_ref(1, j2 - 1, k) + cc_ref(1, j2 - 1, k);
1102 }
1103 }
1104 if (ido != 1) {
1105 if (nbd >= l1) {
1106 for (j = 2; j <= ipph; ++j) {
1107 jc = ipp2 - j;
1108 for (k = 1; k <= l1; ++k) {
1109 for (i = 3; i <= ido; i += 2) {
1110 ic = idp2 - i;
1111 ch_ref(i - 1, k, j) = cc_ref(i - 1, (j << 1) - 1, k) + cc_ref(ic - 1, (j << 1) - 2, k);
1112 ch_ref(i - 1, k, jc) = cc_ref(i - 1, (j << 1) - 1, k) - cc_ref(ic - 1, (j << 1) - 2, k);
1113 ch_ref(i, k, j) = cc_ref(i, (j << 1) - 1, k) - cc_ref(ic, (j << 1) - 2, k);
1114 ch_ref(i, k, jc) = cc_ref(i, (j << 1) - 1, k) + cc_ref(ic, (j << 1) - 2, k);
1115 }
1116 }
1117 }
1118 } else {
1119 for (j = 2; j <= ipph; ++j) {
1120 jc = ipp2 - j;
1121 for (i = 3; i <= ido; i += 2) {
1122 ic = idp2 - i;
1123 for (k = 1; k <= l1; ++k) {
1124 ch_ref(i - 1, k, j) = cc_ref(i - 1, (j << 1) - 1, k) + cc_ref(ic - 1, (j << 1) - 2, k);
1125 ch_ref(i - 1, k, jc) = cc_ref(i - 1, (j << 1) - 1, k) - cc_ref(ic - 1, (j << 1) - 2, k);
1126 ch_ref(i, k, j) = cc_ref(i, (j << 1) - 1, k) - cc_ref(ic, (j << 1) - 2, k);
1127 ch_ref(i, k, jc) = cc_ref(i, (j << 1) - 1, k) + cc_ref(ic, (j << 1) - 2, k);
1128 }
1129 }
1130 }
1131 }
1132 }
1133 ar1 = 1.f;
1134 ai1 = 0.f;
1135 for (l = 2; l <= ipph; ++l) {
1136 lc = ipp2 - l;
1137 ar1h = dcp * ar1 - dsp * ai1;
1138 ai1 = dcp * ai1 + dsp * ar1;
1139 ar1 = ar1h;
1140 for (ik = 1; ik <= idl1; ++ik) {
1141 c2_ref(ik, l) = ch2_ref(ik, 1) + ar1 * ch2_ref(ik, 2);
1142 c2_ref(ik, lc) = ai1 * ch2_ref(ik, ip);
1143 }
1144 dc2 = ar1;
1145 ds2 = ai1;
1146 ar2 = ar1;
1147 ai2 = ai1;
1148 for (j = 3; j <= ipph; ++j) {
1149 jc = ipp2 - j;
1150 ar2h = dc2 * ar2 - ds2 * ai2;
1151 ai2 = dc2 * ai2 + ds2 * ar2;
1152 ar2 = ar2h;
1153 for (ik = 1; ik <= idl1; ++ik) {
1154 c2_ref(ik, l) = c2_ref(ik, l) + ar2 * ch2_ref(ik, j);
1155 c2_ref(ik, lc) = c2_ref(ik, lc) + ai2 * ch2_ref(ik, jc);
1156 }
1157 }
1158 }
1159 for (j = 2; j <= ipph; ++j) {
1160 for (ik = 1; ik <= idl1; ++ik) {
1161 ch2_ref(ik, 1) = ch2_ref(ik, 1) + ch2_ref(ik, j);
1162 }
1163 }
1164 for (j = 2; j <= ipph; ++j) {
1165 jc = ipp2 - j;
1166 for (k = 1; k <= l1; ++k) {
1167 ch_ref(1, k, j) = c1_ref(1, k, j) - c1_ref(1, k, jc);
1168 ch_ref(1, k, jc) = c1_ref(1, k, j) + c1_ref(1, k, jc);
1169 }
1170 }
1171 if (ido != 1) {
1172 if (nbd >= l1) {
1173 for (j = 2; j <= ipph; ++j) {
1174 jc = ipp2 - j;
1175 for (k = 1; k <= l1; ++k) {
1176 for (i = 3; i <= ido; i += 2) {
1177 ch_ref(i - 1, k, j) = c1_ref(i - 1, k, j) - c1_ref(i, k, jc);
1178 ch_ref(i - 1, k, jc) = c1_ref(i - 1, k, j) + c1_ref(i, k, jc);
1179 ch_ref(i, k, j) = c1_ref(i, k, j) + c1_ref(i - 1, k, jc);
1180 ch_ref(i, k, jc) = c1_ref(i, k, j) - c1_ref(i - 1, k, jc);
1181 }
1182 }
1183 }
1184 } else {
1185 for (j = 2; j <= ipph; ++j) {
1186 jc = ipp2 - j;
1187 for (i = 3; i <= ido; i += 2) {
1188 for (k = 1; k <= l1; ++k) {
1189 ch_ref(i - 1, k, j) = c1_ref(i - 1, k, j) - c1_ref(i, k, jc);
1190 ch_ref(i - 1, k, jc) = c1_ref(i - 1, k, j) + c1_ref(i, k, jc);
1191 ch_ref(i, k, j) = c1_ref(i, k, j) + c1_ref(i - 1, k, jc);
1192 ch_ref(i, k, jc) = c1_ref(i, k, j) - c1_ref(i - 1, k, jc);
1193 }
1194 }
1195 }
1196 }
1197 }
1198 if (ido == 1) {
1199 return;
1200 }
1201 for (ik = 1; ik <= idl1; ++ik) {
1202 c2_ref(ik, 1) = ch2_ref(ik, 1);
1203 }
1204 for (j = 2; j <= ip; ++j) {
1205 for (k = 1; k <= l1; ++k) {
1206 c1_ref(1, k, j) = ch_ref(1, k, j);
1207 }
1208 }
1209 if (nbd <= l1) {
1210 is = -(ido);
1211 for (j = 2; j <= ip; ++j) {
1212 is += ido;
1213 idij = is;
1214 for (i = 3; i <= ido; i += 2) {
1215 idij += 2;
1216 for (k = 1; k <= l1; ++k) {
1217 c1_ref(i - 1, k, j) = wa[idij - 1] * ch_ref(i - 1, k, j)
1218 - wa[idij] * ch_ref(i, k, j);
1219 c1_ref(i, k, j) = wa[idij - 1] * ch_ref(i, k, j) + wa[idij] * ch_ref(i - 1, k, j);
1220 }
1221 }
1222 }
1223 } else {
1224 is = -(ido);
1225 for (j = 2; j <= ip; ++j) {
1226 is += ido;
1227 for (k = 1; k <= l1; ++k) {
1228 idij = is;
1229 for (i = 3; i <= ido; i += 2) {
1230 idij += 2;
1231 c1_ref(i - 1, k, j) = wa[idij - 1] * ch_ref(i - 1, k, j)
1232 - wa[idij] * ch_ref(i, k, j);
1233 c1_ref(i, k, j) = wa[idij - 1] * ch_ref(i, k, j) + wa[idij] * ch_ref(i - 1, k, j);
1234 }
1235 }
1236 }
1237 }
1238 } /* radbg */
1239
1240 #undef ch2_ref
1241 #undef ch_ref
1242 #undef cc_ref
1243 #undef c2_ref
1244 #undef c1_ref
1245
1246
radf2(integer ido,integer l1,const real * cc,real * ch,const real * wa1)1247 static void radf2(integer ido, integer l1, const real *cc, real *ch,
1248 const real *wa1)
1249 {
1250 /* System generated locals */
1251 integer ch_offset, cc_offset;
1252
1253 /* Local variables */
1254 integer i, k, ic;
1255 real ti2, tr2;
1256 integer idp2;
1257
1258
1259 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*l1 + (a_2))*ido + a_1]
1260 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*2 + (a_2))*ido + a_1]
1261
1262 /* Parameter adjustments */
1263 ch_offset = 1 + ido * 3;
1264 ch -= ch_offset;
1265 cc_offset = 1 + ido * (1 + l1);
1266 cc -= cc_offset;
1267 --wa1;
1268
1269 /* Function Body */
1270 for (k = 1; k <= l1; ++k) {
1271 ch_ref(1, 1, k) = cc_ref(1, k, 1) + cc_ref(1, k, 2);
1272 ch_ref(ido, 2, k) = cc_ref(1, k, 1) - cc_ref(1, k, 2);
1273 }
1274 if (ido < 2) return;
1275 if (ido != 2) {
1276 idp2 = ido + 2;
1277 for (k = 1; k <= l1; ++k) {
1278 for (i = 3; i <= ido; i += 2) {
1279 ic = idp2 - i;
1280 tr2 = wa1[i - 2] * cc_ref(i - 1, k, 2) + wa1[i - 1] *
1281 cc_ref(i, k, 2);
1282 ti2 = wa1[i - 2] * cc_ref(i, k, 2) - wa1[i - 1] * cc_ref(
1283 i - 1, k, 2);
1284 ch_ref(i, 1, k) = cc_ref(i, k, 1) + ti2;
1285 ch_ref(ic, 2, k) = ti2 - cc_ref(i, k, 1);
1286 ch_ref(i - 1, 1, k) = cc_ref(i - 1, k, 1) + tr2;
1287 ch_ref(ic - 1, 2, k) = cc_ref(i - 1, k, 1) - tr2;
1288 }
1289 }
1290 if (ido % 2 == 1) {
1291 return;
1292 }
1293 }
1294 for (k = 1; k <= l1; ++k) {
1295 ch_ref(1, 2, k) = -cc_ref(ido, k, 2);
1296 ch_ref(ido, 1, k) = cc_ref(ido, k, 1);
1297 }
1298 } /* radf2 */
1299
1300 #undef ch_ref
1301 #undef cc_ref
1302
1303
radf3(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2)1304 static void radf3(integer ido, integer l1, const real *cc, real *ch,
1305 const real *wa1, const real *wa2)
1306 {
1307 static const real taur = -.5f;
1308 static const real taui = .866025403784439f;
1309
1310 /* System generated locals */
1311 integer ch_offset, cc_offset;
1312
1313 /* Local variables */
1314 integer i, k, ic;
1315 real ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
1316 integer idp2;
1317
1318
1319 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*l1 + (a_2))*ido + a_1]
1320 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*3 + (a_2))*ido + a_1]
1321
1322 /* Parameter adjustments */
1323 ch_offset = 1 + (ido << 2);
1324 ch -= ch_offset;
1325 cc_offset = 1 + ido * (1 + l1);
1326 cc -= cc_offset;
1327 --wa1;
1328 --wa2;
1329
1330 /* Function Body */
1331 for (k = 1; k <= l1; ++k) {
1332 cr2 = cc_ref(1, k, 2) + cc_ref(1, k, 3);
1333 ch_ref(1, 1, k) = cc_ref(1, k, 1) + cr2;
1334 ch_ref(1, 3, k) = taui * (cc_ref(1, k, 3) - cc_ref(1, k, 2));
1335 ch_ref(ido, 2, k) = cc_ref(1, k, 1) + taur * cr2;
1336 }
1337 if (ido == 1) {
1338 return;
1339 }
1340 idp2 = ido + 2;
1341 for (k = 1; k <= l1; ++k) {
1342 for (i = 3; i <= ido; i += 2) {
1343 ic = idp2 - i;
1344 dr2 = wa1[i - 2] * cc_ref(i - 1, k, 2) + wa1[i - 1] *
1345 cc_ref(i, k, 2);
1346 di2 = wa1[i - 2] * cc_ref(i, k, 2) - wa1[i - 1] * cc_ref(
1347 i - 1, k, 2);
1348 dr3 = wa2[i - 2] * cc_ref(i - 1, k, 3) + wa2[i - 1] *
1349 cc_ref(i, k, 3);
1350 di3 = wa2[i - 2] * cc_ref(i, k, 3) - wa2[i - 1] * cc_ref(
1351 i - 1, k, 3);
1352 cr2 = dr2 + dr3;
1353 ci2 = di2 + di3;
1354 ch_ref(i - 1, 1, k) = cc_ref(i - 1, k, 1) + cr2;
1355 ch_ref(i, 1, k) = cc_ref(i, k, 1) + ci2;
1356 tr2 = cc_ref(i - 1, k, 1) + taur * cr2;
1357 ti2 = cc_ref(i, k, 1) + taur * ci2;
1358 tr3 = taui * (di2 - di3);
1359 ti3 = taui * (dr3 - dr2);
1360 ch_ref(i - 1, 3, k) = tr2 + tr3;
1361 ch_ref(ic - 1, 2, k) = tr2 - tr3;
1362 ch_ref(i, 3, k) = ti2 + ti3;
1363 ch_ref(ic, 2, k) = ti3 - ti2;
1364 }
1365 }
1366 } /* radf3 */
1367
1368 #undef ch_ref
1369 #undef cc_ref
1370
1371
radf4(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2,const real * wa3)1372 static void radf4(integer ido, integer l1, const real *cc, real *ch,
1373 const real *wa1, const real *wa2, const real *wa3)
1374 {
1375 /* Initialized data */
1376
1377 static const real hsqt2 = .7071067811865475f;
1378
1379 /* System generated locals */
1380 integer cc_offset, ch_offset;
1381
1382 /* Local variables */
1383 integer i, k, ic;
1384 real ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
1385 integer idp2;
1386
1387
1388 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*l1 + (a_2))*ido + a_1]
1389 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*4 + (a_2))*ido + a_1]
1390
1391 /* Parameter adjustments */
1392 ch_offset = 1 + ido * 5;
1393 ch -= ch_offset;
1394 cc_offset = 1 + ido * (1 + l1);
1395 cc -= cc_offset;
1396 --wa1;
1397 --wa2;
1398 --wa3;
1399
1400 /* Function Body */
1401 for (k = 1; k <= l1; ++k) {
1402 tr1 = cc_ref(1, k, 2) + cc_ref(1, k, 4);
1403 tr2 = cc_ref(1, k, 1) + cc_ref(1, k, 3);
1404 ch_ref(1, 1, k) = tr1 + tr2;
1405 ch_ref(ido, 4, k) = tr2 - tr1;
1406 ch_ref(ido, 2, k) = cc_ref(1, k, 1) - cc_ref(1, k, 3);
1407 ch_ref(1, 3, k) = cc_ref(1, k, 4) - cc_ref(1, k, 2);
1408 }
1409 if (ido < 2) return;
1410 if (ido != 2) {
1411 idp2 = ido + 2;
1412 for (k = 1; k <= l1; ++k) {
1413 for (i = 3; i <= ido; i += 2) {
1414 ic = idp2 - i;
1415 cr2 = wa1[i - 2] * cc_ref(i - 1, k, 2) + wa1[i - 1] *
1416 cc_ref(i, k, 2);
1417 ci2 = wa1[i - 2] * cc_ref(i, k, 2) - wa1[i - 1] * cc_ref(
1418 i - 1, k, 2);
1419 cr3 = wa2[i - 2] * cc_ref(i - 1, k, 3) + wa2[i - 1] *
1420 cc_ref(i, k, 3);
1421 ci3 = wa2[i - 2] * cc_ref(i, k, 3) - wa2[i - 1] * cc_ref(
1422 i - 1, k, 3);
1423 cr4 = wa3[i - 2] * cc_ref(i - 1, k, 4) + wa3[i - 1] *
1424 cc_ref(i, k, 4);
1425 ci4 = wa3[i - 2] * cc_ref(i, k, 4) - wa3[i - 1] * cc_ref(
1426 i - 1, k, 4);
1427 tr1 = cr2 + cr4;
1428 tr4 = cr4 - cr2;
1429 ti1 = ci2 + ci4;
1430 ti4 = ci2 - ci4;
1431 ti2 = cc_ref(i, k, 1) + ci3;
1432 ti3 = cc_ref(i, k, 1) - ci3;
1433 tr2 = cc_ref(i - 1, k, 1) + cr3;
1434 tr3 = cc_ref(i - 1, k, 1) - cr3;
1435 ch_ref(i - 1, 1, k) = tr1 + tr2;
1436 ch_ref(ic - 1, 4, k) = tr2 - tr1;
1437 ch_ref(i, 1, k) = ti1 + ti2;
1438 ch_ref(ic, 4, k) = ti1 - ti2;
1439 ch_ref(i - 1, 3, k) = ti4 + tr3;
1440 ch_ref(ic - 1, 2, k) = tr3 - ti4;
1441 ch_ref(i, 3, k) = tr4 + ti3;
1442 ch_ref(ic, 2, k) = tr4 - ti3;
1443 }
1444 }
1445 if (ido % 2 == 1) {
1446 return;
1447 }
1448 }
1449 for (k = 1; k <= l1; ++k) {
1450 ti1 = -hsqt2 * (cc_ref(ido, k, 2) + cc_ref(ido, k, 4));
1451 tr1 = hsqt2 * (cc_ref(ido, k, 2) - cc_ref(ido, k, 4));
1452 ch_ref(ido, 1, k) = tr1 + cc_ref(ido, k, 1);
1453 ch_ref(ido, 3, k) = cc_ref(ido, k, 1) - tr1;
1454 ch_ref(1, 2, k) = ti1 - cc_ref(ido, k, 3);
1455 ch_ref(1, 4, k) = ti1 + cc_ref(ido, k, 3);
1456 }
1457 } /* radf4 */
1458
1459 #undef ch_ref
1460 #undef cc_ref
1461
1462
radf5(integer ido,integer l1,const real * cc,real * ch,const real * wa1,const real * wa2,const real * wa3,const real * wa4)1463 static void radf5(integer ido, integer l1, const real *cc, real *ch,
1464 const real *wa1, const real *wa2, const real *wa3, const real *wa4)
1465 {
1466 /* Initialized data */
1467
1468 static const real tr11 = .309016994374947f;
1469 static const real ti11 = .951056516295154f;
1470 static const real tr12 = -.809016994374947f;
1471 static const real ti12 = .587785252292473f;
1472
1473 /* System generated locals */
1474 integer cc_offset, ch_offset;
1475
1476 /* Local variables */
1477 integer i, k, ic;
1478 real ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
1479 cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
1480 integer idp2;
1481
1482
1483 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*l1 + (a_2))*ido + a_1]
1484 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*5 + (a_2))*ido + a_1]
1485
1486 /* Parameter adjustments */
1487 ch_offset = 1 + ido * 6;
1488 ch -= ch_offset;
1489 cc_offset = 1 + ido * (1 + l1);
1490 cc -= cc_offset;
1491 --wa1;
1492 --wa2;
1493 --wa3;
1494 --wa4;
1495
1496 /* Function Body */
1497 for (k = 1; k <= l1; ++k) {
1498 cr2 = cc_ref(1, k, 5) + cc_ref(1, k, 2);
1499 ci5 = cc_ref(1, k, 5) - cc_ref(1, k, 2);
1500 cr3 = cc_ref(1, k, 4) + cc_ref(1, k, 3);
1501 ci4 = cc_ref(1, k, 4) - cc_ref(1, k, 3);
1502 ch_ref(1, 1, k) = cc_ref(1, k, 1) + cr2 + cr3;
1503 ch_ref(ido, 2, k) = cc_ref(1, k, 1) + tr11 * cr2 + tr12 * cr3;
1504 ch_ref(1, 3, k) = ti11 * ci5 + ti12 * ci4;
1505 ch_ref(ido, 4, k) = cc_ref(1, k, 1) + tr12 * cr2 + tr11 * cr3;
1506 ch_ref(1, 5, k) = ti12 * ci5 - ti11 * ci4;
1507 }
1508 if (ido == 1) {
1509 return;
1510 }
1511 idp2 = ido + 2;
1512 for (k = 1; k <= l1; ++k) {
1513 for (i = 3; i <= ido; i += 2) {
1514 ic = idp2 - i;
1515 dr2 = wa1[i - 2] * cc_ref(i - 1, k, 2) + wa1[i - 1] * cc_ref(i, k, 2);
1516 di2 = wa1[i - 2] * cc_ref(i, k, 2) - wa1[i - 1] * cc_ref(i - 1, k, 2);
1517 dr3 = wa2[i - 2] * cc_ref(i - 1, k, 3) + wa2[i - 1] * cc_ref(i, k, 3);
1518 di3 = wa2[i - 2] * cc_ref(i, k, 3) - wa2[i - 1] * cc_ref(i - 1, k, 3);
1519 dr4 = wa3[i - 2] * cc_ref(i - 1, k, 4) + wa3[i - 1] * cc_ref(i, k, 4);
1520 di4 = wa3[i - 2] * cc_ref(i, k, 4) - wa3[i - 1] * cc_ref(i - 1, k, 4);
1521 dr5 = wa4[i - 2] * cc_ref(i - 1, k, 5) + wa4[i - 1] * cc_ref(i, k, 5);
1522 di5 = wa4[i - 2] * cc_ref(i, k, 5) - wa4[i - 1] * cc_ref(i - 1, k, 5);
1523 cr2 = dr2 + dr5;
1524 ci5 = dr5 - dr2;
1525 cr5 = di2 - di5;
1526 ci2 = di2 + di5;
1527 cr3 = dr3 + dr4;
1528 ci4 = dr4 - dr3;
1529 cr4 = di3 - di4;
1530 ci3 = di3 + di4;
1531 ch_ref(i - 1, 1, k) = cc_ref(i - 1, k, 1) + cr2 + cr3;
1532 ch_ref(i, 1, k) = cc_ref(i, k, 1) + ci2 + ci3;
1533 tr2 = cc_ref(i - 1, k, 1) + tr11 * cr2 + tr12 * cr3;
1534 ti2 = cc_ref(i, k, 1) + tr11 * ci2 + tr12 * ci3;
1535 tr3 = cc_ref(i - 1, k, 1) + tr12 * cr2 + tr11 * cr3;
1536 ti3 = cc_ref(i, k, 1) + tr12 * ci2 + tr11 * ci3;
1537 tr5 = ti11 * cr5 + ti12 * cr4;
1538 ti5 = ti11 * ci5 + ti12 * ci4;
1539 tr4 = ti12 * cr5 - ti11 * cr4;
1540 ti4 = ti12 * ci5 - ti11 * ci4;
1541 ch_ref(i - 1, 3, k) = tr2 + tr5;
1542 ch_ref(ic - 1, 2, k) = tr2 - tr5;
1543 ch_ref(i, 3, k) = ti2 + ti5;
1544 ch_ref(ic, 2, k) = ti5 - ti2;
1545 ch_ref(i - 1, 5, k) = tr3 + tr4;
1546 ch_ref(ic - 1, 4, k) = tr3 - tr4;
1547 ch_ref(i, 5, k) = ti3 + ti4;
1548 ch_ref(ic, 4, k) = ti4 - ti3;
1549 }
1550 }
1551 } /* radf5 */
1552
1553 #undef ch_ref
1554 #undef cc_ref
1555
1556
radfg(integer ido,integer ip,integer l1,integer idl1,real * cc,real * c1,real * c2,real * ch,real * ch2,const real * wa)1557 static void radfg(integer ido, integer ip, integer l1, integer idl1,
1558 real *cc, real *c1, real *c2, real *ch, real *ch2, const real *wa)
1559 {
1560 /* System generated locals */
1561 integer ch_offset, cc_offset,
1562 c1_offset, c2_offset, ch2_offset;
1563
1564 /* Local variables */
1565 integer i, j, k, l, j2, ic, jc, lc, ik, is;
1566 real dc2, ai1, ai2, ar1, ar2, ds2;
1567 integer nbd;
1568 real dcp, arg, dsp, ar1h, ar2h;
1569 integer idp2, ipp2, idij, ipph;
1570
1571
1572 #define c1_ref(a_1,a_2,a_3) c1[((a_3)*l1 + (a_2))*ido + a_1]
1573 #define c2_ref(a_1,a_2) c2[(a_2)*idl1 + a_1]
1574 #define cc_ref(a_1,a_2,a_3) cc[((a_3)*ip + (a_2))*ido + a_1]
1575 #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1]
1576 #define ch2_ref(a_1,a_2) ch2[(a_2)*idl1 + a_1]
1577
1578 /* Parameter adjustments */
1579 ch_offset = 1 + ido * (1 + l1);
1580 ch -= ch_offset;
1581 c1_offset = 1 + ido * (1 + l1);
1582 c1 -= c1_offset;
1583 cc_offset = 1 + ido * (1 + ip);
1584 cc -= cc_offset;
1585 ch2_offset = 1 + idl1;
1586 ch2 -= ch2_offset;
1587 c2_offset = 1 + idl1;
1588 c2 -= c2_offset;
1589 --wa;
1590
1591 /* Function Body */
1592 arg = (2*M_PI) / (real) (ip);
1593 dcp = FFTPACK_COS(arg);
1594 dsp = FFTPACK_SIN(arg);
1595 ipph = (ip + 1) / 2;
1596 ipp2 = ip + 2;
1597 idp2 = ido + 2;
1598 nbd = (ido - 1) / 2;
1599 if (ido == 1) {
1600 for (ik = 1; ik <= idl1; ++ik) {
1601 c2_ref(ik, 1) = ch2_ref(ik, 1);
1602 }
1603 } else {
1604 for (ik = 1; ik <= idl1; ++ik) {
1605 ch2_ref(ik, 1) = c2_ref(ik, 1);
1606 }
1607 for (j = 2; j <= ip; ++j) {
1608 for (k = 1; k <= l1; ++k) {
1609 ch_ref(1, k, j) = c1_ref(1, k, j);
1610 }
1611 }
1612 if (nbd <= l1) {
1613 is = -(ido);
1614 for (j = 2; j <= ip; ++j) {
1615 is += ido;
1616 idij = is;
1617 for (i = 3; i <= ido; i += 2) {
1618 idij += 2;
1619 for (k = 1; k <= l1; ++k) {
1620 ch_ref(i - 1, k, j) = wa[idij - 1] * c1_ref(i - 1, k, j)
1621 + wa[idij] * c1_ref(i, k, j);
1622 ch_ref(i, k, j) = wa[idij - 1] * c1_ref(i, k, j) - wa[
1623 idij] * c1_ref(i - 1, k, j);
1624 }
1625 }
1626 }
1627 } else {
1628 is = -(ido);
1629 for (j = 2; j <= ip; ++j) {
1630 is += ido;
1631 for (k = 1; k <= l1; ++k) {
1632 idij = is;
1633 for (i = 3; i <= ido; i += 2) {
1634 idij += 2;
1635 ch_ref(i - 1, k, j) = wa[idij - 1] * c1_ref(i - 1, k, j)
1636 + wa[idij] * c1_ref(i, k, j);
1637 ch_ref(i, k, j) = wa[idij - 1] * c1_ref(i, k, j) - wa[
1638 idij] * c1_ref(i - 1, k, j);
1639 }
1640 }
1641 }
1642 }
1643 if (nbd >= l1) {
1644 for (j = 2; j <= ipph; ++j) {
1645 jc = ipp2 - j;
1646 for (k = 1; k <= l1; ++k) {
1647 for (i = 3; i <= ido; i += 2) {
1648 c1_ref(i - 1, k, j) = ch_ref(i - 1, k, j) + ch_ref(i -
1649 1, k, jc);
1650 c1_ref(i - 1, k, jc) = ch_ref(i, k, j) - ch_ref(i, k,
1651 jc);
1652 c1_ref(i, k, j) = ch_ref(i, k, j) + ch_ref(i, k, jc);
1653 c1_ref(i, k, jc) = ch_ref(i - 1, k, jc) - ch_ref(i - 1,
1654 k, j);
1655 }
1656 }
1657 }
1658 } else {
1659 for (j = 2; j <= ipph; ++j) {
1660 jc = ipp2 - j;
1661 for (i = 3; i <= ido; i += 2) {
1662 for (k = 1; k <= l1; ++k) {
1663 c1_ref(i - 1, k, j) = ch_ref(i - 1, k, j) + ch_ref(i -
1664 1, k, jc);
1665 c1_ref(i - 1, k, jc) = ch_ref(i, k, j) - ch_ref(i, k,
1666 jc);
1667 c1_ref(i, k, j) = ch_ref(i, k, j) + ch_ref(i, k, jc);
1668 c1_ref(i, k, jc) = ch_ref(i - 1, k, jc) - ch_ref(i - 1,
1669 k, j);
1670 }
1671 }
1672 }
1673 }
1674 }
1675 for (j = 2; j <= ipph; ++j) {
1676 jc = ipp2 - j;
1677 for (k = 1; k <= l1; ++k) {
1678 c1_ref(1, k, j) = ch_ref(1, k, j) + ch_ref(1, k, jc);
1679 c1_ref(1, k, jc) = ch_ref(1, k, jc) - ch_ref(1, k, j);
1680 }
1681 }
1682
1683 ar1 = 1.f;
1684 ai1 = 0.f;
1685 for (l = 2; l <= ipph; ++l) {
1686 lc = ipp2 - l;
1687 ar1h = dcp * ar1 - dsp * ai1;
1688 ai1 = dcp * ai1 + dsp * ar1;
1689 ar1 = ar1h;
1690 for (ik = 1; ik <= idl1; ++ik) {
1691 ch2_ref(ik, l) = c2_ref(ik, 1) + ar1 * c2_ref(ik, 2);
1692 ch2_ref(ik, lc) = ai1 * c2_ref(ik, ip);
1693 }
1694 dc2 = ar1;
1695 ds2 = ai1;
1696 ar2 = ar1;
1697 ai2 = ai1;
1698 for (j = 3; j <= ipph; ++j) {
1699 jc = ipp2 - j;
1700 ar2h = dc2 * ar2 - ds2 * ai2;
1701 ai2 = dc2 * ai2 + ds2 * ar2;
1702 ar2 = ar2h;
1703 for (ik = 1; ik <= idl1; ++ik) {
1704 ch2_ref(ik, l) = ch2_ref(ik, l) + ar2 * c2_ref(ik, j);
1705 ch2_ref(ik, lc) = ch2_ref(ik, lc) + ai2 * c2_ref(ik, jc);
1706 }
1707 }
1708 }
1709 for (j = 2; j <= ipph; ++j) {
1710 for (ik = 1; ik <= idl1; ++ik) {
1711 ch2_ref(ik, 1) = ch2_ref(ik, 1) + c2_ref(ik, j);
1712 }
1713 }
1714
1715 if (ido >= l1) {
1716 for (k = 1; k <= l1; ++k) {
1717 for (i = 1; i <= ido; ++i) {
1718 cc_ref(i, 1, k) = ch_ref(i, k, 1);
1719 }
1720 }
1721 } else {
1722 for (i = 1; i <= ido; ++i) {
1723 for (k = 1; k <= l1; ++k) {
1724 cc_ref(i, 1, k) = ch_ref(i, k, 1);
1725 }
1726 }
1727 }
1728 for (j = 2; j <= ipph; ++j) {
1729 jc = ipp2 - j;
1730 j2 = j + j;
1731 for (k = 1; k <= l1; ++k) {
1732 cc_ref(ido, j2 - 2, k) = ch_ref(1, k, j);
1733 cc_ref(1, j2 - 1, k) = ch_ref(1, k, jc);
1734 }
1735 }
1736 if (ido == 1) {
1737 return;
1738 }
1739 if (nbd >= l1) {
1740 for (j = 2; j <= ipph; ++j) {
1741 jc = ipp2 - j;
1742 j2 = j + j;
1743 for (k = 1; k <= l1; ++k) {
1744 for (i = 3; i <= ido; i += 2) {
1745 ic = idp2 - i;
1746 cc_ref(i - 1, j2 - 1, k) = ch_ref(i - 1, k, j) + ch_ref(
1747 i - 1, k, jc);
1748 cc_ref(ic - 1, j2 - 2, k) = ch_ref(i - 1, k, j) - ch_ref(
1749 i - 1, k, jc);
1750 cc_ref(i, j2 - 1, k) = ch_ref(i, k, j) + ch_ref(i, k,
1751 jc);
1752 cc_ref(ic, j2 - 2, k) = ch_ref(i, k, jc) - ch_ref(i, k, j)
1753 ;
1754 }
1755 }
1756 }
1757 } else {
1758 for (j = 2; j <= ipph; ++j) {
1759 jc = ipp2 - j;
1760 j2 = j + j;
1761 for (i = 3; i <= ido; i += 2) {
1762 ic = idp2 - i;
1763 for (k = 1; k <= l1; ++k) {
1764 cc_ref(i - 1, j2 - 1, k) = ch_ref(i - 1, k, j) + ch_ref(
1765 i - 1, k, jc);
1766 cc_ref(ic - 1, j2 - 2, k) = ch_ref(i - 1, k, j) - ch_ref(
1767 i - 1, k, jc);
1768 cc_ref(i, j2 - 1, k) = ch_ref(i, k, j) + ch_ref(i, k,
1769 jc);
1770 cc_ref(ic, j2 - 2, k) = ch_ref(i, k, jc) - ch_ref(i, k, j)
1771 ;
1772 }
1773 }
1774 }
1775 }
1776 } /* radfg */
1777
1778 #undef ch2_ref
1779 #undef ch_ref
1780 #undef cc_ref
1781 #undef c2_ref
1782 #undef c1_ref
1783
1784
cfftb1(integer n,real * c,real * ch,const real * wa,integer * ifac)1785 static void cfftb1(integer n, real *c, real *ch, const real *wa, integer *ifac)
1786 {
1787 integer i, k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido,
1788 idl1, idot;
1789
1790 /* Function Body */
1791 nf = ifac[1];
1792 na = 0;
1793 l1 = 1;
1794 iw = 0;
1795 for (k1 = 1; k1 <= nf; ++k1) {
1796 ip = ifac[k1 + 1];
1797 l2 = ip * l1;
1798 ido = n / l2;
1799 idot = ido + ido;
1800 idl1 = idot * l1;
1801 switch (ip) {
1802 case 4:
1803 ix2 = iw + idot;
1804 ix3 = ix2 + idot;
1805 passb4(idot, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2], &wa[ix3]);
1806 na = 1 - na;
1807 break;
1808 case 2:
1809 passb2(idot, l1, na?ch:c, na?c:ch, &wa[iw]);
1810 na = 1 - na;
1811 break;
1812 case 3:
1813 ix2 = iw + idot;
1814 passb3(idot, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2]);
1815 na = 1 - na;
1816 break;
1817 case 5:
1818 ix2 = iw + idot;
1819 ix3 = ix2 + idot;
1820 ix4 = ix3 + idot;
1821 passfb5(idot, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], +1);
1822 na = 1 - na;
1823 break;
1824 default:
1825 if (na == 0) {
1826 passfb(&nac, idot, ip, l1, idl1, c, c, c, ch, ch, &wa[iw], +1);
1827 } else {
1828 passfb(&nac, idot, ip, l1, idl1, ch, ch, ch, c, c, &wa[iw], +1);
1829 }
1830 if (nac != 0) {
1831 na = 1 - na;
1832 }
1833 break;
1834 }
1835 l1 = l2;
1836 iw += (ip - 1) * idot;
1837 }
1838 if (na == 0) {
1839 return;
1840 }
1841 for (i = 0; i < 2*n; ++i) {
1842 c[i] = ch[i];
1843 }
1844 } /* cfftb1 */
1845
cfftb(integer n,real * c,real * wsave)1846 void cfftb(integer n, real *c, real *wsave)
1847 {
1848 integer iw1, iw2;
1849
1850 /* Parameter adjustments */
1851 --wsave;
1852 --c;
1853
1854 /* Function Body */
1855 if (n == 1) {
1856 return;
1857 }
1858 iw1 = 2*n + 1;
1859 iw2 = iw1 + 2*n;
1860 cfftb1(n, &c[1], &wsave[1], &wsave[iw1], (int*)&wsave[iw2]);
1861 } /* cfftb */
1862
cfftf1(integer n,real * c,real * ch,const real * wa,integer * ifac)1863 static void cfftf1(integer n, real *c, real *ch, const real *wa, integer *ifac)
1864 {
1865 /* Local variables */
1866 integer i, k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido,
1867 idl1, idot;
1868
1869 /* Function Body */
1870 nf = ifac[1];
1871 na = 0;
1872 l1 = 1;
1873 iw = 0;
1874 for (k1 = 1; k1 <= nf; ++k1) {
1875 ip = ifac[k1 + 1];
1876 l2 = ip * l1;
1877 ido = n / l2;
1878 idot = ido + ido;
1879 idl1 = idot * l1;
1880 switch (ip) {
1881 case 4:
1882 ix2 = iw + idot;
1883 ix3 = ix2 + idot;
1884 passf4(idot, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2], &wa[ix3]);
1885 na = 1 - na;
1886 break;
1887 case 2:
1888 passf2(idot, l1, na?ch:c, na?c:ch, &wa[iw]);
1889 na = 1 - na;
1890 break;
1891 case 3:
1892 ix2 = iw + idot;
1893 passf3(idot, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2]);
1894 na = 1 - na;
1895 break;
1896 case 5:
1897 ix2 = iw + idot;
1898 ix3 = ix2 + idot;
1899 ix4 = ix3 + idot;
1900 passfb5(idot, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], -1);
1901 na = 1 - na;
1902 break;
1903 default:
1904 if (na == 0) {
1905 passfb(&nac, idot, ip, l1, idl1, c, c, c, ch, ch, &wa[iw], -1);
1906 } else {
1907 passfb(&nac, idot, ip, l1, idl1, ch, ch, ch, c, c, &wa[iw], -1);
1908 }
1909 if (nac != 0) {
1910 na = 1 - na;
1911 }
1912 break;
1913 }
1914 l1 = l2;
1915 iw += (ip - 1)*idot;
1916 }
1917 if (na == 0) {
1918 return;
1919 }
1920 for (i = 0; i < 2*n; ++i) {
1921 c[i] = ch[i];
1922 }
1923 } /* cfftf1 */
1924
cfftf(integer n,real * c,real * wsave)1925 void cfftf(integer n, real *c, real *wsave)
1926 {
1927 integer iw1, iw2;
1928
1929 /* Parameter adjustments */
1930 --wsave;
1931 --c;
1932
1933 /* Function Body */
1934 if (n == 1) {
1935 return;
1936 }
1937 iw1 = 2*n + 1;
1938 iw2 = iw1 + 2*n;
1939 cfftf1(n, &c[1], &wsave[1], &wsave[iw1], (int*)&wsave[iw2]);
1940 } /* cfftf */
1941
decompose(integer n,integer * ifac,integer ntryh[4])1942 static int decompose(integer n, integer *ifac, integer ntryh[4]) {
1943 integer ntry=0, nl = n, nf = 0, nq, nr, i, j = 0;
1944 do {
1945 if (j < 4) {
1946 ntry = ntryh[j];
1947 } else {
1948 ntry += 2;
1949 }
1950 ++j;
1951 L104:
1952 nq = nl / ntry;
1953 nr = nl - ntry * nq;
1954 if (nr != 0) continue;
1955 ++nf;
1956 ifac[nf + 2] = ntry;
1957 nl = nq;
1958 if (ntry == 2 && nf != 1) {
1959 for (i = 2; i <= nf; ++i) {
1960 integer ib = nf - i + 2;
1961 ifac[ib + 2] = ifac[ib + 1];
1962 }
1963 ifac[3] = 2;
1964 }
1965 if (nl != 1) {
1966 goto L104;
1967 }
1968 } while (nl != 1);
1969 ifac[1] = n;
1970 ifac[2] = nf;
1971 return nf;
1972 }
1973
cffti1(integer n,real * wa,integer * ifac)1974 static void cffti1(integer n, real *wa, integer *ifac)
1975 {
1976 static integer ntryh[4] = { 3,4,2,5 };
1977
1978 /* Local variables */
1979 integer i, j, i1, k1, l1, l2;
1980 real fi;
1981 integer ld, ii, nf, ip;
1982 real arg;
1983 integer ido, ipm;
1984 real argh;
1985 integer idot;
1986 real argld;
1987
1988 /* Parameter adjustments */
1989 --ifac;
1990 --wa;
1991
1992 nf = decompose(n, ifac, ntryh);
1993
1994 argh = (2*M_PI) / (real) (n);
1995 i = 2;
1996 l1 = 1;
1997 for (k1 = 1; k1 <= nf; ++k1) {
1998 ip = ifac[k1 + 2];
1999 ld = 0;
2000 l2 = l1 * ip;
2001 ido = n / l2;
2002 idot = ido + ido + 2;
2003 ipm = ip - 1;
2004 for (j = 1; j <= ipm; ++j) {
2005 i1 = i;
2006 wa[i - 1] = 1.f;
2007 wa[i] = 0.f;
2008 ld += l1;
2009 fi = 0.f;
2010 argld = (real) ld * argh;
2011 for (ii = 4; ii <= idot; ii += 2) {
2012 i += 2;
2013 fi += 1.f;
2014 arg = fi * argld;
2015 wa[i - 1] = FFTPACK_COS(arg);
2016 wa[i] = FFTPACK_SIN(arg);
2017 }
2018 if (ip > 5) {
2019 wa[i1 - 1] = wa[i - 1];
2020 wa[i1] = wa[i];
2021 };
2022 }
2023 l1 = l2;
2024 }
2025 } /* cffti1 */
2026
cffti(integer n,real * wsave)2027 void cffti(integer n, real *wsave)
2028 {
2029 integer iw1, iw2;
2030 /* Parameter adjustments */
2031 --wsave;
2032
2033 /* Function Body */
2034 if (n == 1) {
2035 return;
2036 }
2037 iw1 = 2*n + 1;
2038 iw2 = iw1 + 2*n;
2039 cffti1(n, &wsave[iw1], (int*)&wsave[iw2]);
2040 return;
2041 } /* cffti */
2042
rfftb1(integer n,real * c,real * ch,const real * wa,integer * ifac)2043 static void rfftb1(integer n, real *c, real *ch, const real *wa, integer *ifac)
2044 {
2045 /* Local variables */
2046 integer i, k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
2047
2048 /* Function Body */
2049 nf = ifac[1];
2050 na = 0;
2051 l1 = 1;
2052 iw = 0;
2053 for (k1 = 1; k1 <= nf; ++k1) {
2054 ip = ifac[k1 + 1];
2055 l2 = ip * l1;
2056 ido = n / l2;
2057 idl1 = ido * l1;
2058 switch (ip) {
2059 case 4:
2060 ix2 = iw + ido;
2061 ix3 = ix2 + ido;
2062 radb4(ido, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2], &wa[ix3]);
2063 na = 1 - na;
2064 break;
2065 case 2:
2066 radb2(ido, l1, na?ch:c, na?c:ch, &wa[iw]);
2067 na = 1 - na;
2068 break;
2069 case 3:
2070 ix2 = iw + ido;
2071 radb3(ido, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2]);
2072 na = 1 - na;
2073 break;
2074 case 5:
2075 ix2 = iw + ido;
2076 ix3 = ix2 + ido;
2077 ix4 = ix3 + ido;
2078 radb5(ido, l1, na?ch:c, na?c:ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
2079 na = 1 - na;
2080 break;
2081 default:
2082 if (na == 0) {
2083 radbg(ido, ip, l1, idl1, c, c, c, ch, ch, &wa[iw]);
2084 } else {
2085 radbg(ido, ip, l1, idl1, ch, ch, ch, c, c, &wa[iw]);
2086 }
2087 if (ido == 1) {
2088 na = 1 - na;
2089 }
2090 break;
2091 }
2092 l1 = l2;
2093 iw += (ip - 1) * ido;
2094 }
2095 if (na == 0) {
2096 return;
2097 }
2098 for (i = 0; i < n; ++i) {
2099 c[i] = ch[i];
2100 }
2101 } /* rfftb1 */
2102
rfftf1(integer n,real * c,real * ch,const real * wa,integer * ifac)2103 static void rfftf1(integer n, real *c, real *ch, const real *wa, integer *ifac)
2104 {
2105 /* Local variables */
2106 integer i, k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
2107
2108 /* Function Body */
2109 nf = ifac[1];
2110 na = 1;
2111 l2 = n;
2112 iw = n-1;
2113 for (k1 = 1; k1 <= nf; ++k1) {
2114 kh = nf - k1;
2115 ip = ifac[kh + 2];
2116 l1 = l2 / ip;
2117 ido = n / l2;
2118 idl1 = ido * l1;
2119 iw -= (ip - 1) * ido;
2120 na = 1 - na;
2121 switch (ip) {
2122 case 4:
2123 ix2 = iw + ido;
2124 ix3 = ix2 + ido;
2125 radf4(ido, l1, na ? ch : c, na ? c : ch, &wa[iw], &wa[ix2], &wa[ix3]);
2126 break;
2127 case 2:
2128 radf2(ido, l1, na ? ch : c, na ? c : ch, &wa[iw]);
2129 break;
2130 case 3:
2131 ix2 = iw + ido;
2132 radf3(ido, l1, na ? ch : c, na ? c : ch, &wa[iw], &wa[ix2]);
2133 break;
2134 case 5:
2135 ix2 = iw + ido;
2136 ix3 = ix2 + ido;
2137 ix4 = ix3 + ido;
2138 radf5(ido, l1, na ? ch : c, na ? c : ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
2139 break;
2140 default:
2141 if (ido == 1) {
2142 na = 1 - na;
2143 }
2144 if (na == 0) {
2145 radfg(ido, ip, l1, idl1, c, c, c, ch, ch, &wa[iw]);
2146 na = 1;
2147 } else {
2148 radfg(ido, ip, l1, idl1, ch, ch, ch, c, c, &wa[iw]);
2149 na = 0;
2150 }
2151 break;
2152 }
2153 l2 = l1;
2154 }
2155 if (na == 1) {
2156 return;
2157 }
2158 for (i = 0; i < n; ++i) {
2159 c[i] = ch[i];
2160 }
2161 }
2162
rfftb(integer n,real * r,real * wsave)2163 void rfftb(integer n, real *r, real *wsave)
2164 {
2165
2166 /* Parameter adjustments */
2167 --wsave;
2168 --r;
2169
2170 /* Function Body */
2171 if (n == 1) {
2172 return;
2173 }
2174 rfftb1(n, &r[1], &wsave[1], &wsave[n + 1], (int*)&wsave[(n << 1) + 1]);
2175 } /* rfftb */
2176
rffti1(integer n,real * wa,integer * ifac)2177 static void rffti1(integer n, real *wa, integer *ifac)
2178 {
2179 static integer ntryh[4] = { 4,2,3,5 };
2180
2181 /* Local variables */
2182 integer i, j, k1, l1, l2;
2183 real fi;
2184 integer ld, ii, nf, ip, is;
2185 real arg;
2186 integer ido, ipm;
2187 integer nfm1;
2188 real argh;
2189 real argld;
2190
2191 /* Parameter adjustments */
2192 --ifac;
2193 --wa;
2194
2195 nf = decompose(n, ifac, ntryh);
2196
2197 argh = (2*M_PI) / (real) (n);
2198 is = 0;
2199 nfm1 = nf - 1;
2200 l1 = 1;
2201 if (nfm1 == 0) {
2202 return;
2203 }
2204 for (k1 = 1; k1 <= nfm1; ++k1) {
2205 ip = ifac[k1 + 2];
2206 ld = 0;
2207 l2 = l1 * ip;
2208 ido = n / l2;
2209 ipm = ip - 1;
2210 for (j = 1; j <= ipm; ++j) {
2211 ld += l1;
2212 i = is;
2213 argld = (real) ld * argh;
2214 fi = 0.f;
2215 for (ii = 3; ii <= ido; ii += 2) {
2216 i += 2;
2217 fi += 1.f;
2218 arg = fi * argld;
2219 wa[i - 1] = FFTPACK_COS(arg);
2220 wa[i] = FFTPACK_SIN(arg);
2221 }
2222 is += ido;
2223 }
2224 l1 = l2;
2225 }
2226 } /* rffti1 */
2227
rfftf(integer n,real * r,real * wsave)2228 void rfftf(integer n, real *r, real *wsave)
2229 {
2230
2231 /* Parameter adjustments */
2232 --wsave;
2233 --r;
2234
2235 /* Function Body */
2236 if (n == 1) {
2237 return;
2238 }
2239 rfftf1(n, &r[1], &wsave[1], &wsave[n + 1], (int*)&wsave[(n << 1) + 1]);
2240 } /* rfftf */
2241
rffti(integer n,real * wsave)2242 void rffti(integer n, real *wsave)
2243 {
2244 /* Parameter adjustments */
2245 --wsave;
2246
2247 /* Function Body */
2248 if (n == 1) {
2249 return;
2250 }
2251 rffti1(n, &wsave[n + 1], (int*)&wsave[(n << 1) + 1]);
2252 return;
2253 } /* rffti */
2254
cosqb1(integer n,real * x,real * w,real * xh)2255 static void cosqb1(integer n, real *x, real *w, real *xh)
2256 {
2257 /* Local variables */
2258 integer i, k, kc, np2, ns2;
2259 real xim1;
2260 integer modn;
2261
2262 /* Parameter adjustments */
2263 --xh;
2264 --w;
2265 --x;
2266
2267 /* Function Body */
2268 ns2 = (n + 1) / 2;
2269 np2 = n + 2;
2270 for (i = 3; i <= n; i += 2) {
2271 xim1 = x[i - 1] + x[i];
2272 x[i] -= x[i - 1];
2273 x[i - 1] = xim1;
2274 }
2275 x[1] += x[1];
2276 modn = n % 2;
2277 if (modn == 0) {
2278 x[n] += x[n];
2279 }
2280 rfftb(n, &x[1], &xh[1]);
2281 for (k = 2; k <= ns2; ++k) {
2282 kc = np2 - k;
2283 xh[k] = w[k - 1] * x[kc] + w[kc - 1] * x[k];
2284 xh[kc] = w[k - 1] * x[k] - w[kc - 1] * x[kc];
2285 }
2286 if (modn == 0) {
2287 x[ns2 + 1] = w[ns2] * (x[ns2 + 1] + x[ns2 + 1]);
2288 }
2289 for (k = 2; k <= ns2; ++k) {
2290 kc = np2 - k;
2291 x[k] = xh[k] + xh[kc];
2292 x[kc] = xh[k] - xh[kc];
2293 }
2294 x[1] += x[1];
2295 } /* cosqb1 */
2296
cosqb(integer n,real * x,real * wsave)2297 void cosqb(integer n, real *x, real *wsave)
2298 {
2299 static const real tsqrt2 = 2.82842712474619f;
2300
2301 /* Local variables */
2302 real x1;
2303
2304 /* Parameter adjustments */
2305 --wsave;
2306 --x;
2307
2308 if (n < 2) {
2309 x[1] *= 4.f;
2310 } else if (n == 2) {
2311 x1 = (x[1] + x[2]) * 4.f;
2312 x[2] = tsqrt2 * (x[1] - x[2]);
2313 x[1] = x1;
2314 } else {
2315 cosqb1(n, &x[1], &wsave[1], &wsave[n + 1]);
2316 }
2317 } /* cosqb */
2318
cosqf1(integer n,real * x,real * w,real * xh)2319 static void cosqf1(integer n, real *x, real *w, real *xh)
2320 {
2321 /* Local variables */
2322 integer i, k, kc, np2, ns2;
2323 real xim1;
2324 integer modn;
2325
2326 /* Parameter adjustments */
2327 --xh;
2328 --w;
2329 --x;
2330
2331 /* Function Body */
2332 ns2 = (n + 1) / 2;
2333 np2 = n + 2;
2334 for (k = 2; k <= ns2; ++k) {
2335 kc = np2 - k;
2336 xh[k] = x[k] + x[kc];
2337 xh[kc] = x[k] - x[kc];
2338 }
2339 modn = n % 2;
2340 if (modn == 0) {
2341 xh[ns2 + 1] = x[ns2 + 1] + x[ns2 + 1];
2342 }
2343 for (k = 2; k <= ns2; ++k) {
2344 kc = np2 - k;
2345 x[k] = w[k - 1] * xh[kc] + w[kc - 1] * xh[k];
2346 x[kc] = w[k - 1] * xh[k] - w[kc - 1] * xh[kc];
2347 }
2348 if (modn == 0) {
2349 x[ns2 + 1] = w[ns2] * xh[ns2 + 1];
2350 }
2351 rfftf(n, &x[1], &xh[1]);
2352 for (i = 3; i <= n; i += 2) {
2353 xim1 = x[i - 1] - x[i];
2354 x[i] = x[i - 1] + x[i];
2355 x[i - 1] = xim1;
2356 }
2357 } /* cosqf1 */
2358
cosqf(integer n,real * x,real * wsave)2359 void cosqf(integer n, real *x, real *wsave)
2360 {
2361 static const real sqrt2 = 1.4142135623731f;
2362
2363 /* Local variables */
2364 real tsqx;
2365
2366 /* Parameter adjustments */
2367 --wsave;
2368 --x;
2369
2370 if (n == 2) {
2371 tsqx = sqrt2 * x[2];
2372 x[2] = x[1] - tsqx;
2373 x[1] += tsqx;
2374 } else if (n > 2) {
2375 cosqf1(n, &x[1], &wsave[1], &wsave[n + 1]);
2376 }
2377 } /* cosqf */
2378
cosqi(integer n,real * wsave)2379 void cosqi(integer n, real *wsave)
2380 {
2381 /* Local variables */
2382 integer k;
2383 real fk, dt;
2384
2385 /* Parameter adjustments */
2386 --wsave;
2387
2388 dt = M_PI/2 / (real) (n);
2389 fk = 0.f;
2390 for (k = 1; k <= n; ++k) {
2391 fk += 1.f;
2392 wsave[k] = FFTPACK_COS(fk * dt);
2393 }
2394 rffti(n, &wsave[n + 1]);
2395 } /* cosqi */
2396
cost(integer n,real * x,real * wsave)2397 void cost(integer n, real *x, real *wsave)
2398 {
2399 /* Local variables */
2400 integer i, k;
2401 real c1, t1, t2;
2402 integer kc;
2403 real xi;
2404 integer nm1, np1;
2405 real x1h;
2406 integer ns2;
2407 real tx2, x1p3, xim2;
2408 integer modn;
2409
2410 /* Parameter adjustments */
2411 --wsave;
2412 --x;
2413
2414 /* Function Body */
2415 nm1 = n - 1;
2416 np1 = n + 1;
2417 ns2 = n / 2;
2418 if (n == 2) {
2419 x1h = x[1] + x[2];
2420 x[2] = x[1] - x[2];
2421 x[1] = x1h;
2422 } else if (n == 3) {
2423 x1p3 = x[1] + x[3];
2424 tx2 = x[2] + x[2];
2425 x[2] = x[1] - x[3];
2426 x[1] = x1p3 + tx2;
2427 x[3] = x1p3 - tx2;
2428 } else if (n > 3) {
2429 c1 = x[1] - x[n];
2430 x[1] += x[n];
2431 for (k = 2; k <= ns2; ++k) {
2432 kc = np1 - k;
2433 t1 = x[k] + x[kc];
2434 t2 = x[k] - x[kc];
2435 c1 += wsave[kc] * t2;
2436 t2 = wsave[k] * t2;
2437 x[k] = t1 - t2;
2438 x[kc] = t1 + t2;
2439 }
2440 modn = n % 2;
2441 if (modn != 0) {
2442 x[ns2 + 1] += x[ns2 + 1];
2443 }
2444 rfftf(nm1, &x[1], &wsave[n + 1]);
2445 xim2 = x[2];
2446 x[2] = c1;
2447 for (i = 4; i <= n; i += 2) {
2448 xi = x[i];
2449 x[i] = x[i - 2] - x[i - 1];
2450 x[i - 1] = xim2;
2451 xim2 = xi;
2452 }
2453 if (modn != 0) {
2454 x[n] = xim2;
2455 }
2456 }
2457 } /* cost */
2458
costi(integer n,real * wsave)2459 void costi(integer n, real *wsave)
2460 {
2461 /* Initialized data */
2462
2463 /* Local variables */
2464 integer k, kc;
2465 real fk, dt;
2466 integer nm1, np1, ns2;
2467
2468 /* Parameter adjustments */
2469 --wsave;
2470
2471 /* Function Body */
2472 if (n <= 3) {
2473 return;
2474 }
2475 nm1 = n - 1;
2476 np1 = n + 1;
2477 ns2 = n / 2;
2478 dt = M_PI / (real) nm1;
2479 fk = 0.f;
2480 for (k = 2; k <= ns2; ++k) {
2481 kc = np1 - k;
2482 fk += 1.f;
2483 wsave[k] = FFTPACK_SIN(fk * dt) * 2.f;
2484 wsave[kc] = FFTPACK_COS(fk * dt) * 2.f;
2485 }
2486 rffti(nm1, &wsave[n + 1]);
2487 } /* costi */
2488
sinqb(integer n,real * x,real * wsave)2489 void sinqb(integer n, real *x, real *wsave)
2490 {
2491 /* Local variables */
2492 integer k, kc, ns2;
2493 real xhold;
2494
2495 /* Parameter adjustments */
2496 --wsave;
2497 --x;
2498
2499 /* Function Body */
2500 if (n <= 1) {
2501 x[1] *= 4.f;
2502 return;
2503 }
2504 ns2 = n / 2;
2505 for (k = 2; k <= n; k += 2) {
2506 x[k] = -x[k];
2507 }
2508 cosqb(n, &x[1], &wsave[1]);
2509 for (k = 1; k <= ns2; ++k) {
2510 kc = n - k;
2511 xhold = x[k];
2512 x[k] = x[kc + 1];
2513 x[kc + 1] = xhold;
2514 }
2515 } /* sinqb */
2516
sinqf(integer n,real * x,real * wsave)2517 void sinqf(integer n, real *x, real *wsave)
2518 {
2519 /* Local variables */
2520 integer k, kc, ns2;
2521 real xhold;
2522
2523 /* Parameter adjustments */
2524 --wsave;
2525 --x;
2526
2527 /* Function Body */
2528 if (n == 1) {
2529 return;
2530 }
2531 ns2 = n / 2;
2532 for (k = 1; k <= ns2; ++k) {
2533 kc = n - k;
2534 xhold = x[k];
2535 x[k] = x[kc + 1];
2536 x[kc + 1] = xhold;
2537 }
2538 cosqf(n, &x[1], &wsave[1]);
2539 for (k = 2; k <= n; k += 2) {
2540 x[k] = -x[k];
2541 }
2542 } /* sinqf */
2543
sinqi(integer n,real * wsave)2544 void sinqi(integer n, real *wsave)
2545 {
2546
2547 /* Parameter adjustments */
2548 --wsave;
2549
2550 /* Function Body */
2551 cosqi(n, &wsave[1]);
2552 } /* sinqi */
2553
sint1(integer n,real * war,real * was,real * xh,real * x,integer * ifac)2554 static void sint1(integer n, real *war, real *was, real *xh, real *
2555 x, integer *ifac)
2556 {
2557 /* Initialized data */
2558
2559 static const real sqrt3 = 1.73205080756888f;
2560
2561 /* Local variables */
2562 integer i, k;
2563 real t1, t2;
2564 integer kc, np1, ns2, modn;
2565 real xhold;
2566
2567 /* Parameter adjustments */
2568 --ifac;
2569 --x;
2570 --xh;
2571 --was;
2572 --war;
2573
2574 /* Function Body */
2575 for (i = 1; i <= n; ++i) {
2576 xh[i] = war[i];
2577 war[i] = x[i];
2578 }
2579
2580 if (n < 2) {
2581 xh[1] += xh[1];
2582 } else if (n == 2) {
2583 xhold = sqrt3 * (xh[1] + xh[2]);
2584 xh[2] = sqrt3 * (xh[1] - xh[2]);
2585 xh[1] = xhold;
2586 } else {
2587 np1 = n + 1;
2588 ns2 = n / 2;
2589 x[1] = 0.f;
2590 for (k = 1; k <= ns2; ++k) {
2591 kc = np1 - k;
2592 t1 = xh[k] - xh[kc];
2593 t2 = was[k] * (xh[k] + xh[kc]);
2594 x[k + 1] = t1 + t2;
2595 x[kc + 1] = t2 - t1;
2596 }
2597 modn = n % 2;
2598 if (modn != 0) {
2599 x[ns2 + 2] = xh[ns2 + 1] * 4.f;
2600 }
2601 rfftf1(np1, &x[1], &xh[1], &war[1], &ifac[1]);
2602 xh[1] = x[1] * .5f;
2603 for (i = 3; i <= n; i += 2) {
2604 xh[i - 1] = -x[i];
2605 xh[i] = xh[i - 2] + x[i - 1];
2606 }
2607 if (modn == 0) {
2608 xh[n] = -x[n + 1];
2609 }
2610 }
2611 for (i = 1; i <= n; ++i) {
2612 x[i] = war[i];
2613 war[i] = xh[i];
2614 }
2615 } /* sint1 */
2616
sinti(integer n,real * wsave)2617 void sinti(integer n, real *wsave)
2618 {
2619 /* Local variables */
2620 integer k;
2621 real dt;
2622 integer np1, ns2;
2623
2624 /* Parameter adjustments */
2625 --wsave;
2626
2627 /* Function Body */
2628 if (n <= 1) {
2629 return;
2630 }
2631 ns2 = n / 2;
2632 np1 = n + 1;
2633 dt = M_PI / (real) np1;
2634 for (k = 1; k <= ns2; ++k) {
2635 wsave[k] = sin(k * dt) * 2.f;
2636 }
2637 rffti(np1, &wsave[ns2 + 1]);
2638 } /* sinti */
2639
sint(integer n,real * x,real * wsave)2640 void sint(integer n, real *x, real *wsave)
2641 {
2642 integer np1, iw1, iw2, iw3;
2643
2644 /* Parameter adjustments */
2645 --wsave;
2646 --x;
2647
2648 /* Function Body */
2649 np1 = n + 1;
2650 iw1 = n / 2 + 1;
2651 iw2 = iw1 + np1;
2652 iw3 = iw2 + np1;
2653 sint1(n, &x[1], &wsave[1], &wsave[iw1], &wsave[iw2], (int*)&wsave[iw3]);
2654 } /* sint */
2655
2656 #ifdef TESTING_FFTPACK
2657 #include <stdio.h>
2658
main(void)2659 int main(void)
2660 {
2661 static integer nd[] = { 120,91,54,49,32,28,24,8,4,3,2 };
2662
2663 /* System generated locals */
2664 real r1, r2, r3;
2665 f77complex q1, q2, q3;
2666
2667 /* Local variables */
2668 integer i, j, k, n;
2669 real w[2000], x[200], y[200], cf, fn, dt;
2670 f77complex cx[200], cy[200];
2671 real xh[200];
2672 integer nz, nm1, np1, ns2;
2673 real arg, tfn;
2674 real sum, arg1, arg2;
2675 real sum1, sum2, dcfb;
2676 integer modn;
2677 real rftb, rftf;
2678 real sqrt2;
2679 real rftfb;
2680 real costt, sintt, dcfftb, dcfftf, cosqfb, costfb;
2681 real sinqfb;
2682 real sintfb;
2683 real cosqbt, cosqft, sinqbt, sinqft;
2684
2685
2686
2687 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2688
2689 /* VERSION 4 APRIL 1985 */
2690
2691 /* A TEST DRIVER FOR */
2692 /* A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE FAST FOURIER */
2693 /* TRANSFORM OF PERIODIC AND OTHER SYMMETRIC SEQUENCES */
2694
2695 /* BY */
2696
2697 /* PAUL N SWARZTRAUBER */
2698
2699 /* NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 */
2700
2701 /* WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION */
2702
2703 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2704
2705
2706 /* THIS PROGRAM TESTS THE PACKAGE OF FAST FOURIER */
2707 /* TRANSFORMS FOR BOTH COMPLEX AND REAL PERIODIC SEQUENCES AND */
2708 /* CERTIAN OTHER SYMMETRIC SEQUENCES THAT ARE LISTED BELOW. */
2709
2710 /* 1. RFFTI INITIALIZE RFFTF AND RFFTB */
2711 /* 2. RFFTF FORWARD TRANSFORM OF A REAL PERIODIC SEQUENCE */
2712 /* 3. RFFTB BACKWARD TRANSFORM OF A REAL COEFFICIENT ARRAY */
2713
2714 /* 4. EZFFTI INITIALIZE EZFFTF AND EZFFTB */
2715 /* 5. EZFFTF A SIMPLIFIED REAL PERIODIC FORWARD TRANSFORM */
2716 /* 6. EZFFTB A SIMPLIFIED REAL PERIODIC BACKWARD TRANSFORM */
2717
2718 /* 7. SINTI INITIALIZE SINT */
2719 /* 8. SINT SINE TRANSFORM OF A REAL ODD SEQUENCE */
2720
2721 /* 9. COSTI INITIALIZE COST */
2722 /* 10. COST COSINE TRANSFORM OF A REAL EVEN SEQUENCE */
2723
2724 /* 11. SINQI INITIALIZE SINQF AND SINQB */
2725 /* 12. SINQF FORWARD SINE TRANSFORM WITH ODD WAVE NUMBERS */
2726 /* 13. SINQB UNNORMALIZED INVERSE OF SINQF */
2727
2728 /* 14. COSQI INITIALIZE COSQF AND COSQB */
2729 /* 15. COSQF FORWARD COSINE TRANSFORM WITH ODD WAVE NUMBERS */
2730 /* 16. COSQB UNNORMALIZED INVERSE OF COSQF */
2731
2732 /* 17. CFFTI INITIALIZE CFFTF AND CFFTB */
2733 /* 18. CFFTF FORWARD TRANSFORM OF A COMPLEX PERIODIC SEQUENCE */
2734 /* 19. CFFTB UNNORMALIZED INVERSE OF CFFTF */
2735
2736
2737 sqrt2 = sqrt(2.f);
2738 int all_ok = 1;
2739 for (nz = 1; nz <= (int)(sizeof nd/sizeof nd[0]); ++nz) {
2740 n = nd[nz - 1];
2741 modn = n % 2;
2742 fn = (real) n;
2743 tfn = fn + fn;
2744 np1 = n + 1;
2745 nm1 = n - 1;
2746 for (j = 1; j <= np1; ++j) {
2747 x[j - 1] = sin((real) j * sqrt2);
2748 y[j - 1] = x[j - 1];
2749 xh[j - 1] = x[j - 1];
2750 }
2751
2752 /* TEST SUBROUTINES RFFTI,RFFTF AND RFFTB */
2753
2754 rffti(n, w);
2755 dt = (2*M_PI) / fn;
2756 ns2 = (n + 1) / 2;
2757 if (ns2 < 2) {
2758 goto L104;
2759 }
2760 for (k = 2; k <= ns2; ++k) {
2761 sum1 = 0.f;
2762 sum2 = 0.f;
2763 arg = (real) (k - 1) * dt;
2764 for (i = 1; i <= n; ++i) {
2765 arg1 = (real) (i - 1) * arg;
2766 sum1 += x[i - 1] * cos(arg1);
2767 sum2 += x[i - 1] * sin(arg1);
2768 }
2769 y[(k << 1) - 3] = sum1;
2770 y[(k << 1) - 2] = -sum2;
2771 }
2772 L104:
2773 sum1 = 0.f;
2774 sum2 = 0.f;
2775 for (i = 1; i <= nm1; i += 2) {
2776 sum1 += x[i - 1];
2777 sum2 += x[i];
2778 }
2779 if (modn == 1) {
2780 sum1 += x[n - 1];
2781 }
2782 y[0] = sum1 + sum2;
2783 if (modn == 0) {
2784 y[n - 1] = sum1 - sum2;
2785 }
2786 rfftf(n, x, w);
2787 rftf = 0.f;
2788 for (i = 1; i <= n; ++i) {
2789 /* Computing MAX */
2790 r2 = rftf, r3 = (r1 = x[i - 1] - y[i - 1], fabs(r1));
2791 rftf = dmax(r2,r3);
2792 x[i - 1] = xh[i - 1];
2793 }
2794 rftf /= fn;
2795 for (i = 1; i <= n; ++i) {
2796 sum = x[0] * .5f;
2797 arg = (real) (i - 1) * dt;
2798 if (ns2 < 2) {
2799 goto L108;
2800 }
2801 for (k = 2; k <= ns2; ++k) {
2802 arg1 = (real) (k - 1) * arg;
2803 sum = sum + x[(k << 1) - 3] * cos(arg1) - x[(k << 1) - 2] *
2804 sin(arg1);
2805 }
2806 L108:
2807 if (modn == 0) {
2808 sum += (real)pow(-1, i-1) * .5f * x[n - 1];
2809 }
2810 y[i - 1] = sum + sum;
2811 }
2812 rfftb(n, x, w);
2813 rftb = 0.f;
2814 for (i = 1; i <= n; ++i) {
2815 /* Computing MAX */
2816 r2 = rftb, r3 = (r1 = x[i - 1] - y[i - 1], fabs(r1));
2817 rftb = dmax(r2,r3);
2818 x[i - 1] = xh[i - 1];
2819 y[i - 1] = xh[i - 1];
2820 }
2821 rfftb(n, y, w);
2822 rfftf(n, y, w);
2823 cf = 1.f / fn;
2824 rftfb = 0.f;
2825 for (i = 1; i <= n; ++i) {
2826 /* Computing MAX */
2827 r2 = rftfb, r3 = (r1 = cf * y[i - 1] - x[i - 1], fabs(
2828 r1));
2829 rftfb = dmax(r2,r3);
2830 }
2831
2832 /* TEST SUBROUTINES SINTI AND SINT */
2833
2834 dt = M_PI / fn;
2835 for (i = 1; i <= nm1; ++i) {
2836 x[i - 1] = xh[i - 1];
2837 }
2838 for (i = 1; i <= nm1; ++i) {
2839 y[i - 1] = 0.f;
2840 arg1 = (real) i * dt;
2841 for (k = 1; k <= nm1; ++k) {
2842 y[i - 1] += x[k - 1] * sin((real) k * arg1);
2843 }
2844 y[i - 1] += y[i - 1];
2845 }
2846 sinti(nm1, w);
2847 sint(nm1, x, w);
2848 cf = .5f / fn;
2849 sintt = 0.f;
2850 for (i = 1; i <= nm1; ++i) {
2851 /* Computing MAX */
2852 r2 = sintt, r3 = (r1 = x[i - 1] - y[i - 1], fabs(r1));
2853 sintt = dmax(r2,r3);
2854 x[i - 1] = xh[i - 1];
2855 y[i - 1] = x[i - 1];
2856 }
2857 sintt = cf * sintt;
2858 sint(nm1, x, w);
2859 sint(nm1, x, w);
2860 sintfb = 0.f;
2861 for (i = 1; i <= nm1; ++i) {
2862 /* Computing MAX */
2863 r2 = sintfb, r3 = (r1 = cf * x[i - 1] - y[i - 1], fabs(
2864 r1));
2865 sintfb = dmax(r2,r3);
2866 }
2867
2868 /* TEST SUBROUTINES COSTI AND COST */
2869
2870 for (i = 1; i <= np1; ++i) {
2871 x[i - 1] = xh[i - 1];
2872 }
2873 for (i = 1; i <= np1; ++i) {
2874 y[i - 1] = (x[0] + (real) pow(-1, i+1) * x[n]) * .5f;
2875 arg = (real) (i - 1) * dt;
2876 for (k = 2; k <= n; ++k) {
2877 y[i - 1] += x[k - 1] * FFTPACK_COS((real) (k - 1) * arg);
2878 }
2879 y[i - 1] += y[i - 1];
2880 }
2881 costi(np1, w);
2882 cost(np1, x, w);
2883 costt = 0.f;
2884 for (i = 1; i <= np1; ++i) {
2885 /* Computing MAX */
2886 r2 = costt, r3 = (r1 = x[i - 1] - y[i - 1], fabs(r1));
2887 costt = dmax(r2,r3);
2888 x[i - 1] = xh[i - 1];
2889 y[i - 1] = xh[i - 1];
2890 }
2891 costt = cf * costt;
2892 cost(np1, x, w);
2893 cost(np1, x, w);
2894 costfb = 0.f;
2895 for (i = 1; i <= np1; ++i) {
2896 /* Computing MAX */
2897 r2 = costfb, r3 = (r1 = cf * x[i - 1] - y[i - 1], fabs(
2898 r1));
2899 costfb = dmax(r2,r3);
2900 }
2901
2902 /* TEST SUBROUTINES SINQI,SINQF AND SINQB */
2903
2904 cf = .25f / fn;
2905 for (i = 1; i <= n; ++i) {
2906 y[i - 1] = xh[i - 1];
2907 }
2908 dt = M_PI / (fn + fn);
2909 for (i = 1; i <= n; ++i) {
2910 x[i - 1] = 0.f;
2911 arg = dt * (real) i;
2912 for (k = 1; k <= n; ++k) {
2913 x[i - 1] += y[k - 1] * sin((real) (k + k - 1) * arg);
2914 }
2915 x[i - 1] *= 4.f;
2916 }
2917 sinqi(n, w);
2918 sinqb(n, y, w);
2919 sinqbt = 0.f;
2920 for (i = 1; i <= n; ++i) {
2921 /* Computing MAX */
2922 r2 = sinqbt, r3 = (r1 = y[i - 1] - x[i - 1], fabs(r1))
2923 ;
2924 sinqbt = dmax(r2,r3);
2925 x[i - 1] = xh[i - 1];
2926 }
2927 sinqbt = cf * sinqbt;
2928 for (i = 1; i <= n; ++i) {
2929 arg = (real) (i + i - 1) * dt;
2930 y[i - 1] = (real) pow(-1, i+1) * .5f * x[n - 1];
2931 for (k = 1; k <= nm1; ++k) {
2932 y[i - 1] += x[k - 1] * sin((real) k * arg);
2933 }
2934 y[i - 1] += y[i - 1];
2935 }
2936 sinqf(n, x, w);
2937 sinqft = 0.f;
2938 for (i = 1; i <= n; ++i) {
2939 /* Computing MAX */
2940 r2 = sinqft, r3 = (r1 = x[i - 1] - y[i - 1], fabs(r1))
2941 ;
2942 sinqft = dmax(r2,r3);
2943 y[i - 1] = xh[i - 1];
2944 x[i - 1] = xh[i - 1];
2945 }
2946 sinqf(n, y, w);
2947 sinqb(n, y, w);
2948 sinqfb = 0.f;
2949 for (i = 1; i <= n; ++i) {
2950 /* Computing MAX */
2951 r2 = sinqfb, r3 = (r1 = cf * y[i - 1] - x[i - 1], fabs(
2952 r1));
2953 sinqfb = dmax(r2,r3);
2954 }
2955
2956 /* TEST SUBROUTINES COSQI,COSQF AND COSQB */
2957
2958 for (i = 1; i <= n; ++i) {
2959 y[i - 1] = xh[i - 1];
2960 }
2961 for (i = 1; i <= n; ++i) {
2962 x[i - 1] = 0.f;
2963 arg = (real) (i - 1) * dt;
2964 for (k = 1; k <= n; ++k) {
2965 x[i - 1] += y[k - 1] * FFTPACK_COS((real) (k + k - 1) * arg);
2966 }
2967 x[i - 1] *= 4.f;
2968 }
2969 cosqi(n, w);
2970 cosqb(n, y, w);
2971 cosqbt = 0.f;
2972 for (i = 1; i <= n; ++i) {
2973 /* Computing MAX */
2974 r2 = cosqbt, r3 = (r1 = x[i - 1] - y[i - 1], fabs(r1))
2975 ;
2976 cosqbt = dmax(r2,r3);
2977 x[i - 1] = xh[i - 1];
2978 }
2979 cosqbt = cf * cosqbt;
2980 for (i = 1; i <= n; ++i) {
2981 y[i - 1] = x[0] * .5f;
2982 arg = (real) (i + i - 1) * dt;
2983 for (k = 2; k <= n; ++k) {
2984 y[i - 1] += x[k - 1] * FFTPACK_COS((real) (k - 1) * arg);
2985 }
2986 y[i - 1] += y[i - 1];
2987 }
2988 cosqf(n, x, w);
2989 cosqft = 0.f;
2990 for (i = 1; i <= n; ++i) {
2991 /* Computing MAX */
2992 r2 = cosqft, r3 = (r1 = y[i - 1] - x[i - 1], fabs(r1))
2993 ;
2994 cosqft = dmax(r2,r3);
2995 x[i - 1] = xh[i - 1];
2996 y[i - 1] = xh[i - 1];
2997 }
2998 cosqft = cf * cosqft;
2999 cosqb(n, x, w);
3000 cosqf(n, x, w);
3001 cosqfb = 0.f;
3002 for (i = 1; i <= n; ++i) {
3003 /* Computing MAX */
3004 r2 = cosqfb, r3 = (r1 = cf * x[i - 1] - y[i - 1], fabs(r1));
3005 cosqfb = dmax(r2,r3);
3006 }
3007
3008 /* TEST CFFTI,CFFTF,CFFTB */
3009
3010 for (i = 1; i <= n; ++i) {
3011 r1 = FFTPACK_COS(sqrt2 * (real) i);
3012 r2 = FFTPACK_SIN(sqrt2 * (real) (i * i));
3013 q1.r = r1, q1.i = r2;
3014 cx[i-1].r = q1.r, cx[i-1].i = q1.i;
3015 }
3016 dt = (2*M_PI) / fn;
3017 for (i = 1; i <= n; ++i) {
3018 arg1 = -((real) (i - 1)) * dt;
3019 cy[i-1].r = 0.f, cy[i-1].i = 0.f;
3020 for (k = 1; k <= n; ++k) {
3021 arg2 = (real) (k - 1) * arg1;
3022 r1 = FFTPACK_COS(arg2);
3023 r2 = FFTPACK_SIN(arg2);
3024 q3.r = r1, q3.i = r2;
3025 q2.r = q3.r * cx[k-1].r - q3.i * cx[k-1].i, q2.i =
3026 q3.r * cx[k-1].i + q3.i * cx[k-1].r;
3027 q1.r = cy[i-1].r + q2.r, q1.i = cy[i-1].i + q2.i;
3028 cy[i-1].r = q1.r, cy[i-1].i = q1.i;
3029 }
3030 }
3031 cffti(n, w);
3032 cfftf(n, (real*)cx, w);
3033 dcfftf = 0.f;
3034 for (i = 1; i <= n; ++i) {
3035 /* Computing MAX */
3036 q1.r = cx[i-1].r - cy[i-1].r, q1.i = cx[i-1].i - cy[i-1]
3037 .i;
3038 r1 = dcfftf, r2 = c_abs(&q1);
3039 dcfftf = dmax(r1,r2);
3040 q1.r = cx[i-1].r / fn, q1.i = cx[i-1].i / fn;
3041 cx[i-1].r = q1.r, cx[i-1].i = q1.i;
3042 }
3043 dcfftf /= fn;
3044 for (i = 1; i <= n; ++i) {
3045 arg1 = (real) (i - 1) * dt;
3046 cy[i-1].r = 0.f, cy[i-1].i = 0.f;
3047 for (k = 1; k <= n; ++k) {
3048 arg2 = (real) (k - 1) * arg1;
3049 r1 = FFTPACK_COS(arg2);
3050 r2 = FFTPACK_SIN(arg2);
3051 q3.r = r1, q3.i = r2;
3052 q2.r = q3.r * cx[k-1].r - q3.i * cx[k-1].i, q2.i =
3053 q3.r * cx[k-1].i + q3.i * cx[k-1].r;
3054 q1.r = cy[i-1].r + q2.r, q1.i = cy[i-1].i + q2.i;
3055 cy[i-1].r = q1.r, cy[i-1].i = q1.i;
3056 }
3057 }
3058 cfftb(n, (real*)cx, w);
3059 dcfftb = 0.f;
3060 for (i = 1; i <= n; ++i) {
3061 /* Computing MAX */
3062 q1.r = cx[i-1].r - cy[i-1].r, q1.i = cx[i-1].i - cy[i-1].i;
3063 r1 = dcfftb, r2 = c_abs(&q1);
3064 dcfftb = dmax(r1,r2);
3065 cx[i-1].r = cy[i-1].r, cx[i-1].i = cy[i-1].i;
3066 }
3067 cf = 1.f / fn;
3068 cfftf(n, (real*)cx, w);
3069 cfftb(n, (real*)cx, w);
3070 dcfb = 0.f;
3071 for (i = 1; i <= n; ++i) {
3072 /* Computing MAX */
3073 q2.r = cf * cx[i-1].r, q2.i = cf * cx[i-1].i;
3074 q1.r = q2.r - cy[i-1].r, q1.i = q2.i - cy[i-1].i;
3075 r1 = dcfb, r2 = c_abs(&q1);
3076 dcfb = dmax(r1,r2);
3077 }
3078 printf("%d\tRFFTF %10.3g\tRFFTB %10.ge\tRFFTFB %10.3g", n, rftf, rftb, rftfb);
3079 printf( "\tSINT %10.3g\tSINTFB %10.ge\tCOST %10.3g\n", sintt, sintfb, costt);
3080 printf( "\tCOSTFB %10.3g\tSINQF %10.ge\tSINQB %10.3g", costfb, sinqft, sinqbt);
3081 printf( "\tSINQFB %10.3g\tCOSQF %10.ge\tCOSQB %10.3g\n", sinqfb, cosqft, cosqbt);
3082 printf( "\tCOSQFB %10.3g\t", cosqfb);
3083 printf( "\tCFFTF %10.ge\tCFFTB %10.3g\n", dcfftf, dcfftb);
3084 printf( "\tCFFTFB %10.3g\n", dcfb);
3085
3086 #define CHECK(x) if (x > 1e-3) { printf(#x " failed: %g\n", x); all_ok = 0; }
3087 CHECK(rftf); CHECK(rftb); CHECK(rftfb); CHECK(sintt); CHECK(sintfb); CHECK(costt);
3088 CHECK(costfb); CHECK(sinqft); CHECK(sinqbt); CHECK(sinqfb); CHECK(cosqft); CHECK(cosqbt);
3089 CHECK(cosqfb); CHECK(dcfftf); CHECK(dcfftb);
3090 }
3091
3092 if (all_ok) printf("Everything looks fine.\n");
3093 else printf("ERRORS WERE DETECTED.\n");
3094 /*
3095 expected:
3096 120 RFFTF 2.786e-06 RFFTB 6.847e-04 RFFTFB 2.795e-07 SINT 1.312e-06 SINTFB 1.237e-06 COST 1.319e-06
3097 COSTFB 4.355e-06 SINQF 3.281e-04 SINQB 1.876e-06 SINQFB 2.198e-07 COSQF 6.199e-07 COSQB 2.193e-06
3098 COSQFB 2.300e-07 DEZF 5.573e-06 DEZB 1.363e-05 DEZFB 1.371e-06 CFFTF 5.590e-06 CFFTB 4.751e-05
3099 CFFTFB 4.215e-07
3100 54 RFFTF 4.708e-07 RFFTB 3.052e-05 RFFTFB 3.439e-07 SINT 3.532e-07 SINTFB 4.145e-07 COST 3.002e-07
3101 COSTFB 6.343e-07 SINQF 4.959e-05 SINQB 4.415e-07 SINQFB 2.882e-07 COSQF 2.826e-07 COSQB 2.472e-07
3102 COSQFB 3.439e-07 DEZF 9.388e-07 DEZB 5.066e-06 DEZFB 5.960e-07 CFFTF 1.426e-06 CFFTB 9.482e-06
3103 CFFTFB 2.980e-07
3104 49 RFFTF 4.476e-07 RFFTB 5.341e-05 RFFTFB 2.574e-07 SINT 9.196e-07 SINTFB 9.401e-07 COST 8.174e-07
3105 COSTFB 1.331e-06 SINQF 4.005e-05 SINQB 9.342e-07 SINQFB 3.057e-07 COSQF 2.530e-07 COSQB 6.228e-07
3106 COSQFB 4.826e-07 DEZF 9.071e-07 DEZB 4.590e-06 DEZFB 5.960e-07 CFFTF 2.095e-06 CFFTB 1.414e-05
3107 CFFTFB 7.398e-07
3108 32 RFFTF 4.619e-07 RFFTB 2.861e-05 RFFTFB 1.192e-07 SINT 3.874e-07 SINTFB 4.172e-07 COST 4.172e-07
3109 COSTFB 1.699e-06 SINQF 2.551e-05 SINQB 6.407e-07 SINQFB 2.980e-07 COSQF 1.639e-07 COSQB 1.714e-07
3110 COSQFB 2.384e-07 DEZF 1.013e-06 DEZB 2.339e-06 DEZFB 7.749e-07 CFFTF 1.127e-06 CFFTB 6.744e-06
3111 CFFTFB 2.666e-07
3112 4 RFFTF 1.490e-08 RFFTB 1.490e-07 RFFTFB 5.960e-08 SINT 7.451e-09 SINTFB 0.000e+00 COST 2.980e-08
3113 COSTFB 1.192e-07 SINQF 4.768e-07 SINQB 2.980e-08 SINQFB 5.960e-08 COSQF 2.608e-08 COSQB 5.960e-08
3114 COSQFB 1.192e-07 DEZF 2.980e-08 DEZB 5.960e-08 DEZFB 0.000e+00 CFFTF 6.664e-08 CFFTB 5.960e-08
3115 CFFTFB 6.144e-08
3116 3 RFFTF 3.974e-08 RFFTB 1.192e-07 RFFTFB 3.303e-08 SINT 1.987e-08 SINTFB 1.069e-08 COST 4.967e-08
3117 COSTFB 5.721e-08 SINQF 8.941e-08 SINQB 2.980e-08 SINQFB 1.259e-07 COSQF 7.451e-09 COSQB 4.967e-08
3118 COSQFB 7.029e-08 DEZF 1.192e-07 DEZB 5.960e-08 DEZFB 5.960e-08 CFFTF 7.947e-08 CFFTB 8.429e-08
3119 CFFTFB 9.064e-08
3120 2 RFFTF 0.000e+00 RFFTB 0.000e+00 RFFTFB 0.000e+00 SINT 0.000e+00 SINTFB 0.000e+00 COST 0.000e+00
3121 COSTFB 0.000e+00 SINQF 1.192e-07 SINQB 2.980e-08 SINQFB 5.960e-08 COSQF 7.451e-09 COSQB 1.490e-08
3122 COSQFB 0.000e+00 DEZF 0.000e+00 DEZB 0.000e+00 DEZFB 0.000e+00 CFFTF 0.000e+00 CFFTB 5.960e-08
3123 CFFTFB 5.960e-08
3124 Everything looks fine.
3125
3126 */
3127
3128 return all_ok ? 0 : 1;
3129 }
3130 #endif /* TESTING_FFTPACK */
3131