Procedure FOURN(var DATA:array of real;
                            NN:array of integer; NDIM,ISIGN:integer);
Label 1,2;
var
    NTOT,IDIM,NPREV,N,I2,I1,IP1,IP2,IP3,IPREW,I3REV,I2REW,NREM:integer;
    II2,II1,II3,I2REV,I3,IBIT,IFP1,IFP2,K1,K2:integer;
    WRR,WII,TEMPR,TEMPI:real;
    WR,WI,WPR,WPI,WTEMP,THETA:double;
begin
    NTOT:=1;
    For IDIM:=1 To NDIM do
        NTOT:=NTOT * NN[IDIM];
    NPREV:=1;
    For IDIM:=1 To NDIM do
    begin
        N:=NN[IDIM];
        NREM:=NTOT div (N * NPREV);
        IP1:=2 * NPREV;
        IP2:=IP1 * N;
        IP3:=IP2 * NREM;
        I2REV:=1;
        For II2:=0 To ((IP2-1) div IP1) do
        begin
            I2:= 1 + II2 * IP1;
            If I2 < I2REV Then
            begin
                For II1:=0 To ((IP1 - 2) div 2) do
                begin
                    I1:= I2 + II1*2;
                    For II3:=0 To ((IP3-I1) div IP2) do
                    begin
                        I3:= I1 + II3 * IP2;
                        I3REV:=I2REV + I3 - I2;
                        TEMPR:=Data[I3];
                        TEMPI:=Data[I3 + 1];
                        DATA[I3]:=Data[I3REV];
                        DATA[I3 + 1]:=Data[I3REV +1];
                        DATA[I3REV]:=TEMPR;
                        DATA[I3REV + 1]:=TEMPI;
                    end;
                end;
            end;
            IBIT:=IP2 div 2;
1:          If (IBIT >= IP1) And (I2REV > IBIT) Then
            begin
                I2REV:=I2REV - IBIT;
                IBIT:=IBIT div 2;
                GoTo 1;
            end;
            I2REV:=I2REV + IBIT;
        end;
        IFP1:=IP1;
2:      If IFP1 < IP2 Then
        begin
            IFP2:=2 * IFP1;
            THETA:=ISIGN * 6.28318530717959 / (IFP2 / IP1);
            WPR:=-2 * Sqr(Sin(0.5 * THETA));
            WPI:=Sin(THETA);
            WR:=1;
            WI:=0;
            For II3:=0 To ((IFP1-1) div IP1) do
            begin
                I3:= 1 + II3*IP1;
                For II1:=0 To ((IP1 - 2) div 2) do
                begin
                    I1:= I3 + II1*2;
                    For II2:=0 To ((IP3-I1) div IFP2) do
                    begin
                        I2:=I1 + II2*IFP2;
                        K1:=I2;
                        K2:=K1 + IFP1;
                        TEMPR:=WR * DATA[K2] - WI * DATA[K2 + 1];
                        TEMPI:=WR * DATA[K2 + 1] + WI * DATA[K2];
                        DATA[K2]:=DATA[K1] - TEMPR;
                        DATA[K2 + 1]:=DATA[K1 + 1] - TEMPI;
                        DATA[K1]:=DATA[K1] + TEMPR;
                        DATA[K1 + 1]:=DATA[K1 + 1] + TEMPI;
                    end;
                end;
                WTEMP:=WR;
                WR:=WR * WPR - WI * WPI + WR;
                WI:=WI * WPR + WTEMP * WPI + WI;
            end;
            IFP1:=IFP2;
            GoTo 2;
        end; 
        NPREV:=N * NPREV;
    end;
end;