+ initial tests based on ALGLIB submitted by Sergey Bochkanov

git-svn-id: trunk@14398 -

.gitattributes

@@ -8276,6 +8276,20 @@ tests/tbs/ub0489.pp svneol=native#text/plain
 tests/tbs/ub0489b.pp svneol=native#text/plain
 tests/tbs/ub0506.pp svneol=native#text/plain
 tests/test/README.txt svneol=native#text/plain
+tests/test/alglib/t_testconvunit.pp svneol=native#text/plain
+tests/test/alglib/t_testcorrunit.pp svneol=native#text/plain
+tests/test/alglib/t_testfftunit.pp svneol=native#text/plain
+tests/test/alglib/t_testfhtunit.pp svneol=native#text/plain
+tests/test/alglib/u_ap.pp svneol=native#text/plain
+tests/test/alglib/u_conv.pp svneol=native#text/plain
+tests/test/alglib/u_corr.pp svneol=native#text/plain
+tests/test/alglib/u_fft.pp svneol=native#text/plain
+tests/test/alglib/u_fht.pp svneol=native#text/plain
+tests/test/alglib/u_ftbase.pp svneol=native#text/plain
+tests/test/alglib/u_testconvunit.pp svneol=native#text/plain
+tests/test/alglib/u_testcorrunit.pp svneol=native#text/plain
+tests/test/alglib/u_testfftunit.pp svneol=native#text/plain
+tests/test/alglib/u_testfhtunit.pp svneol=native#text/plain
 tests/test/cg/cdecl/taoc1.pp svneol=native#text/plain
 tests/test/cg/cdecl/taoc2.pp svneol=native#text/plain
 tests/test/cg/cdecl/taoc3.pp svneol=native#text/plain

tests/Makefile

@@ -1464,7 +1464,7 @@ TESTPACKAGESUBDIRS=cg packages/win-base packages/webtbs packages/hash packages/f
 ifdef QUICKTEST
 export QUICKTEST
 else
-override TESTSUBDIRS+=TESTSUBDIRS $(TESTPACKAGESUBDIRS)
+override TESTSUBDIRS+=$(TESTPACKAGESUBDIRS) alglib
 endif
 TESTDIRS:=test $(addprefix test/,$(TESTSUBDIRS))
 .PHONY: utils units copyfiles testprep

tests/Makefile.fpc

@@ -126,7 +126,7 @@ TESTPACKAGESUBDIRS=cg packages/win-base packages/webtbs packages/hash packages/f
 ifdef QUICKTEST
 export QUICKTEST
 else
-override TESTSUBDIRS+=TESTSUBDIRS $(TESTPACKAGESUBDIRS)
+override TESTSUBDIRS+=$(TESTPACKAGESUBDIRS) alglib
 endif
 
 # All full dirnames in the test/ dir including the subdir self

tests/test/alglib/t_testconvunit.pp

@@ -0,0 +1,19 @@
+{ %opt=-Mdelphi -Sa }
+program t_testconvunit;
+uses Sysutils, u_testconvunit;
+
+begin
+    Randomize();
+    try 
+        if not testconvunit_test_silent() then
+        begin
+            //WriteLn('*');
+            Halt(1);
+        end;
+    except on E: Exception do 
+        begin
+            Halt(2);
+        end;
+    end;
+    Halt(0);
+end.

tests/test/alglib/t_testcorrunit.pp

@@ -0,0 +1,19 @@
+{ %opt=-Mdelphi -Sa }
+program t_testcorrunit;
+uses Sysutils, u_testcorrunit;
+
+begin
+    Randomize();
+    try 
+        if not testcorrunit_test_silent() then
+        begin
+            //WriteLn('*');
+            Halt(1);
+        end;
+    except on E: Exception do 
+        begin
+            Halt(2);
+        end;
+    end;
+    Halt(0);
+end.

tests/test/alglib/t_testfftunit.pp

@@ -0,0 +1,19 @@
+{ %opt=-Mdelphi -Sa }
+program t_testfftunit;
+uses Sysutils, u_testfftunit;
+
+begin
+    Randomize();
+    try 
+        if not testfftunit_test_silent() then
+        begin
+            //WriteLn('*');
+            Halt(1);
+        end;
+    except on E: Exception do 
+        begin
+            Halt(2);
+        end;
+    end;
+    Halt(0);
+end.

tests/test/alglib/t_testfhtunit.pp

@@ -0,0 +1,19 @@
+{ %opt=-Mdelphi -Sa }
+program t_testfhtunit;
+uses Sysutils, u_testfhtunit;
+
+begin
+    Randomize();
+    try 
+        if not testfhtunit_test_silent() then
+        begin
+            //WriteLn('*');
+            Halt(1);
+        end;
+    except on E: Exception do 
+        begin
+            Halt(2);
+        end;
+    end;
+    Halt(0);
+end.

tests/test/alglib/u_ap.pp

@@ -0,0 +1,741 @@
+(*************************************************************************
+AP library
+Copyright (c) 2003-2009 Sergey Bochkanov (ALGLIB project).
+
+>>> LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************)
+unit u_Ap;
+
+interface
+
+uses Math, Sysutils;
+
+/////////////////////////////////////////////////////////////////////////
+// constants
+/////////////////////////////////////////////////////////////////////////
+const
+    MachineEpsilon = 5E-16;
+    MaxRealNumber = 1E300;
+    MinRealNumber = 1E-300;
+
+/////////////////////////////////////////////////////////////////////////
+// arrays
+/////////////////////////////////////////////////////////////////////////
+type
+    PDouble = ^Double;
+
+    Complex = record
+        X, Y: Double;
+    end;
+
+    TInteger1DArray     = array of LongInt;
+    TReal1DArray        = array of Double;
+    TComplex1DArray     = array of Complex;
+    TBoolean1DArray     = array of Boolean;
+
+    TInteger2DArray     = array of array of LongInt;
+    TReal2DArray        = array of array of Double;
+    TComplex2DArray     = array of array of Complex;
+    TBoolean2DArray     = array of array of Boolean;
+
+    RCommState = record
+        Stage:  LongInt;
+        IA:     TInteger1DArray;
+        BA:     TBoolean1DArray;
+        RA:     TReal1DArray;
+        CA:     TComplex1DArray;
+    end;
+
+/////////////////////////////////////////////////////////////////////////
+// Functions/procedures
+/////////////////////////////////////////////////////////////////////////
+function AbsReal(X : Double):Double;
+function AbsInt (I : Integer):Integer;
+function RandomReal():Double;
+function RandomInteger(I : Integer):Integer;
+function Sign(X:Double):Integer;
+function AP_Sqr(X:Double):Double;
+
+function DynamicArrayCopy(const A: TInteger1DArray):TInteger1DArray;overload;
+function DynamicArrayCopy(const A: TReal1DArray):TReal1DArray;overload;
+function DynamicArrayCopy(const A: TComplex1DArray):TComplex1DArray;overload;
+function DynamicArrayCopy(const A: TBoolean1DArray):TBoolean1DArray;overload;
+
+function DynamicArrayCopy(const A: TInteger2DArray):TInteger2DArray;overload;
+function DynamicArrayCopy(const A: TReal2DArray):TReal2DArray;overload;
+function DynamicArrayCopy(const A: TComplex2DArray):TComplex2DArray;overload;
+function DynamicArrayCopy(const A: TBoolean2DArray):TBoolean2DArray;overload;
+
+function AbsComplex(const Z : Complex):Double;
+function Conj(const Z : Complex):Complex;
+function CSqr(const Z : Complex):Complex;
+
+function C_Complex(const X : Double):Complex;
+function C_Opposite(const Z : Complex):Complex;
+function C_Add(const Z1 : Complex; const Z2 : Complex):Complex;
+function C_Mul(const Z1 : Complex; const Z2 : Complex):Complex;
+function C_AddR(const Z1 : Complex; const R : Double):Complex;
+function C_MulR(const Z1 : Complex; const R : Double):Complex;
+function C_Sub(const Z1 : Complex; const Z2 : Complex):Complex;
+function C_SubR(const Z1 : Complex; const R : Double):Complex;
+function C_RSub(const R : Double; const Z1 : Complex):Complex;
+function C_Div(const Z1 : Complex; const Z2 : Complex):Complex;
+function C_DivR(const Z1 : Complex; const R : Double):Complex;
+function C_RDiv(const R : Double; const Z2 : Complex):Complex;
+function C_Equal(const Z1 : Complex; const Z2 : Complex):Boolean;
+function C_NotEqual(const Z1 : Complex; const Z2 : Complex):Boolean;
+function C_EqualR(const Z1 : Complex; const R : Double):Boolean;
+function C_NotEqualR(const Z1 : Complex; const R : Double):Boolean;
+
+/////////////////////////////////////////////////////////////////////////
+// AP BLAS generic interface
+/////////////////////////////////////////////////////////////////////////
+//procedure UseAPBLAS(Flag: Boolean);
+function APVDotProduct(
+   V1: PDouble; I11, I12: Integer;
+   V2: PDouble; I21, I22: Integer):Double;
+procedure APVMove(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);overload;
+procedure APVMove(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer;
+   S: Double);overload;
+procedure APVMoveNeg(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);
+procedure APVAdd(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);overload;
+procedure APVAdd(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer;
+   S: Real);overload;
+procedure APVSub(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);overload;
+procedure APVSub(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer;
+   S: Real);overload;
+procedure APVMul(
+   VOp: PDouble; I1, I2: Integer;
+   S: Real);
+
+/////////////////////////////////////////////////////////////////////////
+// IEEE-compliant functions, placed at the end, under 'non-optimization'
+// compiler switch
+/////////////////////////////////////////////////////////////////////////
+function AP_Double(X:Double):Double;
+function AP_FP_Eq(X:Double; Y:Double):Boolean;
+function AP_FP_NEq(X:Double; Y:Double):Boolean;
+function AP_FP_Less(X:Double; Y:Double):Boolean;
+function AP_FP_Less_Eq(X:Double; Y:Double):Boolean;
+function AP_FP_Greater(X:Double; Y:Double):Boolean;
+function AP_FP_Greater_Eq(X:Double; Y:Double):Boolean;
+
+{var
+    // pointers to AP BLAS functions
+    ASMDotProduct1: function(V1: PDouble; V2: PDouble; N: Integer):Double;cdecl;
+    ASMMove1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;
+    ASMMoveS1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer; S: Double);cdecl;
+    ASMMoveNeg1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;
+    ASMAdd1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;
+    ASMAddS1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer; S: Double);cdecl;
+    ASMSub1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;
+}
+
+implementation
+
+{var
+    // use ablas.dll (ALGLIB BLAS) if found
+    UseAPBLASFlag: Boolean = True;
+}
+    // pointers to AP BLAS functions
+{$IFNDEF NOABLAS}
+{    ASMDotProduct1: function(V1: PDouble; V2: PDouble; N: Integer):Double;cdecl;
+    ASMMove1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;
+    ASMMoveS1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer; S: Double);cdecl;
+    ASMMoveNeg1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;
+    ASMAdd1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;
+    ASMAddS1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer; S: Double);cdecl;
+    ASMSub1: procedure(VDst: PDouble; VSrc: PDouble; N: Integer);cdecl;}
+{$ENDIF}
+
+/////////////////////////////////////////////////////////////////////////
+// Functions/procedures
+/////////////////////////////////////////////////////////////////////////
+function AbsReal(X : Double):Double;
+begin
+    //Result := Abs(X);
+    if X>=0 then
+        Result:=X
+    else
+        Result:=-X;
+end;
+
+function AbsInt (I : Integer):Integer;
+begin
+    //Result := Abs(I);
+    if I>=0 then
+        Result:=I
+    else
+        Result:=-I;
+end;
+
+function RandomReal():Double;
+begin
+    Result := Random;
+end;
+
+function RandomInteger(I : Integer):Integer;
+begin
+    Result := Random(I);
+end;
+
+function Sign(X:Double):Integer;
+begin
+    if X>0 then
+        Result := 1
+    else if X<0 then
+        Result := -1
+    else
+        Result := 0;
+end;
+
+function AP_Sqr(X:Double):Double;
+begin
+    Result := X*X;
+end;
+
+/////////////////////////////////////////////////////////////////////////
+// dynamical arrays copying
+/////////////////////////////////////////////////////////////////////////
+function DynamicArrayCopy(const A: TInteger1DArray):TInteger1DArray;overload;
+var
+    I:  Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+        Result[I]:=A[I];
+end;
+
+function DynamicArrayCopy(const A: TReal1DArray):TReal1DArray;overload;
+var
+    I:  Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+        Result[I]:=A[I];
+end;
+
+function DynamicArrayCopy(const A: TComplex1DArray):TComplex1DArray;overload;
+var
+    I:  Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+        Result[I]:=A[I];
+end;
+
+function DynamicArrayCopy(const A: TBoolean1DArray):TBoolean1DArray;overload;
+var
+    I:  Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+        Result[I]:=A[I];
+end;
+
+function DynamicArrayCopy(const A: TInteger2DArray):TInteger2DArray;overload;
+var
+    I,J:    Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+    begin
+        SetLength(Result[I], High(A[I])+1);
+        for J:=Low(A[I]) to High(A[I]) do
+            Result[I,J]:=A[I,J];
+    end;
+end;
+
+function DynamicArrayCopy(const A: TReal2DArray):TReal2DArray;overload;
+var
+    I,J:    Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+    begin
+        SetLength(Result[I], High(A[I])+1);
+        for J:=Low(A[I]) to High(A[I]) do
+            Result[I,J]:=A[I,J];
+    end;
+end;
+
+function DynamicArrayCopy(const A: TComplex2DArray):TComplex2DArray;overload;
+var
+    I,J:    Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+    begin
+        SetLength(Result[I], High(A[I])+1);
+        for J:=Low(A[I]) to High(A[I]) do
+            Result[I,J]:=A[I,J];
+    end;
+end;
+
+function DynamicArrayCopy(const A: TBoolean2DArray):TBoolean2DArray;overload;
+var
+    I,J:    Integer;
+begin
+    SetLength(Result, High(A)+1);
+    for I:=Low(A) to High(A) do
+    begin
+        SetLength(Result[I], High(A[I])+1);
+        for J:=Low(A[I]) to High(A[I]) do
+            Result[I,J]:=A[I,J];
+    end;
+end;
+
+/////////////////////////////////////////////////////////////////////////
+// complex numbers
+/////////////////////////////////////////////////////////////////////////
+function AbsComplex(const Z : Complex):Double;
+var
+    W : Double;
+    XABS : Double;
+    YABS : Double;
+    V : Double;
+begin
+    XABS := AbsReal(Z.X);
+    YABS := AbsReal(Z.Y);
+    W := Max(XABS, YABS);
+    V := Min(XABS, YABS);
+    if V=0 then
+    begin
+        Result := W;
+    end
+    else
+    begin
+        Result := W*SQRT(1+Sqr(V/W));
+    end;
+end;
+
+
+function Conj(const Z : Complex):Complex;
+begin
+    Result.X := Z.X;
+    Result.Y := -Z.Y;
+end;
+
+
+function CSqr(const Z : Complex):Complex;
+begin
+    Result.X := Sqr(Z.X)-Sqr(Z.Y);
+    Result.Y := 2*Z.X*Z.Y;
+end;
+
+
+function C_Complex(const X : Double):Complex;
+begin
+    Result.X := X;
+    Result.Y := 0;
+end;
+
+
+function C_Opposite(const Z : Complex):Complex;
+begin
+    Result.X := -Z.X;
+    Result.Y := -Z.Y;
+end;
+
+
+function C_Add(const Z1 : Complex; const Z2 : Complex):Complex;
+begin
+    Result.X := Z1.X+Z2.X;
+    Result.Y := Z1.Y+Z2.Y;
+end;
+
+
+function C_Mul(const Z1 : Complex; const Z2 : Complex):Complex;
+begin
+    Result.X := Z1.X*Z2.X-Z1.Y*Z2.Y;
+    Result.Y := Z1.X*Z2.Y+Z1.Y*Z2.X;
+end;
+
+
+function C_AddR(const Z1 : Complex; const R : Double):Complex;
+begin
+    Result.X := Z1.X+R;
+    Result.Y := Z1.Y;
+end;
+
+
+function C_MulR(const Z1 : Complex; const R : Double):Complex;
+begin
+    Result.X := Z1.X*R;
+    Result.Y := Z1.Y*R;
+end;
+
+
+function C_Sub(const Z1 : Complex; const Z2 : Complex):Complex;
+begin
+    Result.X := Z1.X-Z2.X;
+    Result.Y := Z1.Y-Z2.Y;
+end;
+
+
+function C_SubR(const Z1 : Complex; const R : Double):Complex;
+begin
+    Result.X := Z1.X-R;
+    Result.Y := Z1.Y;
+end;
+
+
+function C_RSub(const R : Double; const Z1 : Complex):Complex;
+begin
+    Result.X := R-Z1.X;
+    Result.Y := -Z1.Y;
+end;
+
+
+function C_Div(const Z1 : Complex; const Z2 : Complex):Complex;
+var
+    A : Double;
+    B : Double;
+    C : Double;
+    D : Double;
+    E : Double;
+    F : Double;
+begin
+    A := Z1.X;
+    B := Z1.Y;
+    C := Z2.X;
+    D := Z2.Y;
+    if AbsReal(D)<AbsReal(C) then
+    begin
+        E := D/C;
+        F := C+D*E;
+        Result.X := (A+B*E)/F;
+        Result.Y := (B-A*E)/F;
+    end
+    else
+    begin
+        E := C/D;
+        F := D+C*E;
+        Result.X := (B+A*E)/F;
+        Result.Y := (-A+B*E)/F;
+    end;
+end;
+
+
+function C_DivR(const Z1 : Complex; const R : Double):Complex;
+begin
+    Result.X := Z1.X/R;
+    Result.Y := Z1.Y/R;
+end;
+
+
+function C_RDiv(const R : Double; const Z2 : Complex):Complex;
+var
+    A : Double;
+    C : Double;
+    D : Double;
+    E : Double;
+    F : Double;
+begin
+    A := R;
+    C := Z2.X;
+    D := Z2.Y;
+    if AbsReal(D)<AbsReal(C) then
+    begin
+        E := D/C;
+        F := C+D*E;
+        Result.X := A/F;
+        Result.Y := -A*E/F;
+    end
+    else
+    begin
+        E := C/D;
+        F := D+C*E;
+        Result.X := A*E/F;
+        Result.Y := -A/F;
+    end;
+end;
+
+
+function C_Equal(const Z1 : Complex; const Z2 : Complex):Boolean;
+begin
+    Result := (Z1.X=Z2.X) and (Z1.Y=Z2.Y);
+end;
+
+
+function C_NotEqual(const Z1 : Complex; const Z2 : Complex):Boolean;
+begin
+    Result := (Z1.X<>Z2.X) or (Z1.Y<>Z2.Y);
+end;
+
+function C_EqualR(const Z1 : Complex; const R : Double):Boolean;
+begin
+    Result := (Z1.X=R) and (Z1.Y=0);
+end;
+
+function C_NotEqualR(const Z1 : Complex; const R : Double):Boolean;
+begin
+    Result := (Z1.X<>R) or (Z1.Y<>0);
+end;
+
+
+/////////////////////////////////////////////////////////////////////////
+// AP BLAS generic interface
+/////////////////////////////////////////////////////////////////////////
+{procedure UseAPBLAS(Flag: Boolean);
+begin
+    UseAPBLASFlag:=Flag;
+end;}
+
+function APVDotProduct(
+   V1: PDouble; I11, I12: Integer;
+   V2: PDouble; I21, I22: Integer):Double;
+var
+    I, C: LongInt;
+begin
+    Assert(I12-I11=I22-I21, 'APVDotProduct: arrays of different size!');
+    Inc(V1, I11);
+    Inc(V2, I21);
+
+    //
+    // Generic pascal code
+    //
+    C:=I12-I11;
+    Result:=0;
+    for I:=0 to C do
+    begin
+        Result:=Result+V1^*V2^;
+        Inc(V1);
+        Inc(V2);
+    end;
+end;
+
+
+procedure APVMove(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);overload;
+var
+    I, C: LongInt;
+begin
+    Assert(I12-I11=I22-I21, 'APVMove: arrays of different size!');
+    Inc(VDst, I11);
+    Inc(VSrc, I21);
+
+    //
+    // Generic pascal code
+    //
+    C:=I12-I11;
+    for I:=0 to C do
+    begin
+        VDst^:=VSrc^;
+        Inc(VDst);
+        Inc(VSrc);
+    end;
+end;
+
+
+procedure APVMove(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer;
+   S: Double);overload;
+var
+    I, C: LongInt;
+begin
+    Assert(I12-I11=I22-I21, 'APVMove: arrays of different size!');
+    Inc(VDst, I11);
+    Inc(VSrc, I21);
+
+    //
+    // Generic pascal code
+    //
+    C:=I12-I11;
+    for I:=0 to C do
+    begin
+        VDst^:=S*VSrc^;
+        Inc(VDst);
+        Inc(VSrc);
+    end;
+end;
+
+
+procedure APVMoveNeg(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);
+var
+    I, C: LongInt;
+begin
+    Assert(I12-I11=I22-I21, 'APVMoveNeg: arrays of different size!');
+    Inc(VDst, I11);
+    Inc(VSrc, I21);
+
+    //
+    // Generic pascal code
+    //
+    C:=I12-I11;
+    for I:=0 to C do
+    begin
+        VDst^:=-VSrc^;
+        Inc(VDst);
+        Inc(VSrc);
+    end;
+end;
+
+
+procedure APVAdd(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);overload;
+var
+    I, C: LongInt;
+begin
+    Assert(I12-I11=I22-I21, 'APVAdd: arrays of different size!');
+    Inc(VDst, I11);
+    Inc(VSrc, I21);
+
+    //
+    // Generic pascal code
+    //
+    C:=I12-I11;
+    for I:=0 to C do
+    begin
+        VDst^:=VDst^+VSrc^;
+        Inc(VDst);
+        Inc(VSrc);
+    end;
+end;
+
+
+procedure APVAdd(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer;
+   S: Real);overload;
+var
+    I, C: LongInt;
+begin
+    Assert(I12-I11=I22-I21, 'APVAdd: arrays of different size!');
+    Inc(VDst, I11);
+    Inc(VSrc, I21);
+
+    //
+    // Generic pascal code
+    //
+    C:=I12-I11;
+    for I:=0 to C do
+    begin
+        VDst^:=VDst^+S*VSrc^;
+        Inc(VDst);
+        Inc(VSrc);
+    end;
+end;
+
+
+procedure APVSub(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer);overload;
+var
+    I, C: LongInt;
+begin
+    Assert(I12-I11=I22-I21, 'APVSub arrays of different size!');
+    Inc(VDst, I11);
+    Inc(VSrc, I21);
+
+    //
+    // Generic pascal code
+    //
+    C:=I12-I11;
+    for I:=0 to C do
+    begin
+        VDst^:=VDst^-VSrc^;
+        Inc(VDst);
+        Inc(VSrc);
+    end;
+end;
+
+
+procedure APVSub(
+   VDst: PDouble; I11, I12: Integer;
+   VSrc: PDouble; I21, I22: Integer;
+   S: Real);overload;
+begin
+    Assert(I12-I11=I22-I21, 'APVSub: arrays of different size!');
+    APVAdd(VDst, I11, I12, VSrc, I21, I22, -S);
+end;
+
+
+procedure APVMul(
+   VOp: PDouble; I1, I2: Integer;
+   S: Real);
+var
+    I, C: LongInt;
+begin
+    Inc(VOp, I1);
+    C:=I2-I1;
+    for I:=0 to C do
+    begin
+        VOp^:=S*VOp^;
+        Inc(VOp);
+    end;
+end;
+
+/////////////////////////////////////////////////////////////////////////
+// IEEE-compliant functions
+/////////////////////////////////////////////////////////////////////////
+{$OPTIMIZATION OFF}
+function AP_Double(X:Double):Double;
+begin
+    Result:=X;
+end;
+
+function AP_FP_Eq(X:Double; Y:Double):Boolean;
+begin
+    Result:=X=Y;
+end;
+
+function AP_FP_NEq(X:Double; Y:Double):Boolean;
+begin
+    Result:=X<>Y;
+end;
+
+function AP_FP_Less(X:Double; Y:Double):Boolean;
+begin
+    Result:=X<Y;
+end;
+
+function AP_FP_Less_Eq(X:Double; Y:Double):Boolean;
+begin
+    Result:=X<=Y;
+end;
+
+function AP_FP_Greater(X:Double; Y:Double):Boolean;
+begin
+    Result:=X>Y;
+end;
+
+function AP_FP_Greater_Eq(X:Double; Y:Double):Boolean;
+begin
+    Result:=X>=Y;
+end;
+
+end.

tests/test/alglib/u_conv.pp

@@ -0,0 +1,1776 @@
+(*************************************************************************
+Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************)
+unit u_conv;
+interface
+uses Math, Sysutils, u_ap, u_ftbase, u_fft;
+
+procedure ConvC1D(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+procedure ConvC1DInv(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+procedure ConvC1DCircular(const S : TComplex1DArray;
+     M : Integer;
+     const R : TComplex1DArray;
+     N : Integer;
+     var C : TComplex1DArray);
+procedure ConvC1DCircularInv(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+procedure ConvR1D(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+procedure ConvR1DInv(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+procedure ConvR1DCircular(const S : TReal1DArray;
+     M : Integer;
+     const R : TReal1DArray;
+     N : Integer;
+     var C : TReal1DArray);
+procedure ConvR1DCircularInv(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+procedure ConvC1DX(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     Circular : Boolean;
+     Alg : Integer;
+     Q : Integer;
+     var R : TComplex1DArray);
+procedure ConvR1DX(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     Circular : Boolean;
+     Alg : Integer;
+     Q : Integer;
+     var R : TReal1DArray);
+
+implementation
+
+(*************************************************************************
+1-dimensional complex convolution.
+
+For given A/B returns conv(A,B) (non-circular). Subroutine can automatically
+choose between three implementations: straightforward O(M*N)  formula  for
+very small N (or M), overlap-add algorithm for  cases  where  max(M,N)  is
+significantly larger than min(M,N), but O(M*N) algorithm is too slow,  and
+general FFT-based formula for cases where two previois algorithms are  too
+slow.
+
+Algorithm has max(M,N)*log(max(M,N)) complexity for any M/N.
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - complex function to be transformed
+    M   -   problem size
+    B   -   array[0..N-1] - complex function to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..N+M-2].
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvC1D(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+begin
+    Assert((N>0) and (M>0), 'ConvC1D: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer that B.
+    //
+    if M<N then
+    begin
+        ConvC1D(B, N, A, M, R);
+        Exit;
+    end;
+    ConvC1DX(A, M, B, N, False, -1, 0, R);
+end;
+
+
+(*************************************************************************
+1-dimensional complex non-circular deconvolution (inverse of ConvC1D()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - response
+    N   -   response length, N<=M
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-N].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvC1DInv(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+var
+    I : Integer;
+    P : Integer;
+    Buf : TReal1DArray;
+    Buf2 : TReal1DArray;
+    Plan : FTPlan;
+    C1 : Complex;
+    C2 : Complex;
+    C3 : Complex;
+    T : Double;
+begin
+    Assert((N>0) and (M>0) and (N<=M), 'ConvC1DInv: incorrect N or M!');
+    P := FTBaseFindSmooth(M);
+    FTBaseGenerateComplexFFTPlan(P, Plan);
+    SetLength(Buf, 2*P);
+    I:=0;
+    while I<=M-1 do
+    begin
+        Buf[2*I+0] := A[I].X;
+        Buf[2*I+1] := A[I].Y;
+        Inc(I);
+    end;
+    I:=M;
+    while I<=P-1 do
+    begin
+        Buf[2*I+0] := 0;
+        Buf[2*I+1] := 0;
+        Inc(I);
+    end;
+    SetLength(Buf2, 2*P);
+    I:=0;
+    while I<=N-1 do
+    begin
+        Buf2[2*I+0] := B[I].X;
+        Buf2[2*I+1] := B[I].Y;
+        Inc(I);
+    end;
+    I:=N;
+    while I<=P-1 do
+    begin
+        Buf2[2*I+0] := 0;
+        Buf2[2*I+1] := 0;
+        Inc(I);
+    end;
+    FTBaseExecutePlan(Buf, 0, P, Plan);
+    FTBaseExecutePlan(Buf2, 0, P, Plan);
+    I:=0;
+    while I<=P-1 do
+    begin
+        C1.X := Buf[2*I+0];
+        C1.Y := Buf[2*I+1];
+        C2.X := Buf2[2*I+0];
+        C2.Y := Buf2[2*I+1];
+        C3 := C_Div(C1,C2);
+        Buf[2*I+0] := C3.X;
+        Buf[2*I+1] := -C3.Y;
+        Inc(I);
+    end;
+    FTBaseExecutePlan(Buf, 0, P, Plan);
+    T := AP_Double(1)/P;
+    SetLength(R, M-N+1);
+    I:=0;
+    while I<=M-N do
+    begin
+        R[I].X := +T*Buf[2*I+0];
+        R[I].Y := -T*Buf[2*I+1];
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional circular complex convolution.
+
+For given S/R returns conv(S,R) (circular). Algorithm has linearithmic
+complexity for any M/N.
+
+IMPORTANT:  normal convolution is commutative,  i.e.   it  is symmetric  -
+conv(A,B)=conv(B,A).  Cyclic convolution IS NOT.  One function - S - is  a
+signal,  periodic function, and another - R - is a response,  non-periodic
+function with limited length.
+
+INPUT PARAMETERS
+    S   -   array[0..M-1] - complex periodic signal
+    M   -   problem size
+    B   -   array[0..N-1] - complex non-periodic response
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..M-1].
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvC1DCircular(const S : TComplex1DArray;
+     M : Integer;
+     const R : TComplex1DArray;
+     N : Integer;
+     var C : TComplex1DArray);
+var
+    Buf : TComplex1DArray;
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    Assert((N>0) and (M>0), 'ConvC1DCircular: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer (at least - not shorter) that B.
+    //
+    if M<N then
+    begin
+        SetLength(Buf, M);
+        I1:=0;
+        while I1<=M-1 do
+        begin
+            Buf[I1] := C_Complex(0);
+            Inc(I1);
+        end;
+        I1 := 0;
+        while I1<N do
+        begin
+            I2 := Min(I1+M-1, N-1);
+            J2 := I2-I1;
+            i1_ := (I1) - (0);
+            for i_ := 0 to J2 do
+            begin
+                Buf[i_] := C_Add(Buf[i_], R[i_+i1_]);
+            end;
+            I1 := I1+M;
+        end;
+        ConvC1DCircular(S, M, Buf, M, C);
+        Exit;
+    end;
+    ConvC1DX(S, M, R, N, True, -1, 0, C);
+end;
+
+
+(*************************************************************************
+1-dimensional circular complex deconvolution (inverse of ConvC1DCircular()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved periodic signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - non-periodic response
+    N   -   response length
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-1].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvC1DCircularInv(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+var
+    I : Integer;
+    P : Integer;
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+    Buf : TReal1DArray;
+    Buf2 : TReal1DArray;
+    CBuf : TComplex1DArray;
+    Plan : FTPlan;
+    C1 : Complex;
+    C2 : Complex;
+    C3 : Complex;
+    T : Double;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    Assert((N>0) and (M>0), 'ConvC1DCircularInv: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer (at least - not shorter) that B.
+    //
+    if M<N then
+    begin
+        SetLength(CBuf, M);
+        I:=0;
+        while I<=M-1 do
+        begin
+            CBuf[I] := C_Complex(0);
+            Inc(I);
+        end;
+        I1 := 0;
+        while I1<N do
+        begin
+            I2 := Min(I1+M-1, N-1);
+            J2 := I2-I1;
+            i1_ := (I1) - (0);
+            for i_ := 0 to J2 do
+            begin
+                CBuf[i_] := C_Add(CBuf[i_], B[i_+i1_]);
+            end;
+            I1 := I1+M;
+        end;
+        ConvC1DCircularInv(A, M, CBuf, M, R);
+        Exit;
+    end;
+    
+    //
+    // Task is normalized
+    //
+    FTBaseGenerateComplexFFTPlan(M, Plan);
+    SetLength(Buf, 2*M);
+    I:=0;
+    while I<=M-1 do
+    begin
+        Buf[2*I+0] := A[I].X;
+        Buf[2*I+1] := A[I].Y;
+        Inc(I);
+    end;
+    SetLength(Buf2, 2*M);
+    I:=0;
+    while I<=N-1 do
+    begin
+        Buf2[2*I+0] := B[I].X;
+        Buf2[2*I+1] := B[I].Y;
+        Inc(I);
+    end;
+    I:=N;
+    while I<=M-1 do
+    begin
+        Buf2[2*I+0] := 0;
+        Buf2[2*I+1] := 0;
+        Inc(I);
+    end;
+    FTBaseExecutePlan(Buf, 0, M, Plan);
+    FTBaseExecutePlan(Buf2, 0, M, Plan);
+    I:=0;
+    while I<=M-1 do
+    begin
+        C1.X := Buf[2*I+0];
+        C1.Y := Buf[2*I+1];
+        C2.X := Buf2[2*I+0];
+        C2.Y := Buf2[2*I+1];
+        C3 := C_Div(C1,C2);
+        Buf[2*I+0] := C3.X;
+        Buf[2*I+1] := -C3.Y;
+        Inc(I);
+    end;
+    FTBaseExecutePlan(Buf, 0, M, Plan);
+    T := AP_Double(1)/M;
+    SetLength(R, M);
+    I:=0;
+    while I<=M-1 do
+    begin
+        R[I].X := +T*Buf[2*I+0];
+        R[I].Y := -T*Buf[2*I+1];
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional real convolution.
+
+Analogous to ConvC1D(), see ConvC1D() comments for more details.
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - real function to be transformed
+    M   -   problem size
+    B   -   array[0..N-1] - real function to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..N+M-2].
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvR1D(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+var
+    I : Integer;
+    J : Integer;
+    P : Integer;
+    Q : Integer;
+    ABuf : TComplex1DArray;
+    BBuf : TComplex1DArray;
+    V : Complex;
+    Flop1 : Double;
+    Flop2 : Double;
+begin
+    Assert((N>0) and (M>0), 'ConvR1D: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer that B.
+    //
+    if M<N then
+    begin
+        ConvR1D(B, N, A, M, R);
+        Exit;
+    end;
+    ConvR1DX(A, M, B, N, False, -1, 0, R);
+end;
+
+
+(*************************************************************************
+1-dimensional real deconvolution (inverse of ConvC1D()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - response
+    N   -   response length, N<=M
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-N].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvR1DInv(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+var
+    I : Integer;
+    P : Integer;
+    Buf : TReal1DArray;
+    Buf2 : TReal1DArray;
+    Buf3 : TReal1DArray;
+    Plan : FTPlan;
+    C1 : Complex;
+    C2 : Complex;
+    C3 : Complex;
+    T : Double;
+begin
+    Assert((N>0) and (M>0) and (N<=M), 'ConvR1DInv: incorrect N or M!');
+    P := FTBaseFindSmoothEven(M);
+    SetLength(Buf, P);
+    APVMove(@Buf[0], 0, M-1, @A[0], 0, M-1);
+    I:=M;
+    while I<=P-1 do
+    begin
+        Buf[I] := 0;
+        Inc(I);
+    end;
+    SetLength(Buf2, P);
+    APVMove(@Buf2[0], 0, N-1, @B[0], 0, N-1);
+    I:=N;
+    while I<=P-1 do
+    begin
+        Buf2[I] := 0;
+        Inc(I);
+    end;
+    SetLength(Buf3, P);
+    FTBaseGenerateComplexFFTPlan(P div 2, Plan);
+    FFTR1DInternalEven(Buf, P, Buf3, Plan);
+    FFTR1DInternalEven(Buf2, P, Buf3, Plan);
+    Buf[0] := Buf[0]/Buf2[0];
+    Buf[1] := Buf[1]/Buf2[1];
+    I:=1;
+    while I<=P div 2-1 do
+    begin
+        C1.X := Buf[2*I+0];
+        C1.Y := Buf[2*I+1];
+        C2.X := Buf2[2*I+0];
+        C2.Y := Buf2[2*I+1];
+        C3 := C_Div(C1,C2);
+        Buf[2*I+0] := C3.X;
+        Buf[2*I+1] := C3.Y;
+        Inc(I);
+    end;
+    FFTR1DInvInternalEven(Buf, P, Buf3, Plan);
+    SetLength(R, M-N+1);
+    APVMove(@R[0], 0, M-N, @Buf[0], 0, M-N);
+end;
+
+
+(*************************************************************************
+1-dimensional circular real convolution.
+
+Analogous to ConvC1DCircular(), see ConvC1DCircular() comments for more details.
+
+INPUT PARAMETERS
+    S   -   array[0..M-1] - real signal
+    M   -   problem size
+    B   -   array[0..N-1] - real response
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..M-1].
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvR1DCircular(const S : TReal1DArray;
+     M : Integer;
+     const R : TReal1DArray;
+     N : Integer;
+     var C : TReal1DArray);
+var
+    Buf : TReal1DArray;
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+begin
+    Assert((N>0) and (M>0), 'ConvC1DCircular: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer (at least - not shorter) that B.
+    //
+    if M<N then
+    begin
+        SetLength(Buf, M);
+        I1:=0;
+        while I1<=M-1 do
+        begin
+            Buf[I1] := 0;
+            Inc(I1);
+        end;
+        I1 := 0;
+        while I1<N do
+        begin
+            I2 := Min(I1+M-1, N-1);
+            J2 := I2-I1;
+            APVAdd(@Buf[0], 0, J2, @R[0], I1, I2);
+            I1 := I1+M;
+        end;
+        ConvR1DCircular(S, M, Buf, M, C);
+        Exit;
+    end;
+    
+    //
+    // reduce to usual convolution
+    //
+    ConvR1DX(S, M, R, N, True, -1, 0, C);
+end;
+
+
+(*************************************************************************
+1-dimensional complex deconvolution (inverse of ConvC1D()).
+
+Algorithm has M*log(M)) complexity for any M (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - convolved signal, A = conv(R, B)
+    M   -   convolved signal length
+    B   -   array[0..N-1] - response
+    N   -   response length
+
+OUTPUT PARAMETERS
+    R   -   deconvolved signal. array[0..M-N].
+
+NOTE:
+    deconvolution is unstable process and may result in division  by  zero
+(if your response function is degenerate, i.e. has zero Fourier coefficient).
+
+NOTE:
+    It is assumed that A is zero at T<0, B is zero too.  If  one  or  both
+functions have non-zero values at negative T's, you  can  still  use  this
+subroutine - just shift its result correspondingly.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvR1DCircularInv(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+var
+    I : Integer;
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+    Buf : TReal1DArray;
+    Buf2 : TReal1DArray;
+    Buf3 : TReal1DArray;
+    CBuf : TComplex1DArray;
+    CBuf2 : TComplex1DArray;
+    Plan : FTPlan;
+    C1 : Complex;
+    C2 : Complex;
+    C3 : Complex;
+    T : Double;
+begin
+    Assert((N>0) and (M>0), 'ConvR1DCircularInv: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer (at least - not shorter) that B.
+    //
+    if M<N then
+    begin
+        SetLength(Buf, M);
+        I:=0;
+        while I<=M-1 do
+        begin
+            Buf[I] := 0;
+            Inc(I);
+        end;
+        I1 := 0;
+        while I1<N do
+        begin
+            I2 := Min(I1+M-1, N-1);
+            J2 := I2-I1;
+            APVAdd(@Buf[0], 0, J2, @B[0], I1, I2);
+            I1 := I1+M;
+        end;
+        ConvR1DCircularInv(A, M, Buf, M, R);
+        Exit;
+    end;
+    
+    //
+    // Task is normalized
+    //
+    if M mod 2=0 then
+    begin
+        
+        //
+        // size is even, use fast even-size FFT
+        //
+        SetLength(Buf, M);
+        APVMove(@Buf[0], 0, M-1, @A[0], 0, M-1);
+        SetLength(Buf2, M);
+        APVMove(@Buf2[0], 0, N-1, @B[0], 0, N-1);
+        I:=N;
+        while I<=M-1 do
+        begin
+            Buf2[I] := 0;
+            Inc(I);
+        end;
+        SetLength(Buf3, M);
+        FTBaseGenerateComplexFFTPlan(M div 2, Plan);
+        FFTR1DInternalEven(Buf, M, Buf3, Plan);
+        FFTR1DInternalEven(Buf2, M, Buf3, Plan);
+        Buf[0] := Buf[0]/Buf2[0];
+        Buf[1] := Buf[1]/Buf2[1];
+        I:=1;
+        while I<=M div 2-1 do
+        begin
+            C1.X := Buf[2*I+0];
+            C1.Y := Buf[2*I+1];
+            C2.X := Buf2[2*I+0];
+            C2.Y := Buf2[2*I+1];
+            C3 := C_Div(C1,C2);
+            Buf[2*I+0] := C3.X;
+            Buf[2*I+1] := C3.Y;
+            Inc(I);
+        end;
+        FFTR1DInvInternalEven(Buf, M, Buf3, Plan);
+        SetLength(R, M);
+        APVMove(@R[0], 0, M-1, @Buf[0], 0, M-1);
+    end
+    else
+    begin
+        
+        //
+        // odd-size, use general real FFT
+        //
+        FFTR1D(A, M, CBuf);
+        SetLength(Buf2, M);
+        APVMove(@Buf2[0], 0, N-1, @B[0], 0, N-1);
+        I:=N;
+        while I<=M-1 do
+        begin
+            Buf2[I] := 0;
+            Inc(I);
+        end;
+        FFTR1D(Buf2, M, CBuf2);
+        I:=0;
+        while I<=Floor(AP_Double(M)/2) do
+        begin
+            CBuf[I] := C_Div(CBuf[I],CBuf2[I]);
+            Inc(I);
+        end;
+        FFTR1DInv(CBuf, M, R);
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional complex convolution.
+
+Extended subroutine which allows to choose convolution algorithm.
+Intended for internal use, ALGLIB users should call ConvC1D()/ConvC1DCircular().
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - complex function to be transformed
+    M   -   problem size
+    B   -   array[0..N-1] - complex function to be transformed
+    N   -   problem size, N<=M
+    Alg -   algorithm type:
+            *-2     auto-select Q for overlap-add
+            *-1     auto-select algorithm and parameters
+            * 0     straightforward formula for small N's
+            * 1     general FFT-based code
+            * 2     overlap-add with length Q
+    Q   -   length for overlap-add
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..N+M-1].
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvC1DX(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     Circular : Boolean;
+     Alg : Integer;
+     Q : Integer;
+     var R : TComplex1DArray);
+var
+    I : Integer;
+    J : Integer;
+    P : Integer;
+    PTotal : Integer;
+    I1 : Integer;
+    I2 : Integer;
+    J1 : Integer;
+    J2 : Integer;
+    BBuf : TComplex1DArray;
+    V : Complex;
+    AX : Double;
+    AY : Double;
+    BX : Double;
+    BY : Double;
+    T : Double;
+    TX : Double;
+    TY : Double;
+    FlopCand : Double;
+    FlopBest : Double;
+    AlgBest : Integer;
+    Plan : FTPlan;
+    Buf : TReal1DArray;
+    Buf2 : TReal1DArray;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    Assert((N>0) and (M>0), 'ConvC1DX: incorrect N or M!');
+    Assert(N<=M, 'ConvC1DX: N<M assumption is false!');
+    
+    //
+    // Auto-select
+    //
+    if (Alg=-1) or (Alg=-2) then
+    begin
+        
+        //
+        // Initial candidate: straightforward implementation.
+        //
+        // If we want to use auto-fitted overlap-add,
+        // flop count is initialized by large real number - to force
+        // another algorithm selection
+        //
+        AlgBest := 0;
+        if Alg=-1 then
+        begin
+            FlopBest := 2*M*N;
+        end
+        else
+        begin
+            FlopBest := MaxRealNumber;
+        end;
+        
+        //
+        // Another candidate - generic FFT code
+        //
+        if Alg=-1 then
+        begin
+            if Circular and FTBaseIsSmooth(M) then
+            begin
+                
+                //
+                // special code for circular convolution of a sequence with a smooth length
+                //
+                FlopCand := 3*FTBaseGetFLOPEstimate(M)+6*M;
+                if AP_FP_Less(FlopCand,FlopBest) then
+                begin
+                    AlgBest := 1;
+                    FlopBest := FlopCand;
+                end;
+            end
+            else
+            begin
+                
+                //
+                // general cyclic/non-cyclic convolution
+                //
+                P := FTBaseFindSmooth(M+N-1);
+                FlopCand := 3*FTBaseGetFLOPEstimate(P)+6*P;
+                if AP_FP_Less(FlopCand,FlopBest) then
+                begin
+                    AlgBest := 1;
+                    FlopBest := FlopCand;
+                end;
+            end;
+        end;
+        
+        //
+        // Another candidate - overlap-add
+        //
+        Q := 1;
+        PTotal := 1;
+        while PTotal<N do
+        begin
+            PTotal := PTotal*2;
+        end;
+        while PTotal<=M+N-1 do
+        begin
+            P := PTotal-N+1;
+            FlopCand := Ceil(AP_Double(M)/P)*(2*FTBaseGetFLOPEstimate(PTotal)+8*PTotal);
+            if AP_FP_Less(FlopCand,FlopBest) then
+            begin
+                FlopBest := FlopCand;
+                AlgBest := 2;
+                Q := P;
+            end;
+            PTotal := PTotal*2;
+        end;
+        Alg := AlgBest;
+        ConvC1DX(A, M, B, N, Circular, Alg, Q, R);
+        Exit;
+    end;
+    
+    //
+    // straightforward formula for
+    // circular and non-circular convolutions.
+    //
+    // Very simple code, no further comments needed.
+    //
+    if Alg=0 then
+    begin
+        
+        //
+        // Special case: N=1
+        //
+        if N=1 then
+        begin
+            SetLength(R, M);
+            V := B[0];
+            for i_ := 0 to M-1 do
+            begin
+                R[i_] := C_Mul(V, A[i_]);
+            end;
+            Exit;
+        end;
+        
+        //
+        // use straightforward formula
+        //
+        if Circular then
+        begin
+            
+            //
+            // circular convolution
+            //
+            SetLength(R, M);
+            V := B[0];
+            for i_ := 0 to M-1 do
+            begin
+                R[i_] := C_Mul(V, A[i_]);
+            end;
+            I:=1;
+            while I<=N-1 do
+            begin
+                V := B[I];
+                I1 := 0;
+                I2 := I-1;
+                J1 := M-I;
+                J2 := M-1;
+                i1_ := (J1) - (I1);
+                for i_ := I1 to I2 do
+                begin
+                    R[i_] := C_Add(R[i_], C_Mul(V, A[i_+i1_]));
+                end;
+                I1 := I;
+                I2 := M-1;
+                J1 := 0;
+                J2 := M-I-1;
+                i1_ := (J1) - (I1);
+                for i_ := I1 to I2 do
+                begin
+                    R[i_] := C_Add(R[i_], C_Mul(V, A[i_+i1_]));
+                end;
+                Inc(I);
+            end;
+        end
+        else
+        begin
+            
+            //
+            // non-circular convolution
+            //
+            SetLength(R, M+N-1);
+            I:=0;
+            while I<=M+N-2 do
+            begin
+                R[I] := C_Complex(0);
+                Inc(I);
+            end;
+            I:=0;
+            while I<=N-1 do
+            begin
+                V := B[I];
+                i1_ := (0) - (I);
+                for i_ := I to I+M-1 do
+                begin
+                    R[i_] := C_Add(R[i_], C_Mul(V, A[i_+i1_]));
+                end;
+                Inc(I);
+            end;
+        end;
+        Exit;
+    end;
+    
+    //
+    // general FFT-based code for
+    // circular and non-circular convolutions.
+    //
+    // First, if convolution is circular, we test whether M is smooth or not.
+    // If it is smooth, we just use M-length FFT to calculate convolution.
+    // If it is not, we calculate non-circular convolution and wrap it arount.
+    //
+    // IF convolution is non-circular, we use zero-padding + FFT.
+    //
+    if Alg=1 then
+    begin
+        if Circular and FTBaseIsSmooth(M) then
+        begin
+            
+            //
+            // special code for circular convolution with smooth M
+            //
+            FTBaseGenerateComplexFFTPlan(M, Plan);
+            SetLength(Buf, 2*M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                Buf[2*I+0] := A[I].X;
+                Buf[2*I+1] := A[I].Y;
+                Inc(I);
+            end;
+            SetLength(Buf2, 2*M);
+            I:=0;
+            while I<=N-1 do
+            begin
+                Buf2[2*I+0] := B[I].X;
+                Buf2[2*I+1] := B[I].Y;
+                Inc(I);
+            end;
+            I:=N;
+            while I<=M-1 do
+            begin
+                Buf2[2*I+0] := 0;
+                Buf2[2*I+1] := 0;
+                Inc(I);
+            end;
+            FTBaseExecutePlan(Buf, 0, M, Plan);
+            FTBaseExecutePlan(Buf2, 0, M, Plan);
+            I:=0;
+            while I<=M-1 do
+            begin
+                AX := Buf[2*I+0];
+                AY := Buf[2*I+1];
+                BX := Buf2[2*I+0];
+                BY := Buf2[2*I+1];
+                TX := AX*BX-AY*BY;
+                TY := AX*BY+AY*BX;
+                Buf[2*I+0] := TX;
+                Buf[2*I+1] := -TY;
+                Inc(I);
+            end;
+            FTBaseExecutePlan(Buf, 0, M, Plan);
+            T := AP_Double(1)/M;
+            SetLength(R, M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                R[I].X := +T*Buf[2*I+0];
+                R[I].Y := -T*Buf[2*I+1];
+                Inc(I);
+            end;
+        end
+        else
+        begin
+            
+            //
+            // M is non-smooth, general code (circular/non-circular):
+            // * first part is the same for circular and non-circular
+            //   convolutions. zero padding, FFTs, inverse FFTs
+            // * second part differs:
+            //   * for non-circular convolution we just copy array
+            //   * for circular convolution we add array tail to its head
+            //
+            P := FTBaseFindSmooth(M+N-1);
+            FTBaseGenerateComplexFFTPlan(P, Plan);
+            SetLength(Buf, 2*P);
+            I:=0;
+            while I<=M-1 do
+            begin
+                Buf[2*I+0] := A[I].X;
+                Buf[2*I+1] := A[I].Y;
+                Inc(I);
+            end;
+            I:=M;
+            while I<=P-1 do
+            begin
+                Buf[2*I+0] := 0;
+                Buf[2*I+1] := 0;
+                Inc(I);
+            end;
+            SetLength(Buf2, 2*P);
+            I:=0;
+            while I<=N-1 do
+            begin
+                Buf2[2*I+0] := B[I].X;
+                Buf2[2*I+1] := B[I].Y;
+                Inc(I);
+            end;
+            I:=N;
+            while I<=P-1 do
+            begin
+                Buf2[2*I+0] := 0;
+                Buf2[2*I+1] := 0;
+                Inc(I);
+            end;
+            FTBaseExecutePlan(Buf, 0, P, Plan);
+            FTBaseExecutePlan(Buf2, 0, P, Plan);
+            I:=0;
+            while I<=P-1 do
+            begin
+                AX := Buf[2*I+0];
+                AY := Buf[2*I+1];
+                BX := Buf2[2*I+0];
+                BY := Buf2[2*I+1];
+                TX := AX*BX-AY*BY;
+                TY := AX*BY+AY*BX;
+                Buf[2*I+0] := TX;
+                Buf[2*I+1] := -TY;
+                Inc(I);
+            end;
+            FTBaseExecutePlan(Buf, 0, P, Plan);
+            T := AP_Double(1)/P;
+            if Circular then
+            begin
+                
+                //
+                // circular, add tail to head
+                //
+                SetLength(R, M);
+                I:=0;
+                while I<=M-1 do
+                begin
+                    R[I].X := +T*Buf[2*I+0];
+                    R[I].Y := -T*Buf[2*I+1];
+                    Inc(I);
+                end;
+                I:=M;
+                while I<=M+N-2 do
+                begin
+                    R[I-M].X := R[I-M].X+T*Buf[2*I+0];
+                    R[I-M].Y := R[I-M].Y-T*Buf[2*I+1];
+                    Inc(I);
+                end;
+            end
+            else
+            begin
+                
+                //
+                // non-circular, just copy
+                //
+                SetLength(R, M+N-1);
+                I:=0;
+                while I<=M+N-2 do
+                begin
+                    R[I].X := +T*Buf[2*I+0];
+                    R[I].Y := -T*Buf[2*I+1];
+                    Inc(I);
+                end;
+            end;
+        end;
+        Exit;
+    end;
+    
+    //
+    // overlap-add method for
+    // circular and non-circular convolutions.
+    //
+    // First part of code (separate FFTs of input blocks) is the same
+    // for all types of convolution. Second part (overlapping outputs)
+    // differs for different types of convolution. We just copy output
+    // when convolution is non-circular. We wrap it around, if it is
+    // circular.
+    //
+    if Alg=2 then
+    begin
+        SetLength(Buf, 2*(Q+N-1));
+        
+        //
+        // prepare R
+        //
+        if Circular then
+        begin
+            SetLength(R, M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                R[I] := C_Complex(0);
+                Inc(I);
+            end;
+        end
+        else
+        begin
+            SetLength(R, M+N-1);
+            I:=0;
+            while I<=M+N-2 do
+            begin
+                R[I] := C_Complex(0);
+                Inc(I);
+            end;
+        end;
+        
+        //
+        // pre-calculated FFT(B)
+        //
+        SetLength(BBuf, Q+N-1);
+        for i_ := 0 to N-1 do
+        begin
+            BBuf[i_] := B[i_];
+        end;
+        J:=N;
+        while J<=Q+N-2 do
+        begin
+            BBuf[J] := C_Complex(0);
+            Inc(J);
+        end;
+        FFTC1D(BBuf, Q+N-1);
+        
+        //
+        // prepare FFT plan for chunks of A
+        //
+        FTBaseGenerateComplexFFTPlan(Q+N-1, Plan);
+        
+        //
+        // main overlap-add cycle
+        //
+        I := 0;
+        while I<=M-1 do
+        begin
+            P := Min(Q, M-I);
+            J:=0;
+            while J<=P-1 do
+            begin
+                Buf[2*J+0] := A[I+J].X;
+                Buf[2*J+1] := A[I+J].Y;
+                Inc(J);
+            end;
+            J:=P;
+            while J<=Q+N-2 do
+            begin
+                Buf[2*J+0] := 0;
+                Buf[2*J+1] := 0;
+                Inc(J);
+            end;
+            FTBaseExecutePlan(Buf, 0, Q+N-1, Plan);
+            J:=0;
+            while J<=Q+N-2 do
+            begin
+                AX := Buf[2*J+0];
+                AY := Buf[2*J+1];
+                BX := BBuf[J].X;
+                BY := BBuf[J].Y;
+                TX := AX*BX-AY*BY;
+                TY := AX*BY+AY*BX;
+                Buf[2*J+0] := TX;
+                Buf[2*J+1] := -TY;
+                Inc(J);
+            end;
+            FTBaseExecutePlan(Buf, 0, Q+N-1, Plan);
+            T := AP_Double(1)/(Q+N-1);
+            if Circular then
+            begin
+                J1 := Min(I+P+N-2, M-1)-I;
+                J2 := J1+1;
+            end
+            else
+            begin
+                J1 := P+N-2;
+                J2 := J1+1;
+            end;
+            J:=0;
+            while J<=J1 do
+            begin
+                R[I+J].X := R[I+J].X+Buf[2*J+0]*T;
+                R[I+J].Y := R[I+J].Y-Buf[2*J+1]*T;
+                Inc(J);
+            end;
+            J:=J2;
+            while J<=P+N-2 do
+            begin
+                R[J-J2].X := R[J-J2].X+Buf[2*J+0]*T;
+                R[J-J2].Y := R[J-J2].Y-Buf[2*J+1]*T;
+                Inc(J);
+            end;
+            I := I+P;
+        end;
+        Exit;
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional real convolution.
+
+Extended subroutine which allows to choose convolution algorithm.
+Intended for internal use, ALGLIB users should call ConvR1D().
+
+INPUT PARAMETERS
+    A   -   array[0..M-1] - complex function to be transformed
+    M   -   problem size
+    B   -   array[0..N-1] - complex function to be transformed
+    N   -   problem size, N<=M
+    Alg -   algorithm type:
+            *-2     auto-select Q for overlap-add
+            *-1     auto-select algorithm and parameters
+            * 0     straightforward formula for small N's
+            * 1     general FFT-based code
+            * 2     overlap-add with length Q
+    Q   -   length for overlap-add
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..N+M-1].
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure ConvR1DX(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     Circular : Boolean;
+     Alg : Integer;
+     Q : Integer;
+     var R : TReal1DArray);
+var
+    V : Double;
+    I : Integer;
+    J : Integer;
+    P : Integer;
+    PTotal : Integer;
+    I1 : Integer;
+    I2 : Integer;
+    J1 : Integer;
+    J2 : Integer;
+    AX : Double;
+    AY : Double;
+    BX : Double;
+    BY : Double;
+    T : Double;
+    TX : Double;
+    TY : Double;
+    FlopCand : Double;
+    FlopBest : Double;
+    AlgBest : Integer;
+    Plan : FTPlan;
+    Buf : TReal1DArray;
+    Buf2 : TReal1DArray;
+    Buf3 : TReal1DArray;
+begin
+    Assert((N>0) and (M>0), 'ConvC1DX: incorrect N or M!');
+    Assert(N<=M, 'ConvC1DX: N<M assumption is false!');
+    
+    //
+    // handle special cases
+    //
+    if Min(M, N)<=2 then
+    begin
+        Alg := 0;
+    end;
+    
+    //
+    // Auto-select
+    //
+    if Alg<0 then
+    begin
+        
+        //
+        // Initial candidate: straightforward implementation.
+        //
+        // If we want to use auto-fitted overlap-add,
+        // flop count is initialized by large real number - to force
+        // another algorithm selection
+        //
+        AlgBest := 0;
+        if Alg=-1 then
+        begin
+            FlopBest := Double(0.15)*M*N;
+        end
+        else
+        begin
+            FlopBest := MaxRealNumber;
+        end;
+        
+        //
+        // Another candidate - generic FFT code
+        //
+        if Alg=-1 then
+        begin
+            if Circular and FTBaseIsSmooth(M) and (M mod 2=0) then
+            begin
+                
+                //
+                // special code for circular convolution of a sequence with a smooth length
+                //
+                FlopCand := 3*FTBaseGetFLOPEstimate(M div 2)+AP_Double(6*M)/2;
+                if AP_FP_Less(FlopCand,FlopBest) then
+                begin
+                    AlgBest := 1;
+                    FlopBest := FlopCand;
+                end;
+            end
+            else
+            begin
+                
+                //
+                // general cyclic/non-cyclic convolution
+                //
+                P := FTBaseFindSmoothEven(M+N-1);
+                FlopCand := 3*FTBaseGetFLOPEstimate(P div 2)+AP_Double(6*P)/2;
+                if AP_FP_Less(FlopCand,FlopBest) then
+                begin
+                    AlgBest := 1;
+                    FlopBest := FlopCand;
+                end;
+            end;
+        end;
+        
+        //
+        // Another candidate - overlap-add
+        //
+        Q := 1;
+        PTotal := 1;
+        while PTotal<N do
+        begin
+            PTotal := PTotal*2;
+        end;
+        while PTotal<=M+N-1 do
+        begin
+            P := PTotal-N+1;
+            FlopCand := Ceil(AP_Double(M)/P)*(2*FTBaseGetFLOPEstimate(PTotal div 2)+1*(PTotal div 2));
+            if AP_FP_Less(FlopCand,FlopBest) then
+            begin
+                FlopBest := FlopCand;
+                AlgBest := 2;
+                Q := P;
+            end;
+            PTotal := PTotal*2;
+        end;
+        Alg := AlgBest;
+        ConvR1DX(A, M, B, N, Circular, Alg, Q, R);
+        Exit;
+    end;
+    
+    //
+    // straightforward formula for
+    // circular and non-circular convolutions.
+    //
+    // Very simple code, no further comments needed.
+    //
+    if Alg=0 then
+    begin
+        
+        //
+        // Special case: N=1
+        //
+        if N=1 then
+        begin
+            SetLength(R, M);
+            V := B[0];
+            APVMove(@R[0], 0, M-1, @A[0], 0, M-1, V);
+            Exit;
+        end;
+        
+        //
+        // use straightforward formula
+        //
+        if Circular then
+        begin
+            
+            //
+            // circular convolution
+            //
+            SetLength(R, M);
+            V := B[0];
+            APVMove(@R[0], 0, M-1, @A[0], 0, M-1, V);
+            I:=1;
+            while I<=N-1 do
+            begin
+                V := B[I];
+                I1 := 0;
+                I2 := I-1;
+                J1 := M-I;
+                J2 := M-1;
+                APVAdd(@R[0], I1, I2, @A[0], J1, J2, V);
+                I1 := I;
+                I2 := M-1;
+                J1 := 0;
+                J2 := M-I-1;
+                APVAdd(@R[0], I1, I2, @A[0], J1, J2, V);
+                Inc(I);
+            end;
+        end
+        else
+        begin
+            
+            //
+            // non-circular convolution
+            //
+            SetLength(R, M+N-1);
+            I:=0;
+            while I<=M+N-2 do
+            begin
+                R[I] := 0;
+                Inc(I);
+            end;
+            I:=0;
+            while I<=N-1 do
+            begin
+                V := B[I];
+                APVAdd(@R[0], I, I+M-1, @A[0], 0, M-1, V);
+                Inc(I);
+            end;
+        end;
+        Exit;
+    end;
+    
+    //
+    // general FFT-based code for
+    // circular and non-circular convolutions.
+    //
+    // First, if convolution is circular, we test whether M is smooth or not.
+    // If it is smooth, we just use M-length FFT to calculate convolution.
+    // If it is not, we calculate non-circular convolution and wrap it arount.
+    //
+    // If convolution is non-circular, we use zero-padding + FFT.
+    //
+    // We assume that M+N-1>2 - we should call small case code otherwise
+    //
+    if Alg=1 then
+    begin
+        Assert(M+N-1>2, 'ConvR1DX: internal error!');
+        if Circular and FTBaseIsSmooth(M) and (M mod 2=0) then
+        begin
+            
+            //
+            // special code for circular convolution with smooth even M
+            //
+            SetLength(Buf, M);
+            APVMove(@Buf[0], 0, M-1, @A[0], 0, M-1);
+            SetLength(Buf2, M);
+            APVMove(@Buf2[0], 0, N-1, @B[0], 0, N-1);
+            I:=N;
+            while I<=M-1 do
+            begin
+                Buf2[I] := 0;
+                Inc(I);
+            end;
+            SetLength(Buf3, M);
+            FTBaseGenerateComplexFFTPlan(M div 2, Plan);
+            FFTR1DInternalEven(Buf, M, Buf3, Plan);
+            FFTR1DInternalEven(Buf2, M, Buf3, Plan);
+            Buf[0] := Buf[0]*Buf2[0];
+            Buf[1] := Buf[1]*Buf2[1];
+            I:=1;
+            while I<=M div 2-1 do
+            begin
+                AX := Buf[2*I+0];
+                AY := Buf[2*I+1];
+                BX := Buf2[2*I+0];
+                BY := Buf2[2*I+1];
+                TX := AX*BX-AY*BY;
+                TY := AX*BY+AY*BX;
+                Buf[2*I+0] := TX;
+                Buf[2*I+1] := TY;
+                Inc(I);
+            end;
+            FFTR1DInvInternalEven(Buf, M, Buf3, Plan);
+            SetLength(R, M);
+            APVMove(@R[0], 0, M-1, @Buf[0], 0, M-1);
+        end
+        else
+        begin
+            
+            //
+            // M is non-smooth or non-even, general code (circular/non-circular):
+            // * first part is the same for circular and non-circular
+            //   convolutions. zero padding, FFTs, inverse FFTs
+            // * second part differs:
+            //   * for non-circular convolution we just copy array
+            //   * for circular convolution we add array tail to its head
+            //
+            P := FTBaseFindSmoothEven(M+N-1);
+            SetLength(Buf, P);
+            APVMove(@Buf[0], 0, M-1, @A[0], 0, M-1);
+            I:=M;
+            while I<=P-1 do
+            begin
+                Buf[I] := 0;
+                Inc(I);
+            end;
+            SetLength(Buf2, P);
+            APVMove(@Buf2[0], 0, N-1, @B[0], 0, N-1);
+            I:=N;
+            while I<=P-1 do
+            begin
+                Buf2[I] := 0;
+                Inc(I);
+            end;
+            SetLength(Buf3, P);
+            FTBaseGenerateComplexFFTPlan(P div 2, Plan);
+            FFTR1DInternalEven(Buf, P, Buf3, Plan);
+            FFTR1DInternalEven(Buf2, P, Buf3, Plan);
+            Buf[0] := Buf[0]*Buf2[0];
+            Buf[1] := Buf[1]*Buf2[1];
+            I:=1;
+            while I<=P div 2-1 do
+            begin
+                AX := Buf[2*I+0];
+                AY := Buf[2*I+1];
+                BX := Buf2[2*I+0];
+                BY := Buf2[2*I+1];
+                TX := AX*BX-AY*BY;
+                TY := AX*BY+AY*BX;
+                Buf[2*I+0] := TX;
+                Buf[2*I+1] := TY;
+                Inc(I);
+            end;
+            FFTR1DInvInternalEven(Buf, P, Buf3, Plan);
+            if Circular then
+            begin
+                
+                //
+                // circular, add tail to head
+                //
+                SetLength(R, M);
+                APVMove(@R[0], 0, M-1, @Buf[0], 0, M-1);
+                if N>=2 then
+                begin
+                    APVAdd(@R[0], 0, N-2, @Buf[0], M, M+N-2);
+                end;
+            end
+            else
+            begin
+                
+                //
+                // non-circular, just copy
+                //
+                SetLength(R, M+N-1);
+                APVMove(@R[0], 0, M+N-2, @Buf[0], 0, M+N-2);
+            end;
+        end;
+        Exit;
+    end;
+    
+    //
+    // overlap-add method
+    //
+    if Alg=2 then
+    begin
+        Assert((Q+N-1) mod 2=0, 'ConvR1DX: internal error!');
+        SetLength(Buf, Q+N-1);
+        SetLength(Buf2, Q+N-1);
+        SetLength(Buf3, Q+N-1);
+        FTBaseGenerateComplexFFTPlan((Q+N-1) div 2, Plan);
+        
+        //
+        // prepare R
+        //
+        if Circular then
+        begin
+            SetLength(R, M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                R[I] := 0;
+                Inc(I);
+            end;
+        end
+        else
+        begin
+            SetLength(R, M+N-1);
+            I:=0;
+            while I<=M+N-2 do
+            begin
+                R[I] := 0;
+                Inc(I);
+            end;
+        end;
+        
+        //
+        // pre-calculated FFT(B)
+        //
+        APVMove(@Buf2[0], 0, N-1, @B[0], 0, N-1);
+        J:=N;
+        while J<=Q+N-2 do
+        begin
+            Buf2[J] := 0;
+            Inc(J);
+        end;
+        FFTR1DInternalEven(Buf2, Q+N-1, Buf3, Plan);
+        
+        //
+        // main overlap-add cycle
+        //
+        I := 0;
+        while I<=M-1 do
+        begin
+            P := Min(Q, M-I);
+            APVMove(@Buf[0], 0, P-1, @A[0], I, I+P-1);
+            J:=P;
+            while J<=Q+N-2 do
+            begin
+                Buf[J] := 0;
+                Inc(J);
+            end;
+            FFTR1DInternalEven(Buf, Q+N-1, Buf3, Plan);
+            Buf[0] := Buf[0]*Buf2[0];
+            Buf[1] := Buf[1]*Buf2[1];
+            J:=1;
+            while J<=(Q+N-1) div 2-1 do
+            begin
+                AX := Buf[2*J+0];
+                AY := Buf[2*J+1];
+                BX := Buf2[2*J+0];
+                BY := Buf2[2*J+1];
+                TX := AX*BX-AY*BY;
+                TY := AX*BY+AY*BX;
+                Buf[2*J+0] := TX;
+                Buf[2*J+1] := TY;
+                Inc(J);
+            end;
+            FFTR1DInvInternalEven(Buf, Q+N-1, Buf3, Plan);
+            if Circular then
+            begin
+                J1 := Min(I+P+N-2, M-1)-I;
+                J2 := J1+1;
+            end
+            else
+            begin
+                J1 := P+N-2;
+                J2 := J1+1;
+            end;
+            APVAdd(@R[0], I, I+J1, @Buf[0], 0, J1);
+            if P+N-2>=J2 then
+            begin
+                APVAdd(@R[0], 0, P+N-2-J2, @Buf[0], J2, P+N-2);
+            end;
+            I := I+P;
+        end;
+        Exit;
+    end;
+end;
+
+
+end.

tests/test/alglib/u_corr.pp

@@ -0,0 +1,373 @@
+(*************************************************************************
+Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************)
+unit u_corr;
+interface
+uses Math, Sysutils, u_ap, u_ftbase, u_fft, u_conv;
+
+procedure CorrC1D(const Signal : TComplex1DArray;
+     N : Integer;
+     const Pattern : TComplex1DArray;
+     M : Integer;
+     var R : TComplex1DArray);
+procedure CorrC1DCircular(const Signal : TComplex1DArray;
+     M : Integer;
+     const Pattern : TComplex1DArray;
+     N : Integer;
+     var C : TComplex1DArray);
+procedure CorrR1D(const Signal : TReal1DArray;
+     N : Integer;
+     const Pattern : TReal1DArray;
+     M : Integer;
+     var R : TReal1DArray);
+procedure CorrR1DCircular(const Signal : TReal1DArray;
+     M : Integer;
+     const Pattern : TReal1DArray;
+     N : Integer;
+     var C : TReal1DArray);
+
+implementation
+
+(*************************************************************************
+1-dimensional complex cross-correlation.
+
+For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
+
+Correlation is calculated using reduction to  convolution.  Algorithm with
+max(N,N)*log(max(N,N)) complexity is used (see  ConvC1D()  for  more  info
+about performance).
+
+IMPORTANT:
+    for  historical reasons subroutine accepts its parameters in  reversed
+    order: CorrC1D(Signal, Pattern) = Pattern x Signal (using  traditional
+    definition of cross-correlation, denoting cross-correlation as "x").
+
+INPUT PARAMETERS
+    Signal  -   array[0..N-1] - complex function to be transformed,
+                signal containing pattern
+    N       -   problem size
+    Pattern -   array[0..M-1] - complex function to be transformed,
+                pattern to search withing signal
+    M       -   problem size
+
+OUTPUT PARAMETERS
+    R       -   cross-correlation, array[0..N+M-2]:
+                * positive lags are stored in R[0..N-1],
+                  R[i] = sum(conj(pattern[j])*signal[i+j]
+                * negative lags are stored in R[N..N+M-2],
+                  R[N+M-1-i] = sum(conj(pattern[j])*signal[-i+j]
+
+NOTE:
+    It is assumed that pattern domain is [0..M-1].  If Pattern is non-zero
+on [-K..M-1],  you can still use this subroutine, just shift result by K.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure CorrC1D(const Signal : TComplex1DArray;
+     N : Integer;
+     const Pattern : TComplex1DArray;
+     M : Integer;
+     var R : TComplex1DArray);
+var
+    P : TComplex1DArray;
+    B : TComplex1DArray;
+    I : Integer;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    Assert((N>0) and (M>0), 'CorrC1D: incorrect N or M!');
+    SetLength(P, M);
+    I:=0;
+    while I<=M-1 do
+    begin
+        P[M-1-I] := Conj(Pattern[I]);
+        Inc(I);
+    end;
+    ConvC1D(P, M, Signal, N, B);
+    SetLength(R, M+N-1);
+    i1_ := (M-1) - (0);
+    for i_ := 0 to N-1 do
+    begin
+        R[i_] := B[i_+i1_];
+    end;
+    if M+N-2>=N then
+    begin
+        i1_ := (0) - (N);
+        for i_ := N to M+N-2 do
+        begin
+            R[i_] := B[i_+i1_];
+        end;
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional circular complex cross-correlation.
+
+For given Pattern/Signal returns corr(Pattern,Signal) (circular).
+Algorithm has linearithmic complexity for any M/N.
+
+IMPORTANT:
+    for  historical reasons subroutine accepts its parameters in  reversed
+    order:   CorrC1DCircular(Signal, Pattern) = Pattern x Signal    (using
+    traditional definition of cross-correlation, denoting cross-correlation
+    as "x").
+
+INPUT PARAMETERS
+    Signal  -   array[0..N-1] - complex function to be transformed,
+                periodic signal containing pattern
+    N       -   problem size
+    Pattern -   array[0..M-1] - complex function to be transformed,
+                non-periodic pattern to search withing signal
+    M       -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..M-1].
+
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure CorrC1DCircular(const Signal : TComplex1DArray;
+     M : Integer;
+     const Pattern : TComplex1DArray;
+     N : Integer;
+     var C : TComplex1DArray);
+var
+    P : TComplex1DArray;
+    B : TComplex1DArray;
+    I1 : Integer;
+    I2 : Integer;
+    I : Integer;
+    J2 : Integer;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    Assert((N>0) and (M>0), 'ConvC1DCircular: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer (at least - not shorter) that B.
+    //
+    if M<N then
+    begin
+        SetLength(B, M);
+        I1:=0;
+        while I1<=M-1 do
+        begin
+            B[I1] := C_Complex(0);
+            Inc(I1);
+        end;
+        I1 := 0;
+        while I1<N do
+        begin
+            I2 := Min(I1+M-1, N-1);
+            J2 := I2-I1;
+            i1_ := (I1) - (0);
+            for i_ := 0 to J2 do
+            begin
+                B[i_] := C_Add(B[i_], Pattern[i_+i1_]);
+            end;
+            I1 := I1+M;
+        end;
+        CorrC1DCircular(Signal, M, B, M, C);
+        Exit;
+    end;
+    
+    //
+    // Task is normalized
+    //
+    SetLength(P, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        P[N-1-I] := Conj(Pattern[I]);
+        Inc(I);
+    end;
+    ConvC1DCircular(Signal, M, P, N, B);
+    SetLength(C, M);
+    i1_ := (N-1) - (0);
+    for i_ := 0 to M-N do
+    begin
+        C[i_] := B[i_+i1_];
+    end;
+    if M-N+1<=M-1 then
+    begin
+        i1_ := (0) - (M-N+1);
+        for i_ := M-N+1 to M-1 do
+        begin
+            C[i_] := B[i_+i1_];
+        end;
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional real cross-correlation.
+
+For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
+
+Correlation is calculated using reduction to  convolution.  Algorithm with
+max(N,N)*log(max(N,N)) complexity is used (see  ConvC1D()  for  more  info
+about performance).
+
+IMPORTANT:
+    for  historical reasons subroutine accepts its parameters in  reversed
+    order: CorrR1D(Signal, Pattern) = Pattern x Signal (using  traditional
+    definition of cross-correlation, denoting cross-correlation as "x").
+
+INPUT PARAMETERS
+    Signal  -   array[0..N-1] - real function to be transformed,
+                signal containing pattern
+    N       -   problem size
+    Pattern -   array[0..M-1] - real function to be transformed,
+                pattern to search withing signal
+    M       -   problem size
+
+OUTPUT PARAMETERS
+    R       -   cross-correlation, array[0..N+M-2]:
+                * positive lags are stored in R[0..N-1],
+                  R[i] = sum(pattern[j]*signal[i+j]
+                * negative lags are stored in R[N..N+M-2],
+                  R[N+M-1-i] = sum(pattern[j]*signal[-i+j]
+
+NOTE:
+    It is assumed that pattern domain is [0..M-1].  If Pattern is non-zero
+on [-K..M-1],  you can still use this subroutine, just shift result by K.
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure CorrR1D(const Signal : TReal1DArray;
+     N : Integer;
+     const Pattern : TReal1DArray;
+     M : Integer;
+     var R : TReal1DArray);
+var
+    P : TReal1DArray;
+    B : TReal1DArray;
+    I : Integer;
+begin
+    Assert((N>0) and (M>0), 'CorrR1D: incorrect N or M!');
+    SetLength(P, M);
+    I:=0;
+    while I<=M-1 do
+    begin
+        P[M-1-I] := Pattern[I];
+        Inc(I);
+    end;
+    ConvR1D(P, M, Signal, N, B);
+    SetLength(R, M+N-1);
+    APVMove(@R[0], 0, N-1, @B[0], M-1, M+N-2);
+    if M+N-2>=N then
+    begin
+        APVMove(@R[0], N, M+N-2, @B[0], 0, M-2);
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional circular real cross-correlation.
+
+For given Pattern/Signal returns corr(Pattern,Signal) (circular).
+Algorithm has linearithmic complexity for any M/N.
+
+IMPORTANT:
+    for  historical reasons subroutine accepts its parameters in  reversed
+    order:   CorrR1DCircular(Signal, Pattern) = Pattern x Signal    (using
+    traditional definition of cross-correlation, denoting cross-correlation
+    as "x").
+
+INPUT PARAMETERS
+    Signal  -   array[0..N-1] - real function to be transformed,
+                periodic signal containing pattern
+    N       -   problem size
+    Pattern -   array[0..M-1] - real function to be transformed,
+                non-periodic pattern to search withing signal
+    M       -   problem size
+
+OUTPUT PARAMETERS
+    R   -   convolution: A*B. array[0..M-1].
+
+
+  -- ALGLIB --
+     Copyright 21.07.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure CorrR1DCircular(const Signal : TReal1DArray;
+     M : Integer;
+     const Pattern : TReal1DArray;
+     N : Integer;
+     var C : TReal1DArray);
+var
+    P : TReal1DArray;
+    B : TReal1DArray;
+    I1 : Integer;
+    I2 : Integer;
+    I : Integer;
+    J2 : Integer;
+begin
+    Assert((N>0) and (M>0), 'ConvC1DCircular: incorrect N or M!');
+    
+    //
+    // normalize task: make M>=N,
+    // so A will be longer (at least - not shorter) that B.
+    //
+    if M<N then
+    begin
+        SetLength(B, M);
+        I1:=0;
+        while I1<=M-1 do
+        begin
+            B[I1] := 0;
+            Inc(I1);
+        end;
+        I1 := 0;
+        while I1<N do
+        begin
+            I2 := Min(I1+M-1, N-1);
+            J2 := I2-I1;
+            APVAdd(@B[0], 0, J2, @Pattern[0], I1, I2);
+            I1 := I1+M;
+        end;
+        CorrR1DCircular(Signal, M, B, M, C);
+        Exit;
+    end;
+    
+    //
+    // Task is normalized
+    //
+    SetLength(P, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        P[N-1-I] := Pattern[I];
+        Inc(I);
+    end;
+    ConvR1DCircular(Signal, M, P, N, B);
+    SetLength(C, M);
+    APVMove(@C[0], 0, M-N, @B[0], N-1, M-1);
+    if M-N+1<=M-1 then
+    begin
+        APVMove(@C[0], M-N+1, M-1, @B[0], 0, N-2);
+    end;
+end;
+
+
+end.

tests/test/alglib/u_fft.pp

@@ -0,0 +1,509 @@
+(*************************************************************************
+Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************)
+unit u_fft;
+interface
+uses Math, Sysutils, u_ap, u_ftbase;
+
+procedure FFTC1D(var A : TComplex1DArray; N : Integer);
+procedure FFTC1DInv(var A : TComplex1DArray; N : Integer);
+procedure FFTR1D(const A : TReal1DArray; N : Integer; var F : TComplex1DArray);
+procedure FFTR1DInv(const F : TComplex1DArray;
+     N : Integer;
+     var A : TReal1DArray);
+procedure FFTR1DInternalEven(var A : TReal1DArray;
+     N : Integer;
+     var Buf : TReal1DArray;
+     var Plan : FTPlan);
+procedure FFTR1DInvInternalEven(var A : TReal1DArray;
+     N : Integer;
+     var Buf : TReal1DArray;
+     var Plan : FTPlan);
+
+implementation
+
+(*************************************************************************
+1-dimensional complex FFT.
+
+Array size N may be arbitrary number (composite or prime).  Composite  N's
+are handled with cache-oblivious variation of  a  Cooley-Tukey  algorithm.
+Small prime-factors are transformed using hard coded  codelets (similar to
+FFTW codelets, but without low-level  optimization),  large  prime-factors
+are handled with Bluestein's algorithm.
+
+Fastests transforms are for smooth N's (prime factors are 2, 3,  5  only),
+most fast for powers of 2. When N have prime factors  larger  than  these,
+but orders of magnitude smaller than N, computations will be about 4 times
+slower than for nearby highly composite N's. When N itself is prime, speed
+will be 6 times lower.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - complex function to be transformed
+    N   -   problem size
+    
+OUTPUT PARAMETERS
+    A   -   DFT of a input array, array[0..N-1]
+            A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
+
+
+  -- ALGLIB --
+     Copyright 29.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTC1D(var A : TComplex1DArray; N : Integer);
+var
+    Plan : FTPlan;
+    I : Integer;
+    Buf : TReal1DArray;
+begin
+    Assert(N>0, 'FFTC1D: incorrect N!');
+    
+    //
+    // Special case: N=1, FFT is just identity transform.
+    // After this block we assume that N is strictly greater than 1.
+    //
+    if N=1 then
+    begin
+        Exit;
+    end;
+    
+    //
+    // convert input array to the more convinient format
+    //
+    SetLength(Buf, 2*N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        Buf[2*I+0] := A[I].X;
+        Buf[2*I+1] := A[I].Y;
+        Inc(I);
+    end;
+    
+    //
+    // Generate plan and execute it.
+    //
+    // Plan is a combination of a successive factorizations of N and
+    // precomputed data. It is much like a FFTW plan, but is not stored
+    // between subroutine calls and is much simpler.
+    //
+    FTBaseGenerateComplexFFTPlan(N, Plan);
+    FTBaseExecutePlan(Buf, 0, N, Plan);
+    
+    //
+    // result
+    //
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I].X := Buf[2*I+0];
+        A[I].Y := Buf[2*I+1];
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional complex inverse FFT.
+
+Array size N may be arbitrary number (composite or prime).  Algorithm  has
+O(N*logN) complexity for any N (composite or prime).
+
+See FFTC1D() description for more information about algorithm performance.
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - complex array to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    A   -   inverse DFT of a input array, array[0..N-1]
+            A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
+
+
+  -- ALGLIB --
+     Copyright 29.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTC1DInv(var A : TComplex1DArray; N : Integer);
+var
+    I : Integer;
+begin
+    Assert(N>0, 'FFTC1DInv: incorrect N!');
+    
+    //
+    // Inverse DFT can be expressed in terms of the DFT as
+    //
+    //     invfft(x) = fft(x')'/N
+    //
+    // here x' means conj(x).
+    //
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I].Y := -A[I].Y;
+        Inc(I);
+    end;
+    FFTC1D(A, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I].X := A[I].X/N;
+        A[I].Y := -A[I].Y/N;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional real FFT.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - real function to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    F   -   DFT of a input array, array[0..N-1]
+            F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
+
+NOTE:
+    F[] satisfies symmetry property F[k] = conj(F[N-k]),  so just one half
+of  array  is  usually needed. But for convinience subroutine returns full
+complex array (with frequencies above N/2), so its result may be  used  by
+other FFT-related subroutines.
+
+
+  -- ALGLIB --
+     Copyright 01.06.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTR1D(const A : TReal1DArray; N : Integer; var F : TComplex1DArray);
+var
+    I : Integer;
+    N2 : Integer;
+    Idx : Integer;
+    Hn : Complex;
+    HmnC : Complex;
+    V : Complex;
+    Buf : TReal1DArray;
+    Plan : FTPlan;
+begin
+    Assert(N>0, 'FFTR1D: incorrect N!');
+    
+    //
+    // Special cases:
+    // * N=1, FFT is just identity transform.
+    // * N=2, FFT is simple too
+    //
+    // After this block we assume that N is strictly greater than 2
+    //
+    if N=1 then
+    begin
+        SetLength(F, 1);
+        F[0] := C_Complex(A[0]);
+        Exit;
+    end;
+    if N=2 then
+    begin
+        SetLength(F, 2);
+        F[0].X := A[0]+A[1];
+        F[0].Y := 0;
+        F[1].X := A[0]-A[1];
+        F[1].Y := 0;
+        Exit;
+    end;
+    
+    //
+    // Choose between odd-size and even-size FFTs
+    //
+    if N mod 2=0 then
+    begin
+        
+        //
+        // even-size real FFT, use reduction to the complex task
+        //
+        N2 := N div 2;
+        SetLength(Buf, N);
+        APVMove(@Buf[0], 0, N-1, @A[0], 0, N-1);
+        FTBaseGenerateComplexFFTPlan(N2, Plan);
+        FTBaseExecutePlan(Buf, 0, N2, Plan);
+        SetLength(F, N);
+        I:=0;
+        while I<=N2 do
+        begin
+            Idx := 2*(I mod N2);
+            Hn.X := Buf[Idx+0];
+            Hn.Y := Buf[Idx+1];
+            Idx := 2*((N2-I) mod N2);
+            HmnC.X := Buf[Idx+0];
+            HmnC.Y := -Buf[Idx+1];
+            V.X := -Sin(-2*Pi*I/N);
+            V.Y := Cos(-2*Pi*I/N);
+            F[I] := C_Sub(C_Add(Hn,HmnC),C_Mul(V,C_Sub(Hn,HmnC)));
+            F[I].X := Double(0.5)*F[I].X;
+            F[I].Y := Double(0.5)*F[I].Y;
+            Inc(I);
+        end;
+        I:=N2+1;
+        while I<=N-1 do
+        begin
+            F[I] := Conj(F[N-I]);
+            Inc(I);
+        end;
+        Exit;
+    end
+    else
+    begin
+        
+        //
+        // use complex FFT
+        //
+        SetLength(F, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            F[I] := C_Complex(A[I]);
+            Inc(I);
+        end;
+        FFTC1D(F, N);
+        Exit;
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional real inverse FFT.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    F   -   array[0..floor(N/2)] - frequencies from forward real FFT
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    A   -   inverse DFT of a input array, array[0..N-1]
+
+NOTE:
+    F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
+half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
+is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
+F[floor(N/2)] has no special properties.
+
+Relying on properties noted above, FFTR1DInv subroutine uses only elements
+from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
+N is even it ignores imaginary part of F[floor(N/2)] too.  So you can pass
+either frequencies array with N elements or reduced array with roughly N/2
+elements - subroutine will successfully transform both.
+
+
+  -- ALGLIB --
+     Copyright 01.06.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTR1DInv(const F : TComplex1DArray;
+     N : Integer;
+     var A : TReal1DArray);
+var
+    I : Integer;
+    H : TReal1DArray;
+    FH : TComplex1DArray;
+begin
+    Assert(N>0, 'FFTR1DInv: incorrect N!');
+    
+    //
+    // Special case: N=1, FFT is just identity transform.
+    // After this block we assume that N is strictly greater than 1.
+    //
+    if N=1 then
+    begin
+        SetLength(A, 1);
+        A[0] := F[0].X;
+        Exit;
+    end;
+    
+    //
+    // inverse real FFT is reduced to the inverse real FHT,
+    // which is reduced to the forward real FHT,
+    // which is reduced to the forward real FFT.
+    //
+    // Don't worry, it is really compact and efficient reduction :)
+    //
+    SetLength(H, N);
+    SetLength(A, N);
+    H[0] := F[0].X;
+    I:=1;
+    while I<=Floor(AP_Double(N)/2)-1 do
+    begin
+        H[I] := F[I].X-F[I].Y;
+        H[N-I] := F[I].X+F[I].Y;
+        Inc(I);
+    end;
+    if N mod 2=0 then
+    begin
+        H[Floor(AP_Double(N)/2)] := F[Floor(AP_Double(N)/2)].X;
+    end
+    else
+    begin
+        H[Floor(AP_Double(N)/2)] := F[Floor(AP_Double(N)/2)].X-F[Floor(AP_Double(N)/2)].Y;
+        H[Floor(AP_Double(N)/2)+1] := F[Floor(AP_Double(N)/2)].X+F[Floor(AP_Double(N)/2)].Y;
+    end;
+    FFTR1D(H, N, FH);
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I] := (FH[I].X-FH[I].Y)/N;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Internal subroutine. Never call it directly!
+
+
+  -- ALGLIB --
+     Copyright 01.06.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTR1DInternalEven(var A : TReal1DArray;
+     N : Integer;
+     var Buf : TReal1DArray;
+     var Plan : FTPlan);
+var
+    X : Double;
+    Y : Double;
+    I : Integer;
+    N2 : Integer;
+    Idx : Integer;
+    Hn : Complex;
+    HmnC : Complex;
+    V : Complex;
+begin
+    Assert((N>0) and (N mod 2=0), 'FFTR1DEvenInplace: incorrect N!');
+    
+    //
+    // Special cases:
+    // * N=2
+    //
+    // After this block we assume that N is strictly greater than 2
+    //
+    if N=2 then
+    begin
+        X := A[0]+A[1];
+        Y := A[0]-A[1];
+        A[0] := X;
+        A[1] := Y;
+        Exit;
+    end;
+    
+    //
+    // even-size real FFT, use reduction to the complex task
+    //
+    N2 := N div 2;
+    APVMove(@Buf[0], 0, N-1, @A[0], 0, N-1);
+    FTBaseExecutePlan(Buf, 0, N2, Plan);
+    A[0] := Buf[0]+Buf[1];
+    I:=1;
+    while I<=N2-1 do
+    begin
+        Idx := 2*(I mod N2);
+        Hn.X := Buf[Idx+0];
+        Hn.Y := Buf[Idx+1];
+        Idx := 2*(N2-I);
+        HmnC.X := Buf[Idx+0];
+        HmnC.Y := -Buf[Idx+1];
+        V.X := -Sin(-2*Pi*I/N);
+        V.Y := Cos(-2*Pi*I/N);
+        V := C_Sub(C_Add(Hn,HmnC),C_Mul(V,C_Sub(Hn,HmnC)));
+        A[2*I+0] := Double(0.5)*V.X;
+        A[2*I+1] := Double(0.5)*V.Y;
+        Inc(I);
+    end;
+    A[1] := Buf[0]-Buf[1];
+end;
+
+
+(*************************************************************************
+Internal subroutine. Never call it directly!
+
+
+  -- ALGLIB --
+     Copyright 01.06.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTR1DInvInternalEven(var A : TReal1DArray;
+     N : Integer;
+     var Buf : TReal1DArray;
+     var Plan : FTPlan);
+var
+    X : Double;
+    Y : Double;
+    T : Double;
+    I : Integer;
+    N2 : Integer;
+begin
+    Assert((N>0) and (N mod 2=0), 'FFTR1DInvInternalEven: incorrect N!');
+    
+    //
+    // Special cases:
+    // * N=2
+    //
+    // After this block we assume that N is strictly greater than 2
+    //
+    if N=2 then
+    begin
+        X := Double(0.5)*(A[0]+A[1]);
+        Y := Double(0.5)*(A[0]-A[1]);
+        A[0] := X;
+        A[1] := Y;
+        Exit;
+    end;
+    
+    //
+    // inverse real FFT is reduced to the inverse real FHT,
+    // which is reduced to the forward real FHT,
+    // which is reduced to the forward real FFT.
+    //
+    // Don't worry, it is really compact and efficient reduction :)
+    //
+    N2 := N div 2;
+    Buf[0] := A[0];
+    I:=1;
+    while I<=N2-1 do
+    begin
+        X := A[2*I+0];
+        Y := A[2*I+1];
+        Buf[I] := X-Y;
+        Buf[N-I] := X+Y;
+        Inc(I);
+    end;
+    Buf[N2] := A[1];
+    FFTR1DInternalEven(Buf, N, A, Plan);
+    A[0] := Buf[0]/N;
+    T := AP_Double(1)/N;
+    I:=1;
+    while I<=N2-1 do
+    begin
+        X := Buf[2*I+0];
+        Y := Buf[2*I+1];
+        A[I] := T*(X-Y);
+        A[N-I] := T*(X+Y);
+        Inc(I);
+    end;
+    A[N2] := Buf[1]/N;
+end;
+
+
+end.

tests/test/alglib/u_fht.pp

@@ -0,0 +1,122 @@
+(*************************************************************************
+Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************)
+unit u_fht;
+interface
+uses Math, Sysutils, u_ap, u_ftbase, u_fft;
+
+procedure FHTR1D(var A : TReal1DArray; N : Integer);
+procedure FHTR1DInv(var A : TReal1DArray; N : Integer);
+
+implementation
+
+(*************************************************************************
+1-dimensional Fast Hartley Transform.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - real function to be transformed
+    N   -   problem size
+    
+OUTPUT PARAMETERS
+    A   -   FHT of a input array, array[0..N-1],
+            A_out[k] = sum(A_in[j]*(cos(2*pi*j*k/N)+sin(2*pi*j*k/N)), j=0..N-1)
+
+
+  -- ALGLIB --
+     Copyright 04.06.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FHTR1D(var A : TReal1DArray; N : Integer);
+var
+    Plan : FTPlan;
+    I : Integer;
+    FA : TComplex1DArray;
+begin
+    Assert(N>0, 'FHTR1D: incorrect N!');
+    
+    //
+    // Special case: N=1, FHT is just identity transform.
+    // After this block we assume that N is strictly greater than 1.
+    //
+    if N=1 then
+    begin
+        Exit;
+    end;
+    
+    //
+    // Reduce FHt to real FFT
+    //
+    FFTR1D(A, N, FA);
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I] := FA[I].X-FA[I].Y;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+1-dimensional inverse FHT.
+
+Algorithm has O(N*logN) complexity for any N (composite or prime).
+
+INPUT PARAMETERS
+    A   -   array[0..N-1] - complex array to be transformed
+    N   -   problem size
+
+OUTPUT PARAMETERS
+    A   -   inverse FHT of a input array, array[0..N-1]
+
+
+  -- ALGLIB --
+     Copyright 29.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FHTR1DInv(var A : TReal1DArray; N : Integer);
+var
+    I : Integer;
+begin
+    Assert(N>0, 'FHTR1DInv: incorrect N!');
+    
+    //
+    // Special case: N=1, iFHT is just identity transform.
+    // After this block we assume that N is strictly greater than 1.
+    //
+    if N=1 then
+    begin
+        Exit;
+    end;
+    
+    //
+    // Inverse FHT can be expressed in terms of the FHT as
+    //
+    //     invfht(x) = fht(x)/N
+    //
+    FHTR1D(A, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I] := A[I]/N;
+        Inc(I);
+    end;
+end;
+
+
+end.

tests/test/alglib/u_ftbase.pp

@@ -0,0 +1,1733 @@
+(*************************************************************************
+Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
+
+>>> SOURCE LICENSE >>>
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation (www.fsf.org); either version 2 of the 
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+A copy of the GNU General Public License is available at
+http://www.fsf.org/licensing/licenses
+
+>>> END OF LICENSE >>>
+*************************************************************************)
+unit u_ftbase;
+interface
+uses Math, Sysutils, u_ap;
+
+type
+FTPlan = record
+    Plan : TInteger1DArray;
+    Precomputed : TReal1DArray;
+    TmpBuf : TReal1DArray;
+    StackBuf : TReal1DArray;
+end;
+
+
+
+procedure FTBaseGenerateComplexFFTPlan(N : Integer; var Plan : FTPlan);
+procedure FTBaseGenerateRealFFTPlan(N : Integer; var Plan : FTPlan);
+procedure FTBaseGenerateRealFHTPlan(N : Integer; var Plan : FTPlan);
+procedure FTBaseExecutePlan(var A : TReal1DArray;
+     AOffset : Integer;
+     N : Integer;
+     var Plan : FTPlan);
+procedure FTBaseExecutePlanRec(var A : TReal1DArray;
+     AOffset : Integer;
+     var Plan : FTPlan;
+     EntryOffset : Integer;
+     StackPtr : Integer);
+procedure FTBaseFactorize(N : Integer;
+     TaskType : Integer;
+     var N1 : Integer;
+     var N2 : Integer);
+function FTBaseIsSmooth(N : Integer):Boolean;
+function FTBaseFindSmooth(N : Integer):Integer;
+function FTBaseFindSmoothEven(N : Integer):Integer;
+function FTBaseGetFLOPEstimate(N : Integer):Double;
+
+implementation
+
+const
+    FTBasePlanEntrySize = 8;
+    FTBaseCFFTTask = 0;
+    FTBaseRFHTTask = 1;
+    FTBaseRFFTTask = 2;
+    FFTCooleyTukeyPlan = 0;
+    FFTBluesteinPlan = 1;
+    FFTCodeletPlan = 2;
+    FHTCooleyTukeyPlan = 3;
+    FHTCodeletPlan = 4;
+    FFTRealCooleyTukeyPlan = 5;
+    FFTEmptyPlan = 6;
+    FHTN2Plan = 999;
+    FTBaseUpdateTw = 4;
+    FTBaseCodeletMax = 5;
+    FTBaseCodeletRecommended = 5;
+    FTBaseInefficiencyFactor = Double(1.3);
+    FTBaseMaxSmoothFactor = 5;
+
+procedure FTBaseGeneratePlanRec(N : Integer;
+     TaskType : Integer;
+     var Plan : FTPlan;
+     var PlanSize : Integer;
+     var PrecomputedSize : Integer;
+     var PlanArraySize : Integer;
+     var TmpMemSize : Integer;
+     var StackMemSize : Integer;
+     StackPtr : Integer);forward;
+procedure FTBasePrecomputePlanRec(var Plan : FTPlan;
+     EntryOffset : Integer;
+     StackPtr : Integer);forward;
+procedure FFTTwCalc(var A : TReal1DArray;
+     AOffset : Integer;
+     N1 : Integer;
+     N2 : Integer);forward;
+procedure InternalComplexLinTranspose(var A : TReal1DArray;
+     M : Integer;
+     N : Integer;
+     AStart : Integer;
+     var Buf : TReal1DArray);forward;
+procedure InternalRealLinTranspose(var A : TReal1DArray;
+     M : Integer;
+     N : Integer;
+     AStart : Integer;
+     var Buf : TReal1DArray);forward;
+procedure FFTICLTRec(var A : TReal1DArray;
+     AStart : Integer;
+     AStride : Integer;
+     var B : TReal1DArray;
+     BStart : Integer;
+     BStride : Integer;
+     M : Integer;
+     N : Integer);forward;
+procedure FFTIRLTRec(var A : TReal1DArray;
+     AStart : Integer;
+     AStride : Integer;
+     var B : TReal1DArray;
+     BStart : Integer;
+     BStride : Integer;
+     M : Integer;
+     N : Integer);forward;
+procedure FTBaseFindSmoothRec(N : Integer;
+     Seed : Integer;
+     LeastFactor : Integer;
+     var Best : Integer);forward;
+procedure FFTArrayResize(var A : TInteger1DArray;
+     var ASize : Integer;
+     NewASize : Integer);forward;
+procedure RefFHT(var A : TReal1DArray; N : Integer; Offs : Integer);forward;
+
+
+(*************************************************************************
+This subroutine generates FFT plan - a decomposition of a N-length FFT to
+the more simpler operations. Plan consists of the root entry and the child
+entries.
+
+Subroutine parameters:
+    N               task size
+    
+Output parameters:
+    Plan            plan
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FTBaseGenerateComplexFFTPlan(N : Integer; var Plan : FTPlan);
+var
+    PlanArraySize : Integer;
+    PlanSize : Integer;
+    PrecomputedSize : Integer;
+    TmpMemSize : Integer;
+    StackMemSize : Integer;
+    StackPtr : Integer;
+begin
+    PlanArraySize := 1;
+    PlanSize := 0;
+    PrecomputedSize := 0;
+    StackMemSize := 0;
+    StackPtr := 0;
+    TmpMemSize := 2*N;
+    SetLength(Plan.Plan, PlanArraySize);
+    FTBaseGeneratePlanRec(N, FTBaseCFFTTask, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+    Assert(StackPtr=0, 'Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!');
+    SetLength(Plan.StackBuf, Max(StackMemSize, 1));
+    SetLength(Plan.TmpBuf, Max(TmpMemSize, 1));
+    SetLength(Plan.Precomputed, Max(PrecomputedSize, 1));
+    StackPtr := 0;
+    FTBasePrecomputePlanRec(Plan, 0, StackPtr);
+    Assert(StackPtr=0, 'Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!');
+end;
+
+
+(*************************************************************************
+Generates real FFT plan
+*************************************************************************)
+procedure FTBaseGenerateRealFFTPlan(N : Integer; var Plan : FTPlan);
+var
+    PlanArraySize : Integer;
+    PlanSize : Integer;
+    PrecomputedSize : Integer;
+    TmpMemSize : Integer;
+    StackMemSize : Integer;
+    StackPtr : Integer;
+begin
+    PlanArraySize := 1;
+    PlanSize := 0;
+    PrecomputedSize := 0;
+    StackMemSize := 0;
+    StackPtr := 0;
+    TmpMemSize := 2*N;
+    SetLength(Plan.Plan, PlanArraySize);
+    FTBaseGeneratePlanRec(N, FTBaseRFFTTask, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+    Assert(StackPtr=0, 'Internal error in FTBaseGenerateRealFFTPlan: stack ptr!');
+    SetLength(Plan.StackBuf, Max(StackMemSize, 1));
+    SetLength(Plan.TmpBuf, Max(TmpMemSize, 1));
+    SetLength(Plan.Precomputed, Max(PrecomputedSize, 1));
+    StackPtr := 0;
+    FTBasePrecomputePlanRec(Plan, 0, StackPtr);
+    Assert(StackPtr=0, 'Internal error in FTBaseGenerateRealFFTPlan: stack ptr!');
+end;
+
+
+(*************************************************************************
+Generates real FHT plan
+*************************************************************************)
+procedure FTBaseGenerateRealFHTPlan(N : Integer; var Plan : FTPlan);
+var
+    PlanArraySize : Integer;
+    PlanSize : Integer;
+    PrecomputedSize : Integer;
+    TmpMemSize : Integer;
+    StackMemSize : Integer;
+    StackPtr : Integer;
+begin
+    PlanArraySize := 1;
+    PlanSize := 0;
+    PrecomputedSize := 0;
+    StackMemSize := 0;
+    StackPtr := 0;
+    TmpMemSize := N;
+    SetLength(Plan.Plan, PlanArraySize);
+    FTBaseGeneratePlanRec(N, FTBaseRFHTTask, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+    Assert(StackPtr=0, 'Internal error in FTBaseGenerateRealFHTPlan: stack ptr!');
+    SetLength(Plan.StackBuf, Max(StackMemSize, 1));
+    SetLength(Plan.TmpBuf, Max(TmpMemSize, 1));
+    SetLength(Plan.Precomputed, Max(PrecomputedSize, 1));
+    StackPtr := 0;
+    FTBasePrecomputePlanRec(Plan, 0, StackPtr);
+    Assert(StackPtr=0, 'Internal error in FTBaseGenerateRealFHTPlan: stack ptr!');
+end;
+
+
+(*************************************************************************
+This subroutine executes FFT/FHT plan.
+
+If Plan is a:
+* complex FFT plan  -   sizeof(A)=2*N,
+                        A contains interleaved real/imaginary values
+* real FFT plan     -   sizeof(A)=2*N,
+                        A contains real values interleaved with zeros
+* real FHT plan     -   sizeof(A)=2*N,
+                        A contains real values interleaved with zeros
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FTBaseExecutePlan(var A : TReal1DArray;
+     AOffset : Integer;
+     N : Integer;
+     var Plan : FTPlan);
+var
+    StackPtr : Integer;
+begin
+    StackPtr := 0;
+    FTBaseExecutePlanRec(A, AOffset, Plan, 0, StackPtr);
+end;
+
+
+(*************************************************************************
+Recurrent subroutine for the FTBaseExecutePlan
+
+Parameters:
+    A           FFT'ed array
+    AOffset     offset of the FFT'ed part (distance is measured in doubles)
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FTBaseExecutePlanRec(var A : TReal1DArray;
+     AOffset : Integer;
+     var Plan : FTPlan;
+     EntryOffset : Integer;
+     StackPtr : Integer);
+var
+    I : Integer;
+    J : Integer;
+    K : Integer;
+    Idx : Integer;
+    N1 : Integer;
+    N2 : Integer;
+    N : Integer;
+    M : Integer;
+    Offs : Integer;
+    Offs1 : Integer;
+    Offs2 : Integer;
+    OffsA : Integer;
+    OffsB : Integer;
+    OffsP : Integer;
+    HK : Double;
+    HNK : Double;
+    X : Double;
+    Y : Double;
+    BX : Double;
+    BY : Double;
+    V : Double;
+    EmptyArray : TReal1DArray;
+    A0X : Double;
+    A0Y : Double;
+    A1X : Double;
+    A1Y : Double;
+    A2X : Double;
+    A2Y : Double;
+    A3X : Double;
+    A3Y : Double;
+    A4X : Double;
+    A4Y : Double;
+    V0 : Double;
+    V1 : Double;
+    V2 : Double;
+    V3 : Double;
+    T1X : Double;
+    T1Y : Double;
+    T2X : Double;
+    T2Y : Double;
+    T3X : Double;
+    T3Y : Double;
+    T4X : Double;
+    T4Y : Double;
+    T5X : Double;
+    T5Y : Double;
+    M1X : Double;
+    M1Y : Double;
+    M2X : Double;
+    M2Y : Double;
+    M3X : Double;
+    M3Y : Double;
+    M4X : Double;
+    M4Y : Double;
+    M5X : Double;
+    M5Y : Double;
+    S1X : Double;
+    S1Y : Double;
+    S2X : Double;
+    S2Y : Double;
+    S3X : Double;
+    S3Y : Double;
+    S4X : Double;
+    S4Y : Double;
+    S5X : Double;
+    S5Y : Double;
+    C1 : Double;
+    C2 : Double;
+    C3 : Double;
+    C4 : Double;
+    C5 : Double;
+    Tmp : TReal1DArray;
+begin
+    if Plan.Plan[EntryOffset+3]=FFTEmptyPlan then
+    begin
+        Exit;
+    end;
+    if Plan.Plan[EntryOffset+3]=FFTCooleyTukeyPlan then
+    begin
+        
+        //
+        // Cooley-Tukey plan
+        // * transposition
+        // * row-wise FFT
+        // * twiddle factors:
+        //   - TwBase is a basis twiddle factor for I=1, J=1
+        //   - TwRow is a twiddle factor for a second element in a row (J=1)
+        //   - Tw is a twiddle factor for a current element
+        // * transposition again
+        // * row-wise FFT again
+        //
+        N1 := Plan.Plan[EntryOffset+1];
+        N2 := Plan.Plan[EntryOffset+2];
+        InternalComplexLinTranspose(A, N1, N2, AOffset, Plan.TmpBuf);
+        I:=0;
+        while I<=N2-1 do
+        begin
+            FTBaseExecutePlanRec(A, AOffset+I*N1*2, Plan, Plan.Plan[EntryOffset+5], StackPtr);
+            Inc(I);
+        end;
+        FFTTwCalc(A, AOffset, N1, N2);
+        InternalComplexLinTranspose(A, N2, N1, AOffset, Plan.TmpBuf);
+        I:=0;
+        while I<=N1-1 do
+        begin
+            FTBaseExecutePlanRec(A, AOffset+I*N2*2, Plan, Plan.Plan[EntryOffset+6], StackPtr);
+            Inc(I);
+        end;
+        InternalComplexLinTranspose(A, N1, N2, AOffset, Plan.TmpBuf);
+        Exit;
+    end;
+    if Plan.Plan[EntryOffset+3]=FFTRealCooleyTukeyPlan then
+    begin
+        
+        //
+        // Cooley-Tukey plan
+        // * transposition
+        // * row-wise FFT
+        // * twiddle factors:
+        //   - TwBase is a basis twiddle factor for I=1, J=1
+        //   - TwRow is a twiddle factor for a second element in a row (J=1)
+        //   - Tw is a twiddle factor for a current element
+        // * transposition again
+        // * row-wise FFT again
+        //
+        N1 := Plan.Plan[EntryOffset+1];
+        N2 := Plan.Plan[EntryOffset+2];
+        InternalComplexLinTranspose(A, N2, N1, AOffset, Plan.TmpBuf);
+        I:=0;
+        while I<=N1 div 2-1 do
+        begin
+            
+            //
+            // pack two adjacent smaller real FFT's together,
+            // make one complex FFT,
+            // unpack result
+            //
+            Offs := AOffset+2*I*N2*2;
+            K:=0;
+            while K<=N2-1 do
+            begin
+                A[Offs+2*K+1] := A[Offs+2*N2+2*K+0];
+                Inc(K);
+            end;
+            FTBaseExecutePlanRec(A, Offs, Plan, Plan.Plan[EntryOffset+6], StackPtr);
+            Plan.TmpBuf[0] := A[Offs+0];
+            Plan.TmpBuf[1] := 0;
+            Plan.TmpBuf[2*N2+0] := A[Offs+1];
+            Plan.TmpBuf[2*N2+1] := 0;
+            K:=1;
+            while K<=N2-1 do
+            begin
+                Offs1 := 2*K;
+                Offs2 := 2*N2+2*K;
+                HK := A[Offs+2*K+0];
+                HNK := A[Offs+2*(N2-K)+0];
+                Plan.TmpBuf[Offs1+0] := +Double(0.5)*(HK+HNK);
+                Plan.TmpBuf[Offs2+1] := -Double(0.5)*(HK-HNK);
+                HK := A[Offs+2*K+1];
+                HNK := A[Offs+2*(N2-K)+1];
+                Plan.TmpBuf[Offs2+0] := +Double(0.5)*(HK+HNK);
+                Plan.TmpBuf[Offs1+1] := +Double(0.5)*(HK-HNK);
+                Inc(K);
+            end;
+            APVMove(@A[0], Offs, Offs+2*N2*2-1, @Plan.TmpBuf[0], 0, 2*N2*2-1);
+            Inc(I);
+        end;
+        if N1 mod 2<>0 then
+        begin
+            FTBaseExecutePlanRec(A, AOffset+(N1-1)*N2*2, Plan, Plan.Plan[EntryOffset+6], StackPtr);
+        end;
+        FFTTwCalc(A, AOffset, N2, N1);
+        InternalComplexLinTranspose(A, N1, N2, AOffset, Plan.TmpBuf);
+        I:=0;
+        while I<=N2-1 do
+        begin
+            FTBaseExecutePlanRec(A, AOffset+I*N1*2, Plan, Plan.Plan[EntryOffset+5], StackPtr);
+            Inc(I);
+        end;
+        InternalComplexLinTranspose(A, N2, N1, AOffset, Plan.TmpBuf);
+        Exit;
+    end;
+    if Plan.Plan[EntryOffset+3]=FHTCooleyTukeyPlan then
+    begin
+        
+        //
+        // Cooley-Tukey FHT plan:
+        // * transpose                    \
+        // * smaller FHT's                |
+        // * pre-process                  |
+        // * multiply by twiddle factors  | corresponds to multiplication by H1
+        // * post-process                 |
+        // * transpose again              /
+        // * multiply by H2 (smaller FHT's)
+        // * final transposition
+        //
+        // For more details see Vitezslav Vesely, "Fast algorithms
+        // of Fourier and Hartley transform and their implementation in MATLAB",
+        // page 31.
+        //
+        N1 := Plan.Plan[EntryOffset+1];
+        N2 := Plan.Plan[EntryOffset+2];
+        N := N1*N2;
+        InternalRealLinTranspose(A, N1, N2, AOffset, Plan.TmpBuf);
+        I:=0;
+        while I<=N2-1 do
+        begin
+            FTBaseExecutePlanRec(A, AOffset+I*N1, Plan, Plan.Plan[EntryOffset+5], StackPtr);
+            Inc(I);
+        end;
+        I:=0;
+        while I<=N2-1 do
+        begin
+            J:=0;
+            while J<=N1-1 do
+            begin
+                OffsA := AOffset+I*N1;
+                HK := A[OffsA+J];
+                HNK := A[OffsA+(N1-J) mod N1];
+                Offs := 2*(I*N1+J);
+                Plan.TmpBuf[Offs+0] := -Double(0.5)*(HNK-HK);
+                Plan.TmpBuf[Offs+1] := +Double(0.5)*(HK+HNK);
+                Inc(J);
+            end;
+            Inc(I);
+        end;
+        FFTTwCalc(Plan.TmpBuf, 0, N1, N2);
+        J:=0;
+        while J<=N1-1 do
+        begin
+            A[AOffset+J] := Plan.TmpBuf[2*J+0]+Plan.TmpBuf[2*J+1];
+            Inc(J);
+        end;
+        if N2 mod 2=0 then
+        begin
+            Offs := 2*(N2 div 2)*N1;
+            OffsA := AOffset+N2 div 2*N1;
+            J:=0;
+            while J<=N1-1 do
+            begin
+                A[OffsA+J] := Plan.TmpBuf[Offs+2*J+0]+Plan.TmpBuf[Offs+2*J+1];
+                Inc(J);
+            end;
+        end;
+        I:=1;
+        while I<=(N2+1) div 2-1 do
+        begin
+            Offs := 2*I*N1;
+            Offs2 := 2*(N2-I)*N1;
+            OffsA := AOffset+I*N1;
+            J:=0;
+            while J<=N1-1 do
+            begin
+                A[OffsA+J] := Plan.TmpBuf[Offs+2*J+1]+Plan.TmpBuf[Offs2+2*J+0];
+                Inc(J);
+            end;
+            OffsA := AOffset+(N2-I)*N1;
+            J:=0;
+            while J<=N1-1 do
+            begin
+                A[OffsA+J] := Plan.TmpBuf[Offs+2*J+0]+Plan.TmpBuf[Offs2+2*J+1];
+                Inc(J);
+            end;
+            Inc(I);
+        end;
+        InternalRealLinTranspose(A, N2, N1, AOffset, Plan.TmpBuf);
+        I:=0;
+        while I<=N1-1 do
+        begin
+            FTBaseExecutePlanRec(A, AOffset+I*N2, Plan, Plan.Plan[EntryOffset+6], StackPtr);
+            Inc(I);
+        end;
+        InternalRealLinTranspose(A, N1, N2, AOffset, Plan.TmpBuf);
+        Exit;
+    end;
+    if Plan.Plan[EntryOffset+3]=FHTN2Plan then
+    begin
+        
+        //
+        // Cooley-Tukey FHT plan
+        //
+        N1 := Plan.Plan[EntryOffset+1];
+        N2 := Plan.Plan[EntryOffset+2];
+        N := N1*N2;
+        RefFHT(A, N, AOffset);
+        Exit;
+    end;
+    if Plan.Plan[EntryOffset+3]=FFTCodeletPlan then
+    begin
+        N1 := Plan.Plan[EntryOffset+1];
+        N2 := Plan.Plan[EntryOffset+2];
+        N := N1*N2;
+        if N=2 then
+        begin
+            A0X := A[AOffset+0];
+            A0Y := A[AOffset+1];
+            A1X := A[AOffset+2];
+            A1Y := A[AOffset+3];
+            V0 := A0X+A1X;
+            V1 := A0Y+A1Y;
+            V2 := A0X-A1X;
+            V3 := A0Y-A1Y;
+            A[AOffset+0] := V0;
+            A[AOffset+1] := V1;
+            A[AOffset+2] := V2;
+            A[AOffset+3] := V3;
+            Exit;
+        end;
+        if N=3 then
+        begin
+            Offs := Plan.Plan[EntryOffset+7];
+            C1 := Plan.Precomputed[Offs+0];
+            C2 := Plan.Precomputed[Offs+1];
+            A0X := A[AOffset+0];
+            A0Y := A[AOffset+1];
+            A1X := A[AOffset+2];
+            A1Y := A[AOffset+3];
+            A2X := A[AOffset+4];
+            A2Y := A[AOffset+5];
+            T1X := A1X+A2X;
+            T1Y := A1Y+A2Y;
+            A0X := A0X+T1X;
+            A0Y := A0Y+T1Y;
+            M1X := C1*T1X;
+            M1Y := C1*T1Y;
+            M2X := C2*(A1Y-A2Y);
+            M2Y := C2*(A2X-A1X);
+            S1X := A0X+M1X;
+            S1Y := A0Y+M1Y;
+            A1X := S1X+M2X;
+            A1Y := S1Y+M2Y;
+            A2X := S1X-M2X;
+            A2Y := S1Y-M2Y;
+            A[AOffset+0] := A0X;
+            A[AOffset+1] := A0Y;
+            A[AOffset+2] := A1X;
+            A[AOffset+3] := A1Y;
+            A[AOffset+4] := A2X;
+            A[AOffset+5] := A2Y;
+            Exit;
+        end;
+        if N=4 then
+        begin
+            A0X := A[AOffset+0];
+            A0Y := A[AOffset+1];
+            A1X := A[AOffset+2];
+            A1Y := A[AOffset+3];
+            A2X := A[AOffset+4];
+            A2Y := A[AOffset+5];
+            A3X := A[AOffset+6];
+            A3Y := A[AOffset+7];
+            T1X := A0X+A2X;
+            T1Y := A0Y+A2Y;
+            T2X := A1X+A3X;
+            T2Y := A1Y+A3Y;
+            M2X := A0X-A2X;
+            M2Y := A0Y-A2Y;
+            M3X := A1Y-A3Y;
+            M3Y := A3X-A1X;
+            A[AOffset+0] := T1X+T2X;
+            A[AOffset+1] := T1Y+T2Y;
+            A[AOffset+4] := T1X-T2X;
+            A[AOffset+5] := T1Y-T2Y;
+            A[AOffset+2] := M2X+M3X;
+            A[AOffset+3] := M2Y+M3Y;
+            A[AOffset+6] := M2X-M3X;
+            A[AOffset+7] := M2Y-M3Y;
+            Exit;
+        end;
+        if N=5 then
+        begin
+            Offs := Plan.Plan[EntryOffset+7];
+            C1 := Plan.Precomputed[Offs+0];
+            C2 := Plan.Precomputed[Offs+1];
+            C3 := Plan.Precomputed[Offs+2];
+            C4 := Plan.Precomputed[Offs+3];
+            C5 := Plan.Precomputed[Offs+4];
+            T1X := A[AOffset+2]+A[AOffset+8];
+            T1Y := A[AOffset+3]+A[AOffset+9];
+            T2X := A[AOffset+4]+A[AOffset+6];
+            T2Y := A[AOffset+5]+A[AOffset+7];
+            T3X := A[AOffset+2]-A[AOffset+8];
+            T3Y := A[AOffset+3]-A[AOffset+9];
+            T4X := A[AOffset+6]-A[AOffset+4];
+            T4Y := A[AOffset+7]-A[AOffset+5];
+            T5X := T1X+T2X;
+            T5Y := T1Y+T2Y;
+            A[AOffset+0] := A[AOffset+0]+T5X;
+            A[AOffset+1] := A[AOffset+1]+T5Y;
+            M1X := C1*T5X;
+            M1Y := C1*T5Y;
+            M2X := C2*(T1X-T2X);
+            M2Y := C2*(T1Y-T2Y);
+            M3X := -C3*(T3Y+T4Y);
+            M3Y := C3*(T3X+T4X);
+            M4X := -C4*T4Y;
+            M4Y := C4*T4X;
+            M5X := -C5*T3Y;
+            M5Y := C5*T3X;
+            S3X := M3X-M4X;
+            S3Y := M3Y-M4Y;
+            S5X := M3X+M5X;
+            S5Y := M3Y+M5Y;
+            S1X := A[AOffset+0]+M1X;
+            S1Y := A[AOffset+1]+M1Y;
+            S2X := S1X+M2X;
+            S2Y := S1Y+M2Y;
+            S4X := S1X-M2X;
+            S4Y := S1Y-M2Y;
+            A[AOffset+2] := S2X+S3X;
+            A[AOffset+3] := S2Y+S3Y;
+            A[AOffset+4] := S4X+S5X;
+            A[AOffset+5] := S4Y+S5Y;
+            A[AOffset+6] := S4X-S5X;
+            A[AOffset+7] := S4Y-S5Y;
+            A[AOffset+8] := S2X-S3X;
+            A[AOffset+9] := S2Y-S3Y;
+            Exit;
+        end;
+    end;
+    if Plan.Plan[EntryOffset+3]=FHTCodeletPlan then
+    begin
+        N1 := Plan.Plan[EntryOffset+1];
+        N2 := Plan.Plan[EntryOffset+2];
+        N := N1*N2;
+        if N=2 then
+        begin
+            A0X := A[AOffset+0];
+            A1X := A[AOffset+1];
+            A[AOffset+0] := A0X+A1X;
+            A[AOffset+1] := A0X-A1X;
+            Exit;
+        end;
+        if N=3 then
+        begin
+            Offs := Plan.Plan[EntryOffset+7];
+            C1 := Plan.Precomputed[Offs+0];
+            C2 := Plan.Precomputed[Offs+1];
+            A0X := A[AOffset+0];
+            A1X := A[AOffset+1];
+            A2X := A[AOffset+2];
+            T1X := A1X+A2X;
+            A0X := A0X+T1X;
+            M1X := C1*T1X;
+            M2Y := C2*(A2X-A1X);
+            S1X := A0X+M1X;
+            A[AOffset+0] := A0X;
+            A[AOffset+1] := S1X-M2Y;
+            A[AOffset+2] := S1X+M2Y;
+            Exit;
+        end;
+        if N=4 then
+        begin
+            A0X := A[AOffset+0];
+            A1X := A[AOffset+1];
+            A2X := A[AOffset+2];
+            A3X := A[AOffset+3];
+            T1X := A0X+A2X;
+            T2X := A1X+A3X;
+            M2X := A0X-A2X;
+            M3Y := A3X-A1X;
+            A[AOffset+0] := T1X+T2X;
+            A[AOffset+1] := M2X-M3Y;
+            A[AOffset+2] := T1X-T2X;
+            A[AOffset+3] := M2X+M3Y;
+            Exit;
+        end;
+        if N=5 then
+        begin
+            Offs := Plan.Plan[EntryOffset+7];
+            C1 := Plan.Precomputed[Offs+0];
+            C2 := Plan.Precomputed[Offs+1];
+            C3 := Plan.Precomputed[Offs+2];
+            C4 := Plan.Precomputed[Offs+3];
+            C5 := Plan.Precomputed[Offs+4];
+            T1X := A[AOffset+1]+A[AOffset+4];
+            T2X := A[AOffset+2]+A[AOffset+3];
+            T3X := A[AOffset+1]-A[AOffset+4];
+            T4X := A[AOffset+3]-A[AOffset+2];
+            T5X := T1X+T2X;
+            V0 := A[AOffset+0]+T5X;
+            A[AOffset+0] := V0;
+            M2X := C2*(T1X-T2X);
+            M3Y := C3*(T3X+T4X);
+            S3Y := M3Y-C4*T4X;
+            S5Y := M3Y+C5*T3X;
+            S1X := V0+C1*T5X;
+            S2X := S1X+M2X;
+            S4X := S1X-M2X;
+            A[AOffset+1] := S2X-S3Y;
+            A[AOffset+2] := S4X-S5Y;
+            A[AOffset+3] := S4X+S5Y;
+            A[AOffset+4] := S2X+S3Y;
+            Exit;
+        end;
+    end;
+    if Plan.Plan[EntryOffset+3]=FFTBluesteinPlan then
+    begin
+        
+        //
+        // Bluestein plan:
+        // 1. multiply by precomputed coefficients
+        // 2. make convolution: forward FFT, multiplication by precomputed FFT
+        //    and backward FFT. backward FFT is represented as
+        //
+        //        invfft(x) = fft(x')'/M
+        //
+        //    for performance reasons reduction of inverse FFT to
+        //    forward FFT is merged with multiplication of FFT components
+        //    and last stage of Bluestein's transformation.
+        // 3. post-multiplication by Bluestein factors
+        //
+        N := Plan.Plan[EntryOffset+1];
+        M := Plan.Plan[EntryOffset+4];
+        Offs := Plan.Plan[EntryOffset+7];
+        I:=StackPtr+2*N;
+        while I<=StackPtr+2*M-1 do
+        begin
+            Plan.StackBuf[I] := 0;
+            Inc(I);
+        end;
+        OffsP := Offs+2*M;
+        OffsA := AOffset;
+        OffsB := StackPtr;
+        I:=0;
+        while I<=N-1 do
+        begin
+            BX := Plan.Precomputed[OffsP+0];
+            BY := Plan.Precomputed[OffsP+1];
+            X := A[OffsA+0];
+            Y := A[OffsA+1];
+            Plan.StackBuf[OffsB+0] := X*BX-Y*-BY;
+            Plan.StackBuf[OffsB+1] := X*-BY+Y*BX;
+            OffsP := OffsP+2;
+            OffsA := OffsA+2;
+            OffsB := OffsB+2;
+            Inc(I);
+        end;
+        FTBaseExecutePlanRec(Plan.StackBuf, StackPtr, Plan, Plan.Plan[EntryOffset+5], StackPtr+2*2*M);
+        OffsB := StackPtr;
+        OffsP := Offs;
+        I:=0;
+        while I<=M-1 do
+        begin
+            X := Plan.StackBuf[OffsB+0];
+            Y := Plan.StackBuf[OffsB+1];
+            BX := Plan.Precomputed[OffsP+0];
+            BY := Plan.Precomputed[OffsP+1];
+            Plan.StackBuf[OffsB+0] := X*BX-Y*BY;
+            Plan.StackBuf[OffsB+1] := -(X*BY+Y*BX);
+            OffsB := OffsB+2;
+            OffsP := OffsP+2;
+            Inc(I);
+        end;
+        FTBaseExecutePlanRec(Plan.StackBuf, StackPtr, Plan, Plan.Plan[EntryOffset+5], StackPtr+2*2*M);
+        OffsB := StackPtr;
+        OffsP := Offs+2*M;
+        OffsA := AOffset;
+        I:=0;
+        while I<=N-1 do
+        begin
+            X := +Plan.StackBuf[OffsB+0]/M;
+            Y := -Plan.StackBuf[OffsB+1]/M;
+            BX := Plan.Precomputed[OffsP+0];
+            BY := Plan.Precomputed[OffsP+1];
+            A[OffsA+0] := X*BX-Y*-BY;
+            A[OffsA+1] := X*-BY+Y*BX;
+            OffsP := OffsP+2;
+            OffsA := OffsA+2;
+            OffsB := OffsB+2;
+            Inc(I);
+        end;
+        Exit;
+    end;
+end;
+
+
+(*************************************************************************
+Returns good factorization N=N1*N2.
+
+Usually N1<=N2 (but not always - small N's may be exception).
+if N1<>1 then N2<>1.
+
+Factorization is chosen depending on task type and codelets we have.
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FTBaseFactorize(N : Integer;
+     TaskType : Integer;
+     var N1 : Integer;
+     var N2 : Integer);
+var
+    J : Integer;
+begin
+    N1 := 0;
+    N2 := 0;
+    
+    //
+    // try to find good codelet
+    //
+    if N1*N2<>N then
+    begin
+        J:=FTBaseCodeletRecommended;
+        while J>=2 do
+        begin
+            if N mod J=0 then
+            begin
+                N1 := J;
+                N2 := N div J;
+                Break;
+            end;
+            Dec(J);
+        end;
+    end;
+    
+    //
+    // try to factorize N
+    //
+    if N1*N2<>N then
+    begin
+        J:=FTBaseCodeletRecommended+1;
+        while J<=N-1 do
+        begin
+            if N mod J=0 then
+            begin
+                N1 := J;
+                N2 := N div J;
+                Break;
+            end;
+            Inc(J);
+        end;
+    end;
+    
+    //
+    // looks like N is prime :(
+    //
+    if N1*N2<>N then
+    begin
+        N1 := 1;
+        N2 := N;
+    end;
+    
+    //
+    // normalize
+    //
+    if (N2=1) and (N1<>1) then
+    begin
+        N2 := N1;
+        N1 := 1;
+    end;
+end;
+
+
+(*************************************************************************
+Is number smooth?
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+function FTBaseIsSmooth(N : Integer):Boolean;
+var
+    I : Integer;
+begin
+    I:=2;
+    while I<=FTBaseMaxSmoothFactor do
+    begin
+        while N mod I=0 do
+        begin
+            N := N div I;
+        end;
+        Inc(I);
+    end;
+    Result := N=1;
+end;
+
+
+(*************************************************************************
+Returns smallest smooth (divisible only by 2, 3, 5) number that is greater
+than or equal to max(N,2)
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+function FTBaseFindSmooth(N : Integer):Integer;
+var
+    Best : Integer;
+begin
+    Best := 2;
+    while Best<N do
+    begin
+        Best := 2*Best;
+    end;
+    FTBaseFindSmoothRec(N, 1, 2, Best);
+    Result := Best;
+end;
+
+
+(*************************************************************************
+Returns  smallest  smooth  (divisible only by 2, 3, 5) even number that is
+greater than or equal to max(N,2)
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+function FTBaseFindSmoothEven(N : Integer):Integer;
+var
+    Best : Integer;
+begin
+    Best := 2;
+    while Best<N do
+    begin
+        Best := 2*Best;
+    end;
+    FTBaseFindSmoothRec(N, 2, 2, Best);
+    Result := Best;
+end;
+
+
+(*************************************************************************
+Returns estimate of FLOP count for the FFT.
+
+It is only an estimate based on operations count for the PERFECT FFT
+and relative inefficiency of the algorithm actually used.
+
+N should be power of 2, estimates are badly wrong for non-power-of-2 N's.
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+function FTBaseGetFLOPEstimate(N : Integer):Double;
+begin
+    Result := FTBaseInefficiencyFactor*(4*N*Ln(N)/Ln(2)-6*N+8);
+end;
+
+
+(*************************************************************************
+Recurrent subroutine for the FFTGeneratePlan:
+
+PARAMETERS:
+    N                   plan size
+    IsReal              whether input is real or not.
+                        subroutine MUST NOT ignore this flag because real
+                        inputs comes with non-initialized imaginary parts,
+                        so ignoring this flag will result in corrupted output
+    HalfOut             whether full output or only half of it from 0 to
+                        floor(N/2) is needed. This flag may be ignored if
+                        doing so will simplify calculations
+    Plan                plan array
+    PlanSize            size of used part (in integers)
+    PrecomputedSize     size of precomputed array allocated yet
+    PlanArraySize       plan array size (actual)
+    TmpMemSize          temporary memory required size
+    BluesteinMemSize    temporary memory required size
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FTBaseGeneratePlanRec(N : Integer;
+     TaskType : Integer;
+     var Plan : FTPlan;
+     var PlanSize : Integer;
+     var PrecomputedSize : Integer;
+     var PlanArraySize : Integer;
+     var TmpMemSize : Integer;
+     var StackMemSize : Integer;
+     StackPtr : Integer);
+var
+    J : Integer;
+    K : Integer;
+    M : Integer;
+    N1 : Integer;
+    N2 : Integer;
+    ESize : Integer;
+    EntryOffset : Integer;
+begin
+    
+    //
+    // prepare
+    //
+    if PlanSize+FTBasePlanEntrySize>PlanArraySize then
+    begin
+        FFTArrayResize(Plan.Plan, PlanArraySize, 8*PlanArraySize);
+    end;
+    EntryOffset := PlanSize;
+    ESize := FTBasePlanEntrySize;
+    PlanSize := PlanSize+ESize;
+    
+    //
+    // if N=1, generate empty plan and exit
+    //
+    if N=1 then
+    begin
+        Plan.Plan[EntryOffset+0] := ESize;
+        Plan.Plan[EntryOffset+1] := -1;
+        Plan.Plan[EntryOffset+2] := -1;
+        Plan.Plan[EntryOffset+3] := FFTEmptyPlan;
+        Plan.Plan[EntryOffset+4] := -1;
+        Plan.Plan[EntryOffset+5] := -1;
+        Plan.Plan[EntryOffset+6] := -1;
+        Plan.Plan[EntryOffset+7] := -1;
+        Exit;
+    end;
+    
+    //
+    // generate plans
+    //
+    FTBaseFactorize(N, TaskType, N1, N2);
+    if (TaskType=FTBaseCFFTTask) or (TaskType=FTBaseRFFTTask) then
+    begin
+        
+        //
+        // complex FFT plans
+        //
+        if N1<>1 then
+        begin
+            
+            //
+            // Cooley-Tukey plan (real or complex)
+            //
+            // Note that child plans are COMPLEX
+            // (whether plan itself is complex or not).
+            //
+            TmpMemSize := Max(TmpMemSize, 2*N1*N2);
+            Plan.Plan[EntryOffset+0] := ESize;
+            Plan.Plan[EntryOffset+1] := N1;
+            Plan.Plan[EntryOffset+2] := N2;
+            if TaskType=FTBaseCFFTTask then
+            begin
+                Plan.Plan[EntryOffset+3] := FFTCooleyTukeyPlan;
+            end
+            else
+            begin
+                Plan.Plan[EntryOffset+3] := FFTRealCooleyTukeyPlan;
+            end;
+            Plan.Plan[EntryOffset+4] := 0;
+            Plan.Plan[EntryOffset+5] := PlanSize;
+            FTBaseGeneratePlanRec(N1, FTBaseCFFTTask, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+            Plan.Plan[EntryOffset+6] := PlanSize;
+            FTBaseGeneratePlanRec(N2, FTBaseCFFTTask, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+            Plan.Plan[EntryOffset+7] := -1;
+            Exit;
+        end
+        else
+        begin
+            if (N=2) or (N=3) or (N=4) or (N=5) then
+            begin
+                
+                //
+                // hard-coded plan
+                //
+                Plan.Plan[EntryOffset+0] := ESize;
+                Plan.Plan[EntryOffset+1] := N1;
+                Plan.Plan[EntryOffset+2] := N2;
+                Plan.Plan[EntryOffset+3] := FFTCodeletPlan;
+                Plan.Plan[EntryOffset+4] := 0;
+                Plan.Plan[EntryOffset+5] := -1;
+                Plan.Plan[EntryOffset+6] := -1;
+                Plan.Plan[EntryOffset+7] := PrecomputedSize;
+                if N=3 then
+                begin
+                    PrecomputedSize := PrecomputedSize+2;
+                end;
+                if N=5 then
+                begin
+                    PrecomputedSize := PrecomputedSize+5;
+                end;
+                Exit;
+            end
+            else
+            begin
+                
+                //
+                // Bluestein's plan
+                //
+                // Select such M that M>=2*N-1, M is composite, and M's
+                // factors are 2, 3, 5
+                //
+                K := 2*N2-1;
+                M := FTBaseFindSmooth(K);
+                TmpMemSize := Max(TmpMemSize, 2*M);
+                Plan.Plan[EntryOffset+0] := ESize;
+                Plan.Plan[EntryOffset+1] := N2;
+                Plan.Plan[EntryOffset+2] := -1;
+                Plan.Plan[EntryOffset+3] := FFTBluesteinPlan;
+                Plan.Plan[EntryOffset+4] := M;
+                Plan.Plan[EntryOffset+5] := PlanSize;
+                StackPtr := StackPtr+2*2*M;
+                StackMemSize := Max(StackMemSize, StackPtr);
+                FTBaseGeneratePlanRec(M, FTBaseCFFTTask, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+                StackPtr := StackPtr-2*2*M;
+                Plan.Plan[EntryOffset+6] := -1;
+                Plan.Plan[EntryOffset+7] := PrecomputedSize;
+                PrecomputedSize := PrecomputedSize+2*M+2*N;
+                Exit;
+            end;
+        end;
+    end;
+    if TaskType=FTBaseRFHTTask then
+    begin
+        
+        //
+        // real FHT plans
+        //
+        if N1<>1 then
+        begin
+            
+            //
+            // Cooley-Tukey plan
+            //
+            //
+            TmpMemSize := Max(TmpMemSize, 2*N1*N2);
+            Plan.Plan[EntryOffset+0] := ESize;
+            Plan.Plan[EntryOffset+1] := N1;
+            Plan.Plan[EntryOffset+2] := N2;
+            Plan.Plan[EntryOffset+3] := FHTCooleyTukeyPlan;
+            Plan.Plan[EntryOffset+4] := 0;
+            Plan.Plan[EntryOffset+5] := PlanSize;
+            FTBaseGeneratePlanRec(N1, TaskType, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+            Plan.Plan[EntryOffset+6] := PlanSize;
+            FTBaseGeneratePlanRec(N2, TaskType, Plan, PlanSize, PrecomputedSize, PlanArraySize, TmpMemSize, StackMemSize, StackPtr);
+            Plan.Plan[EntryOffset+7] := -1;
+            Exit;
+        end
+        else
+        begin
+            
+            //
+            // N2 plan
+            //
+            Plan.Plan[EntryOffset+0] := ESize;
+            Plan.Plan[EntryOffset+1] := N1;
+            Plan.Plan[EntryOffset+2] := N2;
+            Plan.Plan[EntryOffset+3] := FHTN2Plan;
+            Plan.Plan[EntryOffset+4] := 0;
+            Plan.Plan[EntryOffset+5] := -1;
+            Plan.Plan[EntryOffset+6] := -1;
+            Plan.Plan[EntryOffset+7] := -1;
+            if (N=2) or (N=3) or (N=4) or (N=5) then
+            begin
+                
+                //
+                // hard-coded plan
+                //
+                Plan.Plan[EntryOffset+0] := ESize;
+                Plan.Plan[EntryOffset+1] := N1;
+                Plan.Plan[EntryOffset+2] := N2;
+                Plan.Plan[EntryOffset+3] := FHTCodeletPlan;
+                Plan.Plan[EntryOffset+4] := 0;
+                Plan.Plan[EntryOffset+5] := -1;
+                Plan.Plan[EntryOffset+6] := -1;
+                Plan.Plan[EntryOffset+7] := PrecomputedSize;
+                if N=3 then
+                begin
+                    PrecomputedSize := PrecomputedSize+2;
+                end;
+                if N=5 then
+                begin
+                    PrecomputedSize := PrecomputedSize+5;
+                end;
+                Exit;
+            end;
+            Exit;
+        end;
+    end;
+end;
+
+
+(*************************************************************************
+Recurrent subroutine for precomputing FFT plans
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FTBasePrecomputePlanRec(var Plan : FTPlan;
+     EntryOffset : Integer;
+     StackPtr : Integer);
+var
+    I : Integer;
+    J : Integer;
+    Idx : Integer;
+    N1 : Integer;
+    N2 : Integer;
+    N : Integer;
+    M : Integer;
+    Offs : Integer;
+    V : Double;
+    EmptyArray : TReal1DArray;
+    BX : Double;
+    BY : Double;
+begin
+    if (Plan.Plan[EntryOffset+3]=FFTCooleyTukeyPlan) or (Plan.Plan[EntryOffset+3]=FFTRealCooleyTukeyPlan) or (Plan.Plan[EntryOffset+3]=FHTCooleyTukeyPlan) then
+    begin
+        FTBasePrecomputePlanRec(Plan, Plan.Plan[EntryOffset+5], StackPtr);
+        FTBasePrecomputePlanRec(Plan, Plan.Plan[EntryOffset+6], StackPtr);
+        Exit;
+    end;
+    if (Plan.Plan[EntryOffset+3]=FFTCodeletPlan) or (Plan.Plan[EntryOffset+3]=FHTCodeletPlan) then
+    begin
+        N1 := Plan.Plan[EntryOffset+1];
+        N2 := Plan.Plan[EntryOffset+2];
+        N := N1*N2;
+        if N=3 then
+        begin
+            Offs := Plan.Plan[EntryOffset+7];
+            Plan.Precomputed[Offs+0] := Cos(2*Pi/3)-1;
+            Plan.Precomputed[Offs+1] := Sin(2*Pi/3);
+            Exit;
+        end;
+        if N=5 then
+        begin
+            Offs := Plan.Plan[EntryOffset+7];
+            V := 2*Pi/5;
+            Plan.Precomputed[Offs+0] := (cos(V)+cos(2*V))/2-1;
+            Plan.Precomputed[Offs+1] := (cos(V)-cos(2*V))/2;
+            Plan.Precomputed[Offs+2] := -sin(V);
+            Plan.Precomputed[Offs+3] := -(sin(V)+sin(2*V));
+            Plan.Precomputed[Offs+4] := sin(V)-sin(2*V);
+            Exit;
+        end;
+    end;
+    if Plan.Plan[EntryOffset+3]=FFTBluesteinPlan then
+    begin
+        FTBasePrecomputePlanRec(Plan, Plan.Plan[EntryOffset+5], StackPtr);
+        N := Plan.Plan[EntryOffset+1];
+        M := Plan.Plan[EntryOffset+4];
+        Offs := Plan.Plan[EntryOffset+7];
+        I:=0;
+        while I<=2*M-1 do
+        begin
+            Plan.Precomputed[Offs+I] := 0;
+            Inc(I);
+        end;
+        I:=0;
+        while I<=N-1 do
+        begin
+            BX := Cos(Pi*AP_Sqr(I)/N);
+            BY := Sin(Pi*AP_Sqr(I)/N);
+            Plan.Precomputed[Offs+2*I+0] := BX;
+            Plan.Precomputed[Offs+2*I+1] := BY;
+            Plan.Precomputed[Offs+2*M+2*I+0] := BX;
+            Plan.Precomputed[Offs+2*M+2*I+1] := BY;
+            if I>0 then
+            begin
+                Plan.Precomputed[Offs+2*(M-I)+0] := BX;
+                Plan.Precomputed[Offs+2*(M-I)+1] := BY;
+            end;
+            Inc(I);
+        end;
+        FTBaseExecutePlanRec(Plan.Precomputed, Offs, Plan, Plan.Plan[EntryOffset+5], StackPtr);
+        Exit;
+    end;
+end;
+
+
+(*************************************************************************
+Twiddle factors calculation
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTTwCalc(var A : TReal1DArray;
+     AOffset : Integer;
+     N1 : Integer;
+     N2 : Integer);
+var
+    I : Integer;
+    J : Integer;
+    N : Integer;
+    Idx : Integer;
+    Offs : Integer;
+    X : Double;
+    Y : Double;
+    TwXM1 : Double;
+    TwY : Double;
+    TwBaseXM1 : Double;
+    TwBaseY : Double;
+    TwRowXM1 : Double;
+    TwRowY : Double;
+    TmpX : Double;
+    TmpY : Double;
+    V : Double;
+begin
+    N := N1*N2;
+    V := -2*Pi/N;
+    TwBaseXM1 := -2*AP_Sqr(Sin(Double(0.5)*V));
+    TwBaseY := Sin(V);
+    TwRowXM1 := 0;
+    TwRowY := 0;
+    I:=0;
+    while I<=N2-1 do
+    begin
+        TwXM1 := 0;
+        TwY := 0;
+        J:=0;
+        while J<=N1-1 do
+        begin
+            Idx := I*N1+J;
+            Offs := AOffset+2*Idx;
+            X := A[Offs+0];
+            Y := A[Offs+1];
+            TmpX := X*TwXM1-Y*TwY;
+            TmpY := X*TwY+Y*TwXM1;
+            A[Offs+0] := X+TmpX;
+            A[Offs+1] := Y+TmpY;
+            
+            //
+            // update Tw: Tw(new) = Tw(old)*TwRow
+            //
+            if J<N1-1 then
+            begin
+                if J mod FTBaseUpdateTw=0 then
+                begin
+                    V := -2*Pi*I*(J+1)/N;
+                    TwXM1 := -2*AP_Sqr(Sin(Double(0.5)*V));
+                    TwY := Sin(V);
+                end
+                else
+                begin
+                    TmpX := TwRowXM1+TwXM1*TwRowXM1-TwY*TwRowY;
+                    TmpY := TwRowY+TwXM1*TwRowY+TwY*TwRowXM1;
+                    TwXM1 := TwXM1+TmpX;
+                    TwY := TwY+TmpY;
+                end;
+            end;
+            Inc(J);
+        end;
+        
+        //
+        // update TwRow: TwRow(new) = TwRow(old)*TwBase
+        //
+        if I<N2-1 then
+        begin
+            if J mod FTBaseUpdateTw=0 then
+            begin
+                V := -2*Pi*(I+1)/N;
+                TwRowXM1 := -2*AP_Sqr(Sin(Double(0.5)*V));
+                TwRowY := Sin(V);
+            end
+            else
+            begin
+                TmpX := TwBaseXM1+TwRowXM1*TwBaseXM1-TwRowY*TwBaseY;
+                TmpY := TwBaseY+TwRowXM1*TwBaseY+TwRowY*TwBaseXM1;
+                TwRowXM1 := TwRowXM1+TmpX;
+                TwRowY := TwRowY+TmpY;
+            end;
+        end;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Linear transpose: transpose complex matrix stored in 1-dimensional array
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure InternalComplexLinTranspose(var A : TReal1DArray;
+     M : Integer;
+     N : Integer;
+     AStart : Integer;
+     var Buf : TReal1DArray);
+begin
+    FFTICLTRec(A, AStart, N, Buf, 0, M, M, N);
+    APVMove(@A[0], AStart, AStart+2*M*N-1, @Buf[0], 0, 2*M*N-1);
+end;
+
+
+(*************************************************************************
+Linear transpose: transpose real matrix stored in 1-dimensional array
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure InternalRealLinTranspose(var A : TReal1DArray;
+     M : Integer;
+     N : Integer;
+     AStart : Integer;
+     var Buf : TReal1DArray);
+begin
+    FFTIRLTRec(A, AStart, N, Buf, 0, M, M, N);
+    APVMove(@A[0], AStart, AStart+M*N-1, @Buf[0], 0, M*N-1);
+end;
+
+
+(*************************************************************************
+Recurrent subroutine for a InternalComplexLinTranspose
+
+Write A^T to B, where:
+* A is m*n complex matrix stored in array A as pairs of real/image values,
+  beginning from AStart position, with AStride stride
+* B is n*m complex matrix stored in array B as pairs of real/image values,
+  beginning from BStart position, with BStride stride
+stride is measured in complex numbers, i.e. in real/image pairs.
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTICLTRec(var A : TReal1DArray;
+     AStart : Integer;
+     AStride : Integer;
+     var B : TReal1DArray;
+     BStart : Integer;
+     BStride : Integer;
+     M : Integer;
+     N : Integer);
+var
+    I : Integer;
+    J : Integer;
+    Idx1 : Integer;
+    Idx2 : Integer;
+    M2 : Integer;
+    M1 : Integer;
+    N1 : Integer;
+begin
+    if (M=0) or (N=0) then
+    begin
+        Exit;
+    end;
+    if Max(M, N)<=8 then
+    begin
+        M2 := 2*BStride;
+        I:=0;
+        while I<=M-1 do
+        begin
+            Idx1 := BStart+2*I;
+            Idx2 := AStart+2*I*AStride;
+            J:=0;
+            while J<=N-1 do
+            begin
+                B[Idx1+0] := A[Idx2+0];
+                B[Idx1+1] := A[Idx2+1];
+                Idx1 := Idx1+M2;
+                Idx2 := Idx2+2;
+                Inc(J);
+            end;
+            Inc(I);
+        end;
+        Exit;
+    end;
+    if N>M then
+    begin
+        
+        //
+        // New partition:
+        //
+        // "A^T -> B" becomes "(A1 A2)^T -> ( B1 )
+        //                                  ( B2 )
+        //
+        N1 := N div 2;
+        if (N-N1>=8) and (N1 mod 8<>0) then
+        begin
+            N1 := N1+(8-N1 mod 8);
+        end;
+        Assert(N-N1>0);
+        FFTICLTRec(A, AStart, AStride, B, BStart, BStride, M, N1);
+        FFTICLTRec(A, AStart+2*N1, AStride, B, BStart+2*N1*BStride, BStride, M, N-N1);
+    end
+    else
+    begin
+        
+        //
+        // New partition:
+        //
+        // "A^T -> B" becomes "( A1 )^T -> ( B1 B2 )
+        //                     ( A2 )
+        //
+        M1 := M div 2;
+        if (M-M1>=8) and (M1 mod 8<>0) then
+        begin
+            M1 := M1+(8-M1 mod 8);
+        end;
+        Assert(M-M1>0);
+        FFTICLTRec(A, AStart, AStride, B, BStart, BStride, M1, N);
+        FFTICLTRec(A, AStart+2*M1*AStride, AStride, B, BStart+2*M1, BStride, M-M1, N);
+    end;
+end;
+
+
+(*************************************************************************
+Recurrent subroutine for a InternalRealLinTranspose
+
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTIRLTRec(var A : TReal1DArray;
+     AStart : Integer;
+     AStride : Integer;
+     var B : TReal1DArray;
+     BStart : Integer;
+     BStride : Integer;
+     M : Integer;
+     N : Integer);
+var
+    I : Integer;
+    J : Integer;
+    Idx1 : Integer;
+    Idx2 : Integer;
+    M1 : Integer;
+    N1 : Integer;
+begin
+    if (M=0) or (N=0) then
+    begin
+        Exit;
+    end;
+    if Max(M, N)<=8 then
+    begin
+        I:=0;
+        while I<=M-1 do
+        begin
+            Idx1 := BStart+I;
+            Idx2 := AStart+I*AStride;
+            J:=0;
+            while J<=N-1 do
+            begin
+                B[Idx1] := A[Idx2];
+                Idx1 := Idx1+BStride;
+                Idx2 := Idx2+1;
+                Inc(J);
+            end;
+            Inc(I);
+        end;
+        Exit;
+    end;
+    if N>M then
+    begin
+        
+        //
+        // New partition:
+        //
+        // "A^T -> B" becomes "(A1 A2)^T -> ( B1 )
+        //                                  ( B2 )
+        //
+        N1 := N div 2;
+        if (N-N1>=8) and (N1 mod 8<>0) then
+        begin
+            N1 := N1+(8-N1 mod 8);
+        end;
+        Assert(N-N1>0);
+        FFTIRLTRec(A, AStart, AStride, B, BStart, BStride, M, N1);
+        FFTIRLTRec(A, AStart+N1, AStride, B, BStart+N1*BStride, BStride, M, N-N1);
+    end
+    else
+    begin
+        
+        //
+        // New partition:
+        //
+        // "A^T -> B" becomes "( A1 )^T -> ( B1 B2 )
+        //                     ( A2 )
+        //
+        M1 := M div 2;
+        if (M-M1>=8) and (M1 mod 8<>0) then
+        begin
+            M1 := M1+(8-M1 mod 8);
+        end;
+        Assert(M-M1>0);
+        FFTIRLTRec(A, AStart, AStride, B, BStart, BStride, M1, N);
+        FFTIRLTRec(A, AStart+M1*AStride, AStride, B, BStart+M1, BStride, M-M1, N);
+    end;
+end;
+
+
+(*************************************************************************
+recurrent subroutine for FFTFindSmoothRec
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FTBaseFindSmoothRec(N : Integer;
+     Seed : Integer;
+     LeastFactor : Integer;
+     var Best : Integer);
+begin
+    Assert(FTBaseMaxSmoothFactor<=5, 'FTBaseFindSmoothRec: internal error!');
+    if Seed>=N then
+    begin
+        Best := Min(Best, Seed);
+        Exit;
+    end;
+    if LeastFactor<=2 then
+    begin
+        FTBaseFindSmoothRec(N, Seed*2, 2, Best);
+    end;
+    if LeastFactor<=3 then
+    begin
+        FTBaseFindSmoothRec(N, Seed*3, 3, Best);
+    end;
+    if LeastFactor<=5 then
+    begin
+        FTBaseFindSmoothRec(N, Seed*5, 5, Best);
+    end;
+end;
+
+
+(*************************************************************************
+Internal subroutine: array resize
+
+  -- ALGLIB --
+     Copyright 01.05.2009 by Bochkanov Sergey
+*************************************************************************)
+procedure FFTArrayResize(var A : TInteger1DArray;
+     var ASize : Integer;
+     NewASize : Integer);
+var
+    Tmp : TInteger1DArray;
+    I : Integer;
+begin
+    SetLength(Tmp, ASize);
+    I:=0;
+    while I<=ASize-1 do
+    begin
+        Tmp[I] := A[I];
+        Inc(I);
+    end;
+    SetLength(A, NewASize);
+    I:=0;
+    while I<=ASize-1 do
+    begin
+        A[I] := Tmp[I];
+        Inc(I);
+    end;
+    ASize := NewASize;
+end;
+
+
+(*************************************************************************
+Reference FHT stub
+*************************************************************************)
+procedure RefFHT(var A : TReal1DArray; N : Integer; Offs : Integer);
+var
+    Buf : TReal1DArray;
+    I : Integer;
+    J : Integer;
+    V : Double;
+begin
+    Assert(N>0, 'RefFHTR1D: incorrect N!');
+    SetLength(Buf, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        V := 0;
+        J:=0;
+        while J<=N-1 do
+        begin
+            V := V+A[Offs+J]*(Cos(2*Pi*I*J/N)+Sin(2*Pi*I*J/N));
+            Inc(J);
+        end;
+        Buf[I] := V;
+        Inc(I);
+    end;
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[Offs+I] := Buf[I];
+        Inc(I);
+    end;
+end;
+
+
+end.

tests/test/alglib/u_testconvunit.pp

@@ -0,0 +1,588 @@
+unit u_testconvunit;
+interface
+uses Math, Sysutils, u_ap, u_ftbase, u_fft, u_conv;
+
+function TestConv(Silent : Boolean):Boolean;
+function testconvunit_test_silent():Boolean;
+function testconvunit_test():Boolean;
+
+implementation
+
+procedure RefConvC1D(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);forward;
+procedure RefConvC1DCircular(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);forward;
+procedure RefConvR1D(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);forward;
+procedure RefConvR1DCircular(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);forward;
+
+
+(*************************************************************************
+Test
+*************************************************************************)
+function TestConv(Silent : Boolean):Boolean;
+var
+    M : Integer;
+    N : Integer;
+    I : Integer;
+    RKind : Integer;
+    CircKind : Integer;
+    RA : TReal1DArray;
+    RB : TReal1DArray;
+    RR1 : TReal1DArray;
+    RR2 : TReal1DArray;
+    CA : TComplex1DArray;
+    CB : TComplex1DArray;
+    CR1 : TComplex1DArray;
+    CR2 : TComplex1DArray;
+    MaxN : Integer;
+    RefErr : Double;
+    RefRErr : Double;
+    InvErr : Double;
+    InvRErr : Double;
+    ErrTol : Double;
+    RefErrors : Boolean;
+    RefRErrors : Boolean;
+    InvErrors : Boolean;
+    InvRErrors : Boolean;
+    WasErrors : Boolean;
+begin
+    MaxN := 32;
+    ErrTol := 100000*Power(MaxN, AP_Double(3)/2)*MachineEpsilon;
+    RefErrors := False;
+    RefRErrors := False;
+    InvErrors := False;
+    InvRErrors := False;
+    WasErrors := False;
+    
+    //
+    // Test against reference O(N^2) implementation.
+    //
+    // Automatic ConvC1D() and different algorithms of ConvC1DX() are tested.
+    //
+    RefErr := 0;
+    RefRErr := 0;
+    M:=1;
+    while M<=MaxN do
+    begin
+        N:=1;
+        while N<=MaxN do
+        begin
+            CircKind:=0;
+            while CircKind<=1 do
+            begin
+                RKind:=-3;
+                while RKind<=1 do
+                begin
+                    
+                    //
+                    // skip impossible combinations of parameters:
+                    // * circular convolution, M<N, RKind<>-3 - internal subroutine does not support M<N.
+                    //
+                    if (CircKind<>0) and (M<N) and (RKind<>-3) then
+                    begin
+                        Inc(RKind);
+                        Continue;
+                    end;
+                    
+                    //
+                    // Complex convolution
+                    //
+                    SetLength(CA, M);
+                    I:=0;
+                    while I<=M-1 do
+                    begin
+                        CA[I].X := 2*RandomReal-1;
+                        CA[I].Y := 2*RandomReal-1;
+                        Inc(I);
+                    end;
+                    SetLength(CB, N);
+                    I:=0;
+                    while I<=N-1 do
+                    begin
+                        CB[I].X := 2*RandomReal-1;
+                        CB[I].Y := 2*RandomReal-1;
+                        Inc(I);
+                    end;
+                    SetLength(CR1, 1);
+                    if RKind=-3 then
+                    begin
+                        
+                        //
+                        // test wrapper subroutine:
+                        // * circular/non-circular
+                        //
+                        if CircKind=0 then
+                        begin
+                            ConvC1D(CA, M, CB, N, CR1);
+                        end
+                        else
+                        begin
+                            ConvC1DCircular(CA, M, CB, N, CR1);
+                        end;
+                    end
+                    else
+                    begin
+                        
+                        //
+                        // test internal subroutine
+                        //
+                        if M>=N then
+                        begin
+                            
+                            //
+                            // test internal subroutine:
+                            // * circular/non-circular mode
+                            //
+                            ConvC1DX(CA, M, CB, N, CircKind<>0, RKind, 0, CR1);
+                        end
+                        else
+                        begin
+                            
+                            //
+                            // test internal subroutine - circular mode only
+                            //
+                            Assert(CircKind=0, 'Convolution test: internal error!');
+                            ConvC1DX(CB, N, CA, M, False, RKind, 0, CR1);
+                        end;
+                    end;
+                    if CircKind=0 then
+                    begin
+                        RefConvC1D(CA, M, CB, N, CR2);
+                    end
+                    else
+                    begin
+                        RefConvC1DCircular(CA, M, CB, N, CR2);
+                    end;
+                    if CircKind=0 then
+                    begin
+                        I:=0;
+                        while I<=M+N-2 do
+                        begin
+                            RefErr := Max(RefErr, AbsComplex(C_Sub(CR1[I],CR2[I])));
+                            Inc(I);
+                        end;
+                    end
+                    else
+                    begin
+                        I:=0;
+                        while I<=M-1 do
+                        begin
+                            RefErr := Max(RefErr, AbsComplex(C_Sub(CR1[I],CR2[I])));
+                            Inc(I);
+                        end;
+                    end;
+                    
+                    //
+                    // Real convolution
+                    //
+                    SetLength(RA, M);
+                    I:=0;
+                    while I<=M-1 do
+                    begin
+                        RA[I] := 2*RandomReal-1;
+                        Inc(I);
+                    end;
+                    SetLength(RB, N);
+                    I:=0;
+                    while I<=N-1 do
+                    begin
+                        RB[I] := 2*RandomReal-1;
+                        Inc(I);
+                    end;
+                    SetLength(RR1, 1);
+                    if RKind=-3 then
+                    begin
+                        
+                        //
+                        // test wrapper subroutine:
+                        // * circular/non-circular
+                        //
+                        if CircKind=0 then
+                        begin
+                            ConvR1D(RA, M, RB, N, RR1);
+                        end
+                        else
+                        begin
+                            ConvR1DCircular(RA, M, RB, N, RR1);
+                        end;
+                    end
+                    else
+                    begin
+                        if M>=N then
+                        begin
+                            
+                            //
+                            // test internal subroutine:
+                            // * circular/non-circular mode
+                            //
+                            ConvR1DX(RA, M, RB, N, CircKind<>0, RKind, 0, RR1);
+                        end
+                        else
+                        begin
+                            
+                            //
+                            // test internal subroutine - non-circular mode only
+                            //
+                            ConvR1DX(RB, N, RA, M, CircKind<>0, RKind, 0, RR1);
+                        end;
+                    end;
+                    if CircKind=0 then
+                    begin
+                        RefConvR1D(RA, M, RB, N, RR2);
+                    end
+                    else
+                    begin
+                        RefConvR1DCircular(RA, M, RB, N, RR2);
+                    end;
+                    if CircKind=0 then
+                    begin
+                        I:=0;
+                        while I<=M+N-2 do
+                        begin
+                            RefRErr := Max(RefRErr, AbsReal(RR1[I]-RR2[I]));
+                            Inc(I);
+                        end;
+                    end
+                    else
+                    begin
+                        I:=0;
+                        while I<=M-1 do
+                        begin
+                            RefRErr := Max(RefRErr, AbsReal(RR1[I]-RR2[I]));
+                            Inc(I);
+                        end;
+                    end;
+                    Inc(RKind);
+                end;
+                Inc(CircKind);
+            end;
+            Inc(N);
+        end;
+        Inc(M);
+    end;
+    RefErrors := RefErrors or AP_FP_Greater(RefErr,ErrTol);
+    RefRErrors := RefRErrors or AP_FP_Greater(RefRErr,ErrTol);
+    
+    //
+    // Test inverse convolution
+    //
+    InvErr := 0;
+    InvRErr := 0;
+    M:=1;
+    while M<=MaxN do
+    begin
+        N:=1;
+        while N<=MaxN do
+        begin
+            
+            //
+            // Complex circilar and non-circular
+            //
+            SetLength(CA, M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                CA[I].X := 2*RandomReal-1;
+                CA[I].Y := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(CB, N);
+            I:=0;
+            while I<=N-1 do
+            begin
+                CB[I].X := 2*RandomReal-1;
+                CB[I].Y := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(CR1, 1);
+            SetLength(CR2, 1);
+            ConvC1D(CA, M, CB, N, CR2);
+            ConvC1DInv(CR2, M+N-1, CB, N, CR1);
+            I:=0;
+            while I<=M-1 do
+            begin
+                InvErr := Max(InvErr, AbsComplex(C_Sub(CR1[I],CA[I])));
+                Inc(I);
+            end;
+            SetLength(CR1, 1);
+            SetLength(CR2, 1);
+            ConvC1DCircular(CA, M, CB, N, CR2);
+            ConvC1DCircularInv(CR2, M, CB, N, CR1);
+            I:=0;
+            while I<=M-1 do
+            begin
+                InvErr := Max(InvErr, AbsComplex(C_Sub(CR1[I],CA[I])));
+                Inc(I);
+            end;
+            
+            //
+            // Real circilar and non-circular
+            //
+            SetLength(RA, M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                RA[I] := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(RB, N);
+            I:=0;
+            while I<=N-1 do
+            begin
+                RB[I] := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(RR1, 1);
+            SetLength(RR2, 1);
+            ConvR1D(RA, M, RB, N, RR2);
+            ConvR1DInv(RR2, M+N-1, RB, N, RR1);
+            I:=0;
+            while I<=M-1 do
+            begin
+                InvRErr := Max(InvRErr, AbsReal(RR1[I]-RA[I]));
+                Inc(I);
+            end;
+            SetLength(RR1, 1);
+            SetLength(RR2, 1);
+            ConvR1DCircular(RA, M, RB, N, RR2);
+            ConvR1DCircularInv(RR2, M, RB, N, RR1);
+            I:=0;
+            while I<=M-1 do
+            begin
+                InvRErr := Max(InvRErr, AbsReal(RR1[I]-RA[I]));
+                Inc(I);
+            end;
+            Inc(N);
+        end;
+        Inc(M);
+    end;
+    InvErrors := InvErrors or AP_FP_Greater(InvErr,ErrTol);
+    InvRErrors := InvRErrors or AP_FP_Greater(InvRErr,ErrTol);
+    
+    //
+    // end
+    //
+    WasErrors := RefErrors or RefRErrors or InvErrors or InvRErrors;
+    if  not Silent then
+    begin
+        Write(Format('TESTING CONVOLUTION'#13#10'',[]));
+        Write(Format('FINAL RESULT:                             ',[]));
+        if WasErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* AGAINST REFERENCE COMPLEX CONV:         ',[]));
+        if RefErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* AGAINST REFERENCE REAL CONV:            ',[]));
+        if RefRErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* COMPLEX INVERSE:                        ',[]));
+        if InvErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* REAL INVERSE:                           ',[]));
+        if InvRErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        if WasErrors then
+        begin
+            Write(Format('TEST FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('TEST PASSED'#13#10'',[]));
+        end;
+    end;
+    Result :=  not WasErrors;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefConvC1D(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+var
+    I : Integer;
+    V : Complex;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    SetLength(R, M+N-1);
+    I:=0;
+    while I<=M+N-2 do
+    begin
+        R[I] := C_Complex(0);
+        Inc(I);
+    end;
+    I:=0;
+    while I<=M-1 do
+    begin
+        V := A[I];
+        i1_ := (0) - (I);
+        for i_ := I to I+N-1 do
+        begin
+            R[i_] := C_Add(R[i_], C_Mul(V, B[i_+i1_]));
+        end;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefConvC1DCircular(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+var
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+    Buf : TComplex1DArray;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    RefConvC1D(A, M, B, N, Buf);
+    SetLength(R, M);
+    for i_ := 0 to M-1 do
+    begin
+        R[i_] := Buf[i_];
+    end;
+    I1 := M;
+    while I1<=M+N-2 do
+    begin
+        I2 := Min(I1+M-1, M+N-2);
+        J2 := I2-I1;
+        i1_ := (I1) - (0);
+        for i_ := 0 to J2 do
+        begin
+            R[i_] := C_Add(R[i_], Buf[i_+i1_]);
+        end;
+        I1 := I1+M;
+    end;
+end;
+
+
+(*************************************************************************
+Reference FFT
+*************************************************************************)
+procedure RefConvR1D(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+var
+    I : Integer;
+    V : Double;
+begin
+    SetLength(R, M+N-1);
+    I:=0;
+    while I<=M+N-2 do
+    begin
+        R[I] := 0;
+        Inc(I);
+    end;
+    I:=0;
+    while I<=M-1 do
+    begin
+        V := A[I];
+        APVAdd(@R[0], I, I+N-1, @B[0], 0, N-1, V);
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefConvR1DCircular(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+var
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+    Buf : TReal1DArray;
+begin
+    RefConvR1D(A, M, B, N, Buf);
+    SetLength(R, M);
+    APVMove(@R[0], 0, M-1, @Buf[0], 0, M-1);
+    I1 := M;
+    while I1<=M+N-2 do
+    begin
+        I2 := Min(I1+M-1, M+N-2);
+        J2 := I2-I1;
+        APVAdd(@R[0], 0, J2, @Buf[0], I1, I2);
+        I1 := I1+M;
+    end;
+end;
+
+
+(*************************************************************************
+Silent unit test
+*************************************************************************)
+function testconvunit_test_silent():Boolean;
+begin
+    Result := TestConv(True);
+end;
+
+
+(*************************************************************************
+Unit test
+*************************************************************************)
+function testconvunit_test():Boolean;
+begin
+    Result := TestConv(False);
+end;
+
+
+end.

tests/test/alglib/u_testcorrunit.pp

@@ -0,0 +1,551 @@
+unit u_testcorrunit;
+interface
+uses Math, Sysutils, u_ap, u_ftbase, u_fft, u_conv, u_corr;
+
+function TestCorr(Silent : Boolean):Boolean;
+function testcorrunit_test_silent():Boolean;
+function testcorrunit_test():Boolean;
+
+implementation
+
+procedure RefCorrC1D(const Signal : TComplex1DArray;
+     N : Integer;
+     const Pattern : TComplex1DArray;
+     M : Integer;
+     var R : TComplex1DArray);forward;
+procedure RefCorrC1DCircular(const Signal : TComplex1DArray;
+     N : Integer;
+     const Pattern : TComplex1DArray;
+     M : Integer;
+     var R : TComplex1DArray);forward;
+procedure RefCorrR1D(const Signal : TReal1DArray;
+     N : Integer;
+     const Pattern : TReal1DArray;
+     M : Integer;
+     var R : TReal1DArray);forward;
+procedure RefCorrR1DCircular(const Signal : TReal1DArray;
+     N : Integer;
+     const Pattern : TReal1DArray;
+     M : Integer;
+     var R : TReal1DArray);forward;
+procedure RefConvC1D(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);forward;
+procedure RefConvC1DCircular(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);forward;
+procedure RefConvR1D(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);forward;
+procedure RefConvR1DCircular(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);forward;
+
+
+(*************************************************************************
+Test
+*************************************************************************)
+function TestCorr(Silent : Boolean):Boolean;
+var
+    M : Integer;
+    N : Integer;
+    I : Integer;
+    RA : TReal1DArray;
+    RB : TReal1DArray;
+    RR1 : TReal1DArray;
+    RR2 : TReal1DArray;
+    CA : TComplex1DArray;
+    CB : TComplex1DArray;
+    CR1 : TComplex1DArray;
+    CR2 : TComplex1DArray;
+    MaxN : Integer;
+    RefErr : Double;
+    RefRErr : Double;
+    InvErr : Double;
+    InvRErr : Double;
+    ErrTol : Double;
+    RefErrors : Boolean;
+    RefRErrors : Boolean;
+    InvErrors : Boolean;
+    InvRErrors : Boolean;
+    WasErrors : Boolean;
+begin
+    MaxN := 32;
+    ErrTol := 100000*Power(MaxN, AP_Double(3)/2)*MachineEpsilon;
+    RefErrors := False;
+    RefRErrors := False;
+    InvErrors := False;
+    InvRErrors := False;
+    WasErrors := False;
+    
+    //
+    // Test against reference O(N^2) implementation.
+    //
+    RefErr := 0;
+    RefRErr := 0;
+    M:=1;
+    while M<=MaxN do
+    begin
+        N:=1;
+        while N<=MaxN do
+        begin
+            
+            //
+            // Complex correlation
+            //
+            SetLength(CA, M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                CA[I].X := 2*RandomReal-1;
+                CA[I].Y := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(CB, N);
+            I:=0;
+            while I<=N-1 do
+            begin
+                CB[I].X := 2*RandomReal-1;
+                CB[I].Y := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(CR1, 1);
+            CorrC1D(CA, M, CB, N, CR1);
+            RefCorrC1D(CA, M, CB, N, CR2);
+            I:=0;
+            while I<=M+N-2 do
+            begin
+                RefErr := Max(RefErr, AbsComplex(C_Sub(CR1[I],CR2[I])));
+                Inc(I);
+            end;
+            SetLength(CR1, 1);
+            CorrC1DCircular(CA, M, CB, N, CR1);
+            RefCorrC1DCircular(CA, M, CB, N, CR2);
+            I:=0;
+            while I<=M-1 do
+            begin
+                RefErr := Max(RefErr, AbsComplex(C_Sub(CR1[I],CR2[I])));
+                Inc(I);
+            end;
+            
+            //
+            // Real correlation
+            //
+            SetLength(RA, M);
+            I:=0;
+            while I<=M-1 do
+            begin
+                RA[I] := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(RB, N);
+            I:=0;
+            while I<=N-1 do
+            begin
+                RB[I] := 2*RandomReal-1;
+                Inc(I);
+            end;
+            SetLength(RR1, 1);
+            CorrR1D(RA, M, RB, N, RR1);
+            RefCorrR1D(RA, M, RB, N, RR2);
+            I:=0;
+            while I<=M+N-2 do
+            begin
+                RefRErr := Max(RefRErr, AbsReal(RR1[I]-RR2[I]));
+                Inc(I);
+            end;
+            SetLength(RR1, 1);
+            CorrR1DCircular(RA, M, RB, N, RR1);
+            RefCorrR1DCircular(RA, M, RB, N, RR2);
+            I:=0;
+            while I<=M-1 do
+            begin
+                RefRErr := Max(RefRErr, AbsReal(RR1[I]-RR2[I]));
+                Inc(I);
+            end;
+            Inc(N);
+        end;
+        Inc(M);
+    end;
+    RefErrors := RefErrors or AP_FP_Greater(RefErr,ErrTol);
+    RefRErrors := RefRErrors or AP_FP_Greater(RefRErr,ErrTol);
+    
+    //
+    // end
+    //
+    WasErrors := RefErrors or RefRErrors;
+    if  not Silent then
+    begin
+        Write(Format('TESTING CORRELATION'#13#10'',[]));
+        Write(Format('FINAL RESULT:                             ',[]));
+        if WasErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* AGAINST REFERENCE COMPLEX CORR:         ',[]));
+        if RefErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* AGAINST REFERENCE REAL CORR:            ',[]));
+        if RefRErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        if WasErrors then
+        begin
+            Write(Format('TEST FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('TEST PASSED'#13#10'',[]));
+        end;
+    end;
+    Result :=  not WasErrors;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefCorrC1D(const Signal : TComplex1DArray;
+     N : Integer;
+     const Pattern : TComplex1DArray;
+     M : Integer;
+     var R : TComplex1DArray);
+var
+    I : Integer;
+    J : Integer;
+    V : Complex;
+    S : TComplex1DArray;
+    i_ : Integer;
+begin
+    SetLength(S, M+N-1);
+    for i_ := 0 to N-1 do
+    begin
+        S[i_] := Signal[i_];
+    end;
+    I:=N;
+    while I<=M+N-2 do
+    begin
+        S[I] := C_Complex(0);
+        Inc(I);
+    end;
+    SetLength(R, M+N-1);
+    I:=0;
+    while I<=N-1 do
+    begin
+        V := C_Complex(0);
+        J:=0;
+        while J<=M-1 do
+        begin
+            if I+J>=N then
+            begin
+                Break;
+            end;
+            V := C_Add(V,C_Mul(Conj(Pattern[J]),S[I+J]));
+            Inc(J);
+        end;
+        R[I] := V;
+        Inc(I);
+    end;
+    I:=1;
+    while I<=M-1 do
+    begin
+        V := C_Complex(0);
+        J:=I;
+        while J<=M-1 do
+        begin
+            V := C_Add(V,C_Mul(Conj(Pattern[J]),S[J-I]));
+            Inc(J);
+        end;
+        R[M+N-1-I] := V;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefCorrC1DCircular(const Signal : TComplex1DArray;
+     N : Integer;
+     const Pattern : TComplex1DArray;
+     M : Integer;
+     var R : TComplex1DArray);
+var
+    I : Integer;
+    J : Integer;
+    V : Complex;
+begin
+    SetLength(R, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        V := C_Complex(0);
+        J:=0;
+        while J<=M-1 do
+        begin
+            V := C_Add(V,C_Mul(Conj(Pattern[J]),Signal[(I+J) mod N]));
+            Inc(J);
+        end;
+        R[I] := V;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefCorrR1D(const Signal : TReal1DArray;
+     N : Integer;
+     const Pattern : TReal1DArray;
+     M : Integer;
+     var R : TReal1DArray);
+var
+    I : Integer;
+    J : Integer;
+    V : Double;
+    S : TReal1DArray;
+begin
+    SetLength(S, M+N-1);
+    APVMove(@S[0], 0, N-1, @Signal[0], 0, N-1);
+    I:=N;
+    while I<=M+N-2 do
+    begin
+        S[I] := 0;
+        Inc(I);
+    end;
+    SetLength(R, M+N-1);
+    I:=0;
+    while I<=N-1 do
+    begin
+        V := 0;
+        J:=0;
+        while J<=M-1 do
+        begin
+            if I+J>=N then
+            begin
+                Break;
+            end;
+            V := V+Pattern[J]*S[I+J];
+            Inc(J);
+        end;
+        R[I] := V;
+        Inc(I);
+    end;
+    I:=1;
+    while I<=M-1 do
+    begin
+        V := 0;
+        J:=I;
+        while J<=M-1 do
+        begin
+            V := V+Pattern[J]*S[-I+J];
+            Inc(J);
+        end;
+        R[M+N-1-I] := V;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefCorrR1DCircular(const Signal : TReal1DArray;
+     N : Integer;
+     const Pattern : TReal1DArray;
+     M : Integer;
+     var R : TReal1DArray);
+var
+    I : Integer;
+    J : Integer;
+    V : Double;
+begin
+    SetLength(R, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        V := 0;
+        J:=0;
+        while J<=M-1 do
+        begin
+            V := V+Pattern[J]*Signal[(I+J) mod N];
+            Inc(J);
+        end;
+        R[I] := V;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefConvC1D(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+var
+    I : Integer;
+    V : Complex;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    SetLength(R, M+N-1);
+    I:=0;
+    while I<=M+N-2 do
+    begin
+        R[I] := C_Complex(0);
+        Inc(I);
+    end;
+    I:=0;
+    while I<=M-1 do
+    begin
+        V := A[I];
+        i1_ := (0) - (I);
+        for i_ := I to I+N-1 do
+        begin
+            R[i_] := C_Add(R[i_], C_Mul(V, B[i_+i1_]));
+        end;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefConvC1DCircular(const A : TComplex1DArray;
+     M : Integer;
+     const B : TComplex1DArray;
+     N : Integer;
+     var R : TComplex1DArray);
+var
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+    Buf : TComplex1DArray;
+    i_ : Integer;
+    i1_ : Integer;
+begin
+    RefConvC1D(A, M, B, N, Buf);
+    SetLength(R, M);
+    for i_ := 0 to M-1 do
+    begin
+        R[i_] := Buf[i_];
+    end;
+    I1 := M;
+    while I1<=M+N-2 do
+    begin
+        I2 := Min(I1+M-1, M+N-2);
+        J2 := I2-I1;
+        i1_ := (I1) - (0);
+        for i_ := 0 to J2 do
+        begin
+            R[i_] := C_Add(R[i_], Buf[i_+i1_]);
+        end;
+        I1 := I1+M;
+    end;
+end;
+
+
+(*************************************************************************
+Reference FFT
+*************************************************************************)
+procedure RefConvR1D(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+var
+    I : Integer;
+    V : Double;
+begin
+    SetLength(R, M+N-1);
+    I:=0;
+    while I<=M+N-2 do
+    begin
+        R[I] := 0;
+        Inc(I);
+    end;
+    I:=0;
+    while I<=M-1 do
+    begin
+        V := A[I];
+        APVAdd(@R[0], I, I+N-1, @B[0], 0, N-1, V);
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference implementation
+*************************************************************************)
+procedure RefConvR1DCircular(const A : TReal1DArray;
+     M : Integer;
+     const B : TReal1DArray;
+     N : Integer;
+     var R : TReal1DArray);
+var
+    I1 : Integer;
+    I2 : Integer;
+    J2 : Integer;
+    Buf : TReal1DArray;
+begin
+    RefConvR1D(A, M, B, N, Buf);
+    SetLength(R, M);
+    APVMove(@R[0], 0, M-1, @Buf[0], 0, M-1);
+    I1 := M;
+    while I1<=M+N-2 do
+    begin
+        I2 := Min(I1+M-1, M+N-2);
+        J2 := I2-I1;
+        APVAdd(@R[0], 0, J2, @Buf[0], I1, I2);
+        I1 := I1+M;
+    end;
+end;
+
+
+(*************************************************************************
+Silent unit test
+*************************************************************************)
+function testcorrunit_test_silent():Boolean;
+begin
+    Result := TestCorr(True);
+end;
+
+
+(*************************************************************************
+Unit test
+*************************************************************************)
+function testcorrunit_test():Boolean;
+begin
+    Result := TestCorr(False);
+end;
+
+
+end.

tests/test/alglib/u_testfftunit.pp

@@ -0,0 +1,557 @@
+unit u_testfftunit;
+interface
+uses Math, Sysutils, u_ap, u_ftbase, u_fft;
+
+function TestFFT(Silent : Boolean):Boolean;
+function testfftunit_test_silent():Boolean;
+function testfftunit_test():Boolean;
+
+implementation
+
+procedure RefFFTC1D(var A : TComplex1DArray; N : Integer);forward;
+procedure RefFFTC1DInv(var A : TComplex1DArray; N : Integer);forward;
+procedure RefInternalCFFT(var a : TReal1DArray;
+     nn : Integer;
+     InverseFFT : Boolean);forward;
+procedure RefInternalRFFT(const a : TReal1DArray;
+     nn : Integer;
+     var F : TComplex1DArray);forward;
+
+
+(*************************************************************************
+Test
+*************************************************************************)
+function TestFFT(Silent : Boolean):Boolean;
+var
+    N : Integer;
+    I : Integer;
+    K : Integer;
+    A1 : TComplex1DArray;
+    A2 : TComplex1DArray;
+    A3 : TComplex1DArray;
+    R1 : TReal1DArray;
+    R2 : TReal1DArray;
+    Buf : TReal1DArray;
+    Plan : FTPlan;
+    MaxN : Integer;
+    BiDiErr : Double;
+    BiDiRErr : Double;
+    RefErr : Double;
+    RefRErr : Double;
+    REIntErr : Double;
+    ErrTol : Double;
+    RefErrors : Boolean;
+    BiDiErrors : Boolean;
+    RefRErrors : Boolean;
+    BiDiRErrors : Boolean;
+    REIntErrors : Boolean;
+    WasErrors : Boolean;
+begin
+    MaxN := 128;
+    ErrTol := 100000*Power(MaxN, AP_Double(3)/2)*MachineEpsilon;
+    BiDiErrors := False;
+    RefErrors := False;
+    BiDiRErrors := False;
+    RefRErrors := False;
+    REIntErrors := False;
+    WasErrors := False;
+    
+    //
+    // Test bi-directional error: norm(x-invFFT(FFT(x)))
+    //
+    BiDiErr := 0;
+    BiDiRErr := 0;
+    N:=1;
+    while N<=MaxN do
+    begin
+        
+        //
+        // Complex FFT/invFFT
+        //
+        SetLength(A1, N);
+        SetLength(A2, N);
+        SetLength(A3, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            A1[I].X := 2*RandomReal-1;
+            A1[I].Y := 2*RandomReal-1;
+            A2[I] := A1[I];
+            A3[I] := A1[I];
+            Inc(I);
+        end;
+        FFTC1D(A2, N);
+        FFTC1DInv(A2, N);
+        FFTC1DInv(A3, N);
+        FFTC1D(A3, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            BiDiErr := Max(BiDiErr, AbsComplex(C_Sub(A1[I],A2[I])));
+            BiDiErr := Max(BiDiErr, AbsComplex(C_Sub(A1[I],A3[I])));
+            Inc(I);
+        end;
+        
+        //
+        // Real
+        //
+        SetLength(R1, N);
+        SetLength(R2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 2*RandomReal-1;
+            R2[I] := R1[I];
+            Inc(I);
+        end;
+        FFTR1D(R2, N, A1);
+        APVMul(@R2[0], 0, N-1, 0);
+        FFTR1DInv(A1, N, R2);
+        I:=0;
+        while I<=N-1 do
+        begin
+            BiDiRErr := Max(BiDiRErr, AbsComplex(C_Complex(R1[I]-R2[I])));
+            Inc(I);
+        end;
+        Inc(N);
+    end;
+    BiDiErrors := BiDiErrors or AP_FP_Greater(BiDiErr,ErrTol);
+    BiDiRErrors := BiDiRErrors or AP_FP_Greater(BiDiRErr,ErrTol);
+    
+    //
+    // Test against reference O(N^2) implementation
+    //
+    RefErr := 0;
+    RefRErr := 0;
+    N:=1;
+    while N<=MaxN do
+    begin
+        
+        //
+        // Complex FFT
+        //
+        SetLength(A1, N);
+        SetLength(A2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            A1[I].X := 2*RandomReal-1;
+            A1[I].Y := 2*RandomReal-1;
+            A2[I] := A1[I];
+            Inc(I);
+        end;
+        FFTC1D(A1, N);
+        RefFFTC1D(A2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            RefErr := Max(RefErr, AbsComplex(C_Sub(A1[I],A2[I])));
+            Inc(I);
+        end;
+        
+        //
+        // Complex inverse FFT
+        //
+        SetLength(A1, N);
+        SetLength(A2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            A1[I].X := 2*RandomReal-1;
+            A1[I].Y := 2*RandomReal-1;
+            A2[I] := A1[I];
+            Inc(I);
+        end;
+        FFTC1DInv(A1, N);
+        RefFFTC1DInv(A2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            RefErr := Max(RefErr, AbsComplex(C_Sub(A1[I],A2[I])));
+            Inc(I);
+        end;
+        
+        //
+        // Real forward/inverse FFT:
+        // * calculate and check forward FFT
+        // * use precalculated FFT to check backward FFT
+        //   fill unused parts of frequencies array with random numbers
+        //   to ensure that they are not really used
+        //
+        SetLength(R1, N);
+        SetLength(R2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 2*RandomReal-1;
+            R2[I] := R1[I];
+            Inc(I);
+        end;
+        FFTR1D(R1, N, A1);
+        RefInternalRFFT(R2, N, A2);
+        I:=0;
+        while I<=N-1 do
+        begin
+            RefRErr := Max(RefRErr, AbsComplex(C_Sub(A1[I],A2[I])));
+            Inc(I);
+        end;
+        SetLength(A3, Floor(AP_Double(N)/2)+1);
+        I:=0;
+        while I<=Floor(AP_Double(N)/2) do
+        begin
+            A3[I] := A2[I];
+            Inc(I);
+        end;
+        A3[0].Y := 2*RandomReal-1;
+        if N mod 2=0 then
+        begin
+            A3[Floor(AP_Double(N)/2)].Y := 2*RandomReal-1;
+        end;
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 0;
+            Inc(I);
+        end;
+        FFTR1DInv(A3, N, R1);
+        I:=0;
+        while I<=N-1 do
+        begin
+            RefRErr := Max(RefRErr, AbsReal(R2[I]-R1[I]));
+            Inc(I);
+        end;
+        Inc(N);
+    end;
+    RefErrors := RefErrors or AP_FP_Greater(RefErr,ErrTol);
+    RefRErrors := RefRErrors or AP_FP_Greater(RefRErr,ErrTol);
+    
+    //
+    // test internal real even FFT
+    //
+    REIntErr := 0;
+    K:=1;
+    while K<=MaxN div 2 do
+    begin
+        N := 2*K;
+        
+        //
+        // Real forward FFT
+        //
+        SetLength(R1, N);
+        SetLength(R2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 2*RandomReal-1;
+            R2[I] := R1[I];
+            Inc(I);
+        end;
+        FTBaseGenerateComplexFFTPlan(N div 2, Plan);
+        SetLength(Buf, N);
+        FFTR1DInternalEven(R1, N, Buf, Plan);
+        RefInternalRFFT(R2, N, A2);
+        REIntErr := Max(REIntErr, AbsReal(R1[0]-A2[0].X));
+        REIntErr := Max(REIntErr, AbsReal(R1[1]-A2[N div 2].X));
+        I:=1;
+        while I<=N div 2-1 do
+        begin
+            REIntErr := Max(REIntErr, AbsReal(R1[2*I+0]-A2[I].X));
+            REIntErr := Max(REIntErr, AbsReal(R1[2*I+1]-A2[I].Y));
+            Inc(I);
+        end;
+        
+        //
+        // Real backward FFT
+        //
+        SetLength(R1, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 2*RandomReal-1;
+            Inc(I);
+        end;
+        SetLength(A2, Floor(AP_Double(N)/2)+1);
+        A2[0] := C_Complex(R1[0]);
+        I:=1;
+        while I<=Floor(AP_Double(N)/2)-1 do
+        begin
+            A2[I].X := R1[2*I+0];
+            A2[I].Y := R1[2*I+1];
+            Inc(I);
+        end;
+        A2[Floor(AP_Double(N)/2)] := C_Complex(R1[1]);
+        FTBaseGenerateComplexFFTPlan(N div 2, Plan);
+        SetLength(Buf, N);
+        FFTR1DInvInternalEven(R1, N, Buf, Plan);
+        FFTR1DInv(A2, N, R2);
+        I:=0;
+        while I<=N-1 do
+        begin
+            REIntErr := Max(REIntErr, AbsReal(R1[I]-R2[I]));
+            Inc(I);
+        end;
+        Inc(K);
+    end;
+    REIntErrors := REIntErrors or AP_FP_Greater(REIntErr,ErrTol);
+    
+    //
+    // end
+    //
+    WasErrors := BiDiErrors or BiDiRErrors or RefErrors or RefRErrors or REIntErrors;
+    if  not Silent then
+    begin
+        Write(Format('TESTING FFT'#13#10'',[]));
+        Write(Format('FINAL RESULT:                             ',[]));
+        if WasErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* BI-DIRECTIONAL COMPLEX TEST:            ',[]));
+        if BiDiErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* AGAINST REFERENCE COMPLEX FFT:          ',[]));
+        if RefErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* BI-DIRECTIONAL REAL TEST:               ',[]));
+        if BiDiRErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* AGAINST REFERENCE REAL FFT:             ',[]));
+        if RefRErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* INTERNAL EVEN FFT:                      ',[]));
+        if REIntErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        if WasErrors then
+        begin
+            Write(Format('TEST FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('TEST PASSED'#13#10'',[]));
+        end;
+    end;
+    Result :=  not WasErrors;
+end;
+
+
+(*************************************************************************
+Reference FFT
+*************************************************************************)
+procedure RefFFTC1D(var A : TComplex1DArray; N : Integer);
+var
+    Buf : TReal1DArray;
+    I : Integer;
+begin
+    Assert(N>0, 'FFTC1D: incorrect N!');
+    SetLength(Buf, 2*N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        Buf[2*I+0] := A[I].X;
+        Buf[2*I+1] := A[I].Y;
+        Inc(I);
+    end;
+    RefInternalCFFT(Buf, N, False);
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I].X := Buf[2*I+0];
+        A[I].Y := Buf[2*I+1];
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference inverse FFT
+*************************************************************************)
+procedure RefFFTC1DInv(var A : TComplex1DArray; N : Integer);
+var
+    Buf : TReal1DArray;
+    I : Integer;
+begin
+    Assert(N>0, 'FFTC1DInv: incorrect N!');
+    SetLength(Buf, 2*N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        Buf[2*I+0] := A[I].X;
+        Buf[2*I+1] := A[I].Y;
+        Inc(I);
+    end;
+    RefInternalCFFT(Buf, N, True);
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I].X := Buf[2*I+0];
+        A[I].Y := Buf[2*I+1];
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Internal complex FFT stub.
+Uses straightforward formula with O(N^2) complexity.
+*************************************************************************)
+procedure RefInternalCFFT(var a : TReal1DArray;
+     nn : Integer;
+     InverseFFT : Boolean);
+var
+    Tmp : TReal1DArray;
+    I : Integer;
+    J : Integer;
+    K : Integer;
+    Hre : Double;
+    Him : Double;
+    C : Double;
+    S : Double;
+    re : Double;
+    im : Double;
+begin
+    SetLength(Tmp, 2*nn-1+1);
+    if  not InverseFFT then
+    begin
+        i:=0;
+        while i<=nn-1 do
+        begin
+            HRe := 0;
+            Him := 0;
+            K:=0;
+            while K<=nn-1 do
+            begin
+                re := A[2*k];
+                im := A[2*k+1];
+                C := Cos(-2*Pi*k*i/nn);
+                S := Sin(-2*Pi*k*i/nn);
+                HRe := HRe+C*Re-S*im;
+                HIm := HIm+C*im+S*re;
+                Inc(K);
+            end;
+            Tmp[2*i] := Hre;
+            Tmp[2*i+1] := Him;
+            Inc(i);
+        end;
+        i:=0;
+        while i<=2*nn-1 do
+        begin
+            A[I] := Tmp[I];
+            Inc(i);
+        end;
+    end
+    else
+    begin
+        k:=0;
+        while k<=nn-1 do
+        begin
+            HRe := 0;
+            Him := 0;
+            i:=0;
+            while i<=nn-1 do
+            begin
+                re := A[2*i];
+                im := A[2*i+1];
+                C := Cos(2*Pi*k*i/nn);
+                S := Sin(2*Pi*k*i/nn);
+                HRe := HRe+C*Re-S*im;
+                HIm := HIm+C*im+S*re;
+                Inc(i);
+            end;
+            Tmp[2*k] := Hre/nn;
+            Tmp[2*k+1] := Him/nn;
+            Inc(k);
+        end;
+        i:=0;
+        while i<=2*nn-1 do
+        begin
+            A[I] := Tmp[I];
+            Inc(i);
+        end;
+    end;
+end;
+
+
+(*************************************************************************
+Internal real FFT stub.
+Uses straightforward formula with O(N^2) complexity.
+*************************************************************************)
+procedure RefInternalRFFT(const a : TReal1DArray;
+     nn : Integer;
+     var F : TComplex1DArray);
+var
+    Tmp : TReal1DArray;
+    I : Integer;
+begin
+    SetLength(Tmp, 2*nn-1+1);
+    I:=0;
+    while I<=nn-1 do
+    begin
+        Tmp[2*I] := A[I];
+        Tmp[2*I+1] := 0;
+        Inc(I);
+    end;
+    RefInternalCFFT(Tmp, nn, False);
+    SetLength(F, nn);
+    I:=0;
+    while I<=nn-1 do
+    begin
+        F[I].X := Tmp[2*I+0];
+        F[I].Y := Tmp[2*I+1];
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Silent unit test
+*************************************************************************)
+function testfftunit_test_silent():Boolean;
+begin
+    Result := TestFFT(True);
+end;
+
+
+(*************************************************************************
+Unit test
+*************************************************************************)
+function testfftunit_test():Boolean;
+begin
+    Result := TestFFT(False);
+end;
+
+
+end.

tests/test/alglib/u_testfhtunit.pp

@@ -0,0 +1,246 @@
+unit u_testfhtunit;
+interface
+uses Math, Sysutils, u_ap, u_ftbase, u_fft, u_fht;
+
+function TestFHT(Silent : Boolean):Boolean;
+function testfhtunit_test_silent():Boolean;
+function testfhtunit_test():Boolean;
+
+implementation
+
+procedure RefFHTR1D(var A : TReal1DArray; N : Integer);forward;
+procedure RefFHTR1DInv(var A : TReal1DArray; N : Integer);forward;
+
+
+(*************************************************************************
+Test
+*************************************************************************)
+function TestFHT(Silent : Boolean):Boolean;
+var
+    N : Integer;
+    I : Integer;
+    R1 : TReal1DArray;
+    R2 : TReal1DArray;
+    R3 : TReal1DArray;
+    MaxN : Integer;
+    BiDiErr : Double;
+    RefErr : Double;
+    ErrTol : Double;
+    RefErrors : Boolean;
+    BiDiErrors : Boolean;
+    WasErrors : Boolean;
+begin
+    MaxN := 128;
+    ErrTol := 100000*Power(MaxN, AP_Double(3)/2)*MachineEpsilon;
+    BiDiErrors := False;
+    RefErrors := False;
+    WasErrors := False;
+    
+    //
+    // Test bi-directional error: norm(x-invFHT(FHT(x)))
+    //
+    BiDiErr := 0;
+    N:=1;
+    while N<=MaxN do
+    begin
+        
+        //
+        // FHT/invFHT
+        //
+        SetLength(R1, N);
+        SetLength(R2, N);
+        SetLength(R3, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 2*RandomReal-1;
+            R2[I] := R1[I];
+            R3[I] := R1[I];
+            Inc(I);
+        end;
+        FHTR1D(R2, N);
+        FHTR1DInv(R2, N);
+        FHTR1DInv(R3, N);
+        FHTR1D(R3, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            BiDiErr := Max(BiDiErr, AbsReal(R1[I]-R2[I]));
+            BiDiErr := Max(BiDiErr, AbsReal(R1[I]-R3[I]));
+            Inc(I);
+        end;
+        Inc(N);
+    end;
+    BiDiErrors := BiDiErrors or AP_FP_Greater(BiDiErr,ErrTol);
+    
+    //
+    // Test against reference O(N^2) implementation
+    //
+    RefErr := 0;
+    N:=1;
+    while N<=MaxN do
+    begin
+        
+        //
+        // FHT
+        //
+        SetLength(R1, N);
+        SetLength(R2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 2*RandomReal-1;
+            R2[I] := R1[I];
+            Inc(I);
+        end;
+        FHTR1D(R1, N);
+        RefFHTR1D(R2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            RefErr := Max(RefErr, AbsReal(R1[I]-R2[I]));
+            Inc(I);
+        end;
+        
+        //
+        // inverse FHT
+        //
+        SetLength(R1, N);
+        SetLength(R2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            R1[I] := 2*RandomReal-1;
+            R2[I] := R1[I];
+            Inc(I);
+        end;
+        FHTR1DInv(R1, N);
+        RefFHTR1DInv(R2, N);
+        I:=0;
+        while I<=N-1 do
+        begin
+            RefErr := Max(RefErr, AbsReal(R1[I]-R2[I]));
+            Inc(I);
+        end;
+        Inc(N);
+    end;
+    RefErrors := RefErrors or AP_FP_Greater(RefErr,ErrTol);
+    
+    //
+    // end
+    //
+    WasErrors := BiDiErrors or RefErrors;
+    if  not Silent then
+    begin
+        Write(Format('TESTING FHT'#13#10'',[]));
+        Write(Format('FINAL RESULT:                             ',[]));
+        if WasErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* BI-DIRECTIONAL TEST:                    ',[]));
+        if BiDiErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        Write(Format('* AGAINST REFERENCE FHT:                  ',[]));
+        if RefErrors then
+        begin
+            Write(Format('FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('OK'#13#10'',[]));
+        end;
+        if WasErrors then
+        begin
+            Write(Format('TEST FAILED'#13#10'',[]));
+        end
+        else
+        begin
+            Write(Format('TEST PASSED'#13#10'',[]));
+        end;
+    end;
+    Result :=  not WasErrors;
+end;
+
+
+(*************************************************************************
+Reference FHT
+*************************************************************************)
+procedure RefFHTR1D(var A : TReal1DArray; N : Integer);
+var
+    Buf : TReal1DArray;
+    I : Integer;
+    J : Integer;
+    V : Double;
+begin
+    Assert(N>0, 'RefFHTR1D: incorrect N!');
+    SetLength(Buf, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        V := 0;
+        J:=0;
+        while J<=N-1 do
+        begin
+            V := V+A[J]*(Cos(2*Pi*I*J/N)+Sin(2*Pi*I*J/N));
+            Inc(J);
+        end;
+        Buf[I] := V;
+        Inc(I);
+    end;
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I] := Buf[I];
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Reference inverse FHT
+*************************************************************************)
+procedure RefFHTR1DInv(var A : TReal1DArray; N : Integer);
+var
+    I : Integer;
+begin
+    Assert(N>0, 'RefFHTR1DInv: incorrect N!');
+    RefFHTR1D(A, N);
+    I:=0;
+    while I<=N-1 do
+    begin
+        A[I] := A[I]/N;
+        Inc(I);
+    end;
+end;
+
+
+(*************************************************************************
+Silent unit test
+*************************************************************************)
+function testfhtunit_test_silent():Boolean;
+begin
+    Result := TestFHT(True);
+end;
+
+
+(*************************************************************************
+Unit test
+*************************************************************************)
+function testfhtunit_test():Boolean;
+begin
+    Result := TestFHT(False);
+end;
+
+
+end.