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