fix packaging (into central file) and copyright dates.
[CommonLispStat.git] / lib / cfft.c
blob653db79de4776f2df63ddc5cc86aed8bf35f5743
1 /* from fftpkg package in cmlib and netlib -- translated by f2c and modified */
3 #include "xmath.h"
5 /*
6 Internal Routines
7 */
9 static int cffti1_(int *n, double *wa, double *ifac)
11 /* Initialized data */
12 static int ntryh[4] = { 3,4,2,5 };
14 /* System generated locals */
15 long i_1, i_2, i_3;
17 /* Local variables */
18 double argh;
19 long ntry;
20 int i, j;
21 long idot;
22 double argld;
23 int i1, ib;
24 long k1, l1, l2;
25 double fi;
26 int ii, nf;
27 long ip, ld, nl, nq, nr;
28 double arg;
29 long ido, ipm;
30 double tpi;
32 /* Parameter adjustments */
33 --wa;
34 --ifac;
36 /* Function Body */
37 nl = *n;
38 nf = 0;
39 j = 0;
41 L101:
42 ++j;
43 if (j - 4 <= 0) ntry = ntryh[j - 1];
44 else ntry += 2;
45 L104:
46 nq = nl / ntry;
47 nr = nl - ntry * nq;
48 if (nr != 0) goto L101;
49 ++nf;
50 ifac[nf + 2] = ntry;
51 nl = nq;
52 if (ntry == 2 && nf != 1) {
53 i_1 = nf;
54 for (i = 2; i <= i_1; ++i) {
55 ib = nf - i + 2;
56 ifac[ib + 2] = ifac[ib + 1];
58 ifac[3] = 2;
60 if (nl != 1) goto L104;
62 ifac[1] = *n;
63 ifac[2] = nf;
64 tpi = 6.28318530717959;
65 argh = tpi / (double) (*n);
66 i = 2;
67 l1 = 1;
68 i_1 = nf;
69 for (k1 = 1; k1 <= i_1; ++k1) {
70 ip = ifac[k1 + 2];
71 ld = 0;
72 l2 = l1 * ip;
73 ido = *n / l2;
74 idot = ido + ido + 2;
75 ipm = ip - 1;
76 i_2 = ipm;
77 for (j = 1; j <= i_2; ++j) {
78 i1 = i;
79 wa[i - 1] = 1.0;
80 wa[i] = 0.0;
81 ld += l1;
82 fi = 0.0;
83 argld = (double) ld * argh;
84 i_3 = idot;
85 for (ii = 4; ii <= i_3; ii += 2) {
86 i += 2;
87 fi += 1.0;
88 arg = fi * argld;
89 wa[i - 1] = cos(arg);
90 wa[i] = sin(arg);
92 if (ip > 5) {
93 wa[i1 - 1] = wa[i - 1];
94 wa[i1] = wa[i];
97 l1 = l2;
99 return 0;
102 static int pass_(int *nac, long *ido, long *ip, long *l1, long *idl1,
103 double *cc, double *c1, double *c2, double *ch,
104 double *ch2, double *wa, int *isign)
106 /* System generated locals */
107 long cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_dim2, ch_offset,
108 ch2_offset, c1_dim1, c1_dim2, c1_offset, ch2_dim1, c2_offset,
109 c2_dim1, i_1, i_2, i_3;
111 /* Local variables */
112 int i, j, k, l, ik;
113 long idij, idlj, idot, ipph, lc, jc, idj, idl, inc, idp;
114 double wai, war;
115 long ipp2;
117 /* Parameter adjustments */
118 cc_dim1 = *ido;
119 cc_dim2 = *ip;
120 cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
121 cc -= cc_offset;
122 c1_dim1 = *ido;
123 c1_dim2 = *l1;
124 c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
125 c1 -= c1_offset;
126 c2_dim1 = *idl1;
127 c2_offset = c2_dim1 + 1;
128 c2 -= c2_offset;
129 ch_dim1 = *ido;
130 ch_dim2 = *l1;
131 ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
132 ch -= ch_offset;
133 ch2_dim1 = *idl1;
134 ch2_offset = ch2_dim1 + 1;
135 ch2 -= ch2_offset;
136 --wa;
138 /* Function Body */
139 idot = *ido / 2;
140 ipp2 = *ip + 2;
141 ipph = (*ip + 1) / 2;
142 idp = *ip * *ido;
144 if (*ido >= *l1) {
145 i_1 = ipph;
146 for (j = 2; j <= i_1; ++j) {
147 jc = ipp2 - j;
148 i_2 = *l1;
149 for (k = 1; k <= i_2; ++k) {
150 i_3 = *ido;
151 for (i = 1; i <= i_3; ++i) {
152 ch[i + (k + j * ch_dim2) * ch_dim1] =
153 cc[i + (j + k * cc_dim2) * cc_dim1]
154 + cc[i + (jc + k * cc_dim2) * cc_dim1];
155 ch[i + (k + jc * ch_dim2) * ch_dim1] =
156 cc[i + (j + k * cc_dim2) * cc_dim1]
157 - cc[i + (jc + k * cc_dim2) * cc_dim1];
161 i_1 = *l1;
162 for (k = 1; k <= i_1; ++k) {
163 i_2 = *ido;
164 for (i = 1; i <= i_2; ++i) {
165 ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
169 else {
170 i_1 = ipph;
171 for (j = 2; j <= i_1; ++j) {
172 jc = ipp2 - j;
173 i_2 = *ido;
174 for (i = 1; i <= i_2; ++i) {
175 i_3 = *l1;
176 for (k = 1; k <= i_3; ++k) {
177 ch[i + (k + j * ch_dim2) * ch_dim1] =
178 cc[i + (j + k * cc_dim2) * cc_dim1]
179 + cc[i + (jc + k * cc_dim2) * cc_dim1];
180 ch[i + (k + jc * ch_dim2) * ch_dim1] =
181 cc[i + (j + k * cc_dim2) * cc_dim1]
182 - cc[i + (jc + k * cc_dim2) * cc_dim1];
186 i_1 = *ido;
187 for (i = 1; i <= i_1; ++i) {
188 i_2 = *l1;
189 for (k = 1; k <= i_2; ++k) {
190 ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
194 idl = 2 - *ido;
195 inc = 0;
196 i_1 = ipph;
197 for (l = 2; l <= i_1; ++l) {
198 lc = ipp2 - l;
199 idl += *ido;
200 i_2 = *idl1;
201 for (ik = 1; ik <= i_2; ++ik) {
202 c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1]
203 * ch2[ik + (ch2_dim1 << 1)];
204 c2[ik + lc * c2_dim1] = - *isign * wa[idl] * ch2[ik + *ip * ch2_dim1];
206 idlj = idl;
207 inc += *ido;
208 i_2 = ipph;
209 for (j = 3; j <= i_2; ++j) {
210 jc = ipp2 - j;
211 idlj += inc;
212 if (idlj > idp) {
213 idlj -= idp;
215 war = wa[idlj - 1];
216 wai = wa[idlj];
217 i_3 = *idl1;
218 for (ik = 1; ik <= i_3; ++ik) {
219 c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
220 c2[ik + lc * c2_dim1] -= *isign * wai * ch2[ik + jc * ch2_dim1];
224 i_1 = ipph;
225 for (j = 2; j <= i_1; ++j) {
226 i_2 = *idl1;
227 for (ik = 1; ik <= i_2; ++ik) {
228 ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
231 i_1 = ipph;
232 for (j = 2; j <= i_1; ++j) {
233 jc = ipp2 - j;
234 i_2 = *idl1;
235 for (ik = 2; ik <= i_2; ik += 2) {
236 ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1]
237 - c2[ik + jc * c2_dim1];
238 ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1]
239 + c2[ik + jc * c2_dim1];
240 ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1]
241 + c2[ik - 1 + jc * c2_dim1];
242 ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1]
243 - c2[ik - 1 + jc * c2_dim1];
246 *nac = 1;
247 if (*ido != 2) {
248 *nac = 0;
249 i_1 = *idl1;
250 for (ik = 1; ik <= i_1; ++ik) {
251 c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
253 i_1 = *ip;
254 for (j = 2; j <= i_1; ++j) {
255 i_2 = *l1;
256 for (k = 1; k <= i_2; ++k) {
257 c1[(k + j * c1_dim2) * c1_dim1 + 1] =
258 ch[(k + j * ch_dim2) * ch_dim1 + 1];
259 c1[(k + j * c1_dim2) * c1_dim1 + 2] =
260 ch[(k + j * ch_dim2) * ch_dim1 + 2];
263 if (idot <= *l1) {
264 idij = 0;
265 i_1 = *ip;
266 for (j = 2; j <= i_1; ++j) {
267 idij += 2;
268 i_2 = *ido;
269 for (i = 4; i <= i_2; i += 2) {
270 idij += 2;
271 i_3 = *l1;
272 for (k = 1; k <= i_3; ++k) {
273 c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
274 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] + *isign * wa[idij]
275 * ch[i + (k + j * ch_dim2) * ch_dim1];
276 c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
277 * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij]
278 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
283 else {
284 idj = 2 - *ido;
285 i_1 = *ip;
286 for (j = 2; j <= i_1; ++j) {
287 idj += *ido;
288 i_2 = *l1;
289 for (k = 1; k <= i_2; ++k) {
290 idij = idj;
291 i_3 = *ido;
292 for (i = 4; i <= i_3; i += 2) {
293 idij += 2;
294 c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
295 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1]
296 + *isign * wa[idij] * ch[i + (k + j * ch_dim2) * ch_dim1];
297 c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
298 * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij]
299 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
305 return 0;
308 static int pass2_(ido, l1, cc, ch, wa1, isign)
309 int *ido, *l1;
310 double *cc, *ch, *wa1;
311 int *isign;
313 /* System generated locals */
314 int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
316 /* Local variables */
317 int i, k;
318 double ti2, tr2;
320 /* Parameter adjustments */
321 cc_dim1 = *ido;
322 cc_offset = cc_dim1 * 3 + 1;
323 cc -= cc_offset;
324 ch_dim1 = *ido;
325 ch_dim2 = *l1;
326 ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
327 ch -= ch_offset;
328 --wa1;
330 /* Function Body */
331 if (*ido <= 2) {
332 i_1 = *l1;
333 for (k = 1; k <= i_1; ++k) {
334 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
335 + cc[((k << 1) + 2) * cc_dim1 + 1];
336 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
337 - cc[((k << 1) + 2) * cc_dim1 + 1];
338 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
339 + cc[((k << 1) + 2) * cc_dim1 + 2];
340 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
341 - cc[((k << 1) + 2) * cc_dim1 + 2];
344 else {
345 i_1 = *l1;
346 for (k = 1; k <= i_1; ++k) {
347 i_2 = *ido;
348 for (i = 2; i <= i_2; i += 2) {
349 ch[i - 1 + (k + ch_dim2) * ch_dim1]
350 = cc[i - 1 + ((k << 1) + 1) * cc_dim1]
351 + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
352 tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1]
353 - cc[i - 1 + ((k << 1) + 2) * cc_dim1];
354 ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1]
355 + cc[i + ((k << 1) + 2) * cc_dim1];
356 ti2 = cc[i + ((k << 1) + 1) * cc_dim1]
357 - cc[i + ((k << 1) + 2) * cc_dim1];
358 ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2
359 - *isign * wa1[i] * tr2;
360 ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2
361 + *isign * wa1[i] * ti2;
365 return 0;
368 static int pass3_(ido, l1, cc, ch, wa1, wa2, isign)
369 int *ido, *l1;
370 double *cc, *ch, *wa1, *wa2;
371 int *isign;
373 /* System generated locals */
374 int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
376 /* Local variables */
377 double taui, taur;
378 int i, k;
379 double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
381 /* Parameter adjustments */
382 cc_dim1 = *ido;
383 cc_offset = (cc_dim1 << 2) + 1;
384 cc -= cc_offset;
385 ch_dim1 = *ido;
386 ch_dim2 = *l1;
387 ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
388 ch -= ch_offset;
389 --wa1;
390 --wa2;
392 /* Function Body */
393 taur = -.5;
394 taui = -(*isign) * .866025403784439;
395 if (*ido == 2) {
396 i_1 = *l1;
397 for (k = 1; k <= i_1; ++k) {
398 tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
399 cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
400 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
401 ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
402 ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
403 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
404 cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1]
405 - cc[(k * 3 + 3) * cc_dim1 + 1]);
406 ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2]
407 - cc[(k * 3 + 3) * cc_dim1 + 2]);
408 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
409 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
410 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
411 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
414 else {
415 i_1 = *l1;
416 for (k = 1; k <= i_1; ++k) {
417 i_2 = *ido;
418 for (i = 2; i <= i_2; i += 2) {
419 tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1]
420 + cc[i - 1 + (k * 3 + 3) * cc_dim1];
421 cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
422 ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1)
423 * cc_dim1] + tr2;
424 ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * cc_dim1];
425 ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
426 ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + ti2;
427 cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1]
428 - cc[i - 1 + (k * 3 + 3) * cc_dim1]);
429 ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1]
430 - cc[i + (k * 3 + 3) * cc_dim1]);
431 dr2 = cr2 - ci3;
432 dr3 = cr2 + ci3;
433 di2 = ci2 + cr3;
434 di3 = ci2 - cr3;
435 ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2
436 - *isign * wa1[i] * dr2;
437 ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2
438 + *isign * wa1[i] * di2;
439 ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3
440 - *isign * wa2[i] * dr3;
441 ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3
442 + *isign * wa2[i] * di3;
446 return 0;
449 static int pass4_(ido, l1, cc, ch, wa1, wa2, wa3, isign)
450 int *ido, *l1;
451 double *cc, *ch, *wa1, *wa2, *wa3;
452 int *isign;
454 /* System generated locals */
455 int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
457 /* Local variables */
458 int i, k;
459 double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
461 /* Parameter adjustments */
462 cc_dim1 = *ido;
463 cc_offset = cc_dim1 * 5 + 1;
464 cc -= cc_offset;
465 ch_dim1 = *ido;
466 ch_dim2 = *l1;
467 ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
468 ch -= ch_offset;
469 --wa1;
470 --wa2;
471 --wa3;
473 /* Function Body */
474 if (*ido == 2) {
475 i_1 = *l1;
476 for (k = 1; k <= i_1; ++k) {
477 ti1 = cc[((k << 2) + 1) * cc_dim1 + 2]
478 - cc[((k << 2) + 3) * cc_dim1 + 2];
479 ti2 = cc[((k << 2) + 1) * cc_dim1 + 2]
480 + cc[((k << 2) + 3) * cc_dim1 + 2];
481 tr4 = *isign * (cc[((k << 2) + 2) * cc_dim1 + 2]
482 - cc[((k << 2) + 4) * cc_dim1 + 2]);
483 ti3 = cc[((k << 2) + 2) * cc_dim1 + 2]
484 + cc[((k << 2) + 4) * cc_dim1 + 2];
485 tr1 = cc[((k << 2) + 1) * cc_dim1 + 1]
486 - cc[((k << 2) + 3) * cc_dim1 + 1];
487 tr2 = cc[((k << 2) + 1) * cc_dim1 + 1]
488 + cc[((k << 2) + 3) * cc_dim1 + 1];
489 ti4 = *isign * (cc[((k << 2) + 4) * cc_dim1 + 1]
490 - cc[((k << 2) + 2) * cc_dim1 + 1]);
491 tr3 = cc[((k << 2) + 2) * cc_dim1 + 1]
492 + cc[((k << 2) + 4) * cc_dim1 + 1];
493 ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
494 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
495 ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
496 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
497 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
498 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
499 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
500 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
503 else {
504 i_1 = *l1;
505 for (k = 1; k <= i_1; ++k) {
506 i_2 = *ido;
507 for (i = 2; i <= i_2; i += 2) {
508 ti1 = cc[i + ((k << 2) + 1) * cc_dim1]
509 - cc[i + ((k << 2) + 3) * cc_dim1];
510 ti2 = cc[i + ((k << 2) + 1) * cc_dim1]
511 + cc[i + ((k << 2) + 3) * cc_dim1];
512 ti3 = cc[i + ((k << 2) + 2) * cc_dim1]
513 + cc[i + ((k << 2) + 4) * cc_dim1];
514 tr4 = *isign * (cc[i + ((k << 2) + 2) * cc_dim1]
515 - cc[i + ((k << 2) + 4) * cc_dim1]);
516 tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1]
517 - cc[i - 1 + ((k << 2) + 3) * cc_dim1];
518 tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1]
519 + cc[i - 1 + ((k << 2) + 3) * cc_dim1];
520 ti4 = *isign * (cc[i - 1 + ((k << 2) + 4) * cc_dim1]
521 - cc[i - 1 + ((k << 2) + 2) * cc_dim1]);
522 tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1]
523 + cc[i - 1 + ((k << 2) + 4) * cc_dim1];
524 ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
525 cr3 = tr2 - tr3;
526 ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
527 ci3 = ti2 - ti3;
528 cr2 = tr1 + tr4;
529 cr4 = tr1 - tr4;
530 ci2 = ti1 + ti4;
531 ci4 = ti1 - ti4;
532 ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2
533 + *isign * wa1[i] * ci2;
534 ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2
535 - *isign * wa1[i] * cr2;
536 ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3
537 + *isign * wa2[i] * ci3;
538 ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3
539 - *isign * wa2[i] * cr3;
540 ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4
541 + *isign * wa3[i] * ci4;
542 ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4
543 - *isign * wa3[i] * cr4;
547 return 0;
550 static int pass5_(ido, l1, cc, ch, wa1, wa2, wa3, wa4, isign)
551 int *ido, *l1;
552 double *cc, *ch, *wa1, *wa2, *wa3, *wa4;
553 int *isign;
555 /* System generated locals */
556 int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
558 /* Local variables */
559 int i, k;
560 double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
561 ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5, ti11, ti12, tr11,
562 tr12;
564 /* Parameter adjustments */
565 cc_dim1 = *ido;
566 cc_offset = cc_dim1 * 6 + 1;
567 cc -= cc_offset;
568 ch_dim1 = *ido;
569 ch_dim2 = *l1;
570 ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
571 ch -= ch_offset;
572 --wa1;
573 --wa2;
574 --wa3;
575 --wa4;
577 /* Function Body */
578 tr11 = .309016994374947;
579 ti11 = -(*isign) * .951056516295154;
580 tr12 = -.809016994374947;
581 ti12 = -(*isign) * .587785252292473;
582 if (*ido == 2) {
583 i_1 = *l1;
584 for (k = 1; k <= i_1; ++k) {
585 ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
586 ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
587 ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
588 ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
589 tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
590 tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
591 tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
592 tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
593 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1]
594 + tr2 + tr3;
595 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2]
596 + ti2 + ti3;
597 cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
598 ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
599 cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
600 ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
601 cr5 = ti11 * tr5 + ti12 * tr4;
602 ci5 = ti11 * ti5 + ti12 * ti4;
603 cr4 = ti12 * tr5 - ti11 * tr4;
604 ci4 = ti12 * ti5 - ti11 * ti4;
605 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
606 ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
607 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
608 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
609 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
610 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
611 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
612 ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
615 else {
616 i_1 = *l1;
617 for (k = 1; k <= i_1; ++k) {
618 i_2 = *ido;
619 for (i = 2; i <= i_2; i += 2) {
620 ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * cc_dim1];
621 ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * cc_dim1];
622 ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * cc_dim1];
623 ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * cc_dim1];
624 tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1]
625 - cc[i - 1 + (k * 5 + 5) * cc_dim1];
626 tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1]
627 + cc[i - 1 + (k * 5 + 5) * cc_dim1];
628 tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1]
629 - cc[i - 1 + (k * 5 + 4) * cc_dim1];
630 tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1]
631 + cc[i - 1 + (k * 5 + 4) * cc_dim1];
632 ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1)
633 * cc_dim1] + tr2 + tr3;
634 ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1)
635 * cc_dim1] + ti2 + ti3;
636 cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
637 ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
639 cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
640 ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
642 cr5 = ti11 * tr5 + ti12 * tr4;
643 ci5 = ti11 * ti5 + ti12 * ti4;
644 cr4 = ti12 * tr5 - ti11 * tr4;
645 ci4 = ti12 * ti5 - ti11 * ti4;
646 dr3 = cr3 - ci4;
647 dr4 = cr3 + ci4;
648 di3 = ci3 + cr4;
649 di4 = ci3 - cr4;
650 dr5 = cr2 + ci5;
651 dr2 = cr2 - ci5;
652 di5 = ci2 - cr5;
653 di2 = ci2 + cr5;
654 ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2
655 + *isign * wa1[i] * di2;
656 ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2
657 - *isign * wa1[i] * dr2;
658 ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3
659 + *isign * wa2[i] * di3;
660 ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3
661 - *isign * wa2[i] * dr3;
662 ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4
663 + *isign * wa3[i] * di4;
664 ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4
665 - *isign * wa3[i] * dr4;
666 ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5
667 + *isign * wa4[i] * di5;
668 ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5
669 - *isign * wa4[i] * dr5;
673 return 0;
677 /*static */
678 int
679 cfft1_(int *n, double *c, double *ch, double *wa,
680 double /* int??*/ *ifac, int *isign)
682 /* System generated locals */
683 long i_1;
685 /* Local variables */
686 long idot;
687 int i, k1, n2, na, nac;
688 long l1, l2, ido, idl1, iw, ix2, ix3, ix4;
689 long nf, ip;
691 /* Parameter adjustments */
692 --c;
693 --ch;
694 --wa;
695 --ifac;
697 /* Function Body */
698 nf = ifac[2];
699 na = 0;
700 l1 = 1;
701 iw = 1;
702 i_1 = nf;
703 for (k1 = 1; k1 <= i_1; ++k1) {
704 ip = ifac[k1 + 2];
705 l2 = ip * l1;
706 ido = *n / l2;
707 idot = ido + ido;
708 idl1 = idot * l1;
709 if (ip == 4) {
710 ix2 = iw + idot;
711 ix3 = ix2 + idot;
712 if (na == 0) {
713 pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
715 else {
716 pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
718 na = 1 - na;
720 else if (ip == 2) {
721 if (na == 0) {
722 pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign);
724 else {
725 pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign);
727 na = 1 - na;
729 else if (ip == 3) {
730 ix2 = iw + idot;
731 if (na == 0) {
732 pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign);
734 else {
735 pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign);
737 na = 1 - na;
739 else if (ip == 5) {
740 ix2 = iw + idot;
741 ix3 = ix2 + idot;
742 ix4 = ix3 + idot;
743 if (na == 0) {
744 pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3],
745 &wa[ix4], isign);
747 else {
748 pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3],
749 &wa[ix4], isign);
751 na = 1 - na;
753 else {
754 if (na == 0) {
755 pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1],
756 &ch[1], &wa[iw], isign);
758 else {
759 pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1],
760 &c[1], &wa[iw], isign);
762 if (nac != 0) {
763 na = 1 - na;
766 l1 = l2;
767 iw += (ip - 1) * idot;
769 if (na != 0) {
770 n2 = *n + *n;
771 i_1 = n2;
772 for (i = 1; i <= i_1; ++i) {
773 c[i] = ch[i];
776 return 0;
781 Public Routine
785 * cfft computes the forward or backward complex discrete fourier transform.
787 * Input Parameters:
789 * n The length of the complex sequence c. The method is
790 * more efficient when n is the product of small primes.
792 * c A complex array of length n which contains the sequence
794 * wsave a real work array which must be dimensioned at least 4n+15
795 * in the program that calls cfft.
796 * isign 1 for transform, -1 for inverse transform.
797 * A call of cfft with isign = 1 followed by a call of cfft with
798 * isign = -1 will multiply the sequence by n.
800 * Output Parameters:
802 * c For j=1,...,n
804 * c(j)=the sum from k=1,...,n of
806 * c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n)
808 * where i=sqrt(-1)
811 int cfft(int n, double *c, double *wsave, int isign)
813 int iw1, iw2;
815 /* Parameter adjustments */
816 --c;
817 --wsave;
819 /* Function Body */
820 if (n != 1) {
821 iw1 = n + n + 1;
822 iw2 = iw1 + n + n;
823 cffti1_(&n, &wsave[iw1], &wsave[iw2]);
824 cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], &wsave[iw2], &isign);
826 return 0;