 | 
|
|
|
{$INCLUDE ..\cDefines.inc}
unit cMaths;
{ }
{ Maths unit v3.61 }
{ }
{ A collection of mathematical routines. }
{ }
{ This unit is copyright © 1995-2003 by David J Butler }
{ }
{ This unit is part of Delphi Fundamentals. }
{ Its original file name is cMaths.pas }
{ The latest version is available from the Fundamentals home page }
{ http://fundementals.sourceforge.net/ }
{ }
{ I invite you to use this unit, free of charge. }
{ I invite you to distibute this unit, but it must be for free. }
{ I also invite you to contribute to its development, }
{ but do not distribute a modified copy of this file. }
{ }
{ A forum is available on SourceForge for general discussion }
{ http://sourceforge.net/forum/forum.php?forum_id=2117 }
{ }
{ Thanks to: }
{ Stephen L. Moshier, Guy Hirson, Ondrej Hrabal, Andrew Driazgov, }
{ Torsten Wasch, Jingyi Peng, Michael Schuemann, Sebastian Schuberth. }
{ }
{ Revision history: }
{ 1995-96 0.01 Wrote statistical and actuarial functions. }
{ 1998/03 0.02 Added Solver and SecantSolver. }
{ 1998/10 0.03 Removed functions now found in Delphi 3's Math unit }
{ Uses Delphi's math exceptions (eg EOverflow, }
{ EInvalidArgument, etc) }
{ 1999/08/29 0.04 Added ASet, TRangeSet, TFlatSet, TSparseFlatSet. }
{ 1999/09/27 0.05 Added TMatrix, TVector. }
{ Added BinomialCoeff. }
{ 1999/10/02 0.06 Added TComplex. }
{ 1999/10/03 0.07 Added DerivateSolvers. }
{ Completed TMatrix. }
{ 1999/10/04 0.08 Added TLifeTable. }
{ 1999/10/05 0.09 T3DPoint }
{ 1999/10/06 0.10 Transform matrices. }
{ 1999/10/13 0.11 TRay, TSphere }
{ 1999/10/14 0.12 TPlane }
{ 1999/10/26 1.13 Upgraded to Delphi 4. Compared the assembly code of the }
{ new dynamic arrays with that of pointers to arrays. }
{ 1999/10/30 1.14 Added TVector.StdDev }
{ Changed some functions to the same name (since Delphi }
{ now supports overloading). }
{ Removed Min and Max functions (now in Math). }
{ 1999/11/04 1.15 Added TVector.Pos, TVector.Append }
{ 1999/11/07 1.16 Added RandomSeed function. }
{ Added assembly bit functions. }
{ 1999/11/10 1.17 Added hashing functions. Coded XOR8 in assembly. }
{ Added MD5 hashing. }
{ 1999/11/11 1.18 Added EncodeBase. }
{ 1999/11/21 1.19 Added TComplex.Power }
{ 1999/11/25 1.20 Moved TRay, TSphere and TPlane to cRayTrace. }
{ Added Primes. }
{ 1999/11/26 1.21 Added Rational numbers (can convert to/from TReal). }
{ Added GCD. }
{ Added RealArray/IntegerArray functions. }
{ 1999/11/27 1.22 Replaced GCD algorithm with Euclid's algorithm. }
{ Added SI constants. }
{ Added TMatrix.Normalise, TMatrix.Multiply (Row, Value) }
{ 1999/12/01 1.23 Added RandomUniform. }
{ 1999/12/03 1.24 Added RandomNormal. }
{ 1999/12/16 1.25 Bug fixes. }
{ 1999/12/23 1.26 Fixed bug in TRational.CalcFrac. }
{ 1999/12/26 1.27 Added Reverse procedures for RealArray/IntegerArray. }
{ 2000/01/22 1.28 Added TStatistic. }
{ 2000/01/23 1.29 Added RandomPseudoWord. }
{ 2000/03/08 1.30 Moved TInteger/TReal, IntegerArray/RealArray to cUtil. }
{ Moved AArray to cDataStructs. }
{ 2000/04/01 1.31 Moved ASet and set implmentations to cDataStructs. }
{ 2000/04/09 1.32 Added SetFPUPrecision. }
{ 2000/05/03 1.33 Added Bit functions (ToggleBit, SetBit, ClearBit, }
{ IsBitSet). }
{ 2000/05/08 1.34 Moved SetFPUPrecision to cUtil. }
{ 2000/05/25 1.35 Started THugeInteger. }
{ 2000/06/06 1.36 Moved bit functions to cUtil. }
{ 2000/06/08 1.37 TVector now inherits from TExtendedArray. }
{ Added TIntegerVector. }
{ Removed TInteger/TReal, TIntegerArray/TRealArray. }
{ 2000/06/16 1.38 Updated documentation for financial functions. }
{ 2000/06/17 1.39 Recalculated all constants to 34 digits. }
{ 2000/06/25 1.40 Fixed bug in CalcCheckSum32 reported by Ondrej Hrabal }
{ <twilight.ia@worldonline.cz> }
{ 2000/07/13 1.41 Fixed bug in RandomUniform reported by Andrew Driazgov }
{ <andrey@asp.tstu.ru> }
{ 2000/07/22 1.42 Replaced Fibonacci function with a much quicker one by }
{ Torsten Wasch (Torsten.Wasch@Post.RWTH-Aachen.de). }
{ 2000/08/04 1.43 Added TMovingStatistic. }
{ 2000/08/22 1.44 Added RandomHex. }
{ 2000/09/20 1.45 More random states to RandomSeed. }
{ 2000/09/23 1.46 Added EncodeBase64. }
{ 2000/10/03 1.47 Fixed bug in TMatrix.Transposed reported by Jingyi Peng }
{ <pengj@rotellacapital.com> }
{ 2000/10/22 1.48 Added CalcKeyedMD5 (RFC 2104 - HMAC) }
{ 2001/05/14 1.49 Moved RandomSeed function to cSysUtils. }
{ 2001/05/19 1.50 Moved Hashing functions to cUtils. }
{ 2001/05/21 1.51 Moved TTRational and TTComplex from cExDataStructs. }
{ 2001/07/31 1.52 Minor Revisions and bug fixes. }
{ 2001/08/18 1.53 Changed Sgn function to return 0 for 0. }
{ Changed Median functions to be average of two middle }
{ values for even number of elements and the geometric }
{ mean to be calculated using logs (as suggested by }
{ Michael Schuemann <schuemann@uke.uni-hamburg.de>) }
{ 2001/11/18 1.54 Added FourierTransform functions. }
{ Added HugeWord functions for Add, Subtract, Invert, }
{ MultiplyLong and Divide. }
{ 2001/11/19 2.55 Added HugeWord function for MultiplyFFT. }
{ Added THugeInteger. }
{ 2002/01/03 2.56 Moved RandomUniform, RandomPseudoWord to cSysUtils. }
{ Moved EncodeBase64, DecodeBase64 to cUtils. }
{ 2002/01/31 2.57 Fixed GCD function to always return positive value. }
{ Thanks to Sebastian Schuberth <s.schuberth@tu-bs.de>. }
{ 2002/02/03 2.58 Added InvMod, ExpMod and Jacobi, contribution by }
{ Sebastian Schuberth <s.schuberth@tu-bs.de>. }
{ 2002/06/01 2.59 Created cHugeInt, cRational, cComplex, cVectors, }
{ cMatrix and c3DMath units from cMaths. }
{ 2003/02/16 3.60 Created cStatistics unit from cMaths. }
{ Revised for Fundamentals 3. }
{ 2003/03/14 3.61 Added FPU functions. }
{ Refactored Prime functions. }
{ }
interface
uses
{ Delphi }
SysUtils,
{ Fundamentals }
cUtils;
const
UnitName = 'cMaths';
UnitVersion = '3.61';
UnitDesc = 'Mathematical functions';
UnitCopyright = '(c) 1995-2003 by David J Butler';
{ }
{ Mathematical constants }
{ }
const
Pi = 3.14159265358979323846 + { Pi (200 digits) }
0.26433832795028841971e-20 + { }
0.69399375105820974944e-40 + { }
0.59230781640628620899e-60 + { }
0.86280348253421170679e-80 + { }
0.82148086513282306647e-100 + { }
0.09384460955058223172e-120 + { }
0.53594081284811174502e-140 + { }
0.84102701938521105559e-160 + { }
0.64462294895493038196e-180; { }
Pi2 = 6.283185307179586476925286766559006; { Pi * 2 }
Pi3 = 9.424777960769379715387930149838509; { Pi * 3 }
Pi4 = 12.56637061435917295385057353311801; { Pi * 4 }
PiOn2 = 1.570796326794896619231321691639751; { Pi / 2 }
PiOn3 = 1.047197551196597746154214461093168; { Pi / 3 }
PiOn4 = 0.785398163397448309615660845819876; { Pi / 4 }
PiSq = 9.869604401089358618834490999876151; { Pi^2 }
PiE = 22.45915771836104547342715220454374; { Pi^e }
LnPi = 1.144729885849400174143427351353059; { Ln (Pi) }
LogPi = 0.497149872694133854351268288290899; { Log (Pi) }
SqrtPi = 1.772453850905516027298167483341145; { Sqrt (Pi) }
Sqrt2Pi = 2.506628274631000502415765284811045; { Sqrt (2 * Pi) }
LnSqrt2Pi = 0.918938533204672741780329736405618; { Ln (Sqrt (2 * Pi)) }
RadPerDeg = 0.017453292519943295769236907684886; { Pi / 180 }
DegPerRad = 57.29577951308232087679815481410517; { 180 / Pi }
E = 2.718281828459045235360287471352663; { e }
E2 = 7.389056098930650227230427460575008; { e^2 }
ExpM2 = 0.135335283236612691893999494972484; { e^-2 }
Ln2 = 0.693147180559945309417232121458177; { Ln (2) }
Ln10 = 2.302585092994045684017991454684364; { Ln (10) }
LogE = 0.434294481903251827651128918916605; { Log (e) }
Log2 = 0.301029995663981195213738894724493; { Log (2) }
Log3 = 0.477121254719662437295027903255115; { Log (3) }
Sqrt2 = 1.414213562373095048801688724209698; { Sqrt (2) }
Sqrt3 = 1.732050807568877293527446341505872; { Sqrt (3) }
Sqrt5 = 2.236067977499789696409173668731276; { Sqrt (5) }
Sqrt7 = 2.645751311064590590501615753639260; { Sqrt (7) }
{ }
{ FPU (Floating Point Unit) functions }
{ }
procedure SetFPUPrecisionSingle;
procedure SetFPUPrecisionDouble;
procedure SetFPUPrecisionExtended;
procedure SetFPURoundingNearest;
procedure SetFPURoundingDown;
procedure SetFPURoundingUp;
procedure SetFPURoundingTruncate;
{ }
{ Polynomial }
{ }
{ Calculates polynomial of degree N. }
{ }
{ Note that coefficients are stored in reverse order, ie Coef[0] = C }
{ N }
{ 2 N }
{ y = C + C x + C x +...+ C x }
{ 0 1 2 N }
{ }
function PolyEval(const X: Extended; const Coef: Array of Extended;
const N: Integer): Extended;
{ }
{ Miscellaneous functions }
{ }
{ Sgn returns the sign of an argument, +1, -1 or 0. }
{ }
{ FloatMod calculates the floating-point value of A mod B. }
{ }
function Sgn(const R: Integer): Integer; overload;
function Sgn(const R: Int64): Integer; overload;
function Sgn(const R: Single): Integer; overload;
function Sgn(const R: Double): Integer; overload;
function Sgn(const R: Extended): Integer; overload;
function FloatMod(const A, B: Extended): Extended;
{ }
{ Trigonometric functions }
{ }
{ Delphi's Math unit includes most commonly used trigonometric functions. }
{ }
function ATan360(const X, Y: Extended): Extended;
function InverseTangentDeg(const X, Y: Extended): Extended;
function InverseTangentRad(const X, Y: Extended): Extended;
function InverseSinDeg(const Y, R: Extended): Extended;
function InverseSinRad(const Y, R: Extended): Extended;
function InverseCosDeg(const X, R: Extended): Extended;
function InverseCosRad(const X, R: Extended): Extended;
function DMSToReal(const Degs, Mins, Secs: Extended): Extended;
procedure RealToDMS(const X: Extended; var Degs, Mins, Secs: Extended);
procedure PolarToRectangular(const R, Theta: Extended; var X, Y: Extended);
procedure RectangularToPolar(const X, Y: Extended; var R, Theta: Extended);
function Distance(const X1, Y1, X2, Y2: Extended): Extended;
{ }
{ Prime numbers }
{ }
{ IsPrime returns True if N is a prime number. }
{ }
{ IsPrimeFactor returns True if F is a prime factor of N. }
{ }
{ PrimeFactors returns an array of the prime factors of N. }
{ }
{ GCD returns the Greatest Common Divisor of N1 and N2. }
{ }
{ IsRelativePrime returns True if the GCD(X, Y) = 1. Relative prime is }
{ also called co-prime. }
{ }
{ InvMod calculates the modular inverse of A (mod N). }
{ If X is the modular inverse of A modulo N then: (X*A) mod N = 1 }
{ [ in mathematical notation it is: X*A = 1 (mod N) ] }
{ }
{ ExpMod calculates A^Z mod N. }
{ }
{ Jacobi calculates the 'Jacobi Symbol (a/n)'. }
{ }
function IsPrime(const N: Int64): Boolean;
function IsPrimeFactor(const N, F: Int64): Boolean;
function PrimeFactors(const N: Int64): Int64Array;
function GCD(const N1, N2: Integer): Integer; overload;
function GCD(const N1, N2: Int64): Int64; overload;
function IsRelativePrime(const X, Y: Int64): Boolean;
function InvMod(const A, N: Integer): Integer; overload;
function InvMod(const A, N: Int64): Int64; overload;
function ExpMod(A, Z: Integer; const N: Integer): Integer; overload;
function ExpMod(A, Z: Int64; const N: Int64): Int64; overload;
function Jacobi(const A, N: Integer): Integer;
{ }
{ Combinatoric functions }
{ }
function Factorial(const N: Integer): Extended;
function Combinations(const N, C: Integer): Extended;
function Permutations(const N, P: Integer): Extended;
function Fibonacci(const N: Integer): Int64;
{ }
{ Statistical functions }
{ }
{ GammaLn returns the natural logarithm of the gamma function. }
{ }
function GammaLn(X: Extended): Extended;
{ }
{ Fourier Transformations }
{ }
procedure FourierTransform(const AngleNumerator: Extended;
const RealIn, ImagIn: Array of Extended;
var RealOut, ImagOut: ExtendedArray);
procedure FFT(const RealIn, ImagIn: Array of Extended;
var RealOut, ImagOut: ExtendedArray);
procedure InverseFFT(const RealIn, ImagIn: Array of Extended;
var RealOut, ImagOut: ExtendedArray);
procedure CalcFrequency(const FrequencyIndex: Integer;
const RealIn, ImagIn: Array of Extended;
var RealOut, ImagOut: Extended);
{ }
{ Numerical routines }
{ }
type
fx = function (const x: Extended): Extended;
function SecantSolver(const f: fx; const y, Guess1, Guess2: Extended): Extended;
{ Uses Secant method to solve for x in f(x) = y }
function NewtonSolver(const f, df: fx; const y, Guess: Extended): Extended;
{ Uses Newton's method to solve for x in f(x) = y. }
{ df = f'(x). }
{ Note: This implementation does not check if the solver goes on a tangent }
{ (which can happen with certain Guess values) }
function FirstDerivative(const f: fx; const x: Extended): Extended;
{ Returns the value of f'(x) }
{ Uses (-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)) / 12h }
function SecondDerivative(const f: fx; const x: Extended): Extended;
{ Returns the value of f''(x) }
{ Uses (-f(x+2h) + 16f(x+h) - 30f(x) + 16f(x-h) - f(x-2h)) / 12h^2 }
function ThirdDerivative(const f: fx; const x: Extended): Extended;
{ Returns the value of f'''(x) }
{ Uses (f(x+2h) - 2f(x+h) + 2f(x-h) - f(x-2*h)) / 2h^3 }
function FourthDerivative(const f: fx; const x: Extended): Extended;
{ Returns the value of f''''(x) }
{ Uses (f(x+2h) - 4f(x+h) + 6f(x) - 4f(x-h) + f(x-2h)) / h^4 }
function SimpsonIntegration(const f: fx; const a, b: Extended; N: Integer): Extended;
{ Returns the area under f from a to b, by applying Simpson's 1/3 Rule with }
{ N subdivisions. }
{ }
{ Self-testing code }
{ }
procedure SelfTest;
implementation
uses
{ Delphi }
Math;
{ }
{ FPU (Floating Point Unit) functions }
{ }
procedure SetFPUPrecisionSingle;
begin
{$IFDEF DELPHI6_UP}
SetPrecisionMode(pmSingle);
{$ENDIF}
end;
procedure SetFPUPrecisionDouble;
begin
{$IFDEF DELPHI6_UP}
SetPrecisionMode(pmDouble);
{$ENDIF}
end;
procedure SetFPUPrecisionExtended;
begin
{$IFDEF DELPHI6_UP}
SetPrecisionMode(pmExtended);
{$ENDIF}
end;
procedure SetFPURoundingNearest;
begin
{$IFDEF DELPHI6_UP}
SetRoundMode(rmNearest);
{$ENDIF}
end;
procedure SetFPURoundingDown;
begin
{$IFDEF DELPHI6_UP}
SetRoundMode(rmDown);
{$ENDIF}
end;
procedure SetFPURoundingUp;
begin
{$IFDEF DELPHI6_UP}
SetRoundMode(rmUp);
{$ENDIF}
end;
procedure SetFPURoundingTruncate;
begin
{$IFDEF DELPHI6_UP}
SetRoundMode(rmTruncate);
{$ENDIF}
end;
{ }
{ Polynomial }
{ }
function PolyEval(const X: Extended; const Coef: Array of Extended; const N: Integer): Extended;
var P : Extended;
L : Integer;
begin
if Length(Coef) <> N + 1 then
raise EInvalidArgument.Create('PolyEval: Invalid number of coefficients');
P := 1.0;
L := N;
Result := 0.0;
While L >= 0 do
begin
Result := Result + Coef[L] * P;
P := P * X;
Dec(L);
end;
end;
{ }
{ Miscellaneous functions }
{ }
{$IFDEF CPU_INTEL386}
function Sgn(const R: Integer): Integer;
asm
TEST EAX, EAX
JL @Negative
JG @Positive
RET
@Negative:
MOV EAX, -1
RET
@Positive:
MOV EAX, 1
end;
{$ELSE}
function Sgn(const R: Integer): Integer;
begin
if R > 0 then
Result := 1 else
if R < 0 then
Result := -1 else
Result := 0;
end;
{$ENDIF}
function Sgn(const R: Int64): Integer;
begin
if R > 0 then
Result := 1 else
if R < 0 then
Result := -1 else
Result := 0;
end;
function Sgn(const R: Single): Integer;
begin
if R > 0.0 then
Result := 1 else
if R < 0.0 then
Result := -1 else
Result := 0;
end;
function Sgn(const R: Double): Integer;
begin
if R > 0.0 then
Result := 1 else
if R < 0.0 then
Result := -1 else
Result := 0;
end;
function Sgn(const R: Extended): Integer;
begin
if R > 0.0 then
Result := 1 else
if R < 0.0 then
Result := -1 else
Result := 0;
end;
function FloatMod(const A, B: Extended): Extended;
begin
Result := A - Floor(A / B) * B;
end;
{ }
{ Trigonometric functions }
{ }
function ATan360(const X, Y: Extended): Extended;
var Angle: Extended;
begin
if FloatZero(X) then
Angle := PiOn2 else
Angle := ArcTan(Y / X);
Angle := Angle * DegPerRad;
if (X <= 0.0) and (Y < 0.0) then
Angle := Angle - 180.0;
if (X < 0.0) and (Y > 0.0) then
Angle := Angle + 180.0;
if Angle < 0.0 then
Angle := Angle + 360.0;
ATan360 := Angle;
end;
function InverseTangentDeg(const X, Y: Extended): Extended;
{ 0 <= Result <= 360 }
var Angle : Extended;
begin
if FloatZero(X) then
Angle := PiOn2 else
Angle := ArcTan (Y / X);
Angle := Angle * 180.0 / Pi;
if (X <= 0.0) and (Y < 0.0) then
Angle := Angle - 180.0 else
if (X < 0.0) and (Y > 0.0) then
Angle := Angle + 180.0;
if Angle < 0.0 then
Angle := Angle + 360.0;
InverseTangentDeg := Angle;
end;
function InverseTangentRad(const X, Y: Extended): Extended;
{ 0 <= result <= 2pi }
var Angle : Extended;
begin
if FloatZero(X) then
Angle := PiOn2 else
Angle := ArcTan(Y / X);
if (X <= 0.0) and (Y < 0) then
Angle := Angle - Pi;
if (X < 0.0) and (Y > 0) then
Angle := Angle + Pi;
If Angle < 0 then
Angle := Angle + Pi2;
InverseTangentRad := Angle;
end;
function InverseSinDeg(const Y, R: Extended): Extended;
{ -90 <= result <= 90 }
var X : Extended;
begin
X := Sqrt(Sqr(R) - Sqr(Y));
Result := InverseTangentDeg(X, Y);
If Result > 90.0 then
Result := Result - 360.0;
end;
function InverseSinRad(const Y, R: Extended): Extended;
{ -90 <= result <= 90 }
var X : Extended;
begin
X := Sqrt(Sqr(R) - Sqr(Y));
Result := InverseTangentRad(X, Y);
if Result > 90.0 then
Result := Result - 360.0;
end;
function InverseCosDeg(const X, R: Extended): Extended;
{ -90 <= result <= 90 }
var Y : Extended;
begin
Y := Sqrt(Sqr(R) - Sqr(X));
Result := InverseTangentDeg(X, Y);
if Result > 90.0 then
Result := Result - 360.0;
end;
function InverseCosRad(const X, R: Extended): Extended;
{ -90 <= result <= 90 }
var Y : Extended;
begin
Y := Sqrt(Sqr(R) - Sqr(X));
Result := InverseTangentRad(X, Y);
if Result > 90.0 then
Result := Result - 360.0;
end;
procedure RealToDMS(const X: Extended; var Degs, Mins, Secs: Extended);
var Y : Extended;
begin
Degs := Int(X);
Y := Frac(X) * 60.0;
Mins := Int(Y);
Secs := Frac(Y) * 60.0;
end;
function DMSToReal(const Degs, Mins, Secs: Extended): Extended;
begin
Result := Degs + Mins / 60.0 + Secs / 3600.0;
end;
function CanonicalForm(const Theta: Extended): Extended; {-PI < theta <= PI}
begin
if Abs(Theta) > Pi then
Result := Round(Theta / Pi2) * Pi2 else
Result := Theta;
end;
procedure PolarToRectangular(const R, Theta: Extended; var X, Y: Extended);
var S, C : Extended;
begin
SinCos(Theta, S, C);
X := R * C;
Y := R * S;
end;
procedure RectangularToPolar(const X, Y: Extended; var R, Theta: Extended);
begin
if FloatZero(X) then
if FloatZero(Y) then
R := 0.0 else
if Y > 0.0 then
R := Y else
R := -Y else
R := Sqrt(Sqr(X) + Sqr(Y));
Theta := ArcTan2(Y, X);
end;
function Distance(const X1, Y1, X2, Y2: Extended): Extended;
begin
Result := Sqrt(Sqr(X1 - X2) + Sqr(Y1 - Y2));
end;
{ }
{ Prime numbers }
{ }
{ The prime functions use the fact that it's only necessary to check up to }
{ Sqrt(N) for factors of N when checking if N is prime. }
{ }
{ Lookups for prime numbers in the Byte range. }
const
BytePrimesCount = 54;
BytePrimesTable: Array[0..BytePrimesCount - 1] of Byte = (
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,
223, 227, 229, 233, 239, 241, 251);
BytePrimesSet: Set of Byte = [
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,
223, 227, 229, 233, 239, 241, 251];
{ Lookup for prime numbers in the Word range. }
const
WordPrimesLookupEntries = $10000; // 65536
WordPrimesLookupLength = WordPrimesLookupEntries div (Sizeof(LongWord) * 8); // 2048
WordPrimesLookupSize = WordPrimesLookupLength * Sizeof(LongWord); // 8192
var
WordPrimesInit : Boolean = False;
WordPrimesSet : Array of LongWord; // 8192 bytes when initialized, 1 bit for
// each number between 0-65535.
procedure ClearWordPrimesBit(const N: Word);
var I : LongWord;
P : PLongWord;
begin
I := N shr 5;
P := @WordPrimesSet[I];
P^ := ClearBit(P^, N and 31);
end;
procedure InitWordPrimesSet;
var I, J, P : Integer;
begin
// Initialize WordPrimesSet bits to True; clear bits 0 and 1
SetLength(WordPrimesSet, WordPrimesLookupLength);
FillChar(Pointer(WordPrimesSet)^, WordPrimesLookupSize, #$FF);
ClearWordPrimesBit(0);
ClearWordPrimesBit(1);
// Clear bits for multiples of Byte primes
For I := 0 to BytePrimesCount - 1 do
begin
P := BytePrimesTable[I];
For J := 2 to $FFFF div P do
ClearWordPrimesBit(P * J);
end;
// Set initialized
WordPrimesInit := True;
end;
{ IsPrime }
function IsPrime(const N: Int64): Boolean;
var V : Int64;
E : Extended;
M : LongWord;
B : LongWord;
I : Longword;
L : Longword;
begin
// Initialize value
V := N;
if V < 0 then
V := -V;
// Check if V is prime
if V <= $FF then
// V is in Byte range, lookup from BytePrimesSet
Result := Byte(V) in BytePrimesSet else
if V <= $FFFF then
begin
// V is in Word range
if not WordPrimesInit then
InitWordPrimesSet;
Result := IsBitSet(WordPrimesSet[Word(V) shr 5], Word(V) and 31);
end else
begin
// V is in LongWord range
// Look for factor from primes in Byte range
For I := 0 to BytePrimesCount - 1 do
begin
B := BytePrimesTable[I];
if V mod B = 0 then
begin
// N is divisable by B, N is not prime
Result := False;
exit;
end;
end;
// Calculate maximum factor
E := V;
M := Ceil(Sqrt(E));
// Look for factor from primes in Word range
if not WordPrimesInit then
InitWordPrimesSet;
if M > $FFFF then
L := $FFFF else
L := M;
For I := $100 to L do
if IsBitSet(WordPrimesSet[I shr 5], I and 31) then
if V mod I = 0 then
begin
// N is divisable by I, N is not prime
Result := False;
exit;
end;
// Look for factor from values in LongWord range
For I := $10000 to M do
if V mod I = 0 then
begin
// N is divisable by I, N is not prime
Result := False;
exit;
end;
// Sqrt(N) reached, N is prime
Result := True;
end;
end;
{ IsPrimeFactor }
function IsPrimeFactor(const N, F: Int64): Boolean;
begin
Result := (N mod F = 0) and IsPrime(F);
end;
{ PrimeFactors }
function PrimeFactors(const N: Int64): Int64Array;
var V : Int64;
E : Extended;
M : LongWord;
I : LongWord;
begin
// Initialize
Result := nil;
V := N;
if V < 0 then
V := -V;
// 0 and 1 has no prime factors
if V <= 1 then
exit;
// Calculate maximum factor
E := V;
M := Ceil(Sqrt(E));
// Find prime factors
For I := 2 to M do
if IsPrime(I) and (V mod I = 0) then
begin
// I is a prime factor
Append(Result, I);
Repeat
V := V div I;
if V = 1 then
// No more factors
exit;
Until V mod I <> 0;
end;
// Check for remaining prime factor
if IsPrime(V) then
Append(Result, V);
end;
{ Find the GCD using Euclid's algorithm }
function GCD(const N1, N2: Integer): Integer;
var X, Y, J : Integer;
begin
X := N1;
Y := N2;
if X < Y then
Swap(X, Y);
While (X <> 1) and (X <> 0) and (Y <> 1) and (Y <> 0) do
begin
J := (X - Y) mod Y;
if J = 0 then
begin
Result := Abs(Y);
exit;
end;
X := Y;
Y := J;
end;
Result := 1;
end;
function GCD(const N1, N2: Int64): Int64;
var X, Y, J : Int64;
begin
X := N1;
Y := N2;
if X < Y then
Swap(X, Y);
While (X <> 1) and (X <> 0) and (Y <> 1) and (Y <> 0) do
begin
J := (X - Y) mod Y;
if J = 0 then
begin
Result := Abs(Y);
exit;
end;
X := Y;
Y := J;
end;
Result := 1;
end;
function IsRelativePrime(const X, Y: Int64): Boolean;
begin
Result := GCD(X, Y) = 1;
end;
{ Find the modular inverse of A modulo N using Euclid's algorithm. }
function InvMod(const A, N: Integer): Integer;
var g0, g1, v0, v1, y, z : Integer;
begin
if N < 2 then
raise EInvalidArgument.CreateFmt('InvMod: n=%d < 2', [N]);
if GCD (A, N) <> 1 then
raise EInvalidArgument.CreateFmt('InvMod: GCD (a=%d, n=%d) <> 1', [A, N]);
g0 := N; g1 := A;
v0 := 0; v1 := 1;
while g1 <> 0 do
begin
y := g0 div g1;
z := g1;
g1 := g0 - y * g1;
g0 := z;
z := v1;
v1 := v0 - y * v1;
v0 := z;
end;
if v0 > 0 then
Result := v0 else
Result := v0 + N;
end;
function InvMod(const A, N: Int64): Int64;
var g0, g1, v0, v1, y, z : Int64;
begin
if N < 2 then
raise EInvalidArgument.CreateFmt ('InvMod: n=%d < 2', [N]);
if GCD(A, N) <> 1 then
raise EInvalidArgument.CreateFmt ('InvMod: GCD(a=%d, n=%d) <> 1', [A, N]);
g0 := N; g1 := A;
v0 := 0; v1 := 1;
while g1 <> 0 do
begin
y := g0 div g1;
z := g1;
g1 := g0 - y * g1;
g0 := z;
z := v1;
v1 := v0 - y * v1;
v0 := z;
end;
if v0 > 0 then
Result := v0 else
Result := v0 + N;
end;
{ Calculates x = a^z mod n }
function ExpMod(A, Z: Integer; const N: Integer): Integer;
var Signed : Boolean;
begin
Signed := Z < 0;
if Signed then
Z := -Z;
Result := 1;
while Z <> 0 do
begin
while not Odd(Z) do
begin
Z := Z shr 1;
A := (A * Int64(A)) mod N;
end;
Dec (Z);
Result := (Result * Int64(A)) mod N;
end;
if Signed then
Result := InvMod(Result, N);
end;
function ExpMod(A, Z: Int64; const N: Int64): Int64;
var Signed : Boolean;
begin
Signed := Z < 0;
if Signed then
Z := -Z;
Result := 1;
while Z <> 0 do
begin
while not Odd(Z) do
begin
Z := Z shr 1;
A := (A * Int64(A)) mod N;
end;
Dec (Z);
Result := (Result * Int64(A)) mod N;
end;
if Signed then
Result := InvMod(Result, N);
end;
{ From: Handbook of Applied Cryptography, Algorithm 2.149 (page 73) }
function Jacobi(const A, N: Integer): Integer;
// Calculates a=2^s+t
procedure GetRepresentationBase2(const a: Integer; var s,t: Integer);
begin
s := 0;
t := a;
while not Odd(T) do
begin
t := t shr 1;
Inc(s);
end;
end;
var e, a1, s : Integer;
begin
// Check for valid input
if a < 0 then
raise EInvalidArgument.CreateFmt('Jacobi: a=%d < 0', [a]);
if a >= n then
raise EInvalidArgument.CreateFmt('Jacobi: a=%d >= n=%d', [a, n]);
if not Odd (n) then
raise EInvalidArgument.CreateFmt('Jacobi: Odd(n=%d) = False', [n]);
if n < 3 then
raise EInvalidArgument.CreateFmt('Jacobi: n=%d < 3', [n]);
if a > 1 then
begin
GetRepresentationBase2(a, e, a1);
s := 1;
if Odd (e) and (((n mod 8) = 3) or ((n mod 8) = 5)) then
s := -s;
if ((n mod 4) = 3) and ((a1 mod 4) = 3) then
s := -s;
if a1=1 then
Result := s else
Result := s * Jacobi(n mod a1, a1);
end else
Result := a;
end;
{ }
{ Numerical solvers }
{ }
function SecantSolver(const f: fx; const y, Guess1, Guess2: Extended): Extended;
var xn, xnm1, xnp1, fxn, fxnm1 : Extended;
begin
xnm1 := Guess1;
xn := Guess2;
fxnm1 := f(xnm1) - y;
Repeat
fxn := f(xn) - y;
xnp1 := xn - fxn * (xn - xnm1) / (fxn - fxnm1);
fxnm1 := fxn;
xnm1 := xn;
xn := xnp1;
Until (f(xn - 0.00000001) - y) * (f(xn + 0.00000001) - y) <= 0.0;
Result := xn;
end;
function NewtonSolver(const f, df: fx; const y, Guess: Extended): Extended;
var xn, xnp1 : Extended;
begin
xnp1 := Guess;
Repeat
xn := xnp1;
xnp1 := xn - f(xn) / df(xn);
Until Abs(xnp1 - xn) < 0.000000000000001;
Result := xn;
end;
const h = 1e-15;
function FirstDerivative(const f: fx; const x: Extended): Extended;
begin
Result := (-f(x + 2 * h) + 8 * f(x + h) - 8 * f(x - h) + f(x - 2 * h)) / (12 * h);
end;
function SecondDerivative(const f: fx; const x: Extended): Extended;
begin
Result := (-f(x + 2 * h) + 16 * f(x + h) - 30 * f(x) +
16 * f(x - h) - f(x - 2 * h)) / (12 * h * h);
end;
function ThirdDerivative(const f: fx; const x: Extended): Extended;
begin
Result := (f(x + 2 * h) - 2 * f(x + h) + 2 * f(x - h) - f(x - 2 * h)) / (2 * h * h * h);
end;
function FourthDerivative(const f: fx; const x: Extended): Extended;
begin
Result := (f(x + 2 * h) - 4 * f(x + h) + 6 * f(x) - 4 * f(x - h) + f(x - 2 * h)) / (h * h * h * h);
end;
function SimpsonIntegration(const f: fx; const a, b: Extended; N: Integer): Extended;
var h : Extended;
I : Integer;
begin
if N mod 2 = 1 then
Inc(N); // N must be multiple of 2
h := (b - a) / N;
Result := 0.0;
For I := 1 to N - 1 do
Result := Result + ((I mod 2) * 2 + 2) * f(a + (I - 0.5) * h);
Result := (Result + f(a) + f(b)) * h / 3.0;
end;
{ }
{ Fast factorial }
{ }
{ For small values of N, calculate using 2*3*..*N, and cache result. }
{ For larger values of N, calculate using Gamma(N+1) }
{ }
const
FactorialSmallLimit = 34;
FactorialCacheLimit = 409;
FactorialMaxLimit = 1754;
var
FactorialCache : ExtendedArray = nil;
function Factorial(const N: Integer): Extended;
var L : Extended;
I : Integer;
begin
// Check that arguments are valid
if N > FactorialMaxLimit then
raise EOverflow.Create('Factorial overflow');
if N < 0 then
raise EInvalidArgument.Create('Factorial: Invalid argument');
// Get result from factorial cache
if (N <= FactorialCacheLimit) and Assigned(FactorialCache) and
(FactorialCache[N] >= 1.0) then
begin
Result := FactorialCache[N];
exit;
end;
// Calculate
if N < FactorialSmallLimit then
begin
L := 1.0;
I := 2;
While I <= N do
begin
L := L * I;
Inc(I);
end;
Result := L;
end
else
Result := Exp(GammaLn(N + 1));
// Save result in cache
if N <= FactorialCacheLimit then
begin
if not Assigned(FactorialCache) then
SetLength(FactorialCache, FactorialCacheLimit + 1);
FactorialCache[N] := Result;
end;
end;
{ }
{ Combinatorics }
{ }
function Combinations(const N, C: Integer): Extended;
begin
Result := Factorial(N) / (Factorial(C) * Factorial(N - C));
end;
function Permutations(const N, P: Integer): Extended;
begin
Result := Factorial(N) / Factorial(N - P);
end;
{ }
{ Fibonacci (N) = Fibonacci (N - 1) + Fibonacci (N - 2) }
{ }
function Fibonacci(const N: Integer): Int64;
var I : Integer;
f1, f2 : Int64;
begin
if (N < 0) or (N > 92) then
raise EInvalidArgument.Create('Fibonacci: Invalid argument');
Result := 1;
if N = 1 then
exit;
f1 := 0; // fib(0) = 0
f2 := 1; // fib(1) = 1
for I := 1 to N do
begin
Result := f1 + f2;
f2 := f1;
f1 := Result;
end;
end;
{ }
{ Statistical functions }
{ }
{ "For arguments greater than 13, the logarithm of the gamma }
{ function is approximated by the logarithmic version of }
{ Stirling's formula using a polynomial approximation of }
{ degree 4. Arguments between -33 and +33 are reduced by }
{ recurrence to the interval [2,3] of a rational approximation. }
{ The cosecant reflection formula is employed for arguments }
{ less than -33. }
{ }
{ Arguments greater than MAXLGM return MAXNUM and an error }
{ message." }
{ }
{ Algorithm translated into Delphi by David Butler from Cephes C }
{ library with permission from Stephen L. Moshier }
{ <moshier@na-net.ornl.gov>. }
function GammaLn (X: Extended): Extended;
const MaxLGM = 2.556348e305;
var P, Q, W, Z : Extended;
{ Stirling's formula expansion of log gamma }
const Stir : Array[0..4] of Extended = (
8.11614167470508450300E-4,
-5.95061904284301438324E-4,
7.93650340457716943945E-4,
-2.77777777730099687205E-3,
8.33333333333331927722E-2);
{ B[], C[]: log gamma function between 2 and 3 }
B : Array[0..5] of Extended = (
-1.37825152569120859100E3,
-3.88016315134637840924E4,
-3.31612992738871184744E5,
-1.16237097492762307383E6,
-1.72173700820839662146E6,
-8.53555664245765465627E5);
C : Array[0..7] of Extended = (
1.00000000000000000000E0,
-3.51815701436523470549E2,
-1.70642106651881159223E4,
-2.20528590553854454839E5,
-1.13933444367982507207E6,
-2.53252307177582951285E6,
-2.01889141433532773231E6,
1.00000000000000000000E0);
begin
if X < -34.0 then
begin
Q := -X;
W := GammaLn(Q);
P := Trunc(Q);
if P = Q then
raise EOverflow.Create('GammaLn') else
begin
Z := Q - P;
if Z > 0.5 then
begin
P := P + 1.0;
Z := P - Q;
end;
Z := Q * Sin(Pi * Z);
if Z = 0.0 then
raise EOverflow.Create('GammaLn') else
GammaLn := LnPi - Ln(Z) - W;
end;
end else
if X <= 13 then
begin
Z := 1.0;
While X >= 3.0 do
begin
X := X - 1.0;
Z := Z * X;
end;
While (X < 2.0) and (X <> 0.0) do
begin
Z := Z / X;
X := X + 1.0;
end;
if X = 0.0 then
raise EOverflow.Create('GammaLn') else
if Z < 0.0 then
Z := -Z;
if X = 2.0 then
GammaLn := Ln(Z)
else
begin
X := X - 2.0;
P := X * PolyEval(X, B, 5) / PolyEval(X, C, 7);
GammaLn := Ln(Z) + P;
end;
end else
if X > MAXLGM then
raise EOverflow.Create('GammaLn')
else
begin
Q := (X - 0.5) * Ln(X) - X + LnSqrt2Pi;
if X > 1.0e8 then
GammaLn := Q else
begin
P := 1.0 / (X * X);
if X >= 1000.0 then
GammaLn := Q + ((7.9365079365079365079365e-4 * P
- 2.7777777777777777777778e-3)
* P + 0.0833333333333333333333) / X else
GammaLn := Q + PolyEval(P, Stir, 4) / X;
end;
end;
end;
{ }
{ Fourier Transformations }
{ FourierTransform is a FFT routine adapted for Delphi from a Pascal }
{ version by Don Cross (see http://www.intersrv.com/~dcross/fft.html) }
{ }
procedure FourierTransform(const AngleNumerator: Extended;
const RealIn, ImagIn: Array of Extended;
var RealOut, ImagOut: ExtendedArray);
var I, J, N, L : LongWord;
NumSamples, NumBits : LongWord;
BlockSize, Blockend : LongWord;
delta_angle, delta_ar : Extended;
alpha, beta : Extended;
tr, ti, ar, ai : Extended;
PRe, PIm, QRe, QIm : ^Extended;
begin
NumSamples := Length(RealIn);
L := Length(ImagIn);
if (L > 0) and (NumSamples <> L) then
raise EInvalidArgument.Create('FourierTransform: RealIn and ImagIn must be of equal length');
if NumSamples < 2 then
raise EInvalidArgument.Create('FourierTransform: At least two samples required');
if not IsPowerOfTwo(NumSamples) then
raise EInvalidArgument.Create('FourierTransform: NumSamples not a power of two');
SetLength(RealOut, NumSamples);
SetLength(ImagOut, NumSamples);
NumBits := SetBitScanForward(NumSamples);
For I := 0 to NumSamples - 1 do
begin
J := ReverseBits (I, NumBits);
RealOut[J] := RealIn[I];
if L > 0 then
ImagOut[J] := ImagIn[I];
end;
Blockend := 1;
BlockSize := 2;
While BlockSize <= NumSamples do
begin
delta_angle := AngleNumerator / BlockSize;
alpha := Sin (0.5 * delta_angle);
alpha := 2.0 * alpha * alpha;
beta := Sin (delta_angle);
I := 0;
While I < NumSamples do
begin
ar := 1.0; { cos(0) }
ai := 0.0; { sin(0) }
PRe := @RealOut[I];
PIm := @ImagOut[I];
J := I + Blockend;
QRe := @RealOut[J];
QIm := @ImagOut[J];
For N := 0 to Blockend - 1 do
begin
tr := ar * QRe^ - ai * QIm^;
ti := ar * QIm^ + ai * QRe^;
QRe^ := PRe^ - tr;
QIm^ := PIm^ - ti;
PRe^ := PRe^ + tr;
PIm^ := PIm^ + ti;
delta_ar := alpha * ar + beta * ai;
ai := ai - (alpha * ai - beta * ar);
ar := ar - delta_ar;
Inc(PRe);
Inc(PIm);
Inc(QRe);
Inc(QIm);
end;
Inc(I, BlockSize);
end;
Blockend := BlockSize;
BlockSize := BlockSize shl 1;
end;
end;
procedure FFT(const RealIn, ImagIn: Array of Extended; var RealOut, ImagOut: ExtendedArray);
begin
FourierTransform(Pi2, RealIn, ImagIn, RealOut, ImagOut);
end;
procedure InverseFFT(const RealIn, ImagIn: Array of Extended; var RealOut, ImagOut: ExtendedArray);
var I, N : Integer;
begin
FourierTransform (-Pi2, RealIn, ImagIn, RealOut, ImagOut);
{ Normalize the resulting time samples }
N := Length(RealOut);
For I := 0 to N - 1 do
begin
RealOut[I] := RealOut[I] / N;
ImagOut[I] := ImagOut[I] / N;
end;
end;
procedure CalcFrequency(const FrequencyIndex: Integer;
const RealIn, ImagIn: Array of Extended;
var RealOut, ImagOut: Extended);
var K, NumSamples : Integer;
cos1, cos2, cos3, theta, beta : Extended;
sin1, sin2, sin3 : Extended;
begin
NumSamples := Length(RealIn);
if NumSamples <> Length(ImagIn) then
raise EInvalidArgument.Create('CalcFrequency: RealIn and ImagIn must be of equal length');
RealOut := 0.0;
ImagOut := 0.0;
theta := Pi2 * FrequencyIndex / NumSamples;
sin1 := sin (-2 * theta);
sin2 := sin (-theta);
cos1 := cos (-2 * theta);
cos2 := cos (-theta);
beta := 2 * cos2;
For K := 0 to NumSamples - 1 do
begin
sin3 := beta * sin2 - sin1;
sin1 := sin2;
sin2 := sin3;
cos3 := beta * cos2 - cos1;
cos1 := cos2;
cos2 := cos3;
RealOut := RealOut + RealIn[K] * cos3 - ImagIn[K] * sin3;
ImagOut := ImagOut + ImagIn[K] * cos3 + RealIn[K] * sin3;
end;
end;
{ }
{ Self-testing code }
{ }
{$ASSERTIONS ON}
procedure SelfTest;
var A : Int64Array;
begin
Assert(Sgn(0) = 0, 'Sgn');
Assert(Sgn(1) = 1, 'Sgn');
Assert(Sgn(-1) = -1, 'Sgn');
Assert(Sgn(2) = 1, 'Sgn');
Assert(Sgn(1.1) = 1, 'Sgn');
Assert(not IsPrime(0), 'IsPrime(0)');
Assert(not IsPrime(1), 'IsPrime(1)');
Assert(IsPrime(2), 'IsPrime(2)');
Assert(IsPrime(3), 'IsPrime(3)');
Assert(IsPrime(-3), 'IsPrime(-3)');
Assert(not IsPrime(4), 'IsPrime(4)');
Assert(not IsPrime(-4), 'IsPrime(-4)');
Assert(IsPrime(257), 'IsPrime(257)');
Assert(not IsPrime($10002), 'IsPrime($10002)');
Assert(IsPrime($10003), 'IsPrime($10003)');
Assert(IsPrimeFactor(17 * 3, 17), 'IsPrimeFactor');
Assert(IsPrimeFactor(17 * 3, 3), 'IsPrimeFactor');
Assert(not IsPrimeFactor(17 * 3, 2), 'IsPrimeFactor');
Assert(not IsPrimeFactor(17, 1), 'IsPrimeFactor');
Assert(IsPrimeFactor(17, 17), 'IsPrimeFactor');
A := PrimeFactors(17 * 3);
Assert(Length(A) = 2, 'PrimeFactors');
Assert(A[0] = 3, 'PrimeFactors');
Assert(A[1] = 17, 'PrimeFactors');
A := PrimeFactors(2 * 2 * 5 * 11 * 19);
Assert(Length(A) = 4, 'PrimeFactors');
Assert(A[0] = 2, 'PrimeFactors');
Assert(A[1] = 5, 'PrimeFactors');
Assert(A[2] = 11, 'PrimeFactors');
Assert(A[3] = 19, 'PrimeFactors');
Assert(GammaLn(2.0) = 0.0, 'GammaLn');
Assert(ApproxEqual(Exp(GammaLn(3.0)), 2.0), 'GammaLn');
Assert(ApproxEqual(Exp(GammaLn(5.0)), 24.0), 'GammaLn');
Assert(ApproxEqual(Exp(GammaLn(6.0)), 120.0), 'GammaLn');
end;
initialization
SetFPUPrecisionExtended;
SetFPURoundingNearest;
finalization
end.