C @(#)sfft.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:01 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C SUBROUTINE FFT VERSION 1.0 JAN. 30, 1983 C IEEE TRANS AU 17 C C FAST FOURIER ANALYSIS SUBROUTINE (FFA) C THIS PROGRAM IS TAKEN FROM IEEE TRANS. AU 17 (1969) 138 C THIS SUBROUTINE REPLACES THE REAL VECTOR B(K), (K=1,2,...,N) C WITH ITS FINITE DISCRETE FOURIER TRANSFORM X(J) WHICH IS RETURNED C SUCH THAT X(0) = B(1), X(N/2) = B(2), AND X(J) = B(2*J+1)+I*B(2*J+2), C (I.E. X(1) = B(3) + I B(4), X(2) = B(5) + I B(6), ETC.) THE C SUBROUTINE IS CALLED AS FFA(M,B) WHERE N = 2**M AND B IS AN C N TERM REAL ARRAY. A REAL-VALUED, RADIX 8 ALGORITHM IS USED C WITH IN-PLACE REORDERING AND THE TRIG FUNCTIONS ARE COMPUTED AS C NEEDED. C C SINGLE PRECISION VERSION C C--------------------------------------------------------------------- SUBROUTINE FFA(M,B) IMPLICIT NONE INTEGER M REAL B(1) C INTEGER N, NN INTEGER N8POW INTEGER INT, IT C N = 2**M N8POW = M/3 C C DO A RADIX 2 OR RADIX 4 ITERATION FIRST IF ONE IS REQUIRED C IF (M-N8POW*3-1) 30,20,10 10 NN = 4 INT = N/NN CALL R4TR(INT,B(1),B(INT+1),B(2*INT+1),B(3*INT+1)) GO TO 40 20 NN = 2 INT = N/NN CALL R2TR(INT,B(1),B(INT+1)) GO TO 40 30 NN = 1 C C PERFORM RADIX 8 ITERATIONS C 40 IF (N8POW) 70,70,45 45 DO 50 IT = 1,N8POW NN = NN*8 INT = N/NN CALL R8TR(INT,NN,B(1),B(INT+1),B(INT*2+1),B(3*INT+1), 1 B(4*INT+1),B(5*INT+1),B(6*INT+1),B(7*INT+1),B(1),B(INT+1), 2 B(2*INT+1),B(3*INT+1),B(4*INT+1),B(5*INT+1),B(6*INT+1), 3 B(7*INT+1)) 50 CONTINUE C C PERFORM IN-PLACE REORDERING C CALL ORD1(M,B) CALL ORD2(M,B) 70 RETURN END C-------------------------------------------------------------------- SUBROUTINE FFS(M,B) C C FAST FOURIER SYNTHESIS SUBROUTINE (FFS) C THIS SUBROUTINE SYNTHESIZES THE REAL VECTOR B(K), WHERE C K = 1,2,...N. THE INITIAL FOURIER COEFFICIENTS ARE PLACED IN C THE B ARRAY WITH THE DC AND HIGHEST HARMONIC COEFFICIENTS C STORED IN B(1) AND B(2) RESPECTIVELY AND ALL OTHERS LOCATED C SUCH THAT THE JTH HARMONIC IS STORED AS B(2*J+1) + I*B(2*J+2). C B IS THE N TERM REAL ARRAY DISCUSSED ABOVE. C------------------------------------------------------------------- IMPLICIT NONE INTEGER M REAL B(1) C INTEGER N, NN INTEGER N8POW INTEGER INT, IT C N = 2**M N8POW = M/3 IF (N8POW) 15,15,5 C C REORDER THE INPUT FOURIER COEFFICIENTS C 5 CALL ORD2(M,B) CALL ORD1(M,B) C C PERFORM THE RADIX 8 ITERATIONS C NN = 8*N DO 10 IT=1,N8POW NN = NN/8 INT = N/NN CALL R8SYN(INT,NN,B(1),B(INT+1),B(2*INT+1),B(3*INT+1), 1 B(4*INT+1),B(5*INT+1),B(6*INT+1),B(7*INT+1),B(1), 2 B(INT+1),B(2*INT+1),B(3*INT+1),B(4*INT+1),B(5*INT+1), 3 B(6*INT+1),B(7*INT+1)) 10 CONTINUE C C DO A RADIX 2 OR RADIX 4 ITERATION IF ONE IS REQUIRED C 15 IF (M-N8POW*3-1) 40,30,20 20 INT = N/4 CALL R4SYN(INT,B(1),B(INT+1),B(2*INT+1),B(INT*3+1)) GO TO 40 30 INT = N/2 CALL R2TR(INT,B(1),B(INT+1)) 40 RETURN END C--------------------------------------------------------------------- SUBROUTINE R2TR(INT,B0,B1) C C RADIX 2 ITERATION SUBROUTINE C IMPLICIT NONE INTEGER INT REAL B0(1),B1(1) C INTEGER K REAL T C DO 100 K=1,INT T = B0(K)+B1(K) B1(K) = B0(K)-B1(K) B0(K) = T 100 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE R4TR(INT,B0,B1,B2,B3) C C RADIX 4 ITERATION SUBROUTINE C IMPLICIT NONE INTEGER INT REAL B0(1),B1(1),B2(1),B3(1) C INTEGER K REAL R0, R1 C DO 200 K=1,INT R0 = B0(K)+B2(K) R1 = B1(K)+B3(K) B2(K) = B0(K)-B2(K) B3(K) = B1(K)-B3(K) B0(K) = R0+R1 B1(K) = R0-R1 200 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE R4SYN(INT,B0,B1,B2,B3) C C RADIX 4 SYNTHESIS SUBROUTINE C IMPLICIT NONE INTEGER INT REAL B0(1),B1(1),B2(1),B3(1) C INTEGER K REAL T0, T1, T2, T3 C DO 200 K=1,INT T0 = B0(K)+B1(K) T1 = B0(K)-B1(K) T2 = B2(K)+B2(K) T3 = B3(K)+B3(K) B0(K) = T0+T2 B2(K) = T0-T2 B1(K) = T1+T3 B3(K) = T1-T3 200 CONTINUE RETURN END C------------------------------------------------------------------- SUBROUTINE R8TR(INT,NN,BR0,BR1,BR2,BR3,BR4,BR5,BR6,BR7 1 ,BI0,BI1,BI2,BI3,BI4,BI5,BI6,BI7) C C RADIX 8 ITERATION SUBROUTINE C IMPLICIT NONE INTEGER INT INTEGER NN REAL BR0(1),BR1(1),BR2(1),BR3(1),BR4(1), 1 BR5(1),BR6(1),BR7(1),BI0(1),BI1(1),BI2(1),BI3(1), 2 BI4(1),BI5(1),BI6(1),BI7(1) INTEGER L(15) C INTEGER INT8 INTEGER J, J0, JLAST INTEGER JI, JL, JR INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9, J10 INTEGER J11, J12, J13, J14 INTEGER JTHET INTEGER K0, KL INTEGER K, KS INTEGER L1, L2, L3, L4, L5, L6, L7, L8, L9, L10 INTEGER L11, L12, L13, L14, L15 REAL ARG REAL C1, C2, C3, C4, C5, C6, C7 REAL C22 REAL PR, PI REAL PIOVN, P7 REAL S1, S2, S3, S4, S5, S6, S7 REAL S22 REAL T0, T1, T2, T3, T4, T5, T6, T7 REAL TH2 REAL TR0, TR1, TR2, TR3, TR4, TR5, TR6, TR7 REAL TI0, TI1, TI2, TI3, TI4, TI5, TI6, TI7 C EQUIVALENCE (L15,L(1)),(L14,L(2)),(L13,L(3)),(L12,L(4)), 1 (L11,L(5)),(L10,L(6)),(L9,L(7)),(L8,L(8)),(L7,L(9)), 2 (L6,L(10)),(L5,L(11)),(L4,L(12)),(L3,L(13)),(L2,L(14)), 3 (L1,L(15)) C L(1) = NN/8 DO 10 K=2,15 IF(L(K-1)-2) 11,12,13 11 L(K-1) = 2 12 L(K) = 2 GO TO 10 13 L(K) = L(K-1)/2 10 CONTINUE PIOVN = 3.14159265358979D0/DBLE(NN) P7 = 0.707106781186548D0 C22 = 0.923879532511287D0 S22 = 0.382683432365090D0 JI = 3 JL = 2 JR = 2 DO 70 J1=2,L1,2 DO 70 J2=J1,L2,L1 DO 70 J3=J2,L3,L2 DO 70 J4=J3,L4,L3 DO 70 J5=J4,L5,L4 DO 70 J6=J5,L6,L5 DO 70 J7=J6,L7,L6 DO 70 J8=J7,L8,L7 DO 70 J9=J8,L9,L8 DO 70 J10=J9,L10,L9 DO 70 J11=J10,L11,L10 DO 70 J12=J11,L12,L11 DO 70 J13=J12,L13,L12 DO 70 J14=J13,L14,L13 DO 70 JTHET=J14,L15,L14 TH2 = JTHET-2 IF(TH2) 71,71,76 71 DO 72 K=1,INT T0 = BR0(K)+BR4(K) T1 = BR1(K)+BR5(K) T2 = BR2(K)+BR6(K) T3 = BR3(K)+BR7(K) T4 = BR0(K)-BR4(K) T5 = BR1(K)-BR5(K) T6 = BR2(K)-BR6(K) T7 = BR3(K)-BR7(K) BR2(K) = T0-T2 BR3(K) = T1-T3 T0 = T0+T2 T1 = T1+T3 BR0(K) = T0+T1 BR1(K) = T0-T1 PR = P7*(T5-T7) PI = P7*(T5+T7) BR4(K) = T4+PR BR7(K) = T6+PI BR6(K) = T4-PR BR5(K) = PI-T6 72 CONTINUE IF(NN-8) 70,70,73 73 K0 = INT*8+1 KL = K0+INT-1 DO 75 K=K0,KL PR = P7*(BI2(K)-BI6(K)) PI = P7*(BI2(K)+BI6(K)) TR0 = BI0(K)+PR TI0 = BI4(K)+PI TR2 = BI0(K)-PR TI2 = BI4(K)-PI PR = P7*(BI3(K)-BI7(K)) PI = P7*(BI3(K)+BI7(K)) TR1 = BI1(K)+PR TI1 = BI5(K)+PI TR3 = BI1(K)-PR TI3 = BI5(K)-PI PR = TR1*C22-TI1*S22 PI = TI1*C22+TR1*S22 BI0(K) = TR0+PR BI6(K) = TR0-PR BI7(K) = TI0+PI BI1(K) = PI-TI0 PR = -TR3*S22-TI3*C22 PI = TR3*C22-TI3*S22 BI2(K) = TR2+PR BI4(K) = TR2-PR BI5(K) = TI2+PI BI3(K) = PI-TI2 75 CONTINUE GO TO 70 76 ARG = TH2*PIOVN C1 = COS(ARG) S1 = SIN(ARG) C2 = C1*C1-S1*S1 S2 = C1*S1+C1*S1 C3 = C1*C2-S1*S2 S3 = C2*S1+S2*C1 C4 = C2*C2-S2*S2 S4 = C2*S2+C2*S2 C5 = C2*C3-S2*S3 S5 = C3*S2+S3*C2 C6 = C3*C3-S3*S3 S6 = C3*S3+C3*S3 C7 = C3*C4-S3*S4 S7 = C4*S3+S4*C3 INT8 = INT*8 J0 = JR*INT8+1 K0 = JI*INT8+1 JLAST = J0+INT-1 DO 77 J=J0,JLAST KS = K0+J-J0 DO 77 K=KS,KS TR1 = BR1(J)*C1-BI1(K)*S1 TI1 = BR1(J)*S1+BI1(K)*C1 TR2 = BR2(J)*C2-BI2(K)*S2 TI2 = BR2(J)*S2+BI2(K)*C2 TR3 = BR3(J)*C3-BI3(K)*S3 TI3 = BR3(J)*S3+BI3(K)*C3 TR4 = BR4(J)*C4-BI4(K)*S4 TI4 = BR4(J)*S4+BI4(K)*C4 TR5 = BR5(J)*C5-BI5(K)*S5 TI5 = BR5(J)*S5+BI5(K)*C5 TR6 = BR6(J)*C6-BI6(K)*S6 TI6 = BR6(J)*S6+BI6(K)*C6 TR7 = BR7(J)*C7-BI7(K)*S7 TI7 = BR7(J)*S7+BI7(K)*C7 T0 = BR0(J)+TR4 T1 = BI0(K)+TI4 TR4 = BR0(J)-TR4 TI4 = BI0(K)-TI4 T2 = TR1+TR5 T3 = TI1+TI5 TR5 = TR1-TR5 TI5 = TI1-TI5 T4 = TR2+TR6 T5 = TI2+TI6 TR6 = TR2-TR6 TI6 = TI2-TI6 T6 = TR3+TR7 T7 = TI3+TI7 TR7 = TR3-TR7 TI7 = TI3-TI7 TR0 = T0+T4 TI0 = T1+T5 TR2 = T0-T4 TI2 = T1-T5 TR1 = T2+T6 TI1 = T3+T7 TR3 = T2-T6 TI3 = T3-T7 T0 = TR4-TI6 T1 = TI4+TR6 T4 = TR4+TI6 T5 = TI4-TR6 T2 = TR5-TI7 T3 = TI5+TR7 T6 = TR5+TI7 T7 = TI5-TR7 BR0(J) = TR0+TR1 BI7(K) = TI0+TI1 BI6(K) = TR0-TR1 BR1(J) = TI1-TI0 BR2(J) = TR2-TI3 BI5(K) = TI2+TR3 BI4(K) = TR2+TI3 BR3(J) = TR3-TI2 PR = P7*(T2-T3) PI = P7*(T2+T3) BR4(J) = T0+PR BI3(K) = T1+PI BI2(K) = T0-PR BR5(J) = PI-T1 PR = -P7*(T6+T7) PI = P7*(T6-T7) BR6(J) = T4+PR BI1(K) = T5+PI BI0(K) = T4-PR BR7(J) = PI-T5 77 CONTINUE JR = JR+2 JI = JI-2 IF (JI-JL) 78,78,70 78 JI = 2*JR-1 JL = JR 70 CONTINUE RETURN END C----------------------------------------------------------------- SUBROUTINE R8SYN(INT,NN,BR0,BR1,BR2,BR3,BR4,BR5,BR6,BR7 1 ,BI0,BI1,BI2,BI3,BI4,BI5,BI6,BI7) C C RADIX 8 SYNTHESIS SUBROUTINE C IMPLICIT NONE INTEGER INT INTEGER NN REAL BR0(1),BR1(1),BR2(1),BR3(1),BR4(1), 1 BR5(1),BR6(1),BR7(1),BI0(1),BI1(1),BI2(1),BI3(1), 2 BI4(1),BI5(1),BI6(1),BI7(1) INTEGER L(15) C INTEGER INT8 INTEGER J, J0, JLAST INTEGER JI, JL, JR INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9, J10 INTEGER J11, J12, J13, J14 INTEGER JTHET INTEGER K0, KL INTEGER K, KS INTEGER L1, L2, L3, L4, L5, L6, L7, L8, L9, L10 INTEGER L11, L12, L13, L14, L15 REAL ARG REAL C1, C2, C3, C4, C5, C6, C7 REAL C22 REAL PR, PI REAL PIOVN, P7 REAL P7TWO REAL RI, RR REAL S1, S2, S3, S4, S5, S6, S7 REAL S22 REAL T0, T1, T2, T3, T4, T5, T6, T7 REAL TH2 REAL TR0, TR1, TR2, TR3, TR4, TR5, TR6, TR7 REAL TI0, TI1, TI2, TI3, TI4, TI5, TI6, TI7 REAL T8, TT0, TT1 REAL TTR6, TTR7 C EQUIVALENCE (L15,L(1)),(L14,L(2)),(L13,L(3)),(L12,L(4)), 1 (L11,L(5)),(L10,L(6)),(L9,L(7)),(L8,L(8)),(L7,L(9)), 2 (L6,L(10)),(L5,L(11)),(L4,L(12)),(L3,L(13)),(L2,L(14)), 3 (L1,L(15)) C L(1) = NN/8 DO 10 K=2,15 IF (L(K-1)-2) 11,12,13 11 L(K-1) = 2 12 L(K) = 2 GO TO 10 13 L(K) = L(K-1)/2 10 CONTINUE PIOVN = 3.14159265358979D0/DBLE(NN) P7 = 0.707106781186548D0 P7TWO = 1.414213562373095D0 C22 = 0.923879532511287D0 S22 = 0.382683432365090D0 JI = 3 JL = 2 JR = 2 DO 70 J1=2,L1,2 DO 70 J2=J1,L2,L1 DO 70 J3=J2,L3,L2 DO 70 J4=J3,L4,L3 DO 70 J5=J4,L5,L4 DO 70 J6=J5,L6,L5 DO 70 J7=J6,L7,L6 DO 70 J8=J7,L8,L7 DO 70 J9=J8,L9,L8 DO 70 J10=J9,L10,L9 DO 70 J11=J10,L11,L10 DO 70 J12=J11,L12,L11 DO 70 J13=J12,L13,L12 DO 70 J14=J13,L14,L13 DO 70 JTHET=J14,L15,L14 TH2 = JTHET-2 IF (TH2) 71,71,76 71 DO 72 K=1,INT T0 = BR0(K)+BR1(K) T1 = BR0(K)-BR1(K) T2 = BR2(K)+BR2(K) T3 = BR3(K)+BR3(K) T4 = BR4(K)+BR6(K) T6 = BR7(K)-BR5(K) T5 = BR4(K)-BR6(K) T7 = BR7(K)+BR5(K) PR = P7*(T7+T5) PI = P7*(T7-T5) TT0 = T0+T2 TT1 = T1+T3 T2 = T0-T2 T3 = T1-T3 T4 = T4+T4 T5 = PR+PR T6 = T6+T6 T7 = PI+PI BR0(K) = TT0+T4 BR1(K) = TT1+T5 BR2(K) = T2+T6 BR3(K) = T3+T7 BR4(K) = TT0-T4 BR5(K) = TT1-T5 BR6(K) = T2-T6 BR7(K) = T3-T7 72 CONTINUE IF (NN-8) 70,70,73 73 K0 = INT*8+1 KL = K0+INT-1 DO 75 K=K0,KL T1 = BI0(K)+BI6(K) T2 = BI7(K)-BI1(K) T3 = BI0(K)-BI6(K) T4 = BI7(K)+BI1(K) PR = T3*C22+T4*S22 PI = T4*C22-T3*S22 T5 = BI2(K)+BI4(K) T6 = BI5(K)-BI3(K) T7 = BI2(K)-BI4(K) T8 = BI5(K)+BI3(K) RR = T8*C22-T7*S22 RI = -T8*S22-T7*C22 BI0(K) = (T1+T5)+(T1+T5) BI4(K) = (T2+T6)+(T2+T6) BI1(K) = (PR+RR)+(PR+RR) BI5(K) = (PI+RI)+(PI+RI) T5 = T1-T5 T6 = T2-T6 BI2(K) = P7TWO*(T6+T5) BI6(K) = P7TWO*(T6-T5) RR = PR-RR RI = PI-RI BI3(K) = P7TWO*(RI+RR) BI7(K) = P7TWO*(RI-RR) 75 CONTINUE GO TO 70 76 ARG = TH2*PIOVN C1 = COS(ARG) S1 = -SIN(ARG) C2 = C1*C1-S1*S1 S2 = C1*S1+C1*S1 C3 = C1*C2-S1*S2 S3 = C2*S1+S2*C1 C4 = C2*C2-S2*S2 S4 = C2*S2+C2*S2 C5 = C2*C3-S2*S3 S5 = C3*S2+S3*C2 C6 = C3*C3-S3*S3 S6 = C3*S3+C3*S3 C7 = C3*C4-S3*S4 S7 = C4*S3+S4*C3 INT8 = INT*8 J0 = JR*INT8+1 K0 = JI*INT8+1 JLAST = J0+INT-1 DO 77 J=J0,JLAST KS = K0+J-J0 DO 77 K=KS,KS TR0 = BR0(J)+BI6(K) TI0 = BI7(K)-BR1(J) TR1 = BR0(J)-BI6(K) TI1 = BI7(K)+BR1(J) TR2 = BR2(J)+BI4(K) TI2 = BI5(K)-BR3(J) TR3 = BI5(K)+BR3(J) TI3 = BI4(K)-BR2(J) TR4 = BR4(J)+BI2(K) TI4 = BI3(K)-BR5(J) T0 = BR4(J)-BI2(K) T1 = BI3(K)+BR5(J) TR5 = P7*(T0+T1) TI5 = P7*(T1-T0) TR6 = BR6(J)+BI0(K) TI6 = BI1(K)-BR7(J) T0 = BR6(J)-BI0(K) T1 = BI1(K)+BR7(J) TR7 = -P7*(T0-T1) TI7 = -P7*(T1+T0) T0 = TR0+TR2 T1 = TI0+TI2 T2 = TR1+TR3 T3 = TI1+TI3 TR2 = TR0-TR2 TI2 = TI0-TI2 TR3 = TR1-TR3 TI3 = TI1-TI3 T4 = TR4+TR6 T5 = TI4+TI6 T6 = TR5+TR7 T7 = TI5+TI7 TTR6 = TI4-TI6 TI6 = TR6-TR4 TTR7 = TI5-TI7 TI7 = TR7-TR5 BR0(J) = T0+T4 BI0(K) = T1+T5 BR1(J) = C1*(T2+T6)-S1*(T3+T7) BI1(K) = C1*(T3+T7)+S1*(T2+T6) BR2(J) = C2*(TR2+TTR6)-S2*(TI2+TI6) BI2(K) = C2*(TI2+TI6)+S2*(TR2+TTR6) BR3(J) = C3*(TR3+TTR7)-S3*(TI3+TI7) BI3(K) = C3*(TI3+TI7)+S3*(TR3+TTR7) BR4(J) = C4*(T0-T4)-S4*(T1-T5) BI4(K) = C4*(T1-T5)+S4*(T0-T4) BR5(J) = C5*(T2-T6)-S5*(T3-T7) BI5(K) = C5*(T3-T7)+S5*(T2-T6) BR6(J) = C6*(TR2-TTR6)-S6*(TI2-TI6) BI6(K) = C6*(TI2-TI6)+S6*(TR2-TTR6) BR7(J) = C7*(TR3-TTR7)-S7*(TI3-TI7) BI7(K) = C7*(TI3-TI7)+S7*(TR3-TTR7) 77 CONTINUE JR = JR+2 JI = JI-2 IF (JI-JL) 78,78,70 78 JI = 2*JR-1 JL = JR 70 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE ORD1(M,B) C C IN-PLACE REORDERING SUBROUTINES C IMPLICIT NONE INTEGER M REAL B(1) C INTEGER K, KL INTEGER N INTEGER J REAL T C K = 4 KL = 2 N = 2**M DO 94 J=4,N,2 IF (K-J) 92,92,91 91 T = B(J) B(J) = B(K) B(K) = T 92 K = K-2 IF (K-KL) 93,93,94 93 K = 2*J KL = J 94 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE ORD2(M,B) C IMPLICIT NONE INTEGER M REAL B(1) C INTEGER IJ INTEGER J1, J2, J3, J4, J5, J6, J7, J8, J9, J10 INTEGER J11, J12, J13, J14, J15, J16, J17, J18 INTEGER L1, L2, L3, L4, L5, L6, L7, L8, L9, L10 INTEGER L11, L12, L13, L14, L15, L16, L17, L18, L19 INTEGER JI INTEGER L(20) INTEGER N INTEGER K REAL T C EQUIVALENCE (L19,L(1)),(L18,L(2)),(L17,L(3)),(L16,L(4)), 1 (L15,L(5)),(L14,L(6)),(L13,L(7)),(L12,L(8)),(L11,L(9)), 2 (L10,L(10)),(L9,L(11)),(L8,L(12)),(L7,L(13)),(L6,L(14)), 3 (L5,L(15)),(L4,L(16)),(L3,L(17)),(L2,L(18)),(L1,L(19)) C N = 2**M L(1) = N DO 101 K=2,M 101 L(K) = L(K-1)/2 DO 102 K=M,19 102 L(K+1) = 2 IJ = 2 DO 103 J1=2,L1,2 DO 103 J2=J1,L2,L1 DO 103 J3=J2,L3,L2 DO 103 J4=J3,L4,L3 DO 103 J5=J4,L5,L4 DO 103 J6=J5,L6,L5 DO 103 J7=J6,L7,L6 DO 103 J8=J7,L8,L7 DO 103 J9=J8,L9,L8 DO 103 J10=J9,L10,L9 DO 103 J11=J10,L11,L10 DO 103 J12=J11,L12,L11 DO 103 J13=J12,L13,L12 DO 103 J14=J13,L14,L13 DO 103 J15=J14,L15,L14 DO 103 J16=J15,L16,L15 DO 103 J17=J16,L17,L16 DO 103 J18=J17,L18,L17 DO 103 JI =J18,L19,L18 IF (IJ-JI) 108,103,103 108 T = B(IJ-1) B(IJ-1) = B(JI-1) B(JI-1) = T T = B(IJ) B(IJ) = B(JI) B(JI) = T 103 IJ = IJ+2 RETURN END