PasMatLib (Delphi):
http://www.shokhirev.com/nikolai/programs/tools/PasMatLib/userguide.html
dlib (C++):
http://dclib.sourceforge.net/compile.html
ALGLIB:
(*************************************************************************
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 Ap;
interface
uses Math, Sysutils;
/////////////////////////////////////////////////////////////////////////
// constants
/////////////////////////////////////////////////////////////////////////
const
MachineEpsilon = 5E-16;
MaxRealNumber = 1E300;
MinRealNumber = 1E-300;
/////////////////////////////////////////////////////////////////////////
// arrays
/////////////////////////////////////////////////////////////////////////
type
AlglibInteger = LongInt;
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 : AlglibInteger):AlglibInteger;
function RandomReal():Double;
function RandomInteger(I : AlglibInteger):AlglibInteger;
function Sign(X:Double):AlglibInteger;
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: AlglibInteger;
V2: PDouble; I21, I22: AlglibInteger):Double;
procedure APVMove(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);overload;
procedure APVMove(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger;
S: Double);overload;
procedure APVMoveNeg(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);
procedure APVAdd(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);overload;
procedure APVAdd(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger;
S: Real);overload;
procedure APVSub(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);overload;
procedure APVSub(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger;
S: Real);overload;
procedure APVMul(
VOp: PDouble; I1, I2: AlglibInteger;
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: AlglibInteger):Double;cdecl;
ASMMove1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);cdecl;
ASMMoveS1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger; S: Double);cdecl;
ASMMoveNeg1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);cdecl;
ASMAdd1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);cdecl;
ASMAddS1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger; S: Double);cdecl;
ASMSub1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);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: AlglibInteger):Double;cdecl;
ASMMove1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);cdecl;
ASMMoveS1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger; S: Double);cdecl;
ASMMoveNeg1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);cdecl;
ASMAdd1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);cdecl;
ASMAddS1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger; S: Double);cdecl;
ASMSub1: procedure(VDst: PDouble; VSrc: PDouble; N: AlglibInteger);cdecl;}
{$ENDIF}
/////////////////////////////////////////////////////////////////////////
// Functions/procedures
/////////////////////////////////////////////////////////////////////////
function AbsReal(X : Double):Double;
begin
//Result := Abs(X);
if X>=0 then
AbsReal:=X
else
AbsReal:=-X;
end;
function AbsInt (I : AlglibInteger):AlglibInteger;
begin
//Result := Abs(I);
if I>=0 then
AbsInt:=I
else
AbsInt:=-I;
end;
function RandomReal():Double;
begin
RandomReal:=Random;
end;
function RandomInteger(I : AlglibInteger):AlglibInteger;
begin
RandomInteger:=Random(I);
end;
function Sign(X:Double):AlglibInteger;
begin
if X>0 then
Sign := 1
else if X<0 then
Sign := -1
else
Sign := 0;
end;
function AP_Sqr(X:Double):Double;
begin
AP_Sqr := X*X;
end;
/////////////////////////////////////////////////////////////////////////
// dynamical arrays copying
/////////////////////////////////////////////////////////////////////////
function DynamicArrayCopy(const A: TInteger1DArray):TInteger1DArray;overload;
var
I: AlglibInteger;
R: TInteger1DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
R[I]:=A[I];
DynamicArrayCopy:=R;
end;
function DynamicArrayCopy(const A: TReal1DArray):TReal1DArray;overload;
var
I: AlglibInteger;
R: TReal1DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
R[I]:=A[I];
DynamicArrayCopy:=R;
end;
function DynamicArrayCopy(const A: TComplex1DArray):TComplex1DArray;overload;
var
I: AlglibInteger;
R: TComplex1DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
R[I]:=A[I];
DynamicArrayCopy:=R;
end;
function DynamicArrayCopy(const A: TBoolean1DArray):TBoolean1DArray;overload;
var
I: AlglibInteger;
R: TBoolean1DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
R[I]:=A[I];
DynamicArrayCopy:=R;
end;
function DynamicArrayCopy(const A: TInteger2DArray):TInteger2DArray;overload;
var
I,J: AlglibInteger;
R: TInteger2DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
begin
SetLength(R[I], High(A[I])+1);
for J:=Low(A[I]) to High(A[I]) do
R[I,J]:=A[I,J];
end;
DynamicArrayCopy:=R;
end;
function DynamicArrayCopy(const A: TReal2DArray):TReal2DArray;overload;
var
I,J: AlglibInteger;
R: TReal2DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
begin
SetLength(R[I], High(A[I])+1);
for J:=Low(A[I]) to High(A[I]) do
R[I,J]:=A[I,J];
end;
DynamicArrayCopy:=R;
end;
function DynamicArrayCopy(const A: TComplex2DArray):TComplex2DArray;overload;
var
I,J: AlglibInteger;
R: TComplex2DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
begin
SetLength(R[I], High(A[I])+1);
for J:=Low(A[I]) to High(A[I]) do
R[I,J]:=A[I,J];
end;
DynamicArrayCopy:=R;
end;
function DynamicArrayCopy(const A: TBoolean2DArray):TBoolean2DArray;overload;
var
I,J: AlglibInteger;
R: TBoolean2DArray;
begin
SetLength(R, High(A)+1);
for I:=Low(A) to High(A) do
begin
SetLength(R[I], High(A[I])+1);
for J:=Low(A[I]) to High(A[I]) do
R[I,J]:=A[I,J];
end;
DynamicArrayCopy:=R;
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
AbsComplex := W;
end
else
begin
AbsComplex := W*SQRT(1+Sqr(V/W));
end;
end;
function Conj(const Z : Complex):Complex;
begin
Conj.X := Z.X;
Conj.Y := -Z.Y;
end;
function CSqr(const Z : Complex):Complex;
begin
CSqr.X := Sqr(Z.X)-Sqr(Z.Y);
CSqr.Y := 2*Z.X*Z.Y;
end;
function C_Complex(const X : Double):Complex;
begin
C_Complex.X := X;
C_Complex.Y := 0;
end;
function C_Opposite(const Z : Complex):Complex;
begin
C_Opposite.X := -Z.X;
C_Opposite.Y := -Z.Y;
end;
function C_Add(const Z1 : Complex; const Z2 : Complex):Complex;
begin
C_Add.X := Z1.X+Z2.X;
C_Add.Y := Z1.Y+Z2.Y;
end;
function C_Mul(const Z1 : Complex; const Z2 : Complex):Complex;
begin
C_Mul.X := Z1.X*Z2.X-Z1.Y*Z2.Y;
C_Mul.Y := Z1.X*Z2.Y+Z1.Y*Z2.X;
end;
function C_AddR(const Z1 : Complex; const R : Double):Complex;
begin
C_AddR.X := Z1.X+R;
C_AddR.Y := Z1.Y;
end;
function C_MulR(const Z1 : Complex; const R : Double):Complex;
begin
C_MulR.X := Z1.X*R;
C_MulR.Y := Z1.Y*R;
end;
function C_Sub(const Z1 : Complex; const Z2 : Complex):Complex;
begin
C_Sub.X := Z1.X-Z2.X;
C_Sub.Y := Z1.Y-Z2.Y;
end;
function C_SubR(const Z1 : Complex; const R : Double):Complex;
begin
C_SubR.X := Z1.X-R;
C_SubR.Y := Z1.Y;
end;
function C_RSub(const R : Double; const Z1 : Complex):Complex;
begin
C_RSub.X := R-Z1.X;
C_RSub.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;
C_Div.X := (A+B*E)/F;
C_Div.Y := (B-A*E)/F;
end
else
begin
E := C/D;
F := D+C*E;
C_Div.X := (B+A*E)/F;
C_Div.Y := (-A+B*E)/F;
end;
end;
function C_DivR(const Z1 : Complex; const R : Double):Complex;
begin
C_DivR.X := Z1.X/R;
C_DivR.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;
C_RDiv.X := A/F;
C_RDiv.Y := -A*E/F;
end
else
begin
E := C/D;
F := D+C*E;
C_RDiv.X := A*E/F;
C_RDiv.Y := -A/F;
end;
end;
function C_Equal(const Z1 : Complex; const Z2 : Complex):Boolean;
begin
C_Equal := (Z1.X=Z2.X) and (Z1.Y=Z2.Y);
end;
function C_NotEqual(const Z1 : Complex; const Z2 : Complex):Boolean;
begin
C_NotEqual := (Z1.X<>Z2.X) or (Z1.Y<>Z2.Y);
end;
function C_EqualR(const Z1 : Complex; const R : Double):Boolean;
begin
C_EqualR := (Z1.X=R) and (Z1.Y=0);
end;
function C_NotEqualR(const Z1 : Complex; const R : Double):Boolean;
begin
C_NotEqualR := (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: AlglibInteger;
V2: PDouble; I21, I22: AlglibInteger):Double;
var
I, C: LongInt;
R: Double;
begin
Assert(I12-I11=I22-I21, 'APVDotProduct: arrays of different size!');
Inc(V1, I11);
Inc(V2, I21);
//
// Generic pascal code
//
C:=I12-I11;
R:=0;
for I:=0 to C do
begin
R:=R+V1^*V2^;
Inc(V1);
Inc(V2);
end;
APVDotProduct:=R;
end;
procedure APVMove(
VDst: PDouble; I11, I12: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);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: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger;
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: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);
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: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);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: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger;
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: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger);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: AlglibInteger;
VSrc: PDouble; I21, I22: AlglibInteger;
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: AlglibInteger;
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
AP_Double:=X;
end;
function AP_FP_Eq(X:Double; Y:Double):Boolean;
begin
AP_FP_Eq:=X=Y;
end;
function AP_FP_NEq(X:Double; Y:Double):Boolean;
begin
AP_FP_NEq:=X<>Y;
end;
function AP_FP_Less(X:Double; Y:Double):Boolean;
begin
AP_FP_Less:=X<Y;
end;
function AP_FP_Less_Eq(X:Double; Y:Double):Boolean;
begin
AP_FP_Less_Eq:=X<=Y;
end;
function AP_FP_Greater(X:Double; Y:Double):Boolean;
begin
AP_FP_Greater:=X>Y;
end;
function AP_FP_Greater_Eq(X:Double; Y:Double):Boolean;
begin
AP_FP_Greater_Eq:=X>=Y;
end;
end.