1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
* SUBROUTINE 'SETDGPFA'
* SETUP ROUTINE FOR SELF-SORTING IN-PLACE
* GENERALIZED PRIME FACTOR (COMPLEX) FFT [DGPFA]
*
* CALL SETDGPFA(TRIGS,N,NPQR,INFO)
*
* INPUT :
* -----
* N IS THE LENGTH OF THE TRANSFORMS. N MUST BE OF THE FORM:
* -----------------------------------
* N = (2**IP) * (3**IQ) * (5**IR)
* -----------------------------------
*
* OUTPUT:
* ------
* TRIGS IS A TABLE OF TWIDDLE FACTORS,
* OF LENGTH 2*IPQR (DOUBLE PRECISION) WORDS, WHERE:
* --------------------------------------
* IPQR = (2**IP) + (3**IQ) + (5**IR)
* --------------------------------------
* NPQR = THREE INTEGERS HOLDING IP, IQ, IR
* INFO = SET TO 0 ON SUCCESS AND -1 ON FAILURE
*
* WRITTEN BY CLIVE TEMPERTON 1990
*
*----------------------------------------------------------------------
*
SUBROUTINE SETDGPFA(TRIGS,N,NPQR,INFO)
*
DOUBLE PRECISION TRIGS(*)
INTEGER N, NPQR(3), INFO
DIMENSION NJ(3)
DOUBLE PRECISION DEL
DOUBLE PRECISION ANGLE, TWOPI
INFO = 0
*
* DECOMPOSE N INTO FACTORS 2,3,5
* ------------------------------
NN = N
IFAC = 2
*
DO 30 LL = 1 , 3
KK = 0
10 CONTINUE
IF (MOD(NN,IFAC).NE.0) GO TO 20
KK = KK + 1
NN = NN / IFAC
GO TO 10
20 CONTINUE
NPQR(LL) = KK
IFAC = IFAC + LL
30 CONTINUE
*
IF (NN.NE.1) THEN
* WRITE(6,40) N
* 40 FORMAT(' *** WARNING!!!',I10,' IS NOT A LEGAL VALUE OF N ***')
INFO = -1
RETURN
ENDIF
*
IP = NPQR(1)
IQ = NPQR(2)
IR = NPQR(3)
*
* COMPUTE LIST OF ROTATED TWIDDLE FACTORS
* ---------------------------------------
NJ(1) = 2**IP
NJ(2) = 3**IQ
NJ(3) = 5**IR
*
TWOPI = 4.0 * ASIN(1.0)
I = 1
*
DO 60 LL = 1 , 3
NI = NJ(LL)
IF (NI.EQ.1) GO TO 60
*
DEL = TWOPI / DFLOAT(NI)
IROT = N / NI
KINK = MOD(IROT,NI)
KK = 0
*
DO 50 K = 1 , NI
ANGLE = DFLOAT(KK) * DEL
TRIGS(I) = COS(ANGLE)
TRIGS(I+1) = SIN(ANGLE)
I = I + 2
KK = KK + KINK
IF (KK.GT.NI) KK = KK - NI
50 CONTINUE
60 CONTINUE
*
RETURN
END
|
