1644
|
1 subroutine cffti1 (n,wa,ifac) |
|
2 implicit double precision (a-h,o-z) |
1645
|
3 dimension wa(*) ,ifac(*) ,ntryh(4) |
1644
|
4 data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/ |
|
5 nl = n |
|
6 nf = 0 |
|
7 j = 0 |
|
8 101 j = j+1 |
|
9 if (j-4) 102,102,103 |
|
10 102 ntry = ntryh(j) |
|
11 go to 104 |
|
12 103 ntry = ntry+2 |
|
13 104 nq = nl/ntry |
|
14 nr = nl-ntry*nq |
|
15 if (nr) 101,105,101 |
|
16 105 nf = nf+1 |
|
17 ifac(nf+2) = ntry |
|
18 nl = nq |
|
19 if (ntry .ne. 2) go to 107 |
|
20 if (nf .eq. 1) go to 107 |
|
21 do 106 i=2,nf |
|
22 ib = nf-i+2 |
|
23 ifac(ib+2) = ifac(ib+1) |
|
24 106 continue |
|
25 ifac(3) = 2 |
|
26 107 if (nl .ne. 1) go to 104 |
|
27 ifac(1) = n |
|
28 ifac(2) = nf |
|
29 tpi = 6.28318530717959d0 |
|
30 argh = tpi/dble(n) |
|
31 i = 2 |
|
32 l1 = 1 |
|
33 do 110 k1=1,nf |
|
34 ip = ifac(k1+2) |
|
35 ld = 0 |
|
36 l2 = l1*ip |
|
37 ido = n/l2 |
|
38 idot = ido+ido+2 |
|
39 ipm = ip-1 |
|
40 do 109 j=1,ipm |
|
41 i1 = i |
|
42 wa(i-1) = 1. |
|
43 wa(i) = 0. |
|
44 ld = ld+l1 |
|
45 fi = 0. |
|
46 argld = dble(ld)*argh |
|
47 do 108 ii=4,idot,2 |
|
48 i = i+2 |
|
49 fi = fi+1. |
|
50 arg = fi*argld |
|
51 wa(i-1) = cos(arg) |
|
52 wa(i) = sin(arg) |
|
53 108 continue |
|
54 if (ip .le. 5) go to 109 |
|
55 wa(i1-1) = wa(i-1) |
|
56 wa(i1) = wa(i) |
|
57 109 continue |
|
58 l1 = l2 |
|
59 110 continue |
|
60 return |
|
61 end |