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