📄 jclmath.pas
字号:
{**************************************************************************************************}
{ }
{ Project JEDI Code Library (JCL) }
{ }
{ The contents of this file are subject to the Mozilla Public License Version 1.1 (the "License"); }
{ you may not use this file except in compliance with the License. You may obtain a copy of the }
{ License at http://www.mozilla.org/MPL/ }
{ }
{ Software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTY OF }
{ ANY KIND, either express or implied. See the License for the specific language governing rights }
{ and limitations under the License. }
{ }
{ The Original Code is JclMath.pas. }
{ }
{ The Initial Developers of the Original Code are Clayton Collie, David Butler, ESB Consultancy, }
{ Jean Debord, Marcel van Brakel and Michael Schnell. }
{ Portions created by these individuals are Copyright (C) of these individuals. }
{ All Rights Reserved. }
{ }
{ Contributors: }
{ Ernesto Benestante }
{ Marcel van Brakel }
{ Aleksei Koudinov }
{ Robert Marquardt (marquardt) }
{ Robert Rossmair (rrossmair) }
{ Matthias Thoma (mthoma) }
{ Mark Vaughan }
{ unknown }
{ }
{**************************************************************************************************}
{ }
{ Various mathematics classes and routines. Includes prime numbers, rational }
{ numbers, generic floating point routines, hyperbolic and transcendenatal }
{ routines, NAN and INF support and more. }
{ }
{ Unit owner: Matthias Thoma }
{ }
{**************************************************************************************************}
// Last modified: $Date: 2005/03/08 08:33:17 $
// For history see end of file
unit JclMath;
{$I jcl.inc}
interface
uses
Classes, SysUtils,
JclBase;
{ Mathematical constants }
const
Bernstein: Float = 0.2801694990238691330364364912307; // Bernstein constant
Cbrt2: Float = 1.2599210498948731647672106072782; // CubeRoot(2)
Cbrt3: Float = 1.4422495703074083823216383107801; // CubeRoot(3)
Cbrt10: Float = 2.1544346900318837217592935665194; // CubeRoot(10)
Cbrt100: Float = 4.6415888336127788924100763509194; // CubeRoot(100)
CbrtPi: Float = 1.4645918875615232630201425272638; // CubeRoot(PI)
Catalan: Float = 0.9159655941772190150546035149324; // Catalan constant
Pi: Float = 3.1415926535897932384626433832795; // PI
PiOn2: Float = 1.5707963267948966192313216916398; // PI / 2
PiOn3: Float = 1.0471975511965977461542144610932; // PI / 3
PiOn4: Float = 0.78539816339744830961566084581988; // PI / 4
Sqrt2: Float = 1.4142135623730950488016887242097; // Sqrt(2)
Sqrt3: Float = 1.7320508075688772935274463415059; // Sqrt(3)
Sqrt5: Float = 2.2360679774997896964091736687313; // Sqrt(5)
Sqrt10: Float = 3.1622776601683793319988935444327; // Sqrt(10)
SqrtPi: Float = 1.7724538509055160272981674833411; // Sqrt(PI)
Sqrt2Pi: Float = 2.506628274631000502415765284811; // Sqrt(2 * PI)
TwoPi: Float = 6.283185307179586476925286766559; // 2 * PI
ThreePi: Float = 9.4247779607693797153879301498385; // 3 * PI
Ln2: Float = 0.69314718055994530941723212145818; // Ln(2)
Ln10: Float = 2.3025850929940456840179914546844; // Ln(10)
LnPi: Float = 1.1447298858494001741434273513531; // Ln(PI)
Log2: Float = 0.30102999566398119521373889472449; // Log10(2)
Log3: Float = 0.47712125471966243729502790325512; // Log10(3)
LogPi: Float = 0.4971498726941338543512682882909; // Log10(PI)
LogE: Float = 0.43429448190325182765112891891661; // Log10(E)
E: Float = 2.7182818284590452353602874713527; // Natural constant
hLn2Pi: Float = 0.91893853320467274178032973640562; // Ln(2*PI)/2
inv2Pi: Float = 0.159154943091895; // 0.5 / Pi
TwoToPower63: Float = 9223372036854775808.0; // 2^63
GoldenMean: Float = 1.618033988749894848204586834365638; // GoldenMean
EulerMascheroni: Float = 0.5772156649015328606065120900824; // Euler GAMMA
const
MaxAngle: Float = 9223372036854775808.0; // 2^63 Rad
{$IFDEF MATH_EXTENDED_PRECISION}
MaxTanH: Float = 5678.2617031470719747459655389854; // Ln(2^16384)/2
MaxFactorial = 1754;
MaxFloatingPoint: Float = 1.189731495357231765085759326628E+4932; // 2^16384
MinFloatingPoint: Float = 3.3621031431120935062626778173218E-4932; // 2^(-16382)
{$ENDIF MATH_EXTENDED_PRECISION}
{$IFDEF MATH_DOUBLE_PRECISION}
MaxTanH: Float = 354.89135644669199842162284618659; // Ln(2^1024)/2
MaxFactorial = 170;
MaxFloatingPoint: Float = 1.797693134862315907729305190789E+308; // 2^1024
MinFloatingPoint: Float = 2.2250738585072013830902327173324E-308; // 2^(-1022)
{$ENDIF MATH_DOUBLE_PRECISION}
{$IFDEF MATH_SINGLE_PRECISION}
MaxTanH: Float = 44.361419555836499802702855773323; // Ln(2^128)/2
MaxFactorial = 33;
MaxFloatingPoint: Float = 3.4028236692093846346337460743177E+38; // 2^128
MinFloatingPoint: Float = 1.1754943508222875079687365372222E-38; // 2^(-126)
{$ENDIF MATH_SINGLE_PRECISION}
var
PrecisionTolerance: Float = 0.0000001;
EpsSingle: Single;
EpsDouble: Double;
EpsExtended: Extended;
Epsilon: Float;
ThreeEpsSingle: Single;
ThreeEpsDouble: Double;
ThreeEpsExtended: Extended;
ThreeEpsilon: Float;
type
TPrimalityTestMethod = (ptTrialDivision, ptRabinMiller);
{ Logarithmic }
function LogBase10(X: Float): Float;
function LogBase2(X: Float): Float;
function LogBaseN(Base, X: Float): Float;
{ Transcendental }
function ArcCos(X: Float): Float;
function ArcCot(X: Float): Float;
function ArcCsc(X: Float): Float;
function ArcSec(X: Float): Float;
function ArcSin(X: Float): Float;
function ArcTan(X: Float): Float;
function ArcTan2(Y, X: Float): Float;
function Cos(X: Float): Float; overload;
function Cot(X: Float): Float; overload;
function Coversine(X: Float): Float;
function Csc(X: Float): Float; overload;
function Exsecans(X: Float): Float;
function Haversine(X: Float): Float;
function Sec(X: Float): Float; overload;
function Sin(X: Float): Float; overload;
procedure SinCos(X: Float; var Sin, Cos: Float);
function Tan(X: Float): Float; overload;
function Versine(X: Float): Float;
{ Hyperbolic }
function ArcCosH(X: Float): Float;
function ArcCotH(X: Float): Float;
function ArcCscH(X: Float): Float;
function ArcSecH(X: Float): Float;
function ArcSinH(X: Float): Float;
function ArcTanH(X: Float): Float;
function CosH(X: Float): Float; overload;
function CotH(X: Float): Float; overload;
function CscH(X: Float): Float; overload;
function SecH(X: Float): Float; overload;
function SinH(X: Float): Float; overload;
function TanH(X: Float): Float; overload;
{ Coordinate conversion }
function DegMinSecToFloat(const Degs, Mins, Secs: Float): Float; // obsolete (see JclUnitConv)
procedure FloatToDegMinSec(const X: Float; var Degs, Mins, Secs: Float); // obsolete (see JclUnitConv)
{ Exponential }
function Exp(const X: Float): Float; overload;
function Power(const Base, Exponent: Float): Float; overload;
function PowerInt(const X: Float; N: Integer): Float; overload;
function TenToY(const Y: Float): Float;
function TruncPower(const Base, Exponent: Float): Float;
function TwoToY(const Y: Float): Float;
{ Floating point numbers support routines }
function IsFloatZero(const X: Float): Boolean;
function FloatsEqual(const X, Y: Float): Boolean;
function MaxFloat(const X, Y: Float): Float;
function MinFloat(const X, Y: Float): Float;
function ModFloat(const X, Y: Float): Float;
function RemainderFloat(const X, Y: Float): Float;
function SetPrecisionTolerance(NewTolerance: Float): Float;
procedure SwapFloats(var X, Y: Float);
procedure CalcMachineEpsSingle;
procedure CalcMachineEpsDouble;
procedure CalcMachineEpsExtended;
procedure CalcMachineEps;
procedure SetPrecisionToleranceToEpsilon;
{ Miscellaneous }
function Ackermann(const A, B: Integer): Integer;
function Ceiling(const X: Float): Integer;
function CommercialRound(const X: Float): Int64;
function Factorial(const N: Integer): Float;
function Fibonacci(const N: Integer): Integer;
function Floor(const X: Float): Integer;
function GCD(const X, Y: Cardinal): Cardinal;
function ISqrt(const I: Smallint): Smallint;
function LCM(const X, Y: Cardinal): Cardinal;
function NormalizeAngle(const Angle: Float): Float;
function Pythagoras(const X, Y: Float): Float;
function Sgn(const X: Float): Integer;
function Signe(const X, Y: Float): Float;
{ Ranges }
function EnsureRange(const AValue, AMin, AMax: Integer): Integer; overload;
function EnsureRange(const AValue, AMin, AMax: Int64): Int64; overload;
function EnsureRange(const AValue, AMin, AMax: Double): Double; overload;
{ Prime numbers }
function IsRelativePrime(const X, Y: Cardinal): Boolean;
function IsPrimeTD(N: Cardinal): Boolean;
function IsPrimeRM(N: Cardinal): Boolean;
function IsPrimeFactor(const F, N: Cardinal): Boolean;
function PrimeFactors(N: Cardinal): TDynCardinalArray;
var
IsPrime: function(N: Cardinal): Boolean = IsPrimeTD;
procedure SetPrimalityTest(const Method: TPrimalityTestMethod);
{ Floating point value classification }
type
TFloatingPointClass =
(
fpZero, // zero
fpNormal, // normal finite <> 0
fpDenormal, // denormalized finite
fpInfinite, // infinite
fpNaN, // not a number
fpInvalid // unsupported floating point format
);
function FloatingPointClass(const Value: Single): TFloatingPointClass; overload;
function FloatingPointClass(const Value: Double): TFloatingPointClass; overload;
function FloatingPointClass(const Value: Extended): TFloatingPointClass; overload;
{ NaN and INF support }
type
TNaNTag = -$3FFFFF..$3FFFFE;
const
Infinity = 1/0; // tricky
NaN = 0/0; // tricky
NegInfinity = -Infinity;
function IsInfinite(const Value: Single): Boolean; overload;
function IsInfinite(const Value: Double): Boolean; overload;
function IsInfinite(const Value: Extended): Boolean; overload;
function IsNaN(const Value: Single): Boolean; overload;
function IsNaN(const Value: Double): Boolean; overload;
function IsNaN(const Value: Extended): Boolean; overload;
function IsSpecialValue(const X: Float): Boolean;
procedure MakeQuietNaN(var X: Single; Tag: TNaNTag = 0); overload;
procedure MakeQuietNaN(var X: Double; Tag: TNaNTag = 0); overload;
procedure MakeQuietNaN(var X: Extended; Tag: TNaNTag = 0); overload;
procedure MakeSignalingNaN(var X: Single; Tag: TNaNTag = 0); overload;
procedure MakeSignalingNaN(var X: Double; Tag: TNaNTag = 0); overload;
procedure MakeSignalingNaN(var X: Extended; Tag: TNaNTag = 0); overload;
{ Mine*Buffer fills "Buffer" with consecutive tagged signaling NaNs.
This allows for real number arrays which enforce initialization: any attempt
to load an uninitialized array element into the FPU will raise an exception
either of class EInvalidOp (Windows 9x/ME) or EJclNaNSignal (Windows NT).
Under Windows NT it is thus possible to derive the violating array index from
the EJclNaNSignal object's Tag property. }
procedure MineSingleBuffer(var Buffer; Count: Integer; StartTag: TNaNTag = 0);
procedure MineDoubleBuffer(var Buffer; Count: Integer; StartTag: TNaNTag = 0);
function MinedSingleArray(Length: Integer): TDynSingleArray;
function MinedDoubleArray(Length: Integer): TDynDoubleArray;
function GetNaNTag(const NaN: Single): TNaNTag; overload;
function GetNaNTag(const NaN: Double): TNaNTag; overload;
function GetNaNTag(const NaN: Extended): TNaNTag; overload;
{ Set support }
type
TJclASet = class(TObject)
protected
function GetBit(const Idx: Integer): Boolean; virtual; abstract;
procedure SetBit(const Idx: Integer; const Value: Boolean); virtual; abstract;
procedure Clear; virtual; abstract;
procedure Invert; virtual; abstract;
function GetRange(const Low, High: Integer; const Value: Boolean): Boolean; virtual; abstract;
procedure SetRange(const Low, High: Integer; const Value: Boolean); virtual; abstract;
end;
type
TJclFlatSet = class(TJclASet)
private
FBits: TBits;
public
constructor Create;
destructor Destroy; override;
procedure Clear; override;
procedure Invert; override;
procedure SetRange(const Low, High: Integer; const Value: Boolean); override;
function GetBit(const Idx: Integer): Boolean; override;
function GetRange(const Low, High: Integer; const Value: Boolean): Boolean; override;
procedure SetBit(const Idx: Integer; const Value: Boolean); override;
end;
type
TPointerArray = array [0..MaxLongint div 256] of Pointer;
PPointerArray = ^TPointerArray;
TDelphiSet = set of Byte; // 256 elements
PDelphiSet = ^TDelphiSet;
const
EmptyDelphiSet: TDelphiSet = [];
CompleteDelphiSet: TDelphiSet = [0..255];
type
TJclSparseFlatSet = class(TJclASet)
private
FSetList: PPointerArray;
FSetListEntries: Integer;
public
destructor Destroy; override;
procedure Clear; override;
procedure Invert; override;
function GetBit(const Idx: Integer): Boolean; override;
procedure SetBit(const Idx: Integer; const Value: Boolean); override;
procedure SetRange(const Low, High: Integer; const Value: Boolean); override;
function GetRange(const Low, High: Integer; const Value: Boolean): Boolean; override;
end;
{ Rational numbers }
type
TJclRational = class(TObject)
private
FT: Integer;
FN: Integer;
function GetAsString: string;
procedure SetAsString(const S: string);
function GetAsFloat: Float;
procedure SetAsFloat(const R: Float);
protected
procedure Simplify;
public
constructor Create; overload;
constructor Create(const R: Float); overload;
constructor Create(const Numerator: Integer; const Denominator: Integer = 1); overload;
property Numerator: Integer read FT;
property Denominator: Integer read FN;
property AsString: string read GetAsString write SetAsString;
property AsFloat: Float read GetAsFloat write SetAsFloat;
procedure Assign(const R: TJclRational); overload;
procedure Assign(const R: Float); overload;
procedure Assign(const Numerator: Integer; const Denominator: Integer = 1); overload;
procedure AssignZero;
procedure AssignOne;
function Duplicate: TJclRational;
function IsEqual(const R: TJclRational): Boolean; reintroduce; overload;
function IsEqual(const Numerator: Integer; const Denominator: Integer = 1) : Boolean; reintroduce; overload;
function IsEqual(const R: Float): Boolean; reintroduce; overload;
function IsZero: Boolean;
function IsOne: Boolean;
procedure Add(const R: TJclRational); overload;
procedure Add(const V: Float); overload;
procedure Add(const V: Integer); overload;
procedure Subtract(const R: TJclRational); overload;
procedure Subtract(const V: Float); overload;
procedure Subtract(const V: Integer); overload;
procedure Negate;
procedure Abs;
function Sgn: Integer;
procedure Multiply(const R: TJclRational); overload;
procedure Multiply(const V: Float); overload;
procedure Multiply(const V: Integer); overload;
procedure Reciprocal;
procedure Divide(const R: TJclRational); overload;
procedure Divide(const V: Float); overload;
procedure Divide(const V: Integer); overload;
procedure Sqrt;
procedure Sqr;
procedure Power(const R: TJclRational); overload;
procedure Power(const V: Integer); overload;
procedure Power(const V: Float); overload;
end;
type
EJclMathError = class(EJclError);
EJclNaNSignal = class(EJclMathError)
private
FTag: TNaNTag;
public
constructor Create(ATag: TNaNTag);
property Tag: TNaNTag read FTag;
end;
procedure DomainCheck(Err: Boolean);
{ Checksums }
function GetParity(Buffer: PByte; Len: Integer): Boolean;
{ CRC }
function Crc16_P(X: PByteArray; N: Integer; Crc: Word = 0): Word;
function Crc16(const X: array of Byte; N: Integer; Crc: Word = 0): Word;
function Crc16_A(const X: array of Byte; Crc: Word = 0): Word;
function CheckCrc16_P(X: PByteArray; N: Integer; Crc: Word): Integer;
function CheckCrc16(var X: array of Byte; N: Integer; Crc: Word): Integer;
function CheckCrc16_A(var X: array of Byte; Crc: Word): Integer;
function Crc32_P(X: PByteArray; N: Integer; Crc: Cardinal = 0): Cardinal;
function Crc32(const X: array of Byte; N: Integer; Crc: Cardinal = 0): Cardinal;
function Crc32_A(const X: array of Byte; Crc: Cardinal = 0): Cardinal;
function CheckCrc32_P(X: PByteArray; N: Integer; Crc: Cardinal): Integer;
function CheckCrc32(var X: array of Byte; N: Integer; Crc: Cardinal): Integer;
function CheckCrc32_A(var X: array of Byte; Crc: Cardinal): Integer;
{$IFDEF CRCINIT}
procedure InitCrc32(Polynom, Start: Cardinal);
procedure InitCrc16(Polynom, Start: Word);
{$ENDIF CRCINIT}
{ Complex numbers }
type
TRectComplex = record
Re: Float;
Im: Float
end;
TPolarComplex = record
Radius: Float;
Angle: Float
end;
function RectComplex(const Re: Float; const Im: Float = 0): TRectComplex; overload;
function RectComplex(const Z: TPolarComplex): TRectComplex; overload;
function PolarComplex(const Radius: Float; const Angle: Float = 0): TPolarComplex; overload;
function PolarComplex(const Z: TRectComplex): TPolarComplex; overload;
function Equal(const Z1, Z2: TRectComplex): Boolean; overload;
function Equal(const Z1, Z2: TPolarComplex): Boolean; overload;
function IsZero(const Z: TRectComplex): Boolean; overload;
function IsZero(const Z: TPolarComplex): Boolean; overload;
function IsInfinite(const Z: TRectComplex): Boolean; overload;
function IsInfinite(const Z: TPolarComplex): Boolean; overload;
function Norm(const Z: TRectComplex): Float; overload;
function Norm(const Z: TPolarComplex): Float; overload;
function AbsSqr(const Z: TRectComplex): Float; overload;
function AbsSqr(const Z: TPolarComplex): Float; overload;
function Conjugate(const Z: TRectComplex): TRectComplex; overload;
function Conjugate(const Z: TPolarComplex): TPolarComplex; overload;
function Inv(const Z: TRectComplex): TRectComplex; overload;
function Inv(const Z: TPolarComplex): TPolarComplex; overload;
function Neg(const Z: TRectComplex): TRectComplex; overload;
function Neg(const Z: TPolarComplex): TPolarComplex; overload;
function Sum(const Z1, Z2: TRectComplex): TRectComplex; overload;
function Sum(const Z: array of TRectComplex): TRectComplex; overload;
function Diff(const Z1, Z2: TRectComplex): TRectComplex;
function Product(const Z1, Z2: TRectComplex): TRectComplex; overload;
function Product(const Z1, Z2: TPolarComplex): TPolarComplex; overload;
function Quotient(const Z1, Z2: TRectComplex): TRectComplex;
function Ln(const Z: TPolarComplex): TRectComplex;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -