unit ASNum;
{$DEFINE REALCHECK}
{$WARN SYMBOL_PLATFORM OFF}
{$WARN DUPLICATE_CTOR_DTOR OFF}
interface
uses
SysUtils, Types, Math, Generics.Defaults, Generics.Collections, GenHelpers;
threadvar
GTYieldProc: TObjProc;
function AltSgn(const N: Integer): Integer; inline;
type
TInt64Guard = record
class function CanUnMin(const A: Int64): Boolean; static; inline;
class function CanAbs(const A: Int64): Boolean; static; inline;
class function CanAdd(const A, B: Int64): Boolean; static; inline;
class function CanSub(const A, B: Int64): Boolean; static; inline;
class function CanMul(const A, B: Int64): Boolean; static; inline;
class function CanDiv(const A, B: Int64): Boolean; static; inline;
class function CanDivEv(const A, B: Int64): Boolean; static; inline;
class function CanSqr(const A: Int64): Boolean; static; inline;
end;
function SameValue2(const A, B: Extended): Boolean; overload;
function SameValue2(const A, B: Double): Boolean; overload;
function SameValueEx(const A, B: Extended; Epsilon: Extended = 0): Boolean; overload;
function SameValueEx(const A, B: Double; Epsilon: Double = 0): Boolean; overload;
function CreateIntSequence(AStart, AEnd: Integer): TArray<Integer>; overload;
function CreateIntSequence(AStart, AEnd, AStep: Integer): TArray<Integer>; overload;
function CreateIntSequence64(AStart, AEnd: Int64): TArray<Int64>; overload;
function CreateIntSequence64(AStart, AEnd, AStep: Int64): TArray<Int64>; overload;
procedure TranslateIntSequence(var ASeq: TArray<Integer>; const Offset: Integer = -1);
function TranslatedIntSequence(const ASeq: array of Integer; const Offset: Integer = -1): TArray<Integer>;
procedure TranslatePoint(var APoint: TPoint; const Offset: Integer = -1);
function TranslatedPoint(const APoint: TPoint; const Offset: Integer = -1): TPoint;
type
TRange = record
From, &To, Step: Integer;
constructor Create(AFrom, ATo: Integer; AStep: Integer = 1); overload;
constructor Create(ASinglePoint: Integer); overload;
end;
function ParseRangeSeq(const ARanges: array of TRange; const ALength: Integer): TArray<Integer>;
type
EMathException = class(Exception);
{$POINTERMATH ON}
PASR = ^TASR;
{$POINTERMATH OFF}
TASR = Extended;
TASRComparer = class
class var StandardOrder, StandardOrderDescending,
AbsoluteValue, AbsoluteValueDescending: IComparer<TASR>;
class constructor ClassCreate;
end;
TASRArray = TArray<TASR>;
TASR2 = record
X, Y: TASR;
constructor Create(AX, AY: TASR);
end;
TASR3 = packed record
constructor Create(AX, AY, AZ: TASR);
case Boolean of
False:
(Elems: array[0..2] of TASR);
True:
(X, Y, Z: TASR);
end;
TASR4 = packed record
constructor Create(AX, AY, AZ, AW: TASR);
case Boolean of
False:
(Elems: array[0..3] of TASR);
True:
(X, Y, Z, W: TASR);
end;
function IntArrToRealArr(const AArray: TArray<Integer>): TArray<TASR>;
function Int64ArrToRealArr(const AArray: TArray<Int64>): TArray<TASR>;
function IntArrToStrArr(const AArray: TArray<Integer>): TArray<string>;
function Int64ArrToStrArr(const AArray: TArray<Int64>): TArray<string>;
function ASR_COUNT(const A, B: TASR): TASR; inline;
function ASR_PLUS(const A, B: TASR): TASR; inline;
function ASR_TIMES(const A, B: TASR): TASR; inline;
type
TFormatStyle = (fsDefault = 0, fsMonth = 1, fsDayOfWeek = 2, fsMidiInstrument = 3);
TFormatStyleHelper = record helper for TFormatStyle
function ToString: string;
class function FromString(const S: string): TFormatStyle; static;
end;
const
FormatStyleNames: array[TFormatStyle] of string =
('default', 'month', 'day of week', 'MIDI instrument');
type
TNumberFormat = (nfDefault, nfFraction, nfFixed, nfExponent, nfDefExp);
TNumberBase = 2..36;
TNumberBases = set of TNumberBase;
TFormatOptions = record
type
TNumberOptions = record
Base: Byte;
NumberFormat: TNumberFormat;
NumDigits: Word;
MinLength: Byte;
DecimalSeparator: Char;
MinusSign: Char;
IntGrouping: Byte;
IntGropingChar: Char;
FracGrouping: Byte;
FracGroupingChar: Char;
PrettyExp: Boolean;
PreferredBases: TNumberBases;
end;
TComplexOptions = record
ImaginaryUnit: Char;
MultiplicationSign: Char;
PlusSign: Char;
MinusSign: Char;
Spaced: Boolean;
function ImaginarySuffix: string; inline;
function PlusStr: string; inline;
function MinusStr: string; inline;
end;
TVectorOptions = record
MaxLen: Integer;
VerticalUntil: Integer;
BasisChar: Char;
LeftDelim,
ComponentSep,
RightDelim: string;
end;
TStringOptions = record
MaxLen: Integer;
NewLineStr: string;
Quoted: Boolean;
end;
TListOptions = record
MaxLen: Integer;
LeftDelim,
ElementSep,
RightDelim: string;
end;
TStructureOptions = record
AccessSep,
ValueSep,
LeftDelim,
RightDelim,
MemberSep: string;
end;
TSetOptions = record
MaxLen: Integer;
LeftDelim,
ElementSep,
RightDelim,
EmptySetSymbol: string;
end;
var
Numbers: TNumberOptions;
Complex: TComplexOptions;
Vectors: TVectorOptions;
Strings: TStringOptions;
Lists: TListOptions;
Structures: TStructureOptions;
Sets: TSetOptions;
end;
const
DOT_OPERATOR = #$22c5;
PLUS_SIGN = '+';
HYPHEN_MINUS = '-';
MINUS_SIGN = #$2212;
SPACE = #$0020;
FIGURE_SPACE = #$2007;
RETURN_SYMBOL_DAWTL = #$21B2;
DefaultFormatOptions: TFormatOptions =
(
Numbers:
(
Base: 10;
NumberFormat: nfDefault;
NumDigits: 12;
DecimalSeparator: '.';
MinusSign: MINUS_SIGN;
IntGrouping: 0;
IntGropingChar: FIGURE_SPACE;
FracGrouping: 0;
FracGroupingChar: FIGURE_SPACE;
PrettyExp: True;
PreferredBases: [10];
);
Complex:
(
ImaginaryUnit: 'i';
MultiplicationSign: DOT_OPERATOR;
PlusSign: '+';
MinusSign: MINUS_SIGN;
Spaced: True;
);
Vectors:
(
MaxLen: 100;
VerticalUntil: 8;
BasisChar: 'e';
LeftDelim: '(';
ComponentSep: ', ';
RightDelim: ')';
);
Strings:
(
MaxLen: 1000;
NewLineStr: RETURN_SYMBOL_DAWTL;
Quoted: False;
);
Lists:
(
MaxLen: 100;
LeftDelim: '(';
ElementSep: ', ';
RightDelim: ')';
);
Structures:
(
AccessSep: '.';
ValueSep: ': ';
LeftDelim: '(';
RightDelim: ')';
MemberSep: ', ';
);
Sets:
(
MaxLen: 100;
LeftDelim: '{';
ElementSep: ', ';
RightDelim: '}';
EmptySetSymbol: #$2205;
)
);
InputFormOptions: TFormatOptions =
(
Numbers:
(
Base: 10;
NumberFormat: nfDefault;
NumDigits: 18;
DecimalSeparator: '.';
MinusSign: MINUS_SIGN;
IntGrouping: 0;
IntGropingChar: FIGURE_SPACE;
FracGrouping: 0;
FracGroupingChar: FIGURE_SPACE;
PrettyExp: False
);
Complex:
(
ImaginaryUnit: 'i';
MultiplicationSign: DOT_OPERATOR;
PlusSign: '+';
MinusSign: MINUS_SIGN;
Spaced: True;
);
Vectors:
(
MaxLen: MaxInt;
VerticalUntil: 1;
BasisChar: 'e';
LeftDelim: '❨';
ComponentSep: ', ';
RightDelim: '❩';
);
Strings:
(
MaxLen: MaxInt;
NewLineStr: RETURN_SYMBOL_DAWTL;
Quoted: True;
);
Lists:
(
MaxLen: MaxInt;
LeftDelim: '''(';
ElementSep: ', ';
RightDelim: ')';
);
Structures:
(
AccessSep: '.';
ValueSep: ': ';
LeftDelim: 'structure(';
RightDelim: ')';
MemberSep: ', ';
);
Sets:
(
MaxLen: MaxInt;
LeftDelim: '{';
ElementSep: ', ';
RightDelim: '}';
EmptySetSymbol: #$2205;
)
);
ExchangeFormOptions: TFormatOptions =
(
Numbers:
(
Base: 10;
NumberFormat: nfFraction;
NumDigits: 18;
DecimalSeparator: '.';
MinusSign: HYPHEN_MINUS;
IntGrouping: 0;
IntGropingChar: FIGURE_SPACE;
FracGrouping: 0;
FracGroupingChar: FIGURE_SPACE;
PrettyExp: False
);
Complex:
(
ImaginaryUnit: 'i';
MultiplicationSign: DOT_OPERATOR;
PlusSign: '+';
MinusSign: HYPHEN_MINUS;
Spaced: True;
);
Vectors:
(
MaxLen: MaxInt;
VerticalUntil: MaxInt;
BasisChar: 'e';
LeftDelim: '(';
ComponentSep: ', ';
RightDelim: ')';
);
Strings:
(
MaxLen: MaxInt;
NewLineStr: RETURN_SYMBOL_DAWTL;
Quoted: True;
);
Lists:
(
MaxLen: MaxInt;
LeftDelim: '''(';
ElementSep: ', ';
RightDelim: ')';
);
Structures:
(
AccessSep: '.';
ValueSep: ': ';
LeftDelim: '(';
RightDelim: ')';
MemberSep: ', ';
);
Sets:
(
MaxLen: MaxInt;
LeftDelim: '{';
ElementSep: ', ';
RightDelim: '}';
EmptySetSymbol: #$2205;
)
);
function FixNumFmt(const AOptions: TFormatOptions; ANumFmt: TNumberFormat): TFormatOptions;
type
TASRFormatter = function(const AOptions: TFormatOptions; const Val: TASR): string;
PASI = ^TASI;
TASI = Int64;
PRationalNumber = ^TRationalNumber;
TRationalNumber = packed record
Numerator,
Denominator: TASI;
class operator Implicit(A: Integer): TRationalNumber;
class operator Implicit(A: TASI): TRationalNumber;
class operator Implicit(const X: TRationalNumber): TASR; inline;
class operator Negative(const X: TRationalNumber): TRationalNumber;
class operator Add(const X, Y: TRationalNumber): TRationalNumber;
class operator Subtract(const X, Y: TRationalNumber): TRationalNumber;
class operator Multiply(const X, Y: TRationalNumber): TRationalNumber;
class operator Divide(const X, Y: TRationalNumber): TRationalNumber;
class operator Equal(const X, Y: TRationalNumber): Boolean; inline;
class operator NotEqual(const X, Y: TRationalNumber): Boolean; inline;
class function Power(const X: TRationalNumber; const N: TASI): TRationalNumber; static;
constructor Create(const ANumerator, ADenominator: TASI);
procedure ToSimplestForm;
function inv: TRationalNumber;
function sqr: TRationalNumber;
function str: string;
function ToString(const AFormatOptions: TFormatOptions): string; overload;
function ToString(const AFormatOptions: TFormatOptions; const ASymbol: string): string; overload;
function ToDecimalString(const AFormatOptions: TFormatOptions): string;
function valid: Boolean; inline;
function Abs: TRationalNumber;
function Sign: TValueSign; inline;
class procedure ErrNoRep; static;
end;
const
InvalidRat: TRationalNumber = (Numerator: 1; Denominator: 0);
type
TSimpleSymbolicForm = record
A, B: TRationalNumber;
Sym: string;
SymVal: TASR;
constructor Create(const AA, AB: TRationalNumber; const ASym: string; ASymVal: TASR); overload;
constructor Create(pA, qA, pB, qB: TASI; const ASym: string; ASymVal: TASR); overload;
constructor Create(const AA: TRationalNumber); overload;
constructor Create(pA, qA: TASI); overload;
constructor Create(const AB: TRationalNumber; const ASym: string; ASymVal: TASR); overload;
constructor Create(pB, qB: TASI; const ASym: string; ASymVal: TASR); overload;
constructor CreateInvalid(const Val: TASR);
class operator Implicit(const X: TRationalNumber): TSimpleSymbolicForm;
class operator Implicit(const X: TSimpleSymbolicForm): TASR;
class operator Add(const X: TRationalNumber; const Y: TSimpleSymbolicForm): TSimpleSymbolicForm;
class operator Multiply(const X: TRationalNumber; const Y: TSimpleSymbolicForm): TSimpleSymbolicForm;
function str: string;
function sstr: string;
function ToString(const AFormatOptions: TFormatOptions): string;
function valid: Boolean; inline;
procedure MakeInvalid; inline;
end;
TCSimpleSymbolicForm = record
Re,
Im: TSimpleSymbolicForm;
end;
const
Sqrt2 = 1.4142135623730950488016887242096980785696718753769480731766797379;
InvSqrt2 = 1 / Sqrt2;
PiDiv2 = Pi / 2;
PiDiv4 = Pi / 4;
TwoDivPi = 2 / Pi;
TwoPi = 2 * Pi;
SqrtPi = 1.77245385090551602729816748334114518;
SqrtPiInv = 1 / SqrtPi;
Sqrt2Pi = 2.506628274631000502415765284811045253;
Sqrt2PiInv = 1 / 2.506628274631000502415765284811045253;
PiSq = Pi * Pi;
PiCb = Pi * Pi * Pi;
PiInv = 1 / Pi;
PiInvSq = 1 / (Pi * Pi);
EulerConstant = 2.718281828459045235360287471352;
SqrtEulerConstant = 1.6487212707001281468486507878141;
EulerConstantSq = EulerConstant * EulerConstant;
EulerConstantCb = EulerConstant * EulerConstant * EulerConstant;
EulerConstantInv = 1 / EulerConstant;
EulerConstantInvSq = 1 / (EulerConstant * EulerConstant);
GoldenRatio = 1.618033988749894848204586834365638;
GoldenRatioConj = 1 - GoldenRatio;
EulerMascheroni = 0.577215664901532860606512090082;
function IntegerToStr(const x: TASI; const AOptions: TFormatOptions): string;
function RealToStr(const x: TASR; const AOptions: TFormatOptions): string;
function pow(const a, b: TASR): TASR; inline;
function intpow(const a, b: TASI): TASI;
function cbrt(const X: TASR): TASR; inline;
function frt(const X: TASR): TASR; inline;
function tanh(const X: TASR): TASR; inline;
function coth(const X: TASR): TASR; inline;
function sech(const X: TASR): TASR; inline;
function csch(const X: TASR): TASR; inline;
function arcsinh(const X: TASR): TASR; inline;
function arccsch(const X: TASR): TASR; inline;
function sinc(const X: TASR): TASR; inline;
type
{$POINTERMATH ON}
PASC = ^TASC;
{$POINTERMATH OFF}
TASC = packed record
Re, Im: TASR;
class operator Implicit(const r: TASR): TASC;
class operator Positive(const z: TASC): TASC; inline;
class operator Negative(const z: TASC): TASC; inline;
class operator Add(const z1: TASC; const z2: TASC): TASC; inline;
class operator Subtract(const z1: TASC; const z2: TASC): TASC; inline;
class operator Multiply(const z1: TASC; const z2: TASC): TASC; inline;
class operator Divide(const z1: TASC; const z2: TASC): TASC;
class operator Equal(const z1: TASC; const z2: TASC): Boolean; inline;
class operator NotEqual(const z1: TASC; const z2: TASC): Boolean; inline;
class operator Round(const z: TASC): TASC;
class operator Trunc(const z: TASC): TASC;
function Modulus: TASR; inline;
function Argument: TASR; inline;
function Conjugate: TASC; inline;
function Sqr: TASC; inline;
function ModSqr: TASR; inline;
function Inverse: TASC; inline;
function IsReal: Boolean; inline;
function Defuzz(const Eps: Double = 1E-8): TASC;
function str: string;
function pstr: string;
end;
TASCComparer = class
class var ReIm, ReImDescending, Modulus, ModulusDescending, Argument,
ArgumentDescending, ModulusArgument, ModulusArgumentDescending: IComparer<TASC>;
class constructor ClassCreate;
end;
TASCArray = TArray<TASC>;
const
ComplexZero: TASC = (Re: 0; Im: 0);
ImaginaryUnit: TASC = (Re: 0; Im: 1);
ImaginaryUnitDiv2: TASC = (Re: 0; Im : 1/2);
NegativeImaginaryUnit: TASC = (Re: 0; Im: -1);
NegativeImaginaryUnitDiv2: TASC = (Re: 0; Im: -1/2);
function ASC(const Re: TASR; const Im: TASR = 0): TASC; overload; inline;
function ASC_COUNT(const A, B: TASC): TASC; inline;
function ASC_PLUS(const A, B: TASC): TASC; inline;
function ASC_TIMES(const A, B: TASC): TASC; inline;
function ComplexToStr(const z: TASC; ApproxEq: Boolean;
const AOptions: TFormatOptions): string;
function TryStringToComplex(AString: string; out Value: TASC;
ASignRequired: Boolean = False): Boolean;
function CSameValue(const z1, z2: TASC; const Epsilon: TASR = 0): Boolean; overload; inline;
function CSameValueEx(const z1, z2: TASC; const Epsilon: TASR = 0): Boolean; overload; inline;
function CIsZero(const z: TASC; const Epsilon: TASR = 0): Boolean; overload; inline;
function IntegerPowerOfImaginaryUnit(const n: Integer): TASC;
function SameValue2(const A, B: TASC): Boolean; overload; inline;
function CompareValue(const A, B: TASC; Epsilon: Extended = 0): TValueRelationship; overload;
function csign(const z: TASC): TASC; inline;
function cexp(const z: TASC): TASC; inline;
function cln(const z: TASC): TASC; inline;
function clog(const z: TASC): TASC; inline;
function cpow(const z, w: TASC): TASC;
function csqrt(const z: TASC): TASC; inline;
function csin(const z: TASC): TASC; inline;
function ccos(const z: TASC): TASC; inline;
function ctan(const z: TASC): TASC; inline;
function ccot(const z: TASC): TASC; inline;
function csec(const z: TASC): TASC; inline;
function ccsc(const z: TASC): TASC; inline;
function carcsin(const z: TASC): TASC; inline;
function carccos(const z: TASC): TASC; inline;
function carctan(const z: TASC): TASC; inline;
function carccot(const z: TASC): TASC; inline;
function carcsec(const z: TASC): TASC; inline;
function carccsc(const z: TASC): TASC; inline;
function csinh(const z: TASC): TASC; inline;
function ccosh(const z: TASC): TASC; inline;
function ctanh(const z: TASC): TASC; inline;
function ccoth(const z: TASC): TASC; inline;
function csech(const z: TASC): TASC; inline;
function ccsch(const z: TASC): TASC; inline;
function carcsinh(const z: TASC): TASC; inline;
function carccosh(const z: TASC): TASC; inline;
function carctanh(const z: TASC): TASC; inline;
function carccoth(const z: TASC): TASC; inline;
function carcsech(const z: TASC): TASC; inline;
function carccsch(const z: TASC): TASC; inline;
function csinc(const z: TASC): TASC; inline;
type
TIntWithMultiplicity = record
Factor: TASI;
Multiplicity: Integer;
end;
function CollapseWithMultiplicity(const ANumbers: array of TASI): TArray<TIntWithMultiplicity>;
const
intfactorials: array[0..20] of TASI =
(1,
1,
2,
6,
24,
120,
720,
5040,
40320,
362880,
3628800,
39916800,
479001600,
6227020800,
87178291200,
1307674368000,
20922789888000,
355687428096000,
6402373705728000,
121645100408832000,
2432902008176640000);
factorials: array[0..100] of TASR =
(1,
1,
2,
6,
24,
120,
720,
5040,
40320,
362880,
3628800,
39916800,
479001600,
6227020800,
87178291200,
1307674368000,
20922789888000,
355687428096000,
6402373705728000,
121645100408832000,
2.43290200817664E18,
5.109094217170944E19,
1.12400072777760768E21,
2.58520167388849766E22,
6.20448401733239439E23,
1.5511210043330986E25,
4.03291461126605636E26,
1.08888694504183522E28,
3.0488834461171386E29,
8.84176199373970196E30,
2.65252859812191059E32,
8.22283865417792282E33,
2.6313083693369353E35,
8.6833176188118865E36,
2.95232799039604141E38,
1.03331479663861449E40,
3.71993326789901217E41,
1.3763753091226345E43,
5.23022617466601112E44,
2.03978820811974434E46,
8.15915283247897734E47,
3.34525266131638071E49,
1.4050061177528799E51,
6.04152630633738356E52,
2.65827157478844877E54,
1.19622220865480195E56,
5.50262215981208895E57,
2.58623241511168181E59,
1.24139155925360727E61,
6.08281864034267561E62,
3.0414093201713378E64,
1.55111875328738228E66,
8.06581751709438786E67,
4.27488328406002556E69,
2.3084369733924138E71,
1.26964033536582759E73,
7.10998587804863452E74,
4.05269195048772168E76,
2.35056133128287857E78,
1.38683118545689836E80,
8.32098711274139014E81,
5.07580213877224799E83,
3.14699732603879375E85,
1.98260831540444006E87,
1.26886932185884164E89,
8.24765059208247066E90,
5.44344939077443064E92,
3.64711109181886853E94,
2.4800355424368306E96,
1.71122452428141311E98,
1.19785716699698918E100,
8.50478588567862317E101,
6.12344583768860868E103,
4.47011546151268434E105,
3.30788544151938641E107,
2.48091408113953981E109,
1.88549470166605025E111,
1.4518309202828587E113,
1.13242811782062978E115,
8.94618213078297528E116,
7.15694570462638023E118,
5.79712602074736798E120,
4.75364333701284175E122,
3.94552396972065865E124,
3.31424013456535327E126,
2.81710411438055028E128,
2.42270953836727324E130,
2.10775729837952772E132,
1.85482642257398439E134,
1.65079551609084611E136,
1.4857159644817615E138,
1.35200152767840296E140,
1.24384140546413073E142,
1.15677250708164157E144,
1.08736615665674308E146,
1.03299784882390593E148,
9.91677934870949689E149,
9.61927596824821198E151,
9.42689044888324774E153,
9.33262154439441526E155,
9.33262154439441526E157);
type
TRatDigitMode = (rdmSingle, rdmFractional, rdmSignificant, rdmFixedSignificant);
function GetDigit(N: TASI; Index: Integer): Integer; overload;
function GetDigits(N: TASI): TArray<Integer>; overload;
function GetDigit(const R: TRationalNumber; Index: Integer): Integer; overload;
function GetDigits(const R: TRationalNumber; Limit: Integer;
ALimitKind: TRatDigitMode; AFullOutput, ARound: Boolean;
AFirstDigitPos, AFirstSigDigitPos: PInteger): TArray<Integer>; overload;
function GetDigits(const R: TRationalNumber; Limit: Integer;
ALimitKind: TRatDigitMode; AFullOutput, ARound: Boolean;
out AFracDigits: TArray<Integer>): TArray<Integer>; overload;
function GetDigit(X: TASR; Index: Integer): Integer; overload;
function IsInteger(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
function IsIntegerEx(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
function IsInteger32(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
function IsInteger32Ex(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
function IsInteger64(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
function IsInteger64Ex(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
function _fact32(const N: Integer): Integer;
function _fact64(const N: Integer): TASI;
function _factf(const N: Integer): TASR;
function factorial(const N: Integer): TASR;
function combinations(n, k: Integer): TASR;
function intcombinations(n, k: Integer): TASI;
function binomial(const n, k: Integer): TASR; inline;
function permutations(n, k: Integer): TASR;
function intpermutations(n, k: Integer): TASI;
function lcm(const A, B: TASI): TASI; overload;
function lcm(const A, B: UInt64): UInt64; overload;
function TryLCM(const A, B: TASI; var LCM: TASI): Boolean;
function lcm(const Values: array of TASI): TASI; overload;
function gcd(const A, B: TASI): TASI; overload;
function gcd(const A, B: UInt64): UInt64; overload;
function gcd(const Values: array of TASI): TASI; overload;
function coprime(const A, B: TASI): Boolean; inline;
function NaiveTotient(const N: Integer): Integer;
function totient(const N: Integer): Integer;
function cototient(const N: Integer): Integer; inline;
function IsPrime(const N: TASI): Boolean;
function NthPrime(N: Integer): Integer;
function NextPrime(const N: TASI): TASI;
function PreviousPrime(const N: TASI): TASI;
function ExpandPrimeCache(N: Integer): Boolean;
procedure ClearPrimeCache;
function GetPrimeCacheMax: Integer;
function GetPrimeCacheSizeBytes: Integer;
function IsCarolNumber(N: TASI): Boolean;
type
TTSFPrio = (tsfpSimplest, tsfpMostExact);
function PrimePi(const N: Integer): Integer;
const
IntPrimorials: array[0..52] of TASI =
(
1,
1,
2,
6,
6,
30,
30,
210,
210,
210,
210,
2310,
2310,
30030,
30030,
30030,
30030,
510510,
510510,
9699690,
9699690,
9699690,
9699690,
223092870,
223092870,
223092870,
223092870,
223092870,
223092870,
6469693230,
6469693230,
200560490130,
200560490130,
200560490130,
200560490130,
200560490130,
200560490130,
7420738134810,
7420738134810,
7420738134810,
7420738134810,
304250263527210,
304250263527210,
13082761331670030,
13082761331670030,
13082761331670030,
13082761331670030,
614889782588491410,
614889782588491410,
614889782588491410,
614889782588491410,
614889782588491410,
614889782588491410
);
function Primorial(const N: Integer): TASR;
const
IntFibonaccis: array[0..92] of TASI =
(
0,
1,
1,
2,
3,
5,
8,
13,
21,
34,
55,
89,
144,
233,
377,
610,
987,
1597,
2584,
4181,
6765,
10946,
17711,
28657,
46368,
75025,
121393,
196418,
317811,
514229,
832040,
1346269,
2178309,
3524578,
5702887,
9227465,
14930352,
24157817,
39088169,
63245986,
102334155,
165580141,
267914296,
433494437,
701408733,
1134903170,
1836311903,
2971215073,
4807526976,
7778742049,
12586269025,
20365011074,
32951280099,
53316291173,
86267571272,
139583862445,
225851433717,
365435296162,
591286729879,
956722026041,
1548008755920,
2504730781961,
4052739537881,
6557470319842,
10610209857723,
17167680177565,
27777890035288,
44945570212853,
72723460248141,
117669030460994,
190392490709135,
308061521170129,
498454011879264,
806515533049393,
1304969544928657,
2111485077978050,
3416454622906707,
5527939700884757,
8944394323791464,
14472334024676221,
23416728348467685,
37889062373143906,
61305790721611591,
99194853094755497,
160500643816367088,
259695496911122585,
420196140727489673,
679891637638612258,
1100087778366101931,
1779979416004714189,
2880067194370816120,
4660046610375530309,
7540113804746346429
);
function Fibonacci(const N: Integer): TASR;
const
IntLucas: array[0..90] of TASI =
(
2,
1,
3,
4,
7,
11,
18,
29,
47,
76,
123,
199,
322,
521,
843,
1364,
2207,
3571,
5778,
9349,
15127,
24476,
39603,
64079,
103682,
167761,
271443,
439204,
710647,
1149851,
1860498,
3010349,
4870847,
7881196,
12752043,
20633239,
33385282,
54018521,
87403803,
141422324,
228826127,
370248451,
599074578,
969323029,
1568397607,
2537720636,
4106118243,
6643838879,
10749957122,
17393796001,
28143753123,
45537549124,
73681302247,
119218851371,
192900153618,
312119004989,
505019158607,
817138163596,
1322157322203,
2139295485799,
3461452808002,
5600748293801,
9062201101803,
14662949395604,
23725150497407,
38388099893011,
62113250390418,
100501350283429,
162614600673847,
263115950957276,
425730551631123,
688846502588399,
1114577054219522,
1803423556807921,
2918000611027443,
4721424167835364,
7639424778862807,
12360848946698171,
20000273725560978,
32361122672259149,
52361396397820127,
84722519070079276,
137083915467899403,
221806434537978679,
358890350005878082,
580696784543856761,
939587134549734843,
1520283919093591604,
2459871053643326447,
3980154972736918051,
6440026026380244498
);
function Lucas(const N: Integer): TASR;
function Floor64(const X: TASR): TASI;
function Ceil64(const X: TASR): TASI;
function imod(const x, y: Integer): Integer; overload; inline;
function imod(const x, y: TASI): TASI; overload; inline;
function rmod(const x, y: TASR): TASR; inline;
function modulo(const x, y: Integer): Integer; overload; inline;
function modulo(const x, y: TASR): TASR; overload; inline;
function PrimeFactors(const N: TASI; AOnlyUnique: Boolean = False): TArray<TASI>;
function PrimeFactorsWithMultiplicity(N: TASI): TArray<TIntWithMultiplicity>;
function Radical(N: TASI): TASI;
function IsSquareFree(N: TASI): Boolean;
function factorize(const N: TASI): TArray<TASI>;
function GetFactorizedString(const N: TASI; AMinusSign: Char = '-';
AMultiplicationSign: Char = '*'; ASpace: Boolean = False): string; overload;
function GetFactorizedString(const N: TASI; AFancy: Boolean): string; overload;
function divisors(const N: TASI): TArray<TASI>;
procedure FactorAsSquareAsPossible(const N: Integer; out A, B: Integer);
function MöbiusMu(const N: TASI): Integer;
function Mertens(const N: TASI): Integer;
function RationalNumber(const ANumerator, ADenominator: TASI): TRationalNumber;
function ToFraction(const X: TASR): TRationalNumber;
function ToSymbolicForm(const X: TASR; APriority: TTSFPrio = tsfpSimplest): TSimpleSymbolicForm;
function CToSymbolicForm(const z: TASC; APriority: TTSFPrio = tsfpSimplest): TCSimpleSymbolicForm;
function IversonBracket(b: Boolean): Integer; inline;
function KroneckerDelta(i, j: Integer): Integer; overload; inline;
function KroneckerDelta(i, j: TASI): Integer; overload; inline;
function LegendreSymbol(a, p: TASI): Integer;
function JacobiSymbol(a, n: TASI): Integer;
function KroneckerSymbol(a, n: TASI): Integer;
function ContinuedFraction(x: TASR; maxlen: Integer = 18): TArray<TASI>; overload;
function ContinuedFraction(const x: TRationalNumber): TArray<TASI>; overload;
function ContinuedFractionToFraction(const AContinuedFraction: TArray<TASI>): TRationalNumber;
function Heaviside(const X: TASR): TASR; inline;
function Ramp(const X: TASR): TASR; inline;
function Rectfcn(const X: TASR): TASR; inline;
function Tri(const X: TASR): TASR; inline;
function SquareWaveUnit(const X: TASR): TASR;
function SquareWave(const X: TASR): TASR; inline;
function TriangleWaveUnit(const X: TASR): TASR;
function TriangleWave(const X: TASR): TASR; inline;
function SawtoothWaveUnit(const X: TASR): TASR;
function SawtoothWave(const X: TASR): TASR; inline;
type
TRealFunction = function(const X: TASR): TASR;
TComplexFunction = function(const z: TASC): TASC;
TRealFunctionRef = reference to function(const X: TASR): TASR;
TComplexFunctionRef = reference to function(const z: TASC): TASC;
function differentiate(AFunction: TRealFunctionRef; const X: TASR;
const ε: TASR = 1E-6): TASR;
type
TIntegrationParams = record
constructor N(const AN: Integer);
constructor Delta(const ADelta: TASR);
case FixedDelta: Boolean of
False: (FN: Integer);
True: (FDelta: TASR);
end;
const
DefaultIntegrationParams: TIntegrationParams =
(FixedDelta: True; FDelta: 8E-7);
function integrate(AFunction: TRealFunctionRef; a, b: TASR;
const AParams: TIntegrationParams): TASR; overload;
function integrate(AFunction: TRealFunctionRef; a, b: TASR): TASR; overload;
function integrate(AFunction: TComplexFunctionRef; a, b: TASR;
const AParams: TIntegrationParams): TASC; overload;
function integrate(AFunction: TComplexFunctionRef; a, b: TASR): TASC; overload;
type
TIntegrationCacheItem = record
X, Val: TASR;
constructor Create(AX: TASR; AVal: TASR = NaN);
end;
PIntegrationCache = ^TIntegrationCache;
TIntegrationCache = array of TIntegrationCacheItem;
function FillIntegrationCache(AFunction: TRealFunctionRef; a, b: TASR;
const CacheDelta: TASR;
var Cache: TIntegrationCache;
const AParams: TIntegrationParams): TASR; overload;
function FillIntegrationCache(AFunction: TRealFunctionRef; a, b: TASR;
const CacheDelta: TASR;
var Cache: TIntegrationCache): TASR; overload;
procedure StepwiseIntegrationCacheFill(AIntegrand: TRealFunctionRef;
var ACache: TIntegrationCache; ACacheDelta: TASR; const Steps: array of TASR;
const x: TASR; const a: TASR = 0); overload;
procedure StepwiseIntegrationCacheFill(AIntegrand: TRealFunctionRef;
var ACache: TIntegrationCache; ACacheDelta: TASR; const Steps: array of TASR;
const x: TASR; const AParams: TIntegrationParams; const a: TASR = 0); overload;
function ClearIntegrationCache(var ACache: TIntegrationCache): Integer;
function GetIntegrationCacheSize(const ACache: TIntegrationCache): Integer;
function CachedIntegration(AFunction: TRealFunctionRef; a, b: TASR;
const ACache: TIntegrationCache;
const AParams: TIntegrationParams): TASR; overload;
function CachedIntegration(AFunction: TRealFunctionRef; a, b: TASR;
const ACache: TIntegrationCache): TASR; overload;
type
TSequence<T> = reference to function(N: Integer): T;
TPredicate<T> = reference to function(const X: T): Boolean;
TAccumulator<T> = reference to function(const AccumulatedValue, NewValue: T): T;
function sum(const Vals: array of TASI): TASI; overload;
function product(const Vals: array of TASI): TASI; overload;
function accumulate(AFunction: TSequence<TASR>; a, b: Integer;
AStart: TASR; AAccumulator: TAccumulator<TASR>): TASR; overload;
function accumulate(AFunction: TSequence<TASC>; a, b: Integer;
AStart: TASC; AAccumulator: TAccumulator<TASC>): TASC; overload;
function sum(AFunction: TSequence<TASR>; a, b: Integer): TASR; overload;
function sum(AFunction: TSequence<TASC>; a, b: Integer): TASC; overload;
function ArithmeticMean(AFunction: TSequence<TASR>; a, b: Integer): TASR; overload;
function ArithmeticMean(AFunction: TSequence<TASC>; a, b: Integer): TASC; overload;
function GeometricMean(AFunction: TSequence<TASR>; a, b: Integer): TASR; overload;
function GeometricMean(AFunction: TSequence<TASC>; a, b: Integer): TASC; overload;
function HarmonicMean(AFunction: TSequence<TASR>; a, b: Integer): TASR; overload;
function HarmonicMean(AFunction: TSequence<TASC>; a, b: Integer): TASC; overload;
function product(AFunction: TSequence<TASR>; a, b: Integer): TASR; overload;
function product(AFunction: TSequence<TASC>; a, b: Integer): TASC; overload;
function max(AFunction: TSequence<TASR>; a, b: Integer): TASR; overload;
function min(AFunction: TSequence<TASR>; a, b: Integer): TASR; overload;
function exists(AFunction: TSequence<TASR>; a, b: Integer;
APredicate: TPredicate<TASR>): Boolean; overload;
function count(AFunction: TSequence<TASR>; a, b: Integer;
APredicate: TPredicate<TASR>): Integer; overload;
function count(AFunction: TSequence<TASR>; a, b: Integer;
AValue: TASR): Integer; overload;
function ForAll(AFunction: TSequence<TASR>; a, b: Integer;
APredicate: TPredicate<TASR>): Boolean; overload;
function contains(AFunction: TSequence<TASR>; a, b: Integer;
const AValue: TASR): Boolean; overload;
function exists(AFunction: TSequence<TASC>; a, b: Integer;
APredicate: TPredicate<TASC>): Boolean; overload;
function count(AFunction: TSequence<TASC>; a, b: Integer;
APredicate: TPredicate<TASC>): Integer; overload;
function count(AFunction: TSequence<TASC>; a, b: Integer;
AValue: TASC): Integer; overload;
function ForAll(AFunction: TSequence<TASC>; a, b: Integer;
APredicate: TPredicate<TASC>): Boolean; overload;
function contains(AFunction: TSequence<TASC>; a, b: Integer;
const AValue: TASC): Boolean; overload;
function GetSpecialFunctionsIntegrationCacheSize: Integer;
procedure ClearSpecialFunctionsIntegrationCaches;
function erf(const X: TASR): TASR;
function erfc(const X: TASR): TASR; inline;
function FresnelC(const X: TASR): TASR;
function FresnelS(const X: TASR): TASR;
function Si(const X: TASR): TASR;
function Ci(const X: TASR): TASR;
function Li(const X: TASR): TASR;
function LiFrom2(const X: TASR): TASR;
function Bessel(N: Integer; const X: TASR;
const AParams: TIntegrationParams): TASR; overload;
function Bessel(N: Integer; const X: TASR): TASR; overload; inline;
function Bessel(N: Integer; const z: TASC;
const AParams: TIntegrationParams): TASC; overload;
function Bessel(N: Integer; const z: TASC): TASC; overload; inline;
function Laguerre(N: Integer; const X: TASR): TASR; overload;
function Laguerre(N: Integer; const z: TASC): TASC; overload;
function Hermite(N: Integer; const X: TASR): TASR; overload;
function Hermite(N: Integer; const z: TASC): TASC; overload;
function Legendre(N: Integer; const X: TASR): TASR; overload;
function Legendre(N: Integer; const z: TASC): TASC; overload;
function GammaFunction(const X: TASR): TASR; overload;
function GammaFunction(const z: TASC): TASC; overload;
function Chebyshev(N: Integer; const X: TASR): TASR; overload;
function Chebyshev(N: Integer; const z: TASC): TASC; overload;
function Bernstein(I, N: Integer; const X: TASR): TASR; overload; inline;
function Bernstein(I, N: Integer; const z: TASC): TASC; overload; inline;
function HarmonicNumber(const N: TASI): TASR;
type
TRealVector = record
strict private
FComponents: TASRArray;
function GetDimension: Integer; inline;
procedure SetDimension(const Value: Integer); inline;
function GetComponent(Index: Integer): TASR; inline;
procedure SetComponent(Index: Integer; const Value: TASR); inline;
public
constructor Create(const Components: array of TASR); overload;
constructor Create(ASize: Integer; const AVal: TASR = 0); overload;
property Components[Index: Integer]: TASR read GetComponent write SetComponent; default;
property Dimension: Integer read GetDimension write SetDimension;
property Data: TASRArray read FComponents write FComponents;
function ToIntegerArray: TArray<Integer>;
class operator Negative(const u: TRealVector): TRealVector;
class operator Add(const u: TRealVector; const v: TRealVector): TRealVector;
class operator Add(const u: TRealVector; const x: TASR): TRealVector;
class operator Subtract(const u: TRealVector; const v: TRealVector): TRealVector;
class operator Subtract(const u: TRealVector; const x: TASR): TRealVector;
class operator Multiply(const u: TRealVector; const v: TRealVector): TASR;
class operator Multiply(const x: TASR; const u: TRealVector): TRealVector;
class operator Multiply(const u: TRealVector; const x: TASR): TRealVector;
class operator Divide(const u: TRealVector; const x: TASR): TRealVector;
class operator Equal(const u: TRealVector; const v: TRealVector): Boolean;
class operator NotEqual(const u: TRealVector; const v: TRealVector): Boolean; inline;
class operator LeftShift(const u: TRealVector; Val: Integer): TRealVector;
class operator RightShift(const u: TRealVector; Val: Integer): TRealVector;
class operator Trunc(const u: TRealVector): TRealVector;
class operator Round(const u: TRealVector): TRealVector;
class operator LessThan(const u, v: TRealVector): Boolean; inline;
class operator LessThanOrEqual(const u, v: TRealVector): Boolean; inline;
class operator GreaterThan(const u, v: TRealVector): Boolean; inline;
class operator GreaterThanOrEqual(const u, v: TRealVector): Boolean; inline;
function IsZeroVector(const Epsilon: Extended = 0): Boolean;
procedure Normalize;
procedure NormalizeIfNonzero;
function Normalized: TRealVector;
function NormalizedIfNonzero: TRealVector; inline;
function Norm: TASR; inline;
function NormSqr: TASR; inline;
function pNorm(const p: TASR): TASR;
function MaxNorm: TASR; inline;
function SumNorm: TASR; inline;
function kNorm(k: Integer): TASR;
function IsPositive(const Epsilon: Extended = 0): Boolean;
function IsNonNegative(const Epsilon: Extended = 0): Boolean;
function IsNegative(const Epsilon: Extended = 0): Boolean;
function IsNonPositive(const Epsilon: Extended = 0): Boolean;
function Abs: TRealVector;
procedure Append(const AValue: TASR); inline;
procedure ExtendWith(const AValue: TRealVector); inline;
procedure Insert(AIndex: Integer; const AValue: TASR); inline;
procedure Remove(const AIndices: array of Integer);
procedure Swap(AIndex1, AIndex2: Integer); inline;
function TruncateAt(ALength: Integer): TRealVector;
function Reduce(AStep: Integer): TRealVector;
function Clone: TRealVector;
function Subvector(AFrom, ATo: Integer): TRealVector; overload;
function Subvector(const AIndices: array of Integer): TRealVector; overload;
function Sort: TRealVector; overload;
function Sort(AComparer: IComparer<TASR>): TRealVector; overload;
function SafeSort(AComparer: IComparer<TASR>): TRealVector; overload;
function ReverseSort: TRealVector; inline;
function Shuffle: TRealVector;
function Unique: TRealVector;
function UniqueAdj: TRealVector;
function UniqueEps(const Epsilon: TASR = 0): TRealVector;
function UniqueAdjEps(const Epsilon: TASR = 0): TRealVector;
function Reverse: TRealVector;
function Apply(AFunction: TRealFunctionRef): TRealVector;
function Filter(APredicate: TPredicate<TASR>): TRealVector;
function Replace(APredicate: TPredicate<TASR>; const ANewValue: TASR): TRealVector; overload;
function Replace(const AOldValue, ANewValue: TASR; const Epsilon: Extended = 0): TRealVector; overload;
function Replace(const ANewValue: TASR): TRealVector; overload;
procedure RemoveFirst(N: Integer);
function Defuzz(const Eps: Double = 1E-8): TRealVector;
function str(const AOptions: TFormatOptions): string;
end;
function ASC(const u: TRealVector): TASC; overload; inline;
function ASR2(const u1, u2: TASR): TRealVector; overload; inline;
function ASR2(const z: TASC): TRealVector; overload; inline;
function ASR3(const u1, u2, u3: TASR): TRealVector; inline;
function ASR4(const u1, u2, u3, u4: TASR): TRealVector; inline;
function ASR5(const u1, u2, u3, u4, u5: TASR): TRealVector; inline;
function ZeroVector(const Dimension: Integer): TRealVector; inline;
function RandomVector(const Dimension: Integer): TRealVector;
function RandomIntVector(const Dimension, A, B: Integer): TRealVector;
function RandomVectorWithSigns(const Dimension: Integer): TRealVector;
function SequenceVector(ALen: Integer; AStart: TASR = 1; AStep: TASR = 1): TRealVector;
function UnitVector(ADim, AIndex: Integer): TRealVector;
function SameVector(const u, v: TRealVector; const Epsilon: Extended = 0): Boolean; overload;
function SameVectorEx(const u, v: TRealVector; const Epsilon: Extended = 0): Boolean; overload;
function AreParallel(const u, v: TRealVector; Epsilon: Extended = 0): Boolean; overload;
function ArePerpendicular(const u, v: TRealVector; const Epsilon: Extended = 0): Boolean; overload; inline;
function accumulate(const u: TRealVector; const AStart: TASR;
AFunc: TAccumulator<TASR>): TASR; overload;
function sum(const u: TRealVector): TASR; overload;
function ArithmeticMean(const u: TRealVector): TASR; overload; inline;
function GeometricMean(const u: TRealVector): TASR; overload;
function HarmonicMean(const u: TRealVector): TASR; overload;
function product(const u: TRealVector): TASR; overload;
function max(const u: TRealVector): TASR; overload;
function min(const u: TRealVector): TASR; overload;
function exists(const u: TRealVector; APredicate: TPredicate<TASR>): Boolean; overload;
function count(const u: TRealVector; APredicate: TPredicate<TASR>): Integer; overload;
function count(const u: TRealVector; const AValue: TASR): Integer; overload;
function ForAll(const u: TRealVector; APredicate: TPredicate<TASR>): Boolean; overload;
function contains(const u: TRealVector; const AValue: TASR;
const AEpsilon: TASR = 0; ALen: Integer = -1): Boolean; overload;
function CrossProduct(const u, v: TRealVector): TRealVector; overload;
function angle(const u, v: TRealVector): TASR; overload;
function VectConcat(const u, v: TRealVector): TRealVector; overload;
type
TComplexVector = record
strict private
FComponents: TASCArray;
function GetDimension: Integer; inline;
procedure SetDimension(const Value: Integer); inline;
function GetComponent(Index: Integer): TASC; inline;
procedure SetComponent(Index: Integer; const Value: TASC); inline;
public
constructor Create(const Components: array of TASC); overload;
constructor Create(ASize: Integer; const AVal: TASC); overload;
property Components[Index: Integer]: TASC read GetComponent write SetComponent; default;
property Dimension: Integer read GetDimension write SetDimension;
property Data: TASCArray read FComponents write FComponents;
class operator Implicit(const u: TRealVector): TComplexVector;
class operator Negative(const u: TComplexVector): TComplexVector;
class operator Add(const u: TComplexVector; const v: TComplexVector): TComplexVector;
class operator Add(const u: TComplexVector; const z: TASC): TComplexVector;
class operator Subtract(const u: TComplexVector; const v: TComplexVector): TComplexVector;
class operator Subtract(const u: TComplexVector; const z: TASC): TComplexVector;
class operator Multiply(const u: TComplexVector; const v: TComplexVector): TASC;
class operator Multiply(const z: TASC; const u: TComplexVector): TComplexVector;
class operator Multiply(const u: TComplexVector; const z: TASC): TComplexVector;
class operator Divide(const u: TComplexVector; const z: TASC): TComplexVector;
class operator Equal(const u: TComplexVector; const v: TComplexVector): Boolean;
class operator NotEqual(const u: TComplexVector; const v: TComplexVector): Boolean; inline;
class operator LeftShift(const u: TComplexVector; Val: Integer): TComplexVector;
class operator RightShift(const u: TComplexVector; Val: Integer): TComplexVector;
class operator Round(const u: TComplexVector): TComplexVector;
function IsZeroVector(const Epsilon: Extended = 0): Boolean;
procedure Normalize;
procedure NormalizeIfNonzero;
function Normalized: TComplexVector;
function NormalizedIfNonzero: TComplexVector;
function Norm: TASR; inline;
function NormSqr: TASR; inline;
function pNorm(const p: TASR): TASR;
function MaxNorm: TASR; inline;
function SumNorm: TASR; inline;
function kNorm(k: Integer): TASR;
function Abs: TRealVector;
procedure Append(const AValue: TASC); inline;
procedure ExtendWith(const AValue: TComplexVector); inline;
procedure Insert(AIndex: Integer; const AValue: TASC); inline;
procedure Remove(const AIndices: array of Integer);
procedure Swap(AIndex1, AIndex2: Integer); inline;
function TruncateAt(ALength: Integer): TComplexVector;
function Reduce(AStep: Integer): TComplexVector;
function Clone: TComplexVector;
function Subvector(AFrom, ATo: Integer): TComplexVector; overload;
function Subvector(const AIndices: array of Integer): TComplexVector; overload;
function Sort(AComparer: IComparer<TASC>): TComplexVector;
function SafeSort(AComparer: IComparer<TASC>): TComplexVector;
function Shuffle: TComplexVector;
function Unique: TComplexVector;
function UniqueAdj: TComplexVector;
function UniqueEps(const Epsilon: TASR = 0): TComplexVector;
function UniqueAdjEps(const Epsilon: TASR = 0): TComplexVector;
function Reverse: TComplexVector;
function Apply(AFunction: TComplexFunctionRef): TComplexVector;
function Filter(APredicate: TPredicate<TASC>): TComplexVector;
function Replace(APredicate: TPredicate<TASC>;
const ANewValue: TASC): TComplexVector; overload;
function Replace(const AOldValue, ANewValue: TASC;
const Epsilon: Extended = 0): TComplexVector; overload;
function Replace(const ANewValue: TASC): TComplexVector; overload;
function RealPart: TRealVector;
function ImaginaryPart: TRealVector;
procedure RemoveFirst(N: Integer);
function Defuzz(const Eps: Double = 1E-8): TComplexVector;
function str(const AOptions: TFormatOptions): string;
function IsReal: Boolean;
end;
function ASC2(const u1, u2: TASC): TComplexVector; inline;
function ASC3(const u1, u2, u3: TASC): TComplexVector; inline;
function ASC4(const u1, u2, u3, u4: TASC): TComplexVector; inline;
function ASC5(const u1, u2, u3, u4, u5: TASC): TComplexVector; inline;
function ComplexZeroVector(const Dimension: Integer): TComplexVector; inline;
function SameVector(const u, v: TComplexVector; const Epsilon: Extended = 0): Boolean; overload;
function SameVectorEx(const u, v: TComplexVector; const Epsilon: Extended = 0): Boolean; overload;
function AreParallel(const u, v: TComplexVector; Epsilon: Extended = 0): Boolean; overload;
function ArePerpendicular(const u, v: TComplexVector; const Epsilon: Extended = 0): Boolean; overload; inline;
function accumulate(const u: TComplexVector; const AStart: TASC;
AFunc: TAccumulator<TASC>): TASC; overload;
function sum(const u: TComplexVector): TASC; overload;
function ArithmeticMean(const u: TComplexVector): TASC; overload; inline;
function GeometricMean(const u: TComplexVector): TASC; overload;
function HarmonicMean(const u: TComplexVector): TASC; overload;
function product(const u: TComplexVector): TASC; overload;
function exists(const u: TComplexVector; APredicate: TPredicate<TASC>): Boolean; overload;
function count(const u: TComplexVector; APredicate: TPredicate<TASC>): Integer; overload;
function count(const u: TComplexVector; const AValue: TASC): Integer; overload;
function ForAll(const u: TComplexVector; APredicate: TPredicate<TASC>): Boolean; overload;
function contains(const u: TComplexVector; const AValue: TASC;
const AEpsilon: TASR = 0; ALen: Integer = -1): Boolean; overload;
function CrossProduct(const u, v: TComplexVector): TComplexVector; overload;
function angle(const u, v: TComplexVector): TASR; overload;
function VectConcat(const u, v: TComplexVector): TComplexVector; overload;
const
INFTY = High(Integer);
type
TDuplicateFinder<T> = record
class function ContainsDuplicates(Arr: array of T): Boolean; static;
class function PresortedContainsDuplicates(const Arr: array of T): Boolean; static;
end;
TMatrixSize = record
strict private
FRows, FCols: Integer;
private
constructor CreateUnsafe(Rows, Cols: Integer); overload;
function LessenedSize: TMatrixSize; inline;
public
constructor Create(Rows, Cols: Integer); overload;
constructor Create(Size: Integer); overload;
property Rows: Integer read FRows;
property Cols: Integer read FCols;
function ElementCount: Integer; inline;
function SmallestDimension: Integer; inline;
function TransposeSize: TMatrixSize; inline;
function IsSquare: Boolean; inline;
class operator Implicit(S: Integer): TMatrixSize;
class operator Implicit(const ASize: TSize): TMatrixSize;
class operator Implicit(const AMatrixSize: TMatrixSize): TSize;
class operator Equal(const S1, S2: TMatrixSize): Boolean; inline;
class operator NotEqual(const S1, S2: TMatrixSize): Boolean; inline;
class operator Add(const S1, S2: TMatrixSize): TMatrixSize;
end;
const
Mat1x1: TMatrixSize = (FRows: 1; FCols: 1);
Mat2x1: TMatrixSize = (FRows: 2; FCols: 1);
Mat3x1: TMatrixSize = (FRows: 3; FCols: 1);
Mat2x2: TMatrixSize = (FRows: 2; FCols: 2);
Mat3x3: TMatrixSize = (FRows: 3; FCols: 3);
Mat4x4: TMatrixSize = (FRows: 4; FCols: 4);
type
TRowOperationType = (roSwap, roScale, roAddMul);
TMatrixIndexFunction<T> = reference to function(Row, Col: Integer): T;
TIndexArray = array of Integer;
function __sign(const P: TIndexArray): Integer;
type
TRealRowOperationRecord = record
RowOperationType: TRowOperationType;
Row1, Row2: Integer;
Factor: TASR;
end;
TRealRowOpSequence = array of TRealRowOperationRecord;
TRealMatrix = record
strict private
const
InitialRowOpSeqSize = 16;
var
FSize: TMatrixSize;
FElements: TASRArray;
FRowOpSeq: TRealRowOpSequence;
FRowOpCount: Integer;
FRowOpFactor: TASR;
_FCollectRowOpData: Boolean;
procedure SetSize(const Value: TMatrixSize);
function GetElement(Row, Col: Integer): TASR; inline;
procedure SetElement(Row, Col: Integer; const Value: TASR); inline;
function GetRow(ARow: Integer): TRealVector;
function GetCol(ACol: Integer): TRealVector;
procedure _DoSetRow(ARowIndex: Integer; const ARow: TRealVector);
procedure _DoSetCol(AColIndex: Integer; const ACol: TRealVector);
function GetFirstCol: TRealVector; inline;
procedure SetFirstCol(const ACol: TRealVector); inline;
function GetLastCol: TRealVector; inline;
procedure SetLastCol(const ACol: TRealVector); inline;
function GetMainDiagonal: TRealVector;
procedure SetMainDiagonal(const ADiagonal: TRealVector);
function GetSuperDiagonal: TRealVector;
procedure SetSuperDiagonal(const ADiagonal: TRealVector);
function GetSubDiagonal: TRealVector;
procedure SetSubDiagonal(const ADiagonal: TRealVector);
function GetAntiDiagonal: TRealVector;
procedure SetAntiDiagonal(const ADiagonal: TRealVector);
procedure Alloc(ARows, ACols: Integer); overload;
procedure Alloc(ASize: Integer); overload;
procedure Alloc(const ASize: TMatrixSize); overload;
function GetRowData(AIndex: Integer): PASR; inline;
function SafeGetRowData(AIndex: Integer): PASR;
function GetMemorySize: Int64;
procedure BeginCollectRowOpData;
procedure AddRowOpRecord(AType: TRowOperationType; ARow1, ARow2: Integer;
AFactor: TASR = 0);
procedure GetHouseholderMap(const AVect: TRealVector;
out tau, gamma: TASR; out u: TRealVector);
function Eigenvalues2x2: TComplexVector;
private
constructor CreateWithoutAllocation(const ASize: TMatrixSize);
function IsEmpty: Boolean; inline;
procedure DoQuickLU(out A: TRealMatrix; out P: TIndexArray);
procedure RequireNonEmpty; inline;
public
constructor CreateUninitialized(const ASize: TMatrixSize);
constructor Create(const AMatrix: TRealMatrix); overload;
constructor Create(const Elements: array of TASRArray); overload;
constructor Create(const Elements: array of TASR; Cols: Integer = 1); overload;
constructor Create(const ASize: TMatrixSize; const AVal: TASR = 0); overload;
constructor CreateFromRows(const Rows: array of TRealVector);
constructor CreateFromColumns(const Columns: array of TRealVector);
constructor Create(const u, v: TRealVector); overload;
constructor Create(const ASize: TMatrixSize;
AFunction: TMatrixIndexFunction<TASR>); overload;
constructor CreateDiagonal(const Elements: array of TASR); overload;
constructor CreateDiagonal(const Elements: TRealVector); overload;
constructor CreateDiagonal(ASize: Integer; AVal: TASR = 1); overload;
constructor Create(const Blocks: array of TRealMatrix; Cols: Integer = 2); overload;
property Size: TMatrixSize read FSize write SetSize;
property Elements[Row, Col: Integer]: TASR read GetElement write SetElement; default;
property Rows[Index: Integer]: TRealVector read GetRow write _DoSetRow;
property Cols[Index: Integer]: TRealVector read GetCol write _DoSetCol;
property FirstColumn: TRealVector read GetFirstCol write SetFirstCol;
property LastColumn: TRealVector read GetLastCol write SetLastCol;
property MainDiagonal: TRealVector read GetMainDiagonal write SetMainDiagonal;
property SuperDiagonal: TRealVector read GetSuperDiagonal write SetSuperDiagonal;
property SubDiagonal: TRealVector read GetSubDiagonal write SetSubDiagonal;
property AntiDiagonal: TRealVector read GetAntiDiagonal write SetAntiDiagonal;
property Data: TASRArray read FElements write FElements;
property RowData[Index: Integer]: PASR read GetRowData;
property SafeRowData[Index: Integer]: PASR read SafeGetRowData;
property MemorySize: Int64 read GetMemorySize;
class operator Implicit(const u: TRealVector): TRealMatrix;
class operator Explicit(const A: TRealMatrix): TRealVector; inline;
class operator Explicit(X: TASR): TRealMatrix;
class operator Negative(const A: TRealMatrix): TRealMatrix;
class operator Add(const A, B: TRealMatrix): TRealMatrix;
class operator Add(const A: TRealMatrix; const X: TASR): TRealMatrix;
class operator Subtract(const A, B: TRealMatrix): TRealMatrix;
class operator Subtract(const A: TRealMatrix; const X: TASR): TRealMatrix;
class operator Multiply(const A, B: TRealMatrix): TRealMatrix;
class operator Multiply(const X: TASR; const A: TRealMatrix): TRealMatrix;
class operator Multiply(const A: TRealMatrix; const X: TASR): TRealMatrix;
class operator Divide(const A: TRealMatrix; const X: TASR): TRealMatrix;
class operator Equal(const A, B: TRealMatrix): Boolean; inline;
class operator NotEqual(const A, B: TRealMatrix): Boolean; inline;
class operator Trunc(const A: TRealMatrix): TRealMatrix;
class operator Round(const A: TRealMatrix): TRealMatrix;
class operator LessThan(const A, B: TRealMatrix): Boolean; inline;
class operator LessThanOrEqual(const A, B: TRealMatrix): Boolean; inline;
class operator GreaterThan(const A, B: TRealMatrix): Boolean; inline;
class operator GreaterThanOrEqual(const A, B: TRealMatrix): Boolean; inline;
function IsRow: Boolean; inline;
function IsColumn: Boolean; inline;
function IsSquare: Boolean; inline;
function IsIdentity(const Epsilon: Extended = 0): Boolean;
function IsZeroMatrix(const Epsilon: Extended = 0): Boolean;
function IsDiagonal(const Epsilon: Extended = 0): Boolean;
function IsAntiDiagonal(const Epsilon: Extended = 0): Boolean;
function IsReversal(const Epsilon: Extended = 0): Boolean;
function IsUpperTriangular(const Epsilon: Extended = 0): Boolean;
function IsLowerTriangular(const Epsilon: Extended = 0): Boolean;
function IsTriangular(const Epsilon: Extended = 0): Boolean; inline;
function PivotPos(ARow: Integer; const Epsilon: Extended = 0): Integer;
function IsZeroRow(ARow: Integer; const Epsilon: Extended = 0): Boolean; inline;
function IsEssentiallyZeroRow(ARow: Integer;
const Epsilon: Extended = 0): Boolean; inline;
function IsRowEchelonForm(const Epsilon: Extended = 0): Boolean;
function IsReducedRowEchelonForm(const Epsilon: Extended = 0): Boolean;
function IsScalar(const Epsilon: Extended = 0): Boolean;
function IsSymmetric(const Epsilon: Extended = 0): Boolean; inline;
function IsSkewSymmetric(const Epsilon: Extended = 0): Boolean; inline;
function IsOrthogonal(const Epsilon: Extended = 0): Boolean; inline;
function IsNormal(const Epsilon: Extended = 0): Boolean; inline;
function IsBinary(const Epsilon: Extended = 0): Boolean;
function IsPermutation(const Epsilon: Extended = 0): Boolean;
function IsCirculant(const Epsilon: Extended = 0): Boolean;
function IsToeplitz(const Epsilon: Extended = 0): Boolean;
function IsHankel(const Epsilon: Extended = 0): Boolean;
function IsUpperHessenberg(const Epsilon: Extended = 0): Boolean;
function IsLowerHessenberg(const Epsilon: Extended = 0): Boolean;
function IsTridiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsUpperBidiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsLowerBidiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsBidiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsCentrosymmetric(const Epsilon: Extended = 0): Boolean; inline;
function IsVandermonde(const Epsilon: Extended = 0): Boolean;
function CommutesWith(const A: TRealMatrix;
const Epsilon: Extended = 0): Boolean; inline;
function IsIdempotent(const Epsilon: Extended = 0): Boolean; inline;
function IsInvolution(const Epsilon: Extended = 0): Boolean; inline;
function IsPositiveDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsPositiveSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsNegativeDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsNegativeSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsIndefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsNilpotent(const Epsilon: Extended = 0): Boolean; inline;
function NilpotencyIndex(const Epsilon: Extended = 0): Integer;
function IsDiagonallyDominant: Boolean;
function IsStrictlyDiagonallyDominant: Boolean;
function IsPositive(const Epsilon: Extended = 0): Boolean; inline;
function IsNonNegative(const Epsilon: Extended = 0): Boolean; inline;
function IsNegative(const Epsilon: Extended = 0): Boolean; inline;
function IsNonPositive(const Epsilon: Extended = 0): Boolean; inline;
procedure MakeLowerTriangular;
procedure MakeUpperTriangular;
procedure MakeUpperHessenberg;
function Sqr: TRealMatrix; inline;
function Transpose: TRealMatrix;
function HermitianSquare: TRealMatrix; inline;
function Modulus: TRealMatrix; inline;
function Determinant: TASR;
function Trace: TASR;
function Inverse: TRealMatrix;
function TryInvert(out AInverse: TRealMatrix): Boolean;
function Rank: Integer; inline;
function Nullity: Integer; inline;
function ConditionNumber(p: Integer = 2): TASR;
function IsSingular: Boolean;
function Norm: TASR; inline;
function NormSqr: TASR; inline;
function pNorm(const p: TASR): TASR; inline;
function MaxNorm: TASR; inline;
function SumNorm: TASR; inline;
function kNorm(const k: Integer): TASR; inline;
function MaxColSumNorm: TASR;
function MaxRowSumNorm: TASR;
function SpectralNorm: TASR; inline;
function DeletedAbsoluteRowSum(ARow: Integer): TASR;
function RowSwap(ARow1, ARow2: Integer): TRealMatrix;
function RowScale(ARow: Integer; AFactor: TASR): TRealMatrix;
function RowAddMul(ATarget, ASource: Integer; AFactor: TASR;
ADefuzz: Boolean = False; AFirstCol: Integer = 0): TRealMatrix;
function RowOp(const ARowOp: TRealRowOperationRecord): TRealMatrix;
function RowEchelonForm(CollectRowOps: Boolean = False): TRealMatrix;
function ReducedRowEchelonForm(CollectRowOps: Boolean = False): TRealMatrix;
function NumZeroRows(const AEpsilon: TASR = 0): Integer;
function NumTrailingZeroRows(const AEpsilon: TASR = 0): Integer;
function GramSchmidt: TRealMatrix;
procedure InplaceGramSchmidt(FirstCol, LastCol: Integer);
function ColumnSpaceBasis: TRealMatrix;
function ColumnSpaceProjection(const AVector: TRealVector): TRealVector;
function DistanceFromColumnSpace(const AVector: TRealVector): TASR;
function SimilarHessenberg(A2x2Bulge: Boolean = False): TRealMatrix;
function UnsortedEigenvalues: TComplexVector;
function eigenvalues: TComplexVector;
function spectrum: TComplexVector; inline;
function eigenvectors(out AEigenvalues: TRealVector;
out AEigenvectors: TRealMatrix; ASkipVerify: Boolean = False): Boolean;
function IsEigenvector(const u: TRealVector;
const Epsilon: Extended = 0): Boolean;
function EigenvalueOf(const u: TRealVector;
const Epsilon: Extended = 0): TASR;
function TryEigenvalueOf(const u: TRealVector; out AEigenvalue: TASR;
const Epsilon: Extended = 0): Boolean;
function IsEigenpair(const lambda: TASR; const u: TRealVector;
const Epsilon: Extended = 0): Boolean;
function EigenvectorOf(const lambda: TASR): TRealVector;
function SpectralRadius: TASR; inline;
function SingularValues: TRealVector;
function Abs: TRealMatrix;
function Defuzz(const Eps: Double = 1E-8): TRealMatrix;
function Clone: TRealMatrix;
function Vectorization: TRealVector; inline;
function AsVector: TRealVector; inline;
function Augment(const A: TRealMatrix): TRealMatrix; overload;
function Augment: TRealMatrix; overload; inline;
procedure SetRow(ARowIndex: Integer; const ARow: array of TASR);
procedure SetCol(AColIndex: Integer; const ACol: array of TASR);
function Submatrix(ARowToRemove: Integer = 0;
AColToRemove: Integer = 0; AAllowEmpty: Boolean = False): TRealMatrix; overload;
function Submatrix(const ARows: array of Integer;
const ACols: array of Integer): TRealMatrix; overload;
function Submatrix(const ARows: array of Integer): TRealMatrix; overload;
function LeadingPrincipalSubmatrix(const ASize: Integer): TRealMatrix;
function Lessened: TRealMatrix;
function Minor(ARow, ACol: Integer): TASR; inline;
function Cofactor(ARow, ACol: Integer): TASR; inline;
function CofactorMatrix: TRealMatrix;
function AdjugateMatrix: TRealMatrix; inline;
function Apply(AFunction: TRealFunctionRef): TRealMatrix;
function Replace(APredicate: TPredicate<TASR>;
const ANewValue: TASR): TRealMatrix; overload;
function Replace(const AOldValue, ANewValue: TASR;
const Epsilon: Extended = 0): TRealMatrix; overload;
function Replace(const ANewValue: TASR): TRealMatrix; overload;
function str(const AOptions: TFormatOptions): string;
procedure AddRow(const ARow: array of TASR);
function Sort(AComparer: IComparer<TASR>): TRealMatrix;
function SafeSort(AComparer: IComparer<TASR>): TRealMatrix;
function Shuffle: TRealMatrix;
function Reverse: TRealMatrix;
procedure LU(out P, L, U: TRealMatrix);
function Cholesky(out R: TRealMatrix): Boolean;
procedure QR(out Q, R: TRealMatrix);
procedure Hessenberg(out A, U: TRealMatrix);
end;
function mpow(const A: TRealMatrix; const N: Integer): TRealMatrix; overload;
function msqrt(const A: TRealMatrix): TRealMatrix; overload;
function SameMatrix(const A, B: TRealMatrix;
const Epsilon: Extended = 0): Boolean; overload; inline;
function SameMatrixEx(const A, B: TRealMatrix;
const Epsilon: Extended = 0): Boolean; overload; inline;
function ZeroMatrix(const ASize: TMatrixSize): TRealMatrix; inline;
function IdentityMatrix(ASize: Integer): TRealMatrix;
function ReversalMatrix(ASize: Integer): TRealMatrix;
function RandomMatrix(const ASize: TMatrixSize): TRealMatrix;
function RandomIntMatrix(const ASize: TMatrixSize; a, b: Integer): TRealMatrix;
function diag(const AElements: array of TASR): TRealMatrix; overload;
function diag(const AElements: TRealVector): TRealMatrix; overload; inline;
function OuterProduct(const u, v: TRealVector): TRealMatrix; overload; inline;
function CirculantMatrix(const AElements: array of TASR): TRealMatrix; overload;
function ToeplitzMatrix(const AFirstRow,
AFirstCol: array of TASR): TRealMatrix; overload;
function HankelMatrix(const AFirstRow,
ALastCol: array of TASR): TRealMatrix; overload;
function BackwardShiftMatrix(ASize: Integer): TRealMatrix;
function ForwardShiftMatrix(ASize: Integer): TRealMatrix;
function VandermondeMatrix(const AElements: array of TASR;
ACols: Integer = 0): TRealMatrix; overload;
function HilbertMatrix(const ASize: TMatrixSize): TRealMatrix;
function RotationMatrix(AAngle: TASR; ADim: Integer = 2; AIndex1: Integer = 0;
AIndex2: Integer = 1): TRealMatrix; overload;
function RotationMatrix(AAngle: TASR; AAxis: TRealVector): TRealMatrix; overload;
function ReflectionMatrix(const u: TRealVector): TRealMatrix; overload;
function QuickReflectionMatrix(const u: TRealVector): TRealMatrix; overload; inline;
function HadamardProduct(const A, B: TRealMatrix): TRealMatrix; overload;
function DirectSum(const A, B: TRealMatrix): TRealMatrix; overload; inline;
function DirectSum(const Blocks: array of TRealMatrix): TRealMatrix; overload;
function Commute(const A, B: TRealMatrix;
const Epsilon: Extended = 0): Boolean; overload; inline;
function accumulate(const A: TRealMatrix; AStart: TASR;
AAccumulator: TAccumulator<TASR>): TASR; overload; inline;
function sum(const A: TRealMatrix): TASR; overload; inline;
function ArithmeticMean(const A: TRealMatrix): TASR; overload; inline;
function GeometricMean(const A: TRealMatrix): TASR; overload; inline;
function HarmonicMean(const A: TRealMatrix): TASR; overload; inline;
function product(const A: TRealMatrix): TASR; overload; inline;
function max(const A: TRealMatrix): TASR; overload; inline;
function min(const A: TRealMatrix): TASR; overload; inline;
function exists(const A: TRealMatrix; APredicate: TPredicate<TASR>): Boolean; overload;
function count(const A: TRealMatrix; APredicate: TPredicate<TASR>): Integer; overload;
function count(const A: TRealMatrix; const AValue: TASR): Integer; overload; inline;
function ForAll(const A: TRealMatrix; APredicate: TPredicate<TASR>): Boolean; overload;
function contains(const A: TRealMatrix; AValue: TASR): Boolean; overload; inline;
function TryForwardSubstitution(const A: TRealMatrix; const Y: TRealVector;
out Solution: TRealVector; IsUnitDiagonal: Boolean = False): Boolean; overload;
function ForwardSubstitution(const A: TRealMatrix; const Y: TRealVector;
IsUnitDiagonal: Boolean = False): TRealVector; overload;
function TryBackSubstitution(const A: TRealMatrix; const Y: TRealVector;
out Solution: TRealVector): Boolean; overload;
function BackSubstitution(const A: TRealMatrix; const Y: TRealVector): TRealVector; overload;
function SysSolve(const AAugmented: TRealMatrix): TRealVector; overload; inline;
function TrySysSolve(const A: TRealMatrix; const Y: TRealVector;
out Solution: TRealVector): Boolean; overload;
function TrySysSolve(const A: TRealMatrix; const Y: TRealMatrix;
out Solution: TRealMatrix): Boolean; overload;
function SysSolve(const A: TRealMatrix; const Y: TRealVector): TRealVector; overload;
function SysSolve(const A: TRealMatrix; const Y: TRealMatrix): TRealMatrix; overload;
function LeastSquaresPolynomialFit(const X, Y: TRealVector;
ADegree: Integer): TRealVector; overload;
procedure VectMove(const ASource: TRealVector; const AFrom, ATo: Integer;
var ATarget: TRealVector; const ATargetFrom: Integer = 0); overload; inline;
procedure VectMoveToMatCol(const ASource: TRealVector; const AFrom, ATo: Integer;
var ATarget: TRealMatrix; const ATargetCol: Integer;
const ATargetFirstRow: Integer = 0); overload;
procedure VectMoveToMatRow(const ASource: TRealVector; const AFrom, ATo: Integer;
var ATarget: TRealMatrix; const ATargetRow: Integer;
const ATargetFirstCol: Integer = 0); overload; inline;
procedure MatMove(const ASource: TRealMatrix; const ARect: TRect;
var ATarget: TRealMatrix; const ATargetTopLeft: TPoint); overload;
procedure MatMove(const ASource: TRealMatrix; var ATarget: TRealMatrix;
const ATargetTopLeft: TPoint); overload; inline;
procedure MatMoveColToVect(const ASource: TRealMatrix; const AColumn: Integer;
const AFrom, ATo: Integer; var ATarget: TRealVector;
const ATargetFrom: Integer = 0); overload;
procedure MatMoveRowToVect(const ASource: TRealMatrix; const ARow: Integer;
const AFrom, ATo: Integer; var ATarget: TRealVector;
const ATargetFrom: Integer = 0); overload; inline;
procedure MatMulBlockInplaceL(var ATarget: TRealMatrix; const ARect: TRect;
const AFactor: TRealMatrix); overload;
procedure MatMulBlockInplaceL(var ATarget: TRealMatrix; const ATopLeft: TPoint;
const AFactor: TRealMatrix); overload; inline;
procedure MatMulBottomRight(var ATarget: TRealMatrix; const AFactor: TRealMatrix); overload;
procedure MatBlockFill(var ATarget: TRealMatrix; const ARect: TRect;
const Value: TASR); overload;
type
TComplexRowOperationRecord = record
RowOperationType: TRowOperationType;
Row1, Row2: Integer;
Factor: TASC;
end;
TComplexRowOpSequence = array of TComplexRowOperationRecord;
TComplexMatrix = record
strict private
const
InitialRowOpSeqSize = 16;
var
FSize: TMatrixSize;
FElements: TASCArray;
FRowOpSeq: TComplexRowOpSequence;
FRowOpCount: Integer;
FRowOpFactor: TASC;
_FCollectRowOpData: Boolean;
procedure SetSize(const Value: TMatrixSize);
function GetElement(Row, Col: Integer): TASC; inline;
procedure SetElement(Row, Col: Integer; const Value: TASC); inline;
function GetRow(ARow: Integer): TComplexVector;
function GetCol(ACol: Integer): TComplexVector;
procedure _DoSetRow(ARowIndex: Integer; const ARow: TComplexVector);
procedure _DoSetCol(AColIndex: Integer; const ACol: TComplexVector);
function GetFirstCol: TComplexVector; inline;
procedure SetFirstCol(const ACol: TComplexVector); inline;
function GetLastCol: TComplexVector; inline;
procedure SetLastCol(const ACol: TComplexVector); inline;
function GetMainDiagonal: TComplexVector;
procedure SetMainDiagonal(const ADiagonal: TComplexVector);
function GetSuperDiagonal: TComplexVector;
procedure SetSuperDiagonal(const ADiagonal: TComplexVector);
function GetSubDiagonal: TComplexVector;
procedure SetSubDiagonal(const ADiagonal: TComplexVector);
function GetAntiDiagonal: TComplexVector;
procedure SetAntiDiagonal(const ADiagonal: TComplexVector);
procedure Alloc(ARows, ACols: Integer); overload;
procedure Alloc(ASize: Integer); overload;
procedure Alloc(const ASize: TMatrixSize); overload;
function GetRowData(AIndex: Integer): PASC; inline;
function SafeGetRowData(AIndex: Integer): PASC;
function GetMemorySize: Int64;
procedure BeginCollectRowOpData;
procedure AddRowOpRecord(AType: TRowOperationType; ARow1, ARow2: Integer;
AFactor: TASC);
procedure GetHouseholderMap(const AVect: TComplexVector;
out tau, gamma: TASC; out u: TComplexVector);
function Eigenvalues2x2: TComplexVector;
private
function IsEmpty: Boolean; inline;
procedure DoQuickLU(out A: TComplexMatrix; out P: TIndexArray);
procedure RequireNonEmpty; inline;
public
constructor CreateUninitialized(const ASize: TMatrixSize);
constructor Create(const AMatrix: TComplexMatrix); overload;
constructor Create(const Elements: array of TASCArray); overload;
constructor Create(const Elements: array of TASC; Cols: Integer = 1); overload;
constructor Create(const ASize: TMatrixSize; const AVal: TASC); overload;
constructor CreateFromRows(const Rows: array of TComplexVector);
constructor CreateFromColumns(const Columns: array of TComplexVector);
constructor Create(const u, v: TComplexVector); overload;
constructor Create(const ASize: TMatrixSize;
AFunction: TMatrixIndexFunction<TASC>); overload;
constructor CreateDiagonal(const Elements: array of TASC); overload;
constructor CreateDiagonal(const Elements: TComplexVector); overload;
constructor CreateDiagonal(ASize: Integer; AVal: TASC); overload;
constructor Create(const Blocks: array of TComplexMatrix; Cols: Integer = 2); overload;
property Size: TMatrixSize read FSize write SetSize;
property Elements[Row, Col: Integer]: TASC read GetElement write SetElement; default;
property Rows[Index: Integer]: TComplexVector read GetRow write _DoSetRow;
property Cols[Index: Integer]: TComplexVector read GetCol write _DoSetCol;
property FirstColumn: TComplexVector read GetFirstCol write SetFirstCol;
property LastColumn: TComplexVector read GetLastCol write SetLastCol;
property MainDiagonal: TComplexVector read GetMainDiagonal write SetMainDiagonal;
property SuperDiagonal: TComplexVector read GetSuperDiagonal write SetSuperDiagonal;
property SubDiagonal: TComplexVector read GetSubDiagonal write SetSubDiagonal;
property AntiDiagonal: TComplexVector read GetAntiDiagonal write SetAntiDiagonal;
property Data: TASCArray read FElements write FElements;
property RowData[Index: Integer]: PASC read GetRowData;
property SafeRowData[Index: Integer]: PASC read SafeGetRowData;
property MemorySize: Int64 read GetMemorySize;
class operator Implicit(const u: TComplexVector): TComplexMatrix;
class operator Implicit(const A: TRealMatrix): TComplexMatrix;
class operator Explicit(const A: TComplexMatrix): TComplexVector; inline;
class operator Explicit(X: TASC): TComplexMatrix;
class operator Negative(const A: TComplexMatrix): TComplexMatrix;
class operator Add(const A, B: TComplexMatrix): TComplexMatrix;
class operator Add(const A: TComplexMatrix; const X: TASC): TComplexMatrix;
class operator Subtract(const A, B: TComplexMatrix): TComplexMatrix;
class operator Subtract(const A: TComplexMatrix; const X: TASC): TComplexMatrix;
class operator Multiply(const A, B: TComplexMatrix): TComplexMatrix;
class operator Multiply(const X: TASC; const A: TComplexMatrix): TComplexMatrix;
class operator Multiply(const A: TComplexMatrix; const X: TASC): TComplexMatrix;
class operator Divide(const A: TComplexMatrix; const X: TASC): TComplexMatrix;
class operator Equal(const A, B: TComplexMatrix): Boolean; inline;
class operator NotEqual(const A, B: TComplexMatrix): Boolean; inline;
class operator Trunc(const A: TComplexMatrix): TComplexMatrix;
class operator Round(const A: TComplexMatrix): TComplexMatrix;
function IsRow: Boolean; inline;
function IsColumn: Boolean; inline;
function IsSquare: Boolean; inline;
function IsIdentity(const Epsilon: Extended = 0): Boolean;
function IsZeroMatrix(const Epsilon: Extended = 0): Boolean;
function IsDiagonal(const Epsilon: Extended = 0): Boolean;
function IsAntiDiagonal(const Epsilon: Extended = 0): Boolean;
function IsReversal(const Epsilon: Extended = 0): Boolean;
function IsUpperTriangular(const Epsilon: Extended = 0): Boolean;
function IsLowerTriangular(const Epsilon: Extended = 0): Boolean;
function IsTriangular(const Epsilon: Extended = 0): Boolean; inline;
function PivotPos(ARow: Integer; const Epsilon: Extended = 0): Integer;
function IsZeroRow(ARow: Integer; const Epsilon: Extended = 0): Boolean; inline;
function IsEssentiallyZeroRow(ARow: Integer;
const Epsilon: Extended = 0): Boolean; inline;
function IsRowEchelonForm(const Epsilon: Extended = 0): Boolean;
function IsReducedRowEchelonForm(const Epsilon: Extended = 0): Boolean;
function IsScalar(const Epsilon: Extended = 0): Boolean;
function IsSymmetric(const Epsilon: Extended = 0): Boolean; inline;
function IsSkewSymmetric(const Epsilon: Extended = 0): Boolean; inline;
function IsHermitian(const Epsilon: Extended = 0): Boolean; inline;
function IsSkewHermitian(const Epsilon: Extended = 0): Boolean; inline;
function IsOrthogonal(const Epsilon: Extended = 0): Boolean; inline;
function IsUnitary(const Epsilon: Extended = 0): Boolean; inline;
function IsNormal(const Epsilon: Extended = 0): Boolean; inline;
function IsBinary(const Epsilon: Extended = 0): Boolean;
function IsPermutation(const Epsilon: Extended = 0): Boolean;
function IsCirculant(const Epsilon: Extended = 0): Boolean;
function IsToeplitz(const Epsilon: Extended = 0): Boolean;
function IsHankel(const Epsilon: Extended = 0): Boolean;
function IsUpperHessenberg(const Epsilon: Extended = 0): Boolean;
function IsLowerHessenberg(const Epsilon: Extended = 0): Boolean;
function IsTridiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsUpperBidiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsLowerBidiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsBidiagonal(const Epsilon: Extended = 0): Boolean; inline;
function IsCentrosymmetric(const Epsilon: Extended = 0): Boolean; inline;
function IsVandermonde(const Epsilon: Extended = 0): Boolean;
function CommutesWith(const A: TComplexMatrix;
const Epsilon: Extended = 0): Boolean; inline;
function IsIdempotent(const Epsilon: Extended = 0): Boolean; inline;
function IsInvolution(const Epsilon: Extended = 0): Boolean; inline;
function IsPositiveDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsPositiveSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsNegativeDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsNegativeSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsIndefinite(const Epsilon: Extended = 0): Boolean; inline;
function IsNilpotent(const Epsilon: Extended = 0): Boolean; inline;
function NilpotencyIndex(const Epsilon: Extended = 0): Integer;
function IsDiagonallyDominant: Boolean;
function IsStrictlyDiagonallyDominant: Boolean;
function RealPart: TRealMatrix;
function ImaginaryPart: TRealMatrix;
procedure MakeLowerTriangular;
procedure MakeUpperTriangular;
procedure MakeUpperHessenberg;
function Sqr: TComplexMatrix; inline;
function Transpose: TComplexMatrix;
function Adjoint: TComplexMatrix;
function HermitianSquare: TComplexMatrix; inline;
function Modulus: TComplexMatrix; inline;
function Determinant: TASC;
function Trace: TASC;
function Inverse: TComplexMatrix;
function TryInvert(out AInverse: TComplexMatrix): Boolean;
function Rank: Integer; inline;
function Nullity: Integer; inline;
function ConditionNumber(p: Integer = 2): TASR;
function IsSingular: Boolean;
function Norm: TASR; inline;
function NormSqr: TASR; inline;
function pNorm(const p: TASR): TASR; inline;
function MaxNorm: TASR; inline;
function SumNorm: TASR; inline;
function kNorm(const k: Integer): TASR; inline;
function MaxColSumNorm: TASR;
function MaxRowSumNorm: TASR;
function SpectralNorm: TASR; inline;
function DeletedAbsoluteRowSum(ARow: Integer): TASR;
function RowSwap(ARow1, ARow2: Integer): TComplexMatrix;
function RowScale(ARow: Integer; AFactor: TASC): TComplexMatrix;
function RowAddMul(ATarget, ASource: Integer; AFactor: TASC;
ADefuzz: Boolean = False; AFirstCol: Integer = 0): TComplexMatrix;
function RowOp(const ARowOp: TComplexRowOperationRecord): TComplexMatrix;
function RowEchelonForm(CollectRowOps: Boolean = False): TComplexMatrix;
function ReducedRowEchelonForm(CollectRowOps: Boolean = False): TComplexMatrix;
function NumZeroRows(const AEpsilon: TASR = 0): Integer;
function NumTrailingZeroRows(const AEpsilon: TASR = 0): Integer;
function GramSchmidt: TComplexMatrix;
procedure InplaceGramSchmidt(FirstCol, LastCol: Integer);
function ColumnSpaceBasis: TComplexMatrix;
function ColumnSpaceProjection(const AVector: TComplexVector): TComplexVector;
function DistanceFromColumnSpace(const AVector: TComplexVector): TASR;
function SimilarHessenberg(A2x2Bulge: Boolean = False): TComplexMatrix;
function UnsortedEigenvalues: TComplexVector;
function eigenvalues: TComplexVector;
function spectrum: TComplexVector; inline;
function eigenvectors(out AEigenvalues: TComplexVector;
out AEigenvectors: TComplexMatrix; ASkipVerify: Boolean = False): Boolean;
function IsEigenvector(const u: TComplexVector;
const Epsilon: Extended = 0): Boolean;
function EigenvalueOf(const u: TComplexVector;
const Epsilon: Extended = 0): TASC;
function TryEigenvalueOf(const u: TComplexVector; out AEigenvalue: TASC;
const Epsilon: Extended = 0): Boolean;
function IsEigenpair(const lambda: TASC; const u: TComplexVector;
const Epsilon: Extended = 0): Boolean;
function EigenvectorOf(const lambda: TASC): TComplexVector;
function SpectralRadius: TASR; inline;
function SingularValues: TRealVector;
function Abs: TRealMatrix;
function Defuzz(const Eps: Double = 1E-8): TComplexMatrix;
function Clone: TComplexMatrix;
function Vectorization: TComplexVector; inline;
function AsVector: TComplexVector; inline;
function Augment(const A: TComplexMatrix): TComplexMatrix; overload;
function Augment: TComplexMatrix; overload; inline;
procedure SetRow(ARowIndex: Integer; const ARow: array of TASC);
procedure SetCol(AColIndex: Integer; const ACol: array of TASC);
function Submatrix(ARowToRemove: Integer = 0;
AColToRemove: Integer = 0; AAllowEmpty: Boolean = False): TComplexMatrix; overload;
function Submatrix(const ARows: array of Integer;
const ACols: array of Integer): TComplexMatrix; overload;
function Submatrix(const ARows: array of Integer): TComplexMatrix; overload;
function LeadingPrincipalSubmatrix(const ASize: Integer): TComplexMatrix;
function Lessened: TComplexMatrix;
function Minor(ARow, ACol: Integer): TASC; inline;
function Cofactor(ARow, ACol: Integer): TASC; inline;
function CofactorMatrix: TComplexMatrix;
function AdjugateMatrix: TComplexMatrix; inline;
function Apply(AFunction: TComplexFunctionRef): TComplexMatrix;
function Replace(APredicate: TPredicate<TASC>;
const ANewValue: TASC): TComplexMatrix; overload;
function Replace(const AOldValue, ANewValue: TASC;
const Epsilon: Extended = 0): TComplexMatrix; overload;
function Replace(const ANewValue: TASC): TComplexMatrix; overload;
function str(const AOptions: TFormatOptions): string;
procedure AddRow(const ARow: array of TASC);
function Sort(AComparer: IComparer<TASC>): TComplexMatrix;
function SafeSort(AComparer: IComparer<TASC>): TComplexMatrix;
function Shuffle: TComplexMatrix;
function Reverse: TComplexMatrix;
procedure LU(out P, L, U: TComplexMatrix);
function Cholesky(out R: TComplexMatrix): Boolean;
procedure QR(out Q, R: TComplexMatrix);
procedure Hessenberg(out A, U: TComplexMatrix);
end;
function mpow(const A: TComplexMatrix; const N: Integer): TComplexMatrix; overload;
function msqrt(const A: TComplexMatrix): TComplexMatrix; overload;
function SameMatrix(const A, B: TComplexMatrix;
const Epsilon: Extended = 0): Boolean; overload; inline;
function SameMatrixEx(const A, B: TComplexMatrix;
const Epsilon: Extended = 0): Boolean; overload; inline;
function ComplexZeroMatrix(const ASize: TMatrixSize): TComplexMatrix; inline;
function ComplexIdentityMatrix(ASize: Integer): TComplexMatrix;
function ComplexReversalMatrix(ASize: Integer): TComplexMatrix;
function diag(const AElements: array of TASC): TComplexMatrix; overload;
function diag(const AElements: TComplexVector): TComplexMatrix; inline; overload;
function OuterProduct(const u, v: TComplexVector): TComplexMatrix; overload; inline;
function CirculantMatrix(const AElements: array of TASC): TComplexMatrix; overload;
function ToeplitzMatrix(const AFirstRow,
AFirstCol: array of TASC): TComplexMatrix; overload;
function HankelMatrix(const AFirstRow,
ALastCol: array of TASC): TComplexMatrix; overload;
function VandermondeMatrix(const AElements: array of TASC;
ACols: Integer = 0): TComplexMatrix; overload;
function ReflectionMatrix(const u: TComplexVector): TComplexMatrix; overload;
function QuickReflectionMatrix(const u: TComplexVector): TComplexMatrix; overload; inline;
function HadamardProduct(const A, B: TComplexMatrix): TComplexMatrix; overload;
function DirectSum(const A, B: TComplexMatrix): TComplexMatrix; overload; inline;
function DirectSum(const Blocks: array of TComplexMatrix): TComplexMatrix; overload;
function Commute(const A, B: TComplexMatrix;
const Epsilon: Extended = 0): Boolean; overload; inline;
function accumulate(const A: TComplexMatrix; const AStart: TASC;
AAccumulator: TAccumulator<TASC>): TASC; overload; inline;
function sum(const A: TComplexMatrix): TASC; overload; inline;
function ArithmeticMean(const A: TComplexMatrix): TASC; overload; inline;
function GeometricMean(const A: TComplexMatrix): TASC; overload; inline;
function HarmonicMean(const A: TComplexMatrix): TASC; overload; inline;
function product(const A: TComplexMatrix): TASC; overload; inline;
function exists(const A: TComplexMatrix; APredicate: TPredicate<TASC>): Boolean; overload;
function count(const A: TComplexMatrix; APredicate: TPredicate<TASC>): Integer; overload;
function count(const A: TComplexMatrix; const AValue: TASC): Integer; overload; inline;
function ForAll(const A: TComplexMatrix; APredicate: TPredicate<TASC>): Boolean; overload;
function contains(const A: TComplexMatrix; AValue: TASC): Boolean; overload; inline;
function TryForwardSubstitution(const A: TComplexMatrix; const Y: TComplexVector;
out Solution: TComplexVector; IsUnitDiagonal: Boolean = False): Boolean; overload;
function ForwardSubstitution(const A: TComplexMatrix; const Y: TComplexVector;
IsUnitDiagonal: Boolean = False): TComplexVector; overload;
function TryBackSubstitution(const A: TComplexMatrix; const Y: TComplexVector;
out Solution: TComplexVector): Boolean; overload;
function BackSubstitution(const A: TComplexMatrix; const Y: TComplexVector): TComplexVector; overload;
function SysSolve(const AAugmented: TComplexMatrix): TComplexVector; overload; inline;
function TrySysSolve(const A: TComplexMatrix; const Y: TComplexVector;
out Solution: TComplexVector): Boolean; overload;
function TrySysSolve(const A: TComplexMatrix; const Y: TComplexMatrix;
out Solution: TComplexMatrix): Boolean; overload;
function SysSolve(const A: TComplexMatrix; const Y: TComplexVector): TComplexVector; overload;
function SysSolve(const A: TComplexMatrix; const Y: TComplexMatrix): TComplexMatrix; overload;
function LeastSquaresPolynomialFit(const X, Y: TComplexVector;
ADegree: Integer): TComplexVector; overload;
procedure VectMove(const ASource: TComplexVector; const AFrom, ATo: Integer;
var ATarget: TComplexVector; const ATargetFrom: Integer = 0); overload; inline;
procedure VectMoveToMatCol(const ASource: TComplexVector; const AFrom, ATo: Integer;
var ATarget: TComplexMatrix; const ATargetCol: Integer;
const ATargetFirstRow: Integer = 0); overload;
procedure VectMoveToMatRow(const ASource: TComplexVector; const AFrom, ATo: Integer;
var ATarget: TComplexMatrix; const ATargetRow: Integer;
const ATargetFirstCol: Integer = 0); overload; inline;
procedure MatMove(const ASource: TComplexMatrix; const ARect: TRect;
var ATarget: TComplexMatrix; const ATargetTopLeft: TPoint); overload;
procedure MatMove(const ASource: TComplexMatrix; var ATarget: TComplexMatrix;
const ATargetTopLeft: TPoint); overload; inline;
procedure MatMoveColToVect(const ASource: TComplexMatrix; const AColumn: Integer;
const AFrom, ATo: Integer; var ATarget: TComplexVector;
const ATargetFrom: Integer = 0); overload;
procedure MatMoveRowToVect(const ASource: TComplexMatrix; const ARow: Integer;
const AFrom, ATo: Integer; var ATarget: TComplexVector;
const ATargetFrom: Integer = 0); overload; inline;
procedure MatMulBlockInplaceL(var ATarget: TComplexMatrix; const ARect: TRect;
const AFactor: TComplexMatrix); overload;
procedure MatMulBlockInplaceL(var ATarget: TComplexMatrix; const ATopLeft: TPoint;
const AFactor: TComplexMatrix); overload; inline;
procedure MatMulBottomRight(var ATarget: TComplexMatrix; const AFactor: TComplexMatrix); overload;
procedure MatBlockFill(var ATarget: TComplexMatrix; const ARect: TRect;
const Value: TASC); overload;
function sqrt(const X: TASR): TASR; overload; inline;
function sqrt(const z: TASC): TASC; overload; inline;
function sqrt(const A: TRealMatrix): TRealMatrix; overload; inline;
function sqrt(const A: TComplexMatrix): TComplexMatrix; overload; inline;
const
RealCheck = {$IFDEF REALCHECK}True{$ELSE}False{$ENDIF};
var
MaxThreadCount: Integer = 8;
function GetTotalCacheSize: Integer;
procedure ClearCaches;
const
InvLn10: TASC = (Re: 0.434294481903251827651128918916605082294397005803666566114; Im: 0);
function DoubleListToASR2s(const LList: TArray<Double>): TArray<TASR2>;
function DoubleListToASR3s(const LList: TArray<Double>): TArray<TASR3>;
implementation
uses
Classes, StrUtils, Windows, Character, WideStrUtils;
var
FS: TFormatSettings;
procedure DoYield; inline;
begin
if Assigned(GTYieldProc) then
GTYieldProc;
end;
class function TInt64Guard.CanUnMin(const A: Int64): Boolean;
begin
Result := A <> A.MinValue;
end;
class function TInt64Guard.CanAbs(const A: Int64): Boolean;
begin
Result := A <> A.MinValue;
end;
class function TInt64Guard.CanAdd(const A, B: Int64): Boolean;
const
SafeMin = Int64.MinValue div 2 + 1;
SafeMax = Int64.MaxValue div 2 - 1;
begin
Result := InRange(A, SafeMin, SafeMax) and InRange(A, SafeMin, SafeMax);
end;
class function TInt64Guard.CanSub(const A, B: Int64): Boolean;
const
SafeMin = Int64.MinValue div 2 + 1;
SafeMax = Int64.MaxValue div 2 - 1;
begin
Result := InRange(A, SafeMin, SafeMax) and InRange(A, SafeMin, SafeMax);
end;
class function TInt64Guard.CanMul(const A, B: Int64): Boolean;
const
SafeMin = Extended(Int64.MinValue) / 2;
SafeMax = Extended(Int64.MaxValue) / 2;
begin
Result := InRange(Extended(A) * Extended(B), SafeMin, SafeMax);
end;
class function TInt64Guard.CanDiv(const A, B: Int64): Boolean;
begin
Result := (A <> A.MinValue) or (B <> -1);
end;
class function TInt64Guard.CanDivEv(const A, B: Int64): Boolean;
begin
Result := CanDiv(A, B) and (A mod B = 0);
end;
class function TInt64Guard.CanSqr(const A: Int64): Boolean;
begin
Result := InRange(A, -3037000499, 3037000499);
end;
function GetThreadCount: Integer; inline;
begin
Result := Min(CPUCount, MaxThreadCount);
end;
const
FuzzFactor = 1000;
ExtendedResolution = 1E-19 * FuzzFactor;
DoubleResolution = 1E-15 * FuzzFactor;
function SameValue2(const A, B: Extended): Boolean;
var
Epsilon: Extended;
begin
Epsilon := Min(Abs(A), Abs(B)) * ExtendedResolution;
if A > B then
Result := (A - B) <= Epsilon
else
Result := (B - A) <= Epsilon;
end;
function SameValue2(const A, B: Double): Boolean;
var
Epsilon: Double;
begin
Epsilon := Min(Abs(A), Abs(B)) * DoubleResolution;
if A > B then
Result := (A - B) <= Epsilon
else
Result := (B - A) <= Epsilon;
end;
function SameValueEx(const A, B: Extended; Epsilon: Extended): Boolean; overload;
begin
if Epsilon = 0 then
Epsilon := Max(Min(Abs(A), Abs(B)) * ExtendedResolution, ExtendedResolution)
else
Epsilon := Max(Min(Abs(A), Abs(B)) * Epsilon, Epsilon);
if A > B then
Result := (A - B) <= Epsilon
else
Result := (B - A) <= Epsilon;
end;
function SameValueEx(const A, B: Double; Epsilon: Double): Boolean; overload;
begin
if Epsilon = 0 then
Epsilon := Max(Min(Abs(A), Abs(B)) * DoubleResolution, DoubleResolution)
else
Epsilon := Max(Min(Abs(A), Abs(B)) * Epsilon, Epsilon);
if A > B then
Result := (A - B) <= Epsilon
else
Result := (B - A) <= Epsilon;
end;
function inv(const X: TASR): TASR; inline;
begin
Result := 1 / X;
end;
function cinv(const X: TASC): TASC; inline;
begin
Result := 1 / X;
end;
function NN(N: Integer): Integer; inline;
begin
if N >= 0 then
Result := N
else
Result := 0;
end;
function AltSgn(const N: Integer): Integer; inline;
begin
if N mod 2 = 0 then
Result := 1
else
Result := -1;
end;
function CreateIntSequence(AStart, AEnd: Integer): TArray<Integer>;
var
i: Integer;
begin
if AEnd >= AStart then
begin
SetLength(Result, AEnd - AStart + 1);
for i := 0 to High(Result) do
Result[i] := AStart + i
end
else
begin
SetLength(Result, AStart - AEnd + 1);
for i := 0 to High(Result) do
Result[i] := AStart - i;
end;
end;
function CreateIntSequence(AStart, AEnd, AStep: Integer): TArray<Integer>;
var
i: Integer;
begin
if AStep = 0 then
raise EMathException.Create('Cannot create an integer sequence with step size 0.');
AStep := Abs(AStep);
if AEnd >= AStart then
begin
SetLength(Result, 1 + (AEnd - AStart) div AStep);
for i := 0 to High(Result) do
Result[i] := AStart + i * AStep;
end
else
begin
SetLength(Result, 1 + (AStart - AEnd) div AStep);
for i := 0 to High(Result) do
Result[i] := AStart - i * AStep;
end;
end;
function CreateIntSequence64(AStart, AEnd: Int64): TArray<Int64>;
var
i: Int64;
begin
if AEnd >= AStart then
begin
SetLength(Result, AEnd - AStart + 1);
for i := 0 to High(Result) do
Result[i] := AStart + i
end
else
begin
SetLength(Result, AStart - AEnd + 1);
for i := 0 to High(Result) do
Result[i] := AStart - i;
end;
end;
function CreateIntSequence64(AStart, AEnd, AStep: Int64): TArray<Int64>;
var
i: Int64;
begin
if AStep = 0 then
raise EMathException.Create('Cannot create an integer sequence with step size 0.');
AStep := Abs(AStep);
if AEnd >= AStart then
begin
SetLength(Result, 1 + (AEnd - AStart) div AStep);
for i := 0 to High(Result) do
Result[i] := AStart + i * AStep;
end
else
begin
SetLength(Result, 1 + (AStart - AEnd) div AStep);
for i := 0 to High(Result) do
Result[i] := AStart - i * AStep;
end;
end;
procedure TranslateIntSequence(var ASeq: TArray<Integer>; const Offset: Integer = -1);
var
i: Integer;
begin
for i := Low(ASeq) to High(ASeq) do
Inc(ASeq[i], Offset);
end;
function TranslatedIntSequence(const ASeq: array of Integer;
const Offset: Integer = -1): TArray<Integer>;
var
i: Integer;
begin
SetLength(Result, Length(ASeq));
for i := 0 to High(ASeq) do
Result[i] := ASeq[i] + Offset;
end;
procedure TranslatePoint(var APoint: TPoint; const Offset: Integer = -1);
begin
APoint.X := APoint.X + Offset;
APoint.Y := APoint.Y + Offset;
end;
function TranslatedPoint(const APoint: TPoint; const Offset: Integer = -1): TPoint;
begin
Result.X := APoint.X + Offset;
Result.Y := APoint.Y + Offset;
end;
function ParseRangeSeq(const ARanges: array of TRange;
const ALength: Integer): TArray<Integer>;
var
List: TList<Integer>;
i, j: Integer;
a, b: Integer;
Step: Integer;
begin
List := TList<Integer>.Create;
try
for i := 0 to High(ARanges) do
begin
a := ARanges[i].From;
if a < 0 then
a := ALength + 1 + a;
b := ARanges[i].&To;
if b < 0 then
b := ALength + 1 + b;
Step := ARanges[i].Step;
case CompareValue(a, b) of
LessThanValue:
begin
j := max(a, 1);
while j <= min(b, ALength) do
begin
List.Add(j);
Inc(j, Step);
end;
end;
EqualsValue:
if InRange(a, 1, ALength) then
List.Add(a);
GreaterThanValue:
begin
j := min(a, ALength);
while j >= max(b, 1) do
begin
List.Add(j);
Dec(j, Step);
end;
end;
end;
end;
Result := List.ToArray;
finally
List.Free;
end;
end;
function IntArrToRealArr(const AArray: TArray<Integer>): TArray<TASR>;
var
i: Integer;
begin
SetLength(Result, Length(AArray));
for i := 0 to High(Result) do
Result[i] := AArray[i];
end;
function Int64ArrToRealArr(const AArray: TArray<Int64>): TArray<TASR>;
var
i: Integer;
begin
SetLength(Result, Length(AArray));
for i := 0 to High(Result) do
Result[i] := AArray[i];
end;
function IntArrToStrArr(const AArray: TArray<Integer>): TArray<string>;
var
i: Integer;
begin
SetLength(Result, Length(AArray));
for i := 0 to High(Result) do
Result[i] := AArray[i].ToString;
end;
function Int64ArrToStrArr(const AArray: TArray<Int64>): TArray<string>;
var
i: Integer;
begin
SetLength(Result, Length(AArray));
for i := 0 to High(Result) do
Result[i] := AArray[i].ToString;
end;
function ASR_COUNT(const A, B: TASR): TASR;
begin
Result := A + 1;
end;
function ASR_PLUS(const A, B: TASR): TASR;
begin
Result := A + B;
end;
function ASR_TIMES(const A, B: TASR): TASR;
begin
Result := A * B;
end;
function DoubleListToASR2s(const LList: TArray<Double>): TArray<TASR2>;
begin
SetLength(Result, Length(LList) div 2);
for var i := 0 to High(Result) do
begin
Result[i].X := LList[2*i ];
Result[i].Y := LList[2*i + 1];
end;
end;
function DoubleListToASR3s(const LList: TArray<Double>): TArray<TASR3>;
begin
SetLength(Result, Length(LList) div 3);
for var i := 0 to High(Result) do
begin
Result[i].X := LList[3*i ];
Result[i].Y := LList[3*i + 1];
Result[i].Z := LList[3*i + 2];
end;
end;
function TFormatStyleHelper.ToString: string;
begin
Result := FormatStyleNames[Self];
end;
class function TFormatStyleHelper.FromString(const S: string): TFormatStyle;
resourcestring
SUnknownNumberFormatName = 'Unknown number format "%s".';
var
i: TFormatStyle;
begin
for i := Low(TFormatStyle) to High(TFormatStyle) do
if i.ToString = S then
Exit(i);
raise Exception.CreateFmt(SUnknownNumberFormatName, [S]);
end;
constructor TRationalNumber.Create(const ANumerator, ADenominator: TASI);
begin
Numerator := ANumerator;
Denominator := ADenominator;
ToSimplestForm;
end;
class procedure TRationalNumber.ErrNoRep;
begin
raise EMathException.Create('Rational number cannot be represented.');
end;
function __DigChr(ADig: Integer): Char; inline;
begin
Result := Chr(Ord('0') + ADig);
end;
function TRationalNumber.ToDecimalString(
const AFormatOptions: TFormatOptions): string;
const
NegExpLimit = 4;
RDMs: array[TNumberFormat] of TRatDigitMode = (rdmFixedSignificant,
rdmFixedSignificant, rdmFractional, rdmSignificant, rdmSignificant);
FO: array[TNumberFormat] of Boolean = (False, False, True, True, False);
var
NF: TNumberFormat;
D: TArray<Integer>;
FirstDigitPos, FirstSigDigitPos: Integer;
NumDigits,
NumIntDigits,
NumFracDigits: Integer;
ResLen: Integer;
c: Integer;
i: Integer;
dgt: Char;
LNumDigits: Integer;
begin
if not valid then
Exit('NaN');
NF := AFormatOptions.Numbers.NumberFormat;
LNumDigits := Max(AFormatOptions.Numbers.NumDigits, IfThen(NF = nfFixed, 0, 1));
if ((Denominator = 1) or (Numerator = 0)) and (NF in [nfDefault, nfFraction]) then
Exit(IntegerToStr(Numerator, AFormatOptions));
try
D := GetDigits(Self, LNumDigits, RDMs[NF], FO[NF], True, @FirstDigitPos, @FirstSigDigitPos);
except
Exit(RealToStr(TASR(Self), AFormatOptions));
end;
NumDigits := Length(D);
NumIntDigits := IfThen(NF in [nfExponent, nfDefExp], 1, Max(FirstDigitPos + 1, 0));
NumFracDigits := NumDigits - NumIntDigits;
Assert(NumDigits >= 1);
Assert(NumIntDigits >= 0);
case NF of
nfDefault, nfFraction, nfFixed:
begin
if NF in [nfDefault, nfFraction] then
if (NumIntDigits > LNumDigits) or (FirstSigDigitPos < -NegExpLimit) then
Exit(ToDecimalString(FixNumFmt(AFormatOptions, nfDefExp)));
ResLen := NumDigits;
if (Numerator < 0) and (FirstSigDigitPos <> FirstSigDigitPos.MaxValue) then
Inc(ResLen, 1);
if (NumIntDigits > 0) and (AFormatOptions.Numbers.IntGrouping <> 0) then
Inc(ResLen, Ceil(NumIntDigits / AFormatOptions.Numbers.IntGrouping) - 1);
if NumFracDigits > 0 then
Inc(ResLen, 1);
if (NumFracDigits > 0) and (AFormatOptions.Numbers.FracGrouping <> 0) then
Inc(ResLen, Ceil(NumFracDigits / AFormatOptions.Numbers.FracGrouping) - 1);
SetLength(Result, ResLen);
c := 0;
if (Numerator < 0) and (FirstSigDigitPos <> FirstSigDigitPos.MaxValue) then
begin
Inc(c);
Result[c] := AFormatOptions.Numbers.MinusSign;
end;
for i := 0 to NumIntDigits - 1 do
begin
if
(AFormatOptions.Numbers.IntGrouping <> 0)
and
(i > 0)
and
((NumIntDigits - i) mod (AFormatOptions.Numbers.IntGrouping) = 0)
then
begin
Inc(c);
Result[c] := AFormatOptions.Numbers.IntGropingChar;
end;
Inc(c);
Result[c] := __DigChr(D[i]);
end;
if NumFracDigits > 0 then
begin
Inc(c);
Result[c] := AFormatOptions.Numbers.DecimalSeparator;
for i := NumIntDigits to NumDigits - 1 do
begin
if
(AFormatOptions.Numbers.FracGrouping <> 0)
and
(i > NumIntDigits)
and
((i - NumIntDigits) mod (AFormatOptions.Numbers.FracGrouping) = 0)
then
begin
Inc(c);
Result[c] := AFormatOptions.Numbers.FracGroupingChar;
end;
Inc(c);
Result[c] := __DigChr(D[i]);
end;
end;
Assert(c = ResLen);
end;
nfExponent, nfDefExp:
begin
ResLen := NumDigits;
if Numerator < 0 then
Inc(ResLen, 1);
if NumFracDigits >= 1 then
Inc(ResLen, 1);
if AFormatOptions.Numbers.PrettyExp then
Inc(ResLen, 4)
else
Inc(ResLen, 1);
if AFormatOptions.Numbers.PrettyExp or (NF = nfDefExp) then
Inc(ResLen, FirstDigitPos.ToString.Length)
else
Inc(ResLen, System.Abs(FirstDigitPos).ToString.Length + 1 );
if (NumFracDigits > 0) and (AFormatOptions.Numbers.FracGrouping <> 0) then
Inc(ResLen, Ceil(NumFracDigits / AFormatOptions.Numbers.FracGrouping) - 1);
SetLength(Result, ResLen);
c := 0;
if Numerator < 0 then
begin
Inc(c);
Result[c] := AFormatOptions.Numbers.MinusSign;
end;
Inc(c);
Result[c] := __DigChr(D[0]);
if NumFracDigits > 0 then
begin
Inc(c);
Result[c] := AFormatOptions.Numbers.DecimalSeparator;
end;
for i := 1 to High(D) do
begin
if
(AFormatOptions.Numbers.FracGrouping <> 0)
and
(i > 1)
and
(Pred(i) mod (AFormatOptions.Numbers.FracGrouping) = 0)
then
begin
Inc(c);
Result[c] := AFormatOptions.Numbers.FracGroupingChar;
end;
Inc(c);
Result[c] := __DigChr(D[i])
end;
if AFormatOptions.Numbers.PrettyExp then
begin
Result[c + 1] := DOT_OPERATOR;
Result[c + 2] := '1';
Result[c + 3] := '0';
Result[c + 4] := '^';
Inc(c, 4);
end
else
begin
Inc(c);
Result[c] := 'E';
end;
if (FirstDigitPos < 0) or not AFormatOptions.Numbers.PrettyExp and (NF = nfExponent) then
Inc(c);
if FirstDigitPos < 0 then
Result[c] := AFormatOptions.Numbers.MinusSign
else if not AFormatOptions.Numbers.PrettyExp and (NF = nfExponent) then
Result[c] := '+';
for dgt in System.Abs(FirstDigitPos).ToString do
begin
Inc(c);
Result[c] := dgt;
end;
Assert(c = ResLen);
end;
else
Result := RealToStr(TASR(Self), AFormatOptions);
end;
end;
procedure TRationalNumber.ToSimplestForm;
var
_gcd: TASI;
begin
if Denominator = 0 then
raise EMathException.Create('Rational number with zero denominator.');
if Numerator = 0 then
begin
Denominator := 1;
Exit;
end;
_gcd := gcd(Numerator, Denominator);
Numerator := Numerator div _gcd;
Denominator := Denominator div _gcd;
if Denominator < 0 then
if TInt64Guard.CanUnMin(Numerator) and TInt64Guard.CanUnMin(Denominator) then
begin
Numerator := -Numerator;
Denominator := -Denominator;
end
else
ErrNoRep;
end;
function TRationalNumber.ToString(const AFormatOptions: TFormatOptions;
const ASymbol: string): string;
begin
ToSimplestForm;
if Denominator = 1 then
begin
if Numerator = 1 then
Result := ASymbol
else if Numerator = -1 then
Result := AFormatOptions.Numbers.MinusSign + ASymbol
else
Result := IntegerToStr(Numerator, AFormatOptions) + DOT_OPERATOR + ASymbol
end
else if Numerator = 1 then
Result := ASymbol + '/' + IntegerToStr(Denominator, AFormatOptions)
else if Numerator = -1 then
Result := AFormatOptions.Numbers.MinusSign + ASymbol + '/' + IntegerToStr(Denominator, AFormatOptions)
else
Result := '(' + ToString(AFormatOptions) + ')' + DOT_OPERATOR + ASymbol;
end;
function TRationalNumber.ToString(const AFormatOptions: TFormatOptions): string;
begin
if AFormatOptions.Numbers.NumberFormat = nfFraction then
begin
if Denominator = 1 then
Result := IntegerToStr(Numerator, AFormatOptions)
else
Result := Format('%s/%s',
[IntegerToStr(Numerator, AFormatOptions), IntegerToStr(Denominator, AFormatOptions)])
end
else
Result := ToDecimalString(AFormatOptions)
end;
class operator TRationalNumber.Implicit(A: Integer): TRationalNumber;
begin
Result.Numerator := A;
Result.Denominator := 1;
end;
class operator TRationalNumber.Implicit(A: TASI): TRationalNumber;
begin
Result.Numerator := A;
Result.Denominator := 1;
end;
class operator TRationalNumber.Implicit(const X: TRationalNumber): TASR;
begin
Result := X.Numerator / X.Denominator;
end;
class operator TRationalNumber.Equal(const X, Y: TRationalNumber): Boolean;
var
tmp1, tmp2: TRationalNumber;
begin
tmp1 := X;
tmp2 := Y;
tmp1.ToSimplestForm;
tmp2.ToSimplestForm;
Result := (tmp1.Numerator = tmp2.Numerator) and (tmp1.Denominator = tmp2.Denominator);
end;
class operator TRationalNumber.NotEqual(const X, Y: TRationalNumber): Boolean;
begin
Result := not (X = Y);
end;
class function TRationalNumber.Power(const X: TRationalNumber;
const N: TASI): TRationalNumber;
var
c: Integer;
begin
if not X.valid then
Exit(InvalidRat);
if (X.Numerator = 0) and (N = 0) then
Exit(InvalidRat);
if InRange(N, -1000, 1000) then
begin
c := System.Abs(N);
Result := 1;
while (c > 0) and Result.valid do
begin
Result := Result * X;
Dec(c);
end;
if (N < 0) and Result.valid then
Result := Result.inv;
end
else
Result := InvalidRat;
end;
function TRationalNumber.Sign: TValueSign;
begin
Result := Math.Sign(Numerator) * Math.Sign(Denominator);
end;
class operator TRationalNumber.Negative(const X: TRationalNumber): TRationalNumber;
begin
if TInt64Guard.CanUnMin(X.Numerator) then
begin
Result.Numerator := -X.Numerator;
Result.Denominator := X.Denominator;
end
else
Result := InvalidRat;
end;
class operator TRationalNumber.Add(const X, Y: TRationalNumber): TRationalNumber;
var
_lcm, XFactor, YFactor, XPart, YPart: TASI;
begin
if TryLCM(X.Denominator, Y.Denominator, _lcm) then
begin
XFactor := _lcm div X.Denominator;
YFactor := _lcm div Y.Denominator;
if TInt64Guard.CanMul(X.Numerator, XFactor) and TInt64Guard.CanMul(Y.Numerator, YFactor) then
begin
XPart := X.Numerator * XFactor;
YPart := Y.Numerator * YFactor;
if TInt64Guard.CanAdd(XPart, YPart) then
begin
Result.Numerator := XPart + YPart;
Result.Denominator := _lcm;
Result.ToSimplestForm;
Exit;
end;
end;
end;
Result := InvalidRat;
end;
class operator TRationalNumber.Subtract(const X, Y: TRationalNumber): TRationalNumber;
var
_lcm, XFactor, YFactor, XPart, YPart: TASI;
begin
if TryLCM(X.Denominator, Y.Denominator, _lcm) then
begin
XFactor := _lcm div X.Denominator;
YFactor := _lcm div Y.Denominator;
if TInt64Guard.CanMul(X.Numerator, XFactor) and TInt64Guard.CanMul(Y.Numerator, YFactor) then
begin
XPart := X.Numerator * XFactor;
YPart := Y.Numerator * YFactor;
if TInt64Guard.CanSub(XPart, YPart) then
begin
Result.Numerator := XPart - YPart;
Result.Denominator := _lcm;
Result.ToSimplestForm;
Exit;
end;
end;
end;
Result := InvalidRat;
end;
class operator TRationalNumber.Multiply(const X, Y: TRationalNumber): TRationalNumber;
var
tmp1, tmp2: TRationalNumber;
begin
if TInt64Guard.CanMul(X.Numerator, Y.Numerator) and TInt64Guard.CanMul(X.Denominator, Y.Denominator) then
begin
Result.Numerator := X.Numerator * Y.Numerator;
Result.Denominator := X.Denominator * Y.Denominator;
Result.ToSimplestForm;
end
else
begin
tmp1 := RationalNumber(X.Numerator, Y.Denominator);
tmp2 := RationalNumber(Y.Numerator, X.Denominator);
if TInt64Guard.CanMul(tmp1.Numerator, tmp2.Numerator) and TInt64Guard.CanMul(tmp1.Denominator, tmp2.Denominator) then
begin
Result.Numerator := tmp1.Numerator * tmp2.Numerator;
Result.Denominator := tmp1.Denominator * tmp2.Denominator;
Result.ToSimplestForm;
end
else
Result := InvalidRat;
end
end;
class operator TRationalNumber.Divide(const X, Y: TRationalNumber): TRationalNumber;
begin
Result := X * Y.inv;
end;
function TRationalNumber.inv: TRationalNumber;
begin
if Numerator = 0 then
raise EMathException.Create('Division by zero.');
Result.Denominator := Numerator;
Result.Numerator := Denominator;
if Result.Denominator < 0 then
begin
if TInt64Guard.CanUnMin(Result.Numerator) and TInt64Guard.CanUnMin(Result.Denominator) then
begin
Result.Numerator := -Result.Numerator;
Result.Denominator := -Result.Denominator;
end
else
Result := InvalidRat;
end
end;
function TRationalNumber.sqr: TRationalNumber;
begin
Result := Self * Self;
end;
function TRationalNumber.str: string;
begin
if Denominator = 1 then
Result := Format('%d', [Numerator])
else
Result := Format('%d/%d', [Numerator, Denominator]);
end;
function TRationalNumber.valid: Boolean;
begin
Result := Denominator <> 0;
end;
function TRationalNumber.Abs: TRationalNumber;
begin
if TInt64Guard.CanAbs(Numerator) then
begin
Result.Numerator := System.Abs(Numerator);
Result.Denominator := Denominator;
end
else
Result := InvalidRat;
end;
constructor TSimpleSymbolicForm.Create(const AA, AB: TRationalNumber;
const ASym: string; ASymVal: TASR);
begin
A := AA;
B := AB;
Sym := ASym;
SymVal := ASymVal;
end;
constructor TSimpleSymbolicForm.Create(pA, qA, pB, qB: TASI; const ASym: string;
ASymVal: TASR);
begin
A := RationalNumber(pA, qA);
B := RationalNumber(pB, qB);
Sym := ASym;
SymVal := ASymVal;
end;
constructor TSimpleSymbolicForm.Create(const AA: TRationalNumber);
begin
A := AA;
B := 0;
Sym := '';
SymVal := 0;
end;
constructor TSimpleSymbolicForm.Create(pA, qA: TASI);
begin
A := RationalNumber(pA, qA);
B := 0;
Sym := '';
SymVal := 0;
end;
constructor TSimpleSymbolicForm.Create(const AB: TRationalNumber;
const ASym: string; ASymVal: TASR);
begin
A := 0;
B := AB;
Sym := ASym;
SymVal := ASymVal;
end;
constructor TSimpleSymbolicForm.Create(pB, qB: TASI; const ASym: string;
ASymVal: TASR);
begin
A := 0;
B := RationalNumber(pB, qB);
Sym := ASym;
SymVal := ASymVal;
end;
constructor TSimpleSymbolicForm.CreateInvalid(const Val: TASR);
begin
SymVal := Val;
MakeInvalid;
end;
class operator TSimpleSymbolicForm.Implicit(const X: TRationalNumber): TSimpleSymbolicForm;
begin
Result := TSimpleSymbolicForm.Create(X);
end;
class operator TSimpleSymbolicForm.Implicit(const X: TSimpleSymbolicForm): TASR;
begin
if X.valid then
Result := TASR(X.A) + TASR(X.B) * X.SymVal
else
Result := X.SymVal;
end;
class operator TSimpleSymbolicForm.Add(const X: TRationalNumber;
const Y: TSimpleSymbolicForm): TSimpleSymbolicForm;
begin
Result := Y;
Result.A := Result.A + X;
end;
class operator TSimpleSymbolicForm.Multiply(const X: TRationalNumber;
const Y: TSimpleSymbolicForm): TSimpleSymbolicForm;
begin
Result := TSimpleSymbolicForm.Create(X * Y.A, X * Y.B, Y.Sym, Y.SymVal);
end;
function TSimpleSymbolicForm.str: string;
begin
Result := Format('(%s, %s, %s = %g)', [A.str, B.str, Sym, SymVal]);
end;
function TSimpleSymbolicForm.ToString(
const AFormatOptions: TFormatOptions): string;
var
LFormatOptions: TFormatOptions;
begin
Result := '';
if not valid then
Exit;
LFormatOptions := FixNumFmt(AFormatOptions, nfFraction);
if (A = 0) and (B = 0) then
Result := IntegerToStr(0, LFormatOptions)
else if (A <> 0) and (B = 0) then
Result := A.ToString(LFormatOptions)
else if (A = 0) and (B <> 0) then
Result := B.ToString(LFormatOptions, Sym)
else
begin
if B.Numerator > 0 then
Result := A.ToString(LFormatOptions) + ' + ' + B.ToString(LFormatOptions, Sym)
else
Result := A.ToString(LFormatOptions) + #32 + MINUS_SIGN + #32 + (-B).ToString(LFormatOptions, Sym)
end;
end;
function TSimpleSymbolicForm.sstr: string;
begin
Result := Format('(%s, %s, %s)', [A.str, B.str, Sym]);
end;
function TSimpleSymbolicForm.valid: Boolean;
begin
Result := A.valid and B.valid;
end;
procedure TSimpleSymbolicForm.MakeInvalid;
begin
A.Denominator := 0;
end;
function FixNumFmt(const AOptions: TFormatOptions; ANumFmt: TNumberFormat): TFormatOptions;
begin
Result := AOptions;
Result.Numbers.NumberFormat := ANumFmt;
end;
function IntegerToStr(const x: TASI; const AOptions: TFormatOptions): string;
var
y, Base, Rem: UInt64;
i, c, pStart: Integer;
NumDigits: array[0..255] of Char;
const
Digits: array[0..35] of Char = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ';
begin
if (AOptions.Numbers.Base = 10) or (AOptions.Numbers.Base < 2) or (AOptions.Numbers.Base > 36) or not TInt64Guard.CanAbs(x) then
begin
Result := x.ToString;
if AOptions.Numbers.MinLength > 0 then
begin
if x < 0 then
begin
Assert((Result.Length > 0) and (Result[1] = '-'));
Delete(Result, 1, 1);
end;
Result := Result.PadLeft(AOptions.Numbers.MinLength, '0');
if x < 0 then
Result := AOptions.Numbers.MinusSign + Result;
end;
end
else
begin
if x = 0 then
Result := StringOfChar('0', Max(1, AOptions.Numbers.MinLength))
else
begin
Assert(x.Size = y.Size);
Assert(8*x.Size <= Length(NumDigits));
{$WARN COMPARISON_TRUE OFF}
Assert(AOptions.Numbers.MinLength.MaxValue < Length(NumDigits));
{$WARN COMPARISON_TRUE DEFAULT}
y := UInt64(Abs(x));
Base := AOptions.Numbers.Base;
FillChar(NumDigits, SizeOf(NumDigits), 0);
c := 0;
while y <> 0 do
begin
DivMod(y, Base, y, Rem);
Assert(InRange(Rem, Low(Digits), High(Digits)));
Assert(c <= High(NumDigits));
NumDigits[c] := Digits[Rem];
Inc(c);
end;
while c < AOptions.Numbers.MinLength do
begin
Assert(c <= High(NumDigits));
NumDigits[c] := '0';
Inc(c);
end;
if x < 0 then
begin
Assert(c <= High(NumDigits));
NumDigits[c] := AOptions.Numbers.MinusSign;
Inc(c);
end;
SetLength(Result, c);
for i := 1 to c do
Result[i] := NumDigits[c - i];
end;
end;
if (Result.Length > 0) and (Result[1] = '-') then
Result[1] := AOptions.Numbers.MinusSign;
if AOptions.Numbers.IntGrouping > 0 then
begin
i := Result.Length;
if (Result.Length > 0) and (Result[1] = AOptions.Numbers.MinusSign) then
pStart := 2
else
pStart := 1;
c := 1;
while i > pStart do
begin
if c mod AOptions.Numbers.IntGrouping = 0 then
Insert(AOptions.Numbers.IntGropingChar, Result, i);
Dec(i);
Inc(c);
end;
end;
end;
function RealToStr(const x: TASR; const AOptions: TFormatOptions): string;
var
i, c: Integer;
pPeriod, pExp, pStart, pEnd: Integer;
begin
if IsInfinite(x) then
case Sign(x) of
+1:
Exit('∞');
-1:
Exit('−∞');
end;
case AOptions.Numbers.NumberFormat of
nfFixed:
Result := FloatToStrF(x, ffFixed, 18, AOptions.Numbers.NumDigits, FS);
nfExponent:
Result := FloatToStrF(x, ffExponent, AOptions.Numbers.NumDigits, 0, FS);
else
Result := FloatToStrF(x, ffGeneral, AOptions.Numbers.NumDigits, 0, FS);
end;
for i := 1 to Result.Length - 1 do
if Result[i] = '-' then
Result[i] := AOptions.Numbers.MinusSign;
if AOptions.Numbers.PrettyExp then
for i := 1 to Result.Length - 1 do
if Result[i] = 'E' then
begin
Delete(Result, i, 1);
if Result[i] = '+' then
Delete(Result, i, 1);
Insert('⋅10^', Result, i);
Break;
end;
if AOptions.Numbers.IntGrouping + AOptions.Numbers.FracGrouping > 0 then
begin
pPeriod := 0;
for i := 1 to Result.Length do
if Result[i] = '.' then
begin
pPeriod := i;
Break;
end;
pExp := 0;
for i := 1 to Result.Length do
if (Result[i] = 'E') or (Result[i] = DOT_OPERATOR) then
begin
pExp := i;
Break;
end;
if AOptions.Numbers.IntGrouping > 0 then
begin
if pPeriod > 0 then
i := pPeriod - 1
else if pExp > 0 then
i := pExp - 1
else
i := Result.Length;
if (Result.Length > 0) and (Result[1] = AOptions.Numbers.MinusSign) then
pStart := 2
else
pStart := 1;
c := 1;
while i > pStart do
begin
if c mod AOptions.Numbers.IntGrouping = 0 then
begin
Insert(AOptions.Numbers.IntGropingChar, Result, i);
if pPeriod > 0 then
Inc(pPeriod);
if pExp > 0 then
Inc(pExp);
end;
Dec(i);
Inc(c);
end;
end;
if (AOptions.Numbers.FracGrouping > 0) and (pPeriod > 0) then
begin
if pExp > 0 then
pEnd := pExp - 1
else
pEnd := Result.Length;
i := pPeriod + 2;
c := 1;
while i <= pEnd do
begin
if c mod AOptions.Numbers.FracGrouping = 0 then
begin
Insert(AOptions.Numbers.FracGroupingChar, Result, i);
Inc(i);
Inc(pEnd);
end;
Inc(i);
Inc(c);
end;
end;
end;
end;
type
TStandardOrderASRComparer = class(TComparer<TASR>)
function Compare(const Left, Right: TASR): Integer; override;
end;
TDescendingStandardOrderASRComparer = class(TComparer<TASR>)
function Compare(const Left, Right: TASR): Integer; override;
end;
TAbsoluteValueASRComparer = class(TComparer<TASR>)
function Compare(const Left, Right: TASR): Integer; override;
end;
TDescendingAbsoluteValueASRComparer = class(TComparer<TASR>)
function Compare(const Left, Right: TASR): Integer; override;
end;
function TStandardOrderASRComparer.Compare(const Left, Right: TASR): Integer;
begin
Result := CompareValue(Left, Right);
end;
function TDescendingStandardOrderASRComparer.Compare(const Left, Right: TASR): Integer;
begin
Result := -CompareValue(Left, Right);
end;
function TAbsoluteValueASRComparer.Compare(const Left, Right: TASR): Integer;
begin
Result := CompareValue(Abs(Left), Abs(Right));
if Result = 0 then
Result := CompareValue(Left, Right);
end;
function TDescendingAbsoluteValueASRComparer.Compare(const Left, Right: TASR): Integer;
begin
Result := CompareValue(Abs(Left), Abs(Right));
if Result = 0 then
Result := CompareValue(Left, Right);
Result := -Result;
end;
class constructor TASRComparer.ClassCreate;
begin
TASRComparer.StandardOrder := TStandardOrderASRComparer.Create;
TASRComparer.StandardOrderDescending := TDescendingStandardOrderASRComparer.Create;
TASRComparer.AbsoluteValue := TAbsoluteValueASRComparer.Create;
TASRComparer.AbsoluteValueDescending := TDescendingAbsoluteValueASRComparer.Create;
end;
function pow(const a, b: TASR): TASR;
begin
if (a = 0) and (b = 0) then
raise EMathException.Create('Zero raised to the power of zero.');
Result := Math.Power(a, b);
end;
function intpow(const a, b: TASI): TASI;
var
i: TASI;
begin
if b < 0 then
raise EMathException.Create('Integer power requires a non-negative exponent.');
Result := 1;
for i := 1 to b do
Result := Result * a;
end;
function cbrt(const X: TASR): TASR; inline;
begin
Result := Math.Power(X, 1/3);
end;
function frt(const X: TASR): TASR; inline;
begin
Result := Math.Power(X, 1/4);
end;
function tanh(const X: TASR): TASR; inline;
var
yp, ym: TASR;
begin
yp := exp(x);
ym := exp(-x);
Result := (yp - ym) / (yp + ym);
end;
function coth(const X: TASR): TASR; inline;
var
yp, ym: TASR;
begin
yp := exp(x);
ym := exp(-x);
Result := (yp + ym) / (yp - ym);
end;
function sech(const X: TASR): TASR; inline;
begin
{$IF sizeof(TASR) = 10}
if Abs(X) > 11300 then
Exit(0.0);
{$ELSE}
if Abs(X) > 700 then
Exit(0.0);
{$ENDIF}
Result := 2 / (exp(X) + exp(-X));
end;
function csch(const X: TASR): TASR; inline;
begin
{$IF sizeof(TASR) = 10}
if Abs(X) > 11300 then
Exit(0.0);
{$ELSE}
if Abs(X) > 700 then
Exit(0.0);
{$ENDIF}
Result := 2 / (exp(X) - exp(-X));
end;
function arcsinh(const X: TASR): TASR; inline;
begin
if X >= 0 then
Result := Math.ArcSinh(X)
else
Result := -Math.ArcSinh(-X);
end;
function arccsch(const X: TASR): TASR; inline;
var
y: TASR;
begin
if X >= 0 then
begin
y := 1/x;
Result := ln(y + sqrt(y*y + 1));
end
else
begin
y := 1/-x;
Result := -ln(y + sqrt(y*y + 1));
end;
end;
function sinc(const X: TASR): TASR; inline;
begin
if IsZero(X) then
Result := 1
else
Result := sin(x) / x;
end;
class operator TASC.Implicit(const r: TASR): TASC;
begin
Result.Re := r;
Result.Im := 0;
end;
class operator TASC.Positive(const z: TASC): TASC;
begin
Result := z;
end;
function TASC.pstr: string;
begin
Result := ComplexToStr(Self, False, DefaultFormatOptions);
end;
class operator TASC.Negative(const z: TASC): TASC;
begin
Result.Re := -z.Re;
Result.Im := -z.Im;
end;
class operator TASC.Add(const z1: TASC; const z2: TASC): TASC;
begin
Result.Re := z1.Re + z2.Re;
Result.Im := z1.Im + z2.Im;
end;
class operator TASC.Subtract(const z1: TASC; const z2: TASC): TASC;
begin
Result.Re := z1.Re - z2.Re;
Result.Im := z1.Im - z2.Im;
end;
class operator TASC.Multiply(const z1: TASC; const z2: TASC): TASC;
begin
Result.Re := z1.Re * z2.Re - z1.Im * z2.Im;
Result.Im := z1.Re * z2.Im + z1.Im * z2.Re;
end;
class operator TASC.Divide(const z1: TASC; const z2: TASC): TASC;
var
denom: TASR;
begin
denom := z2.Re * z2.Re + z2.Im * z2.Im;
Result.Re := (z1.Re * z2.Re + z1.Im * z2.Im) / denom;
Result.Im := (z1.Im * z2.Re - z1.Re * z2.Im) / denom;
end;
class operator TASC.Equal(const z1: TASC; const z2: TASC): Boolean;
begin
Result := (z1.Re = z2.Re) and (z1.Im = z2.Im);
end;
class operator TASC.NotEqual(const z1: TASC; const z2: TASC): Boolean;
begin
Result := not (z1 = z2);
end;
class operator TASC.Round(const z: TASC): TASC;
begin
Result.Re := Round(z.Re);
Result.Im := Round(z.Im);
end;
class operator TASC.Trunc(const z: TASC): TASC;
begin
Result.Re := Trunc(z.Re);
Result.Im := Trunc(z.Im);
end;
function TASC.Modulus: TASR;
begin
Result := Hypot(Re, Im);
end;
function TASC.Argument: TASR;
begin
Result := ArcTan2(Im, Re);
end;
function TASC.Conjugate: TASC;
begin
Result.Re := Re;
Result.Im := -Im;
end;
function TASC.Sqr: TASC;
begin
Result.Re := System.Sqr(Re) - System.Sqr(Im);
Result.Im := 2*Re*Im;
end;
function TASC.ModSqr: TASR;
begin
Result := System.Sqr(Re) + System.Sqr(Im);
end;
function TASC.Inverse: TASC;
var
h: TASR;
begin
h := System.Sqr(Re) + System.Sqr(Im);
Result.Re := Re/h;
Result.Im := -Im/h;
end;
function TASC.IsReal: Boolean;
begin
Result := IsZero(Im);
end;
function TASC.Defuzz(const Eps: Double): TASC;
begin
if IsInteger(Re, Eps) then
Result.Re := Round(Re)
else
Result.Re := Re;
if IsInteger(Im, Eps) then
Result.Im := Round(Im)
else
Result.Im := Im;
end;
function TASC.str: string;
begin
Result := Format('(%.18g, %.18g)', [Re, Im]);
end;
function ASC(const Re, Im: TASR): TASC; inline;
begin
Result.Re := Re;
Result.Im := Im;
end;
function ASC_COUNT(const A, B: TASC): TASC;
begin
Result := A + 1;
end;
function ASC_PLUS(const A, B: TASC): TASC;
begin
Result := A + B;
end;
function ASC_TIMES(const A, B: TASC): TASC;
begin
Result := A * B;
end;
function ComplexToStr(const z: TASC; ApproxEq: Boolean;
const AOptions: TFormatOptions): string;
var
Re0, Im0, Im1, ImN1: Boolean;
begin
if ApproxEq then
begin
Re0 := IsZero(z.Re);
Im0 := IsZero(z.Im);
Im1 := IsZero(z.Im - 1);
ImN1 := IsZero(z.Im + 1);
end
else
begin
Re0 := z.Re = 0;
Im0 := z.Im = 0;
Im1 := z.Im = 1;
ImN1 := z.Im = -1;
end;
if Re0 and Im0 then
Result := RealToStr(0, AOptions)
else if Im0 then
Result := RealToStr(z.Re, AOptions)
else if Re0 then
if Im1 then
Result := 'i'
else if ImN1 then
Result := AOptions.Complex.MinusSign + 'i'
else
Result := RealToStr(z.Im, AOptions) + AOptions.Complex.ImaginarySuffix
else
if Im1 then
Result := RealToStr(z.Re, AOptions) + AOptions.Complex.PlusStr + 'i'
else if ImN1 then
Result := RealToStr(z.Re, AOptions) + AOptions.Complex.MinusStr + 'i'
else if z.Im > 0 then
Result := RealToStr(z.Re, AOptions) + AOptions.Complex.PlusStr + RealToStr(z.Im, AOptions) + AOptions.Complex.ImaginarySuffix
else
Result := RealToStr(z.Re, AOptions) + AOptions.Complex.MinusStr + RealToStr(-z.Im, AOptions) + AOptions.Complex.ImaginarySuffix;
end;
function TryStringToComplex(AString: string; out Value: TASC;
ASignRequired: Boolean = False): Boolean;
var
num: TASR;
part1, part2: TASC;
len: Integer;
i, p, p2, pp: Integer;
signed, period, exp, neg, times, imag: Boolean;
procedure SkipWhitespace;
begin
while (i <= len) and AString[i].IsWhiteSpace do Inc(i);
end;
procedure SkipTrailingWhitespace;
var
j: Integer;
begin
j := i;
while (j <= len) and AString[j].IsWhiteSpace do Inc(j);
if j > len then
begin
p2 := i;
i := j;
end;
end;
procedure SkipDigits;
begin
while (i <= len) and AString[i].IsDigit do Inc(i);
end;
function EOF: Boolean;
begin
Result := i > len;
end;
procedure SkipSign;
begin
if (i <= len) and inOpArray(AString[i], ['+', '-']) then Inc(i);
end;
procedure SkipSignEx;
begin
if i <= len then
if AString[i] = '+' then
begin
signed := True;
Inc(i);
end
else if AString[i] = '-' then
begin
signed := True;
neg := not neg;
Inc(i);
end;
end;
procedure SkipPeriod;
begin
if (i <= len) and (AString[i] = '.') and not period then
begin
period := True;
Inc(i);
end;
end;
function sign: TASR;
begin
Result := IfThen(neg, -1, 1);
end;
function basis: TASC;
begin
if imag then
Result := ImaginaryUnit
else
Result := 1;
end;
function TryDone: Boolean;
var
res: Boolean;
begin
if i > len then
begin
res := TryStrToFloat(Copy(AString, p, p2 - p), num, FS);
TryStringToComplex := res;
if res then Value := sign * num;
Exit(True);
end
else
Result := False;
end;
procedure SkipExp;
begin
if (i <= len) and (AString[i] = 'e') and not exp then
begin
exp := True;
Inc(i);
end;
end;
procedure SkipTimes;
begin
if (i <= len) and (AString[i] = '*') then
begin
times := True;
Inc(i);
end;
end;
procedure SkipImUnit;
begin
if (i <= len) and (AString[i] = 'i') and not imag then
begin
imag := True;
Inc(i);
end;
end;
begin
len := Length(AString);
signed := False;
period := False;
exp := False;
times := False;
imag := False;
neg := False;
p2 := 128;
for i := 1 to len do
case AString[i] of
MINUS_SIGN:
AString[i] := '-';
'I', 'j', 'J':
AString[i] := 'i';
DOT_OPERATOR:
AString[i] := '*';
'E':
AString[i] := 'e';
end;
i := 1;
SkipWhitespace;
if EOF then Exit(False);
SkipSignEx;
if ASignRequired and not signed then Exit(False);
SkipWhitespace;
p := i;
SkipPeriod;
if TryDone then Exit;
SkipDigits;
SkipTrailingWhitespace;
if TryDone then Exit;
SkipPeriod;
SkipTrailingWhitespace;
if TryDone then Exit;
SkipDigits;
SkipTrailingWhitespace;
if TryDone then Exit;
SkipExp;
if EOF then Exit(False);
if exp then SkipSign;
if EOF then Exit(False);
pp := i;
SkipDigits;
if exp and (i = pp) then Exit(False);
p2 := i;
SkipWhitespace;
if TryDone then Exit;
SkipTimes;
SkipWhitespace;
if EOF then Exit(False);
SkipImUnit;
SkipWhitespace;
if times and not imag then Exit(False);
if i > len then
begin
Assert(imag);
num := 1;
Result := (p = p2) or TryStrToFloat(Copy(AString, p, p2 - p), num, FS);
if Result then Value := sign * num * ImaginaryUnit;
Exit;
end;
num := 1;
if (imag and (p = p2)) or TryStrToFloat(Copy(AString, p, p2 - p), num, FS) then
begin
part1 := sign * num * basis;
if i = 1 then Exit(False);
if TryStringToComplex(Copy(AString, i), part2, True) then
begin
Value := part1 + part2;
Exit(True);
end;
end;
Result := False;
end;
function CSameValue(const z1, z2: TASC; const Epsilon: TASR = 0): Boolean; inline;
begin
Result := SameValue(z1.Re, z2.Re, Epsilon) and SameValue(z1.Im, z2.Im, Epsilon);
end;
function CSameValueEx(const z1, z2: TASC; const Epsilon: TASR = 0): Boolean; inline;
begin
Result := SameValueEx(z1.Re, z2.Re, Epsilon) and SameValueEx(z1.Im, z2.Im, Epsilon);
end;
function CIsZero(const z: TASC; const Epsilon: TASR): Boolean; overload; inline;
begin
Result := IsZero(z.Re, Epsilon) and IsZero(z.Im, Epsilon);
end;
function IntegerPowerOfImaginaryUnit(const n: Integer): TASC;
const
vals: array[0..3] of TASC = ((Re: 1; Im: 0), (Re: 0; Im: 1), (Re: -1; Im: 0), (Re: 0; Im: -1));
begin
Result := vals[imod(n, Length(vals))];
end;
function SameValue2(const A, B: TASC): Boolean; overload;
begin
Result := SameValue2(A.Re, B.Re) and SameValue2(A.Im, B.Im);
end;
function CompareValue(const A, B: TASC; Epsilon: Extended): TValueRelationship;
begin
if CSameValue(A, B, Epsilon) then
Result := EqualsValue
else if (A.Re < B.Re) or ((A.Re = B.Re) and (A.Im < B.Im)) then
Result := LessThanValue
else
Result := GreaterThanValue;
end;
type
TReImASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
TDescendingReImASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
TModulusASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
TDescendingModulusASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
TArgASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
TDescendingArgASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
TModulusArgumentASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
TDescendingModulusArgumentASCComparer = class(TComparer<TASC>)
function Compare(const Left, Right: TASC): Integer; override;
end;
function TReImASCComparer.Compare(const Left, Right: TASC): Integer;
begin
Result := CompareValue(Left.Re, Right.Re);
if Result = 0 then
Result := CompareValue(Left.Im, Right.Im);
end;
function TDescendingReImASCComparer.Compare(const Left, Right: TASC): Integer;
begin
Result := CompareValue(Left.Re, Right.Re);
if Result = 0 then
Result := CompareValue(Left.Im, Right.Im);
Result := -Result;
end;
function TModulusASCComparer.Compare(const Left, Right: TASC): Integer;
begin
Result := CompareValue(Left.Modulus, Right.Modulus);
if Result = 0 then
Result := CompareValue(Left.Re, Right.Re);
if Result = 0 then
Result := CompareValue(Left.Im, Right.Im);
end;
function TDescendingModulusASCComparer.Compare(const Left, Right: TASC): Integer;
begin
Result := CompareValue(Left.Modulus, Right.Modulus);
if Result = 0 then
Result := CompareValue(Left.Re, Right.Re);
if Result = 0 then
Result := CompareValue(Left.Im, Right.Im);
Result := -Result;
end;
function TArgASCComparer.Compare(const Left, Right: TASC): Integer;
var
LA, RA: TASR;
begin
LA := Left.Argument;
if SameValue(LA, -Pi) then
LA := Pi;
RA := Right.Argument;
if SameValue(RA, -Pi) then
RA := Pi;
Result := CompareValue(LA, RA);
if Result = 0 then
Result := CompareValue(Left.Modulus, Right.Modulus);
end;
function TDescendingArgASCComparer.Compare(const Left, Right: TASC): Integer;
var
LA, RA: TASR;
begin
LA := Left.Argument;
if SameValue(LA, -Pi) then
LA := Pi;
RA := Right.Argument;
if SameValue(RA, -Pi) then
RA := Pi;
Result := CompareValue(LA, RA);
if Result = 0 then
Result := CompareValue(Left.Modulus, Right.Modulus);
Result := -Result;
end;
function TModulusArgumentASCComparer.Compare(const Left, Right: TASC): Integer;
var
LA, RA: TASR;
begin
Result := CompareValue(Left.Modulus, Right.Modulus);
if Result = 0 then
begin
LA := Left.Argument;
if SameValue(LA, -Pi) then
LA := Pi;
RA := Right.Argument;
if SameValue(RA, -Pi) then
RA := Pi;
Result := CompareValue(LA, RA);
end;
end;
function TDescendingModulusArgumentASCComparer.Compare(const Left, Right: TASC): Integer;
var
LA, RA: TASR;
begin
Result := CompareValue(Left.Modulus, Right.Modulus);
if Result = 0 then
begin
LA := Left.Argument;
if SameValue(LA, -Pi) then
LA := Pi;
RA := Right.Argument;
if SameValue(RA, -Pi) then
RA := Pi;
Result := CompareValue(LA, RA);
end;
Result := -Result;
end;
class constructor TASCComparer.ClassCreate;
begin
TASCComparer.ReIm := TReImASCComparer.Create;
TASCComparer.ReImDescending := TDescendingReImASCComparer.Create;
TASCComparer.Modulus := TModulusASCComparer.Create;
TASCComparer.ModulusDescending := TDescendingModulusASCComparer.Create;
TASCComparer.Argument := TArgASCComparer.Create;
TASCComparer.ArgumentDescending := TDescendingArgASCComparer.Create;
TASCComparer.ModulusArgument := TModulusArgumentASCComparer.Create;
TASCComparer.ModulusArgumentDescending := TDescendingModulusArgumentASCComparer.Create;
end;
function csign(const z: TASC): TASC; inline;
begin
if CIsZero(z) then
Result := 0
else
Result := z / z.Modulus;
end;
function cexp(const z: TASC): TASC; inline;
var
e, s, c: TASR;
begin
e := exp(z.Re);
SinCos(z.Im, s, c);
Result.Re := e * c;
Result.Im := e * s;
end;
function cln(const z: TASC): TASC; inline;
begin
Result.Re := ln(z.Modulus);
Result.Im := z.Argument;
end;
function clog(const z: TASC): TASC; inline;
begin
Result := cln(z) * InvLn10;
end;
function cpow(const z, w: TASC): TASC;
var
n: Integer;
i: Integer;
begin
if z = 0 then
if w.Re > 0 then
Exit(0)
else
raise EMathException.Create('Zero raised to a complex number with non-positive real part.');
if IsZero(w.Im) and IsZero(frac(w.Re)) and (abs(w.Re) < 13) then
begin
n := round(w.Re);
if n > 0 then
begin
Result := z;
for i := 2 to n do
Result := Result * z;
Exit;
end
else if n < 0 then
begin
Result := z;
for i := 2 to -n do
Result := Result * z;
Exit(1/Result);
end
else if (n = 0) and (z <> 0) then
Exit(1)
else
raise EMathException.Create('Zero raised to the power of zero.');
end
else if ((z.Im = 0) and (w.Im = 0)) and
((z.Re > 0) or ((z.Re = 0) and (w.Re > 0)) or ((z.Re < 0) and (frac(w.Re) = 0))) then
Exit(pow(z.Re, w.Re))
else if (z.Re = 0) and (w.Im = 0) and (frac(w.Re) = 0) then
Exit(IntPower(z.Im, round(w.Re)) * IntegerPowerOfImaginaryUnit(round(w.Re)));
Result := cexp(w * cln(z));
end;
function csqrt(const z: TASC): TASC; inline;
begin
if (z.Im = 0) and (z.Re >= 0) then
Exit(sqrt(z.Re));
if (z.Im = 0) and (z.Re < 0) then
Exit(ASC(0, sqrt(-z.Re)));
Result := cexp(cln(z)/2);
end;
function csin(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(sin(z.Re));
{$ENDIF}
Result := (cexp(ImaginaryUnit*z) - cexp(NegativeImaginaryUnit*z)) / (2*ImaginaryUnit);
end;
function ccos(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(cos(z.Re));
{$ENDIF}
Result := (cexp(ImaginaryUnit*z) + cexp(NegativeImaginaryUnit*z)) / 2;
end;
function ctan(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(tan(z.Re));
{$ENDIF}
w := cexp(2*ImaginaryUnit*z);
Result := NegativeImaginaryUnit * (w - 1)/(w + 1);
end;
function ccot(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(cot(z.Re));
{$ENDIF}
w := cexp(2*ImaginaryUnit*z);
Result := ImaginaryUnit * (w + 1)/(w - 1);
end;
function csec(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(sec(z.Re));
{$ENDIF}
Result := 2 / (cexp(ImaginaryUnit*z) + cexp(NegativeImaginaryUnit*z));
end;
function ccsc(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(csc(z.Re));
{$ENDIF}
Result := (2*ImaginaryUnit) / (cexp(ImaginaryUnit*z) - cexp(NegativeImaginaryUnit*z));
end;
function carcsin(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and InRange(z.Re, -1, 1) then
Exit(arcsin(z.Re));
{$ENDIF}
Result := NegativeImaginaryUnit * cln(ImaginaryUnit * z + csqrt(1 - z.Sqr));
end;
function carccos(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and InRange(z.Re, -1, 1) then
Exit(arccos(z.Re));
{$ENDIF}
Result := PiDiv2 + ImaginaryUnit * cln(ImaginaryUnit * z + csqrt(1 - z.Sqr));
end;
function carctan(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(arctan(z.Re));
{$ENDIF}
w := ImaginaryUnit * z;
Result := ImaginaryUnitDiv2 * (cln(1 - w) - cln(1 + w));
end;
function carccot(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(arccot(z.Re));
{$ENDIF}
if CIsZero(z) then
Result := PiDiv2
else
begin
w := ImaginaryUnit / z;
Result := ImaginaryUnitDiv2 * (cln(1 - w) - cln(1 + w))
end;
end;
function carcsec(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and (abs(z.Re) >= 1) then
Exit(arcsec(z.Re));
{$ENDIF}
Result := PiDiv2 + ImaginaryUnit * cln(ImaginaryUnit / z + csqrt(1 - z.Inverse.Sqr));
end;
function carccsc(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and (abs(z.Re) >= 1) then
Exit(arccsc(z.Re));
{$ENDIF}
Result := NegativeImaginaryUnit * cln(ImaginaryUnit / z + csqrt(1 - z.Inverse.Sqr));
end;
function csinh(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(sinh(z.Re));
{$ENDIF}
Result := (cexp(z) - cexp(-z)) / 2;
end;
function ccosh(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(cosh(z.Re));
{$ENDIF}
Result := (cexp(z) + cexp(-z)) / 2;
end;
function ctanh(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(tanh(z.Re));
{$ENDIF}
w := cexp(2*z);
Result := (w - 1)/(w + 1);
end;
function ccoth(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(coth(z.Re));
{$ENDIF}
w := cexp(2*z);
Result := (w + 1)/(w - 1);
end;
function csech(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(sech(z.Re));
{$ENDIF}
Result := 2 / (cexp(z) + cexp(-z));
end;
function ccsch(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(csch(z.Re));
{$ENDIF}
Result := 2 / (cexp(z) - cexp(-z));
end;
function carcsinh(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(arcsinh(z.Re));
{$ENDIF}
Result := cln(z + csqrt(1 + z.Sqr));
end;
function carccosh(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and (z.Re >= 1) then
Exit(arccosh(z.Re));
{$ENDIF}
Result := cln(z + csqrt(z + 1) * csqrt(z - 1));
end;
function carctanh(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and (z.Re > -1) and (z.Re < 1) then
Exit(arctanh(z.Re));
{$ENDIF}
Result := (cln(1 + z) - cln(1 - z)) / 2;
end;
function carccoth(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and ((z.Re < -1) or (z.Re > 1)) then
Exit(arccoth(z.Re));
{$ENDIF}
if CIsZero(z) then
Result := PiDiv2 * ImaginaryUnit
else
begin
w := z.Inverse;
Result := (cln(1 + w) - cln(1 - w)) / 2;
end;
end;
function carcsech(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if (z.Im = 0) and (z.Re > 0) and (z.Re <= 1) then
Exit(arcsech(z.Re));
{$ENDIF}
w := z.Inverse;
Result := cln(w + csqrt(w + 1) * csqrt(w - 1));
end;
function carccsch(const z: TASC): TASC; inline;
var
w: TASC;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(arccsch(z.Re));
{$ENDIF}
w := z.Inverse;
Result := cln(w + csqrt(1 + w.Sqr));
end;
function csinc(const z: TASC): TASC; inline;
begin
{$IFDEF REALCHECK}
if z.Im = 0 then
Exit(sinc(z.Re));
{$ENDIF}
if CIsZero(z) then
Result := 1
else
Result := csin(z) / z;
end;
function CollapseWithMultiplicity(const ANumbers: array of TASI): TArray<TIntWithMultiplicity>;
var
ActualLength: Integer;
i, c: Integer;
PrevNumber: TASI;
procedure New(ANumber: TASI);
begin
PrevNumber := ANumber;
c := 1;
end;
procedure Add;
begin
Result[ActualLength].Factor := PrevNumber;
Result[ActualLength].Multiplicity := c;
Inc(ActualLength);
end;
begin
SetLength(Result, Length(ANumbers));
ActualLength := 0;
for i := 0 to High(ANumbers) do
begin
if i = 0 then
New(ANumbers[i])
else if ANumbers[i] <> PrevNumber then
begin
Add;
New(ANumbers[i]);
end
else
Inc(c);
end;
if Length(ANumbers) > 0 then
Add;
SetLength(Result, ActualLength);
end;
function GetDigit(N: TASI; Index: Integer): Integer;
begin
if Index < 0 then
Exit(0);
while (Index > 0) and (N <> 0) do
begin
N := N div 10;
Dec(Index);
end;
Result := Abs(N mod 10);
end;
function GetDigits(N: TASI): TArray<Integer>;
begin
if N = 0 then
Result := [0]
else
while N <> 0 do
begin
TArrBuilder<Integer>.Add(Result, Abs(N mod 10));
N := N div 10;
end;
end;
procedure _RatDigits(const R: TRationalNumber; Index: Integer;
const Mode: TRatDigitMode; AFullOutput, ARound: Boolean; out Digit: Integer;
out Digits: TArray<Integer>; AFirstDigitPos, AFirstSigDigitPos: PInteger);
function makenum(source: TArray<Integer>; start, len: Integer): TASI;
var
i: Integer;
j: TASI;
begin
Result := 0;
j := 1;
for i := start to start + len - 1 do
begin
if i >= 0 then
Inc(Result, source[i] * j);
j := 10 * j;
end;
end;
var
LFirstDigitPos: Integer;
Numerator, Denominator: TArray<Integer>;
Pos: Integer;
CurNumerator, Quotient, Remainder: TASI;
Sig: Integer;
Rounding: Boolean;
EffectRound: Boolean;
i: Integer;
begin
Rounding := False;
EffectRound := False;
Numerator := GetDigits(R.Numerator);
Denominator := GetDigits(R.Denominator);
if Length(Denominator) >= 18 then
raise EMathException.Create('Denominator too large.');
Pos := Length(Numerator) - Length(Denominator);
CurNumerator := makenum(Numerator, Pos, Length(Denominator));
LFirstDigitPos := LFirstDigitPos.MaxValue;
if Mode <> rdmSingle then
begin
if Assigned(AFirstDigitPos) then
AFirstDigitPos^ := AFirstDigitPos^.MaxValue;
if Assigned(AFirstSigDigitPos) then
AFirstSigDigitPos^ := AFirstSigDigitPos^.MaxValue;
end;
case Mode of
rdmSingle:
begin
Digit := 0;
if (R.Numerator = 0) or (Index > Pos) then
Exit;
end;
rdmFractional:
begin
if R.Numerator = 0 then
begin
SetLength(Digits, Index + 1);
if Assigned(AFirstDigitPos) then
AFirstDigitPos^ := 0;
Exit;
end;
if Pos < 0 then
begin
SetLength(Digits, Min(-Pos, 1 + Index));
if Assigned(AFirstDigitPos) then
AFirstDigitPos^ := 0;
LFirstDigitPos := 0;
if Length(Digits) - 1 >= Index then
begin
if ARound and (-Pos = Index + 1) then
Rounding := True
else
Exit;
end;
end;
end;
rdmSignificant:
begin
if R.Numerator = 0 then
begin
SetLength(Digits, Index);
if Assigned(AFirstDigitPos) then
AFirstDigitPos^ := 0;
Exit;
end;
end;
rdmFixedSignificant:
begin
if R.Numerator = 0 then
begin
SetLength(Digits, Index);
if Assigned(AFirstDigitPos) then
AFirstDigitPos^ := 0;
Exit;
end;
if Pos < 0 then
begin
SetLength(Digits, -Pos);
if Assigned(AFirstDigitPos) then
AFirstDigitPos^ := 0;
LFirstDigitPos := 0;
end;
end;
else
raise Exception.CreateFmt('Unsupported rational digit extraction mode: %d', [Ord(Mode)]);
end;
Sig := 0;
while True do
begin
Quotient := CurNumerator div R.Denominator;
Remainder := CurNumerator mod R.Denominator;
case Mode of
rdmSingle:
begin
if Pos = Index then
begin
Digit := Quotient;
Exit;
end;
if (Pos < 0) and (Remainder = 0) then
Exit;
end;
rdmFractional:
begin
if Rounding then
begin
if Quotient >= 5 then
begin
EffectRound := True;
if Assigned(AFirstSigDigitPos) and (AFirstSigDigitPos^ = AFirstSigDigitPos^.MaxValue) then
AFirstSigDigitPos^ := Pos;
end;
Break;
end;
if (Pos <= 0) or (Length(Digits) > 0) or (Quotient <> 0) then
begin
if Assigned(AFirstDigitPos) then
if AFirstDigitPos^ = AFirstDigitPos^.MaxValue then
AFirstDigitPos^ := Pos;
if LFirstDigitPos = LFirstDigitPos.MaxValue then
LFirstDigitPos := Pos;
TArrBuilder<Integer>.Add(Digits, Quotient);
if (Quotient <> 0) or (Sig > 0) then
begin
if (Sig = 0) and Assigned(AFirstSigDigitPos) then
AFirstSigDigitPos^ := Pos;
Inc(Sig);
end;
end;
if not AFullOutput and (Pos <= 0) and (Remainder = 0) then
Break;
if Pos <= -Index then
begin
if ARound then
Rounding := True
else
Break;
end;
end;
rdmSignificant, rdmFixedSignificant:
begin
if Rounding then
begin
if Quotient >= 5 then
EffectRound := True;
Break;
end;
if (Length(Digits) > 0) or (Quotient <> 0) or (Pos <= 0) and (Mode = rdmFixedSignificant) then
begin
if Assigned(AFirstDigitPos) then
if AFirstDigitPos^ = AFirstDigitPos^.MaxValue then
AFirstDigitPos^ := Pos;
if LFirstDigitPos = LFirstDigitPos.MaxValue then
LFirstDigitPos := Pos;
TArrBuilder<Integer>.Add(Digits, Quotient);
if (Quotient <> 0) or (Sig > 0) then
begin
if (Sig = 0) and Assigned(AFirstSigDigitPos) then
AFirstSigDigitPos^ := Pos;
Inc(Sig);
end;
end;
if not AFullOutput and (Pos <= 0) and (Remainder = 0) then
Break;
if Sig >= Index then
begin
if ARound then
Rounding := True
else
Break;
end;
end;
else
raise Exception.CreateFmt('Unsupported rational digit extraction mode: %d', [Ord(Mode)]);
end;
CurNumerator := 10 * Remainder;
if Pos > 0 then
Inc(CurNumerator, Numerator[Pos - 1]);
Dec(Pos);
end;
if EffectRound then
begin
for i := High(Digits) downto 0 do
if Digits[i] <> 9 then
begin
Inc(Digits[i]);
Break;
end
else
begin
Digits[i] := 0;
if i = 0 then
begin
Insert([1], Digits, 0);
if Assigned(AFirstDigitPos) and (AFirstDigitPos^ <> AFirstDigitPos^.MaxValue) then
Inc(AFirstDigitPos^);
if LFirstDigitPos <> LFirstDigitPos.MaxValue then
Inc(LFirstDigitPos);
if Assigned(AFirstSigDigitPos) and (AFirstSigDigitPos^ <> AFirstSigDigitPos^.MaxValue) then
Inc(AFirstSigDigitPos^);
if Mode in [rdmSignificant, rdmFixedSignificant] then
SetLength(Digits, Pred(Length(Digits)));
Break;
end;
end;
end;
if not AFullOutput then
begin
i := High(Digits);
while (i > 0) and (Digits[i] = 0) and ((i > LFirstDigitPos) or (Mode <> rdmFixedSignificant)) do
Dec(i);
SetLength(Digits, i + 1);
end;
end;
function GetDigit(const R: TRationalNumber; Index: Integer): Integer; overload;
var
dummy: TArray<Integer>;
begin
_RatDigits(R, Index, rdmSingle, False, False, Result, dummy, nil, nil);
end;
function GetDigits(const R: TRationalNumber; Limit: Integer;
ALimitKind: TRatDigitMode; AFullOutput, ARound: Boolean;
AFirstDigitPos, AFirstSigDigitPos: PInteger): TArray<Integer>; overload;
var
dummy: Integer;
begin
_RatDigits(R, Limit, ALimitKind, AFullOutput, ARound, dummy, Result,
AFirstDigitPos, AFirstSigDigitPos);
end;
function GetDigits(const R: TRationalNumber; Limit: Integer;
ALimitKind: TRatDigitMode; AFullOutput, ARound: Boolean;
out AFracDigits: TArray<Integer>): TArray<Integer>; overload;
var
AIntDigits: TArray<Integer> absolute Result;
LDigits: TArray<Integer>;
LFirstDigitPos: Integer;
begin
LDigits := GetDigits(R, Limit, ALimitKind, AFullOutput, ARound, @LFirstDigitPos, nil);
if LFirstDigitPos >= 0 then
AIntDigits := Copy(LDigits, 0, LFirstDigitPos + 1)
else
AIntDigits := [0];
if LFirstDigitPos < 0 then
SetLength(AFracDigits, -LFirstDigitPos - 1);
TReverser<Integer>.Reverse(AIntDigits);
AFracDigits := AFracDigits + Copy(LDigits, LFirstDigitPos + 1);
end;
function GetDigit(X: TASR; Index: Integer): Integer; overload;
procedure ErrPrec;
begin
raise EMathException.Create('Cannot extract required digit from floating-point number due to insufficient precision.');
end;
var
IndexHigh: Integer;
IntVal: Int64;
const
PRECISION_DIGITS = 15;
begin
X := Abs(X);
if X = 0 then Exit(0);
IndexHigh := Floor(Log10(X));
if Index > IndexHigh then
Exit(0);
if IndexHigh - Index + 1 > PRECISION_DIGITS then
ErrPrec;
IntVal := Round(X * Math.IntPower(10, PRECISION_DIGITS - IndexHigh - 1));
Result := GetDigit(IntVal, PRECISION_DIGITS - 1 + (Index - IndexHigh));
end;
function IsInteger(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
begin
Result := InRange(X, TASI.MinValue + 1, TASI.MaxValue - 1) and SameValue(X, Round(X), Epsilon);
end;
function IsIntegerEx(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
begin
Result := InRange(X, TASI.MinValue + 1, TASI.MaxValue - 1) and SameValueEx(X, Round(X), Epsilon);
end;
function IsInteger32(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
begin
Result := InRange(X, Int32.MinValue + 1, Int32.MaxValue - 1) and SameValue(X, Round(X), Epsilon);
end;
function IsInteger32Ex(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
begin
Result := InRange(X, Int32.MinValue + 1, Int32.MaxValue - 1) and SameValueEx(X, Round(X), Epsilon);
end;
function IsInteger64(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
begin
Result := InRange(X, Int64.MinValue + 1, Int64.MaxValue - 1) and SameValue(X, Round(X), Epsilon);
end;
function IsInteger64Ex(const X: TASR; const Epsilon: Extended = 0): Boolean; inline;
begin
Result := InRange(X, Int64.MinValue + 1, Int64.MaxValue - 1) and SameValueEx(X, Round(X), Epsilon);
end;
function _fact32(const N: Integer): Integer;
var
i: Integer;
begin
if N < 0 then
raise EMathException.Create('Factorial of a negative number.');
if N > 12 then
raise EIntOverflow.Create('Integer overflow when computing 32-bit integer factorial.');
Result := 1;
for i := 2 to N do
Result := Result * i;
end;
function _fact64(const N: Integer): TASI;
var
i: Integer;
begin
if N < 0 then
raise EMathException.Create('Factorial of a negative number.');
if N > 20 then
raise EIntOverflow.Create('Integer overflow when computing 64-bit integer factorial.');
Result := 1;
for i := 2 to N do
Result := Result * i;
end;
function _factf(const N: Integer): TASR;
var
i: Integer;
begin
if N < 0 then
raise EMathException.Create('Factorial of a negative number.');
Result := 1;
for i := 2 to N do
Result := Result * i;
end;
function factorial(const N: Integer): TASR;
var
i: Integer;
begin
if N < 0 then
raise EMathException.Create('Factorial of a negative number.');
if N <= High(factorials) then
Exit(factorials[N]);
Result := factorials[High(factorials)];
for i := High(factorials) + 1 to N do
Result := Result * i;
end;
function combinations(n, k: Integer): TASR;
var
i: Integer;
begin
if not ((k >= 0) and (n >= k)) then
raise EMathException.Create('Binomial coefficient C(n, k) only defined for integers n and k satisfying 0 ≤ k ≤ n.');
if k > n - k then
k := n - k;
Result := 1.0;
for i := 0 to k - 1 do
Result := Result * (n - i) / (i + 1);
end;
function intcombinations(n, k: Integer): TASI;
var
i: Integer;
begin
if not ((k >= 0) and (n >= k)) then
raise EMathException.Create('Binomial coefficient C(n, k) only defined for integers n and k satisfying 0 ≤ k ≤ n.');
if k > n - k then
k := n - k;
Result := 1;
for i := 1 to k do
begin
if Result div i > Result.MaxValue div n then
Exit(0);
Result := (Result div i) * n + ((Result mod i) * n) div i;
Dec(n);
end;
end;
function binomial(const n, k: Integer): TASR; inline;
begin
Result := combinations(n, k);
end;
function permutations(n, k: Integer): TASR;
var
i: Integer;
begin
if not ((k >= 0) and (n >= k)) then
raise EMathException.Create('Binomial coefficient C(n, k) only defined for integers n and k satisfying 0 ≤ k ≤ n.');
Result := 1.0;
for i := n - k + 1 to n do
Result := Result * i;
end;
function intpermutations(n, k: Integer): TASI;
var
i: Integer;
begin
if not ((k >= 0) and (n >= k)) then
raise EMathException.Create('Binomial coefficient C(n, k) only defined for integers n and k satisfying 0 ≤ k ≤ n.');
Result := 1;
for i := n - k + 1 to n do
if Result < Result.MaxValue div i then
Result := Result * i
else
Exit(0);
end;
function lcm(const A, B: TASI): TASI; inline;
begin
if (a = 0) or (b = 0) then
Result := 0
else
Result := Abs((A div gcd(A, B)) * B);
end;
function lcm(const A, B: UInt64): UInt64; inline;
begin
if (a = 0) or (b = 0) then
Result := 0
else
Result := (A div gcd(A, B)) * B;
end;
function TryLCM(const A, B: TASI; var LCM: TASI): Boolean;
var
tmp: TASI;
begin
if (a = 0) or (b = 0) then
begin
Result := True;
LCM := 0;
end
else
begin
tmp := A div gcd(A, B);
Result := TInt64Guard.CanMul(tmp, B);
if Result then
LCM := tmp * B;
end;
end;
function lcm(const Values: array of TASI): TASI;
var
i: Integer;
begin
if Length(Values) = 0 then
raise EMathException.Create('Cannot compute the LCM of an empty list.');
Result := Abs(Values[0]);
for i := 1 to High(Values) do
Result := lcm(Result, Values[i]);
for i := 0 to High(Values) do
if (Values[i] <> 0) and (Result mod Values[i] <> 0) then
raise Exception.Create('LCM too large for the integer type.');
end;
function gcd(const A, B: TASI): TASI;
begin
if B = 0 then
Result := Abs(A)
else
Result := gcd(B, A mod B);
end;
function gcd(const A, B: UInt64): UInt64;
begin
if B = 0 then
Result := A
else
Result := gcd(B, A mod B);
end;
function gcd(const Values: array of TASI): TASI;
var
i: Integer;
begin
if Length(Values) = 0 then
raise EMathException.Create('Cannot compute the GCD of an empty list.');
Result := Abs(Values[0]);
for i := 1 to High(Values) do
Result := gcd(Result, Values[i]);
end;
function coprime(const A, B: TASI): Boolean; inline;
begin
Result := gcd(A, B) = 1;
end;
function NaiveTotient(const N: Integer): Integer;
var
i: Integer;
begin
if n <= 0 then
raise EMathException.Create('Totient only defined for positive integers.');
Result := 1;
for i := 2 to N - 1 do
if coprime(i, N) then
Inc(Result);
end;
function totient(const N: Integer): Integer;
var
prevfactor, factor: TASI;
res: Double;
begin
if n <= 0 then
raise EMathException.Create('Totient only defined for positive integers.');
res := n;
prevfactor := 1;
for factor in PrimeFactors(N) do
begin
if factor <> prevfactor then
res := res * (1 - 1 / factor);
prevfactor := factor;
end;
Result := Round(res);
Assert(IsIntegerEx(res));
end;
function cototient(const N: Integer): Integer; inline;
begin
Result := N - totient(N);
end;
var
_IsPrime: TBits;
_Primes: TList<Integer>;
_MöbiusCache: array of Int8;
_MertensCache: TList<Integer>;
const
PRIMES_ALLOC_BY = 65536;
MAX_ISPRIME_INDEX = 1000000000;
MAX_ISPRIME_INDEX_SQUARE = TASI(MAX_ISPRIME_INDEX) * MAX_ISPRIME_INDEX;
MAX_PRIME = 999999937;
INDEX_OF_MAX_PRIME = 50847534;
INVALID_MÖBIUS_VALUE = 2;
function ExpandPrimeCache(N: Integer): Boolean;
var
i, j: Integer;
begin
if Assigned(_IsPrime) and (_IsPrime.Size >= N + 1) then
Exit(False);
N := EnsureRange(Round(1.25*N), 1000, MAX_ISPRIME_INDEX);
if Assigned(_IsPrime) and (_IsPrime.Size >= N + 1) then
Exit(False);
Result := True;
FreeAndNil(_IsPrime);
FreeAndNil(_Primes);
_IsPrime := TBits.Create;
_Primes := TList<Integer>.Create;
try
_IsPrime.Size := N + 1;
_IsPrime[0] := False;
_IsPrime[1] := False;
for i := 2 to N do
_IsPrime[i] := True;
for i := 2 to Trunc(Sqrt(N)) do
if _IsPrime[i] then
begin
j := Sqr(i);
while j <= N do
begin
_IsPrime[j] := False;
Inc(j, i);
end;
end;
for i := 0 to _IsPrime.Size - 1 do
if _IsPrime[i] then
_Primes.Add(i);
except
FreeAndNil(_Primes);
FreeAndNil(_IsPrime);
raise;
end;
end;
procedure ClearPrimeCache;
begin
FreeAndNil(_Primes);
FreeAndNil(_IsPrime);
end;
function GetPrimeCacheMax: Integer;
begin
if _IsPrime = nil then
Result := -1
else
Result := _IsPrime.Size - 1;
end;
function GetPrimeCacheSizeBytes: Integer;
begin
Result := 0;
if Assigned(_IsPrime) then
Inc(Result, Ceil(_IsPrime.Size / 8));
if Assigned(_Primes) then
Inc(Result, SizeOf(Integer) * _Primes.Count);
end;
function IsPrime(const N: TASI): Boolean;
var
i: Integer;
attempt: Integer;
begin
if N < 2 then
Exit(False);
if N > MAX_ISPRIME_INDEX_SQUARE then
raise EMathException.Create('Cannot test integers this large for primality.');
ExpandPrimeCache(10000);
for attempt := 0 to 1 do
begin
if N <= _IsPrime.Size - 1 then
Exit(_IsPrime[N]);
for i := 0 to _Primes.Count - 1 do
if _Primes[i] > Sqrt(N) then
Exit(True)
else if N mod _Primes[i] = 0 then
Exit(False);
ExpandPrimeCache(Round(Sqrt(N)));
end;
raise EMathException.Create('Couldn''t determine if the integer is prime.');
end;
function NthPrime(N: Integer): Integer;
begin
if N < 1 then
raise EMathException.Create('To obtain the Nth prime number, N must be a positive integer.');
if N > INDEX_OF_MAX_PRIME then
raise EMathException.Create('Cannot test integers this large for primality.');
ExpandPrimeCache(10000);
if N > 6 then
ExpandPrimeCache(Ceil(N * Ln(N*Ln(N))));
Dec(N);
if N <= _Primes.Count - 1 then
Exit(_Primes[N]);
raise Exception.Create('Couldn''t obtain the Nth prime number.');
end;
function NextPrime(const N: TASI): TASI;
var
i: TASI;
begin
if N < 2 then
Exit(2);
for i := N + 1 to TASI.MaxValue do
if IsPrime(i) then
Exit(i);
raise EMathException.Create('Cannot test integers this large for primality.');
end;
function PreviousPrime(const N: TASI): TASI;
var
i: TASI;
begin
if N <= 2 then
raise EMathException.Create('There are no prime numbers smaller than 2.');
for i := N - 1 downto 2 do
if IsPrime(i) then
Exit(i);
raise Exception.Create('Internal error.');
end;
function PrimePi(const N: Integer): Integer;
begin
if N > MAX_ISPRIME_INDEX then
raise EMathException.Create('Cannot test integers this large for primality.');
ExpandPrimeCache(N);
if _Primes.BinarySearch(N, Result) then
Inc(Result);
end;
function Primorial(const N: Integer): TASR;
var
p: TASI;
begin
if N < 0 then
Result := 1
else if InRange(N, Low(IntPrimorials), High(IntPrimorials)) then
Result := IntPrimorials[N]
else
begin
Result := IntPrimorials[High(IntPrimorials)];
p := NextPrime(High(IntPrimorials));
while p <= N do
begin
Result := p * Result;
p := NextPrime(p);
end;
end;
end;
function Fibonacci(const N: Integer): TASR;
var
i: Integer;
prev, prevprev: TASR;
begin
if N < 0 then
raise EInvalidArgument.Create('To obtain the Nth Fibonacci number, N must be zero or a positive integer.');
if N = 0 then
Exit(0);
if N = 1 then
Exit(1);
prev := 1;
prevprev := 0;
Result := 1;
for i := 2 to N do
begin
Result := prev + prevprev;
prevprev := prev;
prev := Result;
end;
end;
function Lucas(const N: Integer): TASR;
var
i: Integer;
prev, prevprev: TASR;
begin
if N < 0 then
raise EInvalidArgument.Create('To obtain the Nth Lucas number, N must be zero or a positive integer.');
if N = 0 then
Exit(2);
if N = 1 then
Exit(1);
prev := 1;
prevprev := 2;
Result := 3;
for i := 2 to N do
begin
Result := prev + prevprev;
prevprev := prev;
prev := Result;
end;
end;
function Floor64(const X: TASR): TASI;
begin
Result := TASI(Trunc(X));
if Frac(X) < 0 then
Dec(Result);
end;
function Ceil64(const X: TASR): TASI;
begin
Result := TASI(Trunc(X));
if Frac(X) > 0 then
Inc(Result);
end;
function imod(const x, y: Integer): Integer; inline;
begin
Result := x mod y;
if Result < 0 then
Inc(Result, y);
end;
function imod(const x, y: TASI): TASI; inline;
begin
Result := x mod y;
if Result < 0 then
Inc(Result, y);
end;
function rmod(const x, y: TASR): TASR; inline;
begin
Result := x - Floor64(x / y) * y;
end;
function modulo(const x, y: Integer): Integer; overload; inline;
begin
Result := imod(x, y);
end;
function modulo(const x, y: TASR): TASR; overload; inline;
begin
Result := rmod(x, y);
end;
function PrimeFactors(const N: TASI; AOnlyUnique: Boolean = False): TArray<TASI>;
var
i: Integer;
a: TASI;
begin
if N < 1 then
raise EMathException.Create('Cannot find the prime factors of a non-positive integer.');
SetLength(Result, 0);
if N = 1 then
Exit;
if Assigned(_IsPrime) and (N <= _IsPrime.Size - 1) and IsPrime(N) then
begin
SetLength(Result, 1);
Result[0] := N;
Exit;
end;
a := N;
ExpandPrimeCache(10000);
for i := 0 to _Primes.Count - 1 do
begin
if _Primes[i] > Sqrt(a) then
begin
if AOnlyUnique then
TArrBuilder<TASI>.AddUnique(Result, a)
else
TArrBuilder<TASI>.Add(Result, a);
Exit;
end;
while a mod _Primes[i] = 0 do
begin
if AOnlyUnique then
TArrBuilder<TASI>.AddUnique(Result, _Primes[i])
else
TArrBuilder<TASI>.Add(Result, _Primes[i]);
a := a div _Primes[i];
if a = 1 then Exit;
end
end;
if ExpandPrimeCache(Round(Sqrt(N))) then
Result := PrimeFactors(N, AOnlyUnique)
else
raise EMathException.Create('Couldn''t factor integer.');
end;
function PrimeFactorsWithMultiplicity(N: TASI): TArray<TIntWithMultiplicity>;
begin
Result := CollapseWithMultiplicity(PrimeFactors(N));
end;
function Radical(N: TASI): TASI;
var
p: TASI;
begin
Result := 1;
for p in PrimeFactors(N, True) do
Result := Result * p;
end;
function IsSquareFree(N: TASI): Boolean;
var
i: Integer;
a: TASI;
begin