ASNum.pas

File name
C:\Users\Andreas Rejbrand\Documents\Dev\AlgoSim\ASNum.pas
Date exported
Time exported
Formatting processor
TPascalFormattingProcessor
unit ASNum;

{ **************************************************************************** }
{ Rejbrand AlgoSim numerics library                                            }
{ Copyright © 2017-2018 Andreas Rejbrand                                       }
{ https://english.rejbrand.se/                                                 }
{ **************************************************************************** }

{$DEFINE REALCHECK}
{$WARN SYMBOL_PLATFORM OFF}
{$WARN DUPLICATE_CTOR_DTOR OFF}

interface

uses
  SysUtils, Types, Math, Generics.Defaults, Generics.Collections, GenHelpers;

threadvar
  GTYieldProc: TObjProc;

/// <summary>Returns (-1)^N.</summary>
/// <returns>(-1)^N</returns>
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;

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; // mainly for debugging
    function sstr: string; // mainly for debugging
    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;

// Elementary functions for real numbers
// exp, ln from lib
function pow(const a, b: TASR): TASR; inline;
function intpow(const a, b: TASI): TASI;
// sqrt from Math.pas
function cbrt(const X: TASR): TASR; inline;
function frt(const X: TASR): TASR; inline;
// sin, cos, tan, cot, sec, csc from lib
// arcsin, arccos, arctan, arccot, arcsec, arccsc from lib
// sinh, cosh from lib
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;
// arccosh, arctanh, arccoth, arcsech from lib
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; // mainly for debugging
    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;

// Elementary functions for complex numbers
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;

// Integer functions

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; // alias
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;

// Differentiation
function differentiate(AFunction: TRealFunctionRef; const X: TASR;
  const ε: TASR = 1E-6): TASR;

// Numerical integration
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;

// Sums, products, max, and min
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;


// Special functions
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;

// Vectors
type
  /// <summary>The AlgoSim real vector type.</summary>
  /// <remarks>A TASRArray can be cast to a TRealVector. That is, if <c>v</c> is
  ///  a TASRArray, then <c>TRealVector(v)</c> is legal. Similarly, a TRealVector
  ///  can be cast to a TASRArray (or the Data member can be used).</remarks>
  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; {copy}
    function Reduce(AStep: Integer): TRealVector; {copy}
    function Clone: TRealVector;
    function Subvector(AFrom, ATo: Integer): TRealVector; overload;
    function Subvector(const AIndices: array of Integer): TRealVector; overload;
    function Sort: TRealVector; overload; {in place, no copy}
    function Sort(AComparer: IComparer<TASR>): TRealVector; overload; {in place, no copy}
    function ReverseSort: TRealVector; inline; {in place, no copy}
    function Shuffle: TRealVector; {in place, no copy}
    function Unique: TRealVector; {copy}
    function UniqueAdj: TRealVector; {copy}
    function UniqueEps(const Epsilon: TASR = 0): TRealVector; {copy}
    function UniqueAdjEps(const Epsilon: TASR = 0): TRealVector; {copy}
    function Reverse: TRealVector; {in place, no copy}
    function Apply(AFunction: TRealFunctionRef): TRealVector; {copy}
    function Filter(APredicate: TPredicate<TASR>): TRealVector; {copy}
    function Replace(APredicate: TPredicate<TASR>; const ANewValue: TASR): TRealVector; overload; {copy}
    function Replace(const AOldValue, ANewValue: TASR; const Epsilon: Extended = 0): TRealVector; overload; {copy}
    function Replace(const ANewValue: TASR): TRealVector; overload; {copy}
    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
  /// <summary>The AlgoSim complex vector type.</summary>
  /// <remarks>A TASCArray can be cast to a TComplexVector. That is, if <c>v</c>
  ///  is a TASCArray, then <c>TComplexVector(v)</c> is legal.</remarks>
  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; {copy}
    function Reduce(AStep: Integer): TComplexVector; {copy}
    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; {in place, no copy}
    function Shuffle: TComplexVector; {in place, no copy}
    function Unique: TComplexVector; {copy}
    function UniqueAdj: TComplexVector; {copy}
    function UniqueEps(const Epsilon: TASR = 0): TComplexVector; {copy}
    function UniqueAdjEps(const Epsilon: TASR = 0): TComplexVector; {copy}
    function Reverse: TComplexVector; {in place, no copy}
    function Apply(AFunction: TComplexFunctionRef): TComplexVector; {copy}
    function Filter(APredicate: TPredicate<TASC>): TComplexVector; {copy}
    function Replace(APredicate: TPredicate<TASC>;
      const ANewValue: TASC): TComplexVector; overload; {copy}
    function Replace(const AOldValue, ANewValue: TASC;
      const Epsilon: Extended = 0): TComplexVector; overload; {copy}
    function Replace(const ANewValue: TASC): TComplexVector; overload; {copy}
    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);

// Matrices
type
  TDuplicateFinder<T> = record
    class function ContainsDuplicates(Arr: array of T): Boolean; static;
    class function PresortedContainsDuplicates(const Arr: array of T): Boolean; static;
  end;

  /// <summary>A record that holds the size of a matrix.</summary>
  /// <remarks><para>If cast from an integer N, the matrix size is N×N. Such a
  ///  cast can be implicit.</para>
  ///  <para>TMatrixSize records can be compared using = and <>. Records can be
  ///  added using +; the result is the size of the direct sum of two matrices
  ///  of the given sizes.</para></remarks>
  TMatrixSize = record
  strict private
    FRows, FCols: Integer;
  private
    constructor CreateUnsafe(Rows, Cols: Integer); overload; // not child safe
    function LessenedSize: TMatrixSize; inline; // not child safe
  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; // size of direct sum
  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 // Due to a D2009 bug, a generic type cannot be used
    RowOperationType: TRowOperationType;
    Row1, Row2: Integer;
    Factor: TASR;
  end;
  TRealRowOpSequence = array of TRealRowOperationRecord;

  /// <summary>The AlgoSim real matrix type.</summary>
  TRealMatrix = record
  strict private
    const
      InitialRowOpSeqSize = 16; // strategy: then double
    var
      FSize: TMatrixSize;
      FElements: TASRArray;
      FRowOpSeq: TRealRowOpSequence;
      FRowOpCount: Integer;
      FRowOpFactor: TASR;
      _FCollectRowOpData: Boolean; // never set directly -- use BeginCollectRowOpData and EndCollectRowOpData

    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;

    /// <summary>Clears any existing row op data, initializes the row op buffers
    /// and sets the collect row op data flag to true, so that row op data will
    /// be collected after this call.</summary>
    procedure BeginCollectRowOpData;

    /// <summary>Clears any existing row op data and sets the collect row op data
    ///  flag to false.</summary>
//    procedure EndCollectRowOpData;

    /// <summary>Adds a row operation record to the buffer.</summary>
    /// <remarks>Warning! This method MUST NOT be called if FCollectRowOpData
    ///  is set to false. Doing so will result in undefined behaviour.</remarks>
    procedure AddRowOpRecord(AType: TRowOperationType; ARow1, ARow2: Integer;
      AFactor: TASR = 0);

    procedure GetHouseholderMap(const AVect: TRealVector;
      out tau, gamma: TASR; out u: TRealVector);

    /// <summary>Returns the two complex eigenvalues of the second leading
    ///  principal submatrix of the matrix.</summary>
    /// <remarks>The matrix MUST has at least two columns and it MUST have
    ///  at least two rows.</remarks>
    function Eigenvalues2x2: TComplexVector;

  private
    /// <summary>Creates a new real matrix record but does not allocate data
    ///  for the matrix entries. The <c>Data</c> property is set to <c>nil</c>.
    ///  </summary>
    constructor CreateWithoutAllocation(const ASize: TMatrixSize);

    /// <summary>Returns true iff the matrix is an empty matrix, that is, if
    ///  it has exactly zero elements.</summary>
    function IsEmpty: Boolean; inline;

    /// <summary>Performs a fast LU decomposition of the matrix, which has to be
    ///  an n×n square.</summary>
    /// <param name="A">An n×n matrix that receives the elements of L and U.
    ///  The upper triangular part, including the main diagonal, receives the
    ///  elments of U, and the lower triangular part, excluding the diagonal,
    ///  receives the elements of L that are below the main diagonal. The main
    ///  diagonal elements of L are all equal to 1, so A contains all
    ///  information about L and U.</param>
    /// <param name="P">An n-dimensional array of integers that contain all
    ///  information about the row exchanges made during the LU decomposition.
    ///  If the kth element of the array is m, then at the kth step, rows k and
    ///  m (>= k) were exchanged. If m = k, no row exchange was made at this
    ///  step. If m = -1, no row exchange was made at this step, and, in
    ///  addition, it was found that the matrix was singular at this step.
    ///  </param>
    /// <remarks>This method is as fast as possible. The public method
    ///  LU uses this method to construct the actual matrixes L, U, and P in
    ///  the LU decomposition Self = P.Transpose * L * U.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    procedure DoQuickLU(out A: TRealMatrix; out P: TIndexArray);

    procedure RequireNonEmpty; inline;

  public

    // ***
    // Constructors
    // ***

    /// <summary>Allocates memory for a new real matrix of a given size.
    ///  Equivalently, it creates a new real matrix of a given size with
    ///  elements with undefined initial values.</summary>
    /// <remarks>Creating an uninitialized matrix of a given size is faster
    ///  compared to creating an initialized matrix of the same size.</remarks>
    constructor CreateUninitialized(const ASize: TMatrixSize);

    /// <summary>Copy constructor for <c>TRealMatrix</c>.</summary>
    constructor Create(const AMatrix: TRealMatrix); overload;

    /// <summary>Creates a new real matrix from an array of real-number arrays.
    ///  Each real-number array will form a row in the matrix.</summary>
    /// <param name="Elements">An open array of real-number arrays that will
    ///  be used to create the matrix.</param>
    /// <remarks>The real-number arrays have to have the same length; this
    ///  common length becomes the number of columns of the matrix.</remarks>
    /// <see cref="CreateFromRows" />
    /// <see cref="CreateFromColumns" />
    constructor Create(const Elements: array of TASRArray); overload;

    /// <summary>Creates a new real matrix from a single array of real numbers
    ///  and a given matrix format; the real numbers will populate the matrix
    ///  top to bottom and left to right. The real number with index i + 1 will
    ///  be to the right of the element with index i unless the latter is the
    ///  right-most element of a row, in which case element i + 1 will become
    ///  the first element of the next row.</summary>
    /// <param name="Elements">An array of real numbers.</param>
    /// <param name="Cols">The number of columns in the matrix.</param>
    /// <remarks>The number of rows in the matrix will be the length of
    ///  <c>Elements</c> divided by <c>Cols</c>, which has to be an integer.
    ///  The default value of <c>Cols</c> is 1, which will create a single-column
    ///  matrix (a column vector) with the length of <c>Elements</c>.</remarks>
    constructor Create(const Elements: array of TASR; Cols: Integer = 1); overload;

    /// <summary>Creates a new real matrix with a given format and fills it with
    ///  a given real number.</summary>
    /// <param name="ASize">The format of the matrix.</param>
    /// <param name="AVal">The real number to fill the matrix with.</param>
    /// <remarks>The default value of <c>AVal</c> is zero, which will create the
    ///  zero matrix of the given size.</remarks>
    constructor Create(const ASize: TMatrixSize; const AVal: TASR = 0); overload;

    /// <summary>Creates a new real matrix from a sequence of real vectors, which
    ///  will form the rows of the matrix.</summary>
    /// <param name="Rows">An open array of real vectors; each real vector will
    ///  form a row in the matrix. The number of real vectors will become the
    ///  number of rows in the matrix.</param>
    /// <remarks>The vectors need to have the same length, which will become the
    ///  number of columns of the matrix.</remarks>
    constructor CreateFromRows(const Rows: array of TRealVector);

    /// <summary>Creates a new real matrix from a sequence of real vectors, which
    ///  will form the columns of the matrix.</summary>
    /// <param name="Columns">An open array of real vectors; each real vector will
    ///  form a column in the matrix. The number of real vectors will become the
    ///  number of columns in the matrix.</param>
    /// <remarks>The vectors need to have the same length, which will become the
    ///  number of rows of the matrix.</remarks>
    constructor CreateFromColumns(const Columns: array of TRealVector);

    /// <summary>Creates a new real matrix A from two real vectors u and v; the
    ///  matrix is the outer product of the vectors, that is, A_ij = u_i * v_j.
    ///  </summary>
    /// <param name="u">The first real vector.</param>
    /// <param name="v">The second real vector.</param>
    /// <remarks>The vectors need not be of the same dimension.</remarks>
    constructor Create(const u, v: TRealVector); overload;

    /// <summary>Creates a new real matrix using given format and a real-valued
    ///  function of valid matrix indices.</summary>
    /// <param name="ASize">The format of the matrix.</param>
    /// <param name="AFunction">A reference to a real-valued function of matrix
    ///  indices (integers). This function is used to populate the matrix. The
    ///  element (i, j) in the matrix is the value AFunction(i, j). Here,
    ///  row and column indices are 1-based.</param>
    constructor Create(const ASize: TMatrixSize;
      AFunction: TMatrixIndexFunction<TASR>); overload;

    /// <summary>Creates a new diagonal real matrix using a given array of real
    ///  numbers.</summary>
    /// <param name="Elements">An open array of real numbers that will be used as
    ///  the diagonal entries of the matrix.</param>
    constructor CreateDiagonal(const Elements: array of TASR); overload;

    /// <summary>Creates a new diagonal real matrix using a given real vector.
    ///  </summary>
    /// <param name="Elements">A real vector that will be used to populate the
    ///  diagonal of the matrix.</param>
    constructor CreateDiagonal(const Elements: TRealVector); overload;

    /// <summary>Creates a new diagonal real matrix of a given size and a (fixed)
    ///  given value on the diagonal.</summary>
    /// <param name="ASize">The format of the matrix.</param>
    /// <param name="AVal">The value on the main diagonal.</param>
    constructor CreateDiagonal(ASize: Integer; AVal: TASR = 1); overload;

    /// <summary>Creates a new real matrix using pre-defined blocks.</summary>
    /// <param name="Blocks">An open array of the blocks of the matrix in row-
    ///  major order.</param>
    /// <param name="Cols">The number of blocks per row of the block matrix;
    ///  this value must be at least 1.</param>
    /// <remarks>Each block must contain at least a single element. In each row
    ///  of the block matrix, the number of rows in each block must be the same.
    ///  In each column of the block matrix, the number of columns in each block
    ///  must be the same. If the number of elements of <c>Blocks</c> is not
    ///  divisible by <c>Cols</c>, it will be padded with zero matrices until
    ///  the last row is filled.</remarks>
    /// <exception cref="EMathException">Raised if any of the constraints
    ///  specified above are violated.</exception>
    constructor Create(const Blocks: array of TRealMatrix; Cols: Integer = 2); overload;

    // ***
    // Properties
    // ***

    /// <summary>The size of the matrix.</summary>
    property Size: TMatrixSize read FSize write SetSize;

    /// <summary>Gives access to the matrix elements using row and column indices.
    ///  </summary>
    /// <remarks>Row and column indices are zero-based.</remarks>
    property Elements[Row, Col: Integer]: TASR read GetElement write SetElement; default;

    /// <summary>Gives access to the rows of the matrix.</summary>
    /// <remarks>Row indices are zero-based.</remarks>
    property Rows[Index: Integer]: TRealVector read GetRow write _DoSetRow;

    /// <summary>Gives access to the columns of the matrix.</summary>
    /// <remarks>Column indices are zero-based.</remarks>
    property Cols[Index: Integer]: TRealVector read GetCol write _DoSetCol;

    /// <summary>Gives access to the first column.</summary>
    property FirstColumn: TRealVector read GetFirstCol write SetFirstCol;

    /// <summary>Gives access to the last column.</summary>
    property LastColumn: TRealVector read GetLastCol write SetLastCol;

    /// <summary>Gives access to the main diagonal of the matrix.</summary>
    /// <remarks>Diagonal indices are zero-based.</remarks>
    property MainDiagonal: TRealVector read GetMainDiagonal write SetMainDiagonal;

    /// <summary>Gives access to the superdiagonal of the matrix.</summary>
    /// <remarks>Superdiagonal indices are zero-based.</remarks>
    property SuperDiagonal: TRealVector read GetSuperDiagonal write SetSuperDiagonal;

    /// <summary>Gives access to the subdiagonal of the matrix.</summary>
    /// <remarks>Subdiagonal indices are zero-based.</remarks>
    property SubDiagonal: TRealVector read GetSubDiagonal write SetSubDiagonal;

    /// <summary>Gives access to the antidiagonal of the matrix.</summary>
    /// <remarks>Antidiagonal indices are zero-based.</remarks>
    property AntiDiagonal: TRealVector read GetAntiDiagonal write SetAntiDiagonal;

    /// <summary>Gives access to the underlying data as an array of real numbers
    ///  in row-major order.</summary>
    /// <remarks>This may be cast to a TRealVector. Notice that all changes to the
    ///  data will affect the original matrix. The AsVector member function instead
    ///  makes a copy of the data.</remarks>
    /// <see cref="AsVector" />
    property Data: TASRArray read FElements write FElements;

    /// <summary>Returns a pointer to the first element of a specific row.
    ///  The row index MUST be valid.</summary>
    property RowData[Index: Integer]: PASR read GetRowData;

    /// <summary>Returns a pointer to the first element of a specific row.
    ///  The function validates that the row index is valid, and raises an
    ///  exception if it is not.</summary>
    property SafeRowData[Index: Integer]: PASR read SafeGetRowData;

    /// <summary>The total size in bytes of this matrix.</summary>
    property MemorySize: Int64 read GetMemorySize;

    // ***
    // Operators
    // ***

    /// <summary>Creates a single-column matrix from a vector.</summary>
    class operator Implicit(const u: TRealVector): TRealMatrix;

    /// <summary>Treats a single-column matrix as a vector.</summary>
    class operator Explicit(const A: TRealMatrix): TRealVector; inline;

    /// <summary>Creates a single-entry matrix from a real number.</summary>
    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;

    // ***
    // Tests
    // ***

    /// <summary>Returns true iff the matrix consists of a single row.</summary>
    function IsRow: Boolean; inline;

    /// <summary>Returns true iff the matrix consists of a single column.</summary>
    function IsColumn: Boolean; inline;

    /// <summary>Returns true iff the matrix is square.</summary>
    function IsSquare: Boolean; inline;

    /// <summary>Returns true iff the matrix is the N-dimensional identity
    ///  matrix.</summary>
    function IsIdentity(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is the m×n zero matrix (that is, if
    ///  it consists of only zeros).</summary>
    function IsZeroMatrix(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff all non-zero entries are located on the main
    ///  diagonal.</summary>
    function IsDiagonal(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is square and all non-zero entries
    ///  are located on the antidiagonal.</summary>
    function IsAntiDiagonal(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is the N-dimensional reversal
    ///  matrix.</summary>
    function IsReversal(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff all non-zero entries are located on or above
    ///  the main diagonal.</summary>
    /// <remarks>The matrix need not be square in order to be upper diagonal.
    ///  </remarks>
    function IsUpperTriangular(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff all non-zero entries are located on or below
    ///  the main diagonal.</summary>
    /// <remarks>The matrix need not be square in order to be lower diagonal.
    ///  </remarks>
    function IsLowerTriangular(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is upper triangular or lower
    ///  triangular.</summary>
    function IsTriangular(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns the column index of the first non-zero element on a
    ///  given row.</summary>
    /// <returns>The column index of the first non-zero element on the given
    ///  row or -1 if the row consists of only zeros.</returns>
    function PivotPos(ARow: Integer; const Epsilon: Extended = 0): Integer;

    /// <summary>Returns whether or not a given row consists of only zeros.
    ///  </summary>
    function IsZeroRow(ARow: Integer; const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns whether or not a given row consists of only zeros except
    ///  for the right-most entry on the row.</summary>
    function IsEssentiallyZeroRow(ARow: Integer;
      const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is in row echelon form.</summary>
    function IsRowEchelonForm(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is in reduced row echelon form.
    ///  A matrix that is in reduced row echelon form is also in row echelon
    ///  form.</summary>
    function IsReducedRowEchelonForm(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determins whether the matrix is scalar, that is, if it is a
    ///  multiple of the identity matrix.</summary>
    function IsScalar(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is symmetric.</summary>
    function IsSymmetric(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is skew-symmetric.</summary>
    function IsSkewSymmetric(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is orthogonal.</summary>
    function IsOrthogonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is normal.</summary>
    function IsNormal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix consists only of ones and zeros.
    ///  </summary>
    function IsBinary(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is a permutation matrix.</summary>
    function IsPermutation(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is circulant.</summary>
    /// <remarks>A circulant matrix need not be square.</remarks>
    function IsCirculant(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is a Toeplitz matrix.</summary>
    /// <remarks>A Toeplitz matrix need not be square.</remarks>
    function IsToeplitz(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is a Hankel matrix.</summary>
    /// <remarks>A Hankel matrix need not be square.</remarks>
    function IsHankel(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is upper Hessenberg.</summary>
    /// <remarks>An upper Hessenberg matrix is always square.</remarks>
    function IsUpperHessenberg(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is lower Hessenberg.</summary>
    /// <remarks>A lower Hessenberg matrix is always square.</remarks>
    function IsLowerHessenberg(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is tridiagonal.</summary>
    /// <remarks>A tridiagonal matrix is always square.</remarks>
    function IsTridiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is upper bidiagonal.</summary>
    /// <remarks>An upper bidiagonal matrix is always square.</remarks>
    function IsUpperBidiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is lower bidiagonal.</summary>
    /// <remarks>A lower bidiagonal matrix is always square.</remarks>
    function IsLowerBidiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is bidiagonal.</summary>
    /// <remarks>A bidiagonal matrix is always square.</remarks>
    function IsBidiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is centrosymmetric.</summary>
    function IsCentrosymmetric(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is a Vandermonde matrix.</summary>
    /// <remarks>A Vandermonde matrix need not be square.</remarks>
    function IsVandermonde(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix commutes with a given matrix.
    ///  </summary>
    function CommutesWith(const A: TRealMatrix;
      const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is idempotent, that is, if A^2 =
    ///  A.</summary>
    /// <remarks>An idempotent matrix is always square.</remarks>
    function IsIdempotent(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is an involution, that is, if A^2 =
    ///  I.</summary>
    /// <remarks>An involutory matrix is always square.</remarks>
    function IsInvolution(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is positive definite.</summary>
    /// <remarks>A positive definite matrix is always square and symmetric.
    ///  </remarks>
    function IsPositiveDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is positive semidefinite.</summary>
    /// <remarks>A positive semidefinite matrix is always square and symmetric.
    ///  </remarks>
    function IsPositiveSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is negative definite.</summary>
    /// <remarks>A negative definite matrix is always square and symmetric.
    ///  </remarks>
    function IsNegativeDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is negative semidefinite.</summary>
    /// <remarks>A negative semidefinite matrix is always square and symmetric.
    ///  </remarks>
    function IsNegativeSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is indefinite.</summary>
    /// <remarks>An indefinite matrix is always square and symmetric.</remarks>
    function IsIndefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the square matrix is nilpotent.</summary>
    /// <exception cref="EMathException">Raised if the matrix isn't square or if
    ///  the method failed to determine nilpotency.</exception>
    function IsNilpotent(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns the nilpotency index of the square matrix.</summary>
    /// <returns>The nilpotency index if the matrix is nilpotent or -1 otherwise.
    ///  </returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square or if
    ///  the method failed to determine nilpotency.</exception>
    function NilpotencyIndex(const Epsilon: Extended = 0): Integer;

    /// <summary>Returns true iff the matrix is diagonally dominant, that is,
    ///  iff for every row the absolute value of the diagonal entry is greater
    ///  than or equal to the deleted absolute row sum.</summary>
    /// <remarks>Only defined for square matrices.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function IsDiagonallyDominant: Boolean;

    /// <summary>Returns true iff the matrix is strictly diagonally dominant,
    ///  that is, iff for every row the absolute value of the diagonal entry is
    ///  greater than the deleted absolute row sum.</summary>
    /// <remarks>Only defined for square matrices.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function IsStrictlyDiagonallyDominant: Boolean;

    /// <summary>Returns true iff the matrix is positive, that is, if all
    ///  elements are positive.</summary>
    function IsPositive(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns true iff the matrix is non-negative, that is, if all
    ///  elements are non-negative.</summary>
    function IsNonNegative(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns true iff the matrix is negative, that is, if all
    ///  elements are negative.</summary>
    function IsNegative(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns true iff the matrix is non-positive, that is, if all
    ///  elements are non-positive.</summary>
    function IsNonPositive(const Epsilon: Extended = 0): Boolean; inline;

    // ***
    // Operations and computations
    // ***

    /// <summary>Sets all elements to zero above the main diagonal.</summary>
    procedure MakeLowerTriangular;

    /// <summary>Sets all elements to zero below the main diagonal.</summary>
    procedure MakeUpperTriangular;

    /// <summary>Sets all elements to zero below the main subdiagonal.</summary>
    procedure MakeUpperHessenberg;

    /// <summary>Returns the square of the matrix.</summary>
    function Sqr: TRealMatrix; inline;

    /// <summary>Returns the transpose of the matrix.</summary>
    function Transpose: TRealMatrix;

    /// <summary>Returns the Hermitian square of the matrix.</summary>
    function HermitianSquare: TRealMatrix; inline;

    /// <summary>Returns the modulus of the matrix.</summary>
    /// <remarks>The modulus is the unique positive-semidefinite square root of
    ///  the Hermitian square.</remarks>
    function Modulus: TRealMatrix; inline;

    /// <summary>Returns the determinant of the matrix.</summary>
    function Determinant: TASR;

    /// <summary>Returns the trace of the matrix.</summary>
    function Trace: TASR;

    /// <summary>Returns the inverse of the matrix, if it exists.</summary>
    /// <exception cref="EMathException">Thrown if the matrix doesn't have an
    ///  inverse. This happens if the matrix isn't square or if its determinant
    ///  is zero (iff it doesn't have full rank).</exception>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    /// <see cref="TryInvert" />
    function Inverse: TRealMatrix;

    /// <summary>Attempts to invert the matrix, which has to be square.</summary>
    /// <param name="AInverse">Receives the inverse of the matrix, if it can
    ///  be computed.</param>
    /// <returns>True iff the matrix can be inverted, that is, if and only if
    ///  it is non-singular.</returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    /// <see cref="Inverse" />
    function TryInvert(out AInverse: TRealMatrix): Boolean;

    /// <summary>Returns the Moore-Penrose pseudoinverse of the matrix.</summary>
//    function Pseudoinverse: TRealMatrix; experimental;

    /// <summary>Returns the rank of the matrix.</summary>
    function Rank: Integer; inline;

    /// <summary>Returns the nullity of the matrix.</summary>
    function Nullity: Integer; inline;

    /// <summary>Returns the condition number of the matrix using a specified
    ///  matrix norm.</summary>
    /// <exception cref="EMathException">Raised if the matrix isn't square and
    ///  invertible.</exception>
    function ConditionNumber(p: Integer = 2): TASR;

    /// <summary>Determines whether the matrix, which has to be square, is
    ///  singular.</summary>
    /// <returns>True iff the matrix is singular.</returns>
    /// <exception cref="EMathException">Thrown if the matrix isn't square.
    ///  </exception>
    function IsSingular: Boolean;

    /// <summary>Returns the Frobenius norm of the matrix, that is, the square
    ///  root of the sum of the squares of all its elements.</summary>
    /// <remarks>This norm is submultiplicative.</remarks>
    function Norm: TASR; inline;

    /// <summary>Returns the square of the Frobenius norm of the matrix.</summary>
    function NormSqr: TASR; inline;

    /// <summary>Returns the (vector) p-norm of the matrix.</summary>
    function pNorm(const p: TASR): TASR; inline;

    /// <summary>Returns the max norm of the matrix, that is, the largest
    ///  absolute value of its elements.</summary>
    /// <remarks>This norm is NOT submultiplicative.</remarks>
    function MaxNorm: TASR; inline;

    /// <summary>Returns the sum norm of the matrix, that is, the sum of all
    ///  absolute values of the entries.</summary>
    /// <remarks>This norm is submultiplicative.</remarks>
    function SumNorm: TASR; inline;

    /// <summary>Returns the (vector) k-norm of the matrix.</summary>
    function kNorm(const k: Integer): TASR; inline;

    /// <summary>Returns the maximum column sum norm of the matrix.</summary>
    /// <remarks>This norm is submultiplicative. It is the operator norm induced
    ///  by the vector 1-norm (l_1).</remarks>
    function MaxColSumNorm: TASR;

    /// <summary>Returns the maximum row sum norm of the matrix.</summary>
    /// <remarks>This norm is submultiplicative. It is the operator norm induced
    ///  by the vector maximum norm (l_∞).</remarks>
    function MaxRowSumNorm: TASR;

    /// <summary>Returns the spectral norm of the matrix, that is, the largest
    ///  singular value.</summary>
    /// <remarks>This norm is submultiplicative. It is the operator norm induced
    ///  by the vector 2-norm (l_2).</remarks>
    function SpectralNorm: TASR; inline;

    /// <summary>Returns a deleted absolute row sum, that is, the sum of all
    ///  absolute-value entries on a given row except the row's diagonal entry
    ///  (if it has a diagonal entry).</summary>
    function DeletedAbsoluteRowSum(ARow: Integer): TASR;

    /// <summary>Swaps two rows in a matrix (in place).</summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowSwap(ARow1, ARow2: Integer): TRealMatrix;

    /// <summary>Multiplies a single row in a matrix with a constant real number
    ///  (in place).</summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowScale(ARow: Integer; AFactor: TASR): TRealMatrix;

    /// <summary>Adds a multiple of one row to a second row in a matrix (in place).
    ///  </summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowAddMul(ATarget, ASource: Integer; AFactor: TASR;
      ADefuzz: Boolean = False; AFirstCol: Integer = 0): TRealMatrix;

    /// <summary>Performs a given row operation to the matrix. The row operation
    ///  is given by a TRowOperationRecord.</summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowOp(const ARowOp: TRealRowOperationRecord): TRealMatrix;

    /// <summary>Uses row operations to find a row echelon form of the matrix.
    ///  </summary>
    /// <param name="CollectRowOps">If set to true, the row operations performed
    ///  will be saved to the row op buffer.</param>
    /// <returns>Returns a row echelon form of the matrix.</returns>
    /// <remarks>The original matrix is not altered by this method; all row
    ///  operations are performed on a copy.</remarks>
    function RowEchelonForm(CollectRowOps: Boolean = False): TRealMatrix;

    /// <summary>Uses row operations to find the reduced row echelon form of the
    ///  matrix.</summary>
    /// <param name="CollectRowOps">If set to true, the row operations performed
    ///  will be saved to the row op buffer.</param>
    /// <returns>Returns the reduced row echelon form of the matrix.</returns>
    /// <remarks>The original matrix is not altered by this method; all row
    ///  operations are performed on a copy.</remarks>
    function ReducedRowEchelonForm(CollectRowOps: Boolean = False): TRealMatrix;

    /// <summary>Returns the number of zero rows in the matrix.</summary>
    function NumZeroRows(const AEpsilon: TASR = 0): Integer;

    /// <summary>Returns the number of zero rows at the bottom of the matrix.
    ///  </summary>
    function NumTrailingZeroRows(const AEpsilon: TASR = 0): Integer;

    /// <summary>Performs the Gram-Schmidt orthonormalization procedure on the
    ///  columns of the matrix; the result is returned.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function GramSchmidt: TRealMatrix;

    procedure InplaceGramSchmidt(FirstCol, LastCol: Integer);

    /// <summary>Returns the matrix with zero or more columns removed so that
    ///  the remaining columns form a basis for the column space (range) of the
    ///  matrix.</summary>
    function ColumnSpaceBasis: TRealMatrix;

    function ColumnSpaceProjection(const AVector: TRealVector): TRealVector;

    /// <summary>Returns the Euclidean distance between a vector and the column
    ///  space of the matrix.</summary>
    function DistanceFromColumnSpace(const AVector: TRealVector): TASR; // empty matrix aware

    /// <summary>Computes and returns an upper Hessenberg matrix that is similar
    ///  to the original matrix, which must be square.</summary>
    /// <param name="A2x2Bulge">If true, the original matrix MUST be a PROPER
    ///  upper Hessenberg matrix except for a 2×2 bulge at the upper-left
    ///  position. In this case, a different algorithm suitable for this
    ///  particular case might be used, which is faster than the general-case
    ///  algorithm.</param>
    /// <remarks>If the original matrix is symmetric, the result is tridiagonal.
    ///  </remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function SimilarHessenberg(A2x2Bulge: Boolean = False): TRealMatrix;

    /// <summary>Returns the array of eigenvalues of the matrix.</summary>
    /// <returns>A complex vector of the n eigenvalues.</returns>
    /// <exception cref="EMathException">Raised if the eigenvalues couldn't be
    ///  computed, for any reason.</exception>
    function UnsortedEigenvalues: TComplexVector;

    function eigenvalues: TComplexVector;

    /// <summary>An alias for <c>eigenvalues</c>.</summary>
    function spectrum: TComplexVector; inline; // alias

    /// <summary>Computes the eigenvalues and eigenvectors of the matrix.
    ///  </summary>
    /// <param name="AEigenvalues">Receives the eigenvalues.</param>
    /// <param name="AEigenvectors">Receives a real n×n matrix, the columns
    ///  of which are linearly independent eigenvectors of the matrix in the same
    ///  order as the eigenvalues.</param>
    /// <param name="ASkipVerify">If false (the default), the method verifies
    ///  that the obtained vectors really are linearly independent eigenvectors
    ///  of the matrix, in the same order as the eigenvalues. If true, this
    ///  verification is skipped. In any case, the eigenvalues are verified.
    ///  </param>
    /// <returns>True iff the eigenvalues and eigenvectors were successfully
    ///  computed and all checks complete successfully.</returns>
    /// <remarks>Currently only implemented for symmetric matrices.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square
    ///  and symmetric.</exception>
    function eigenvectors(out AEigenvalues: TRealVector;
      out AEigenvectors: TRealMatrix; ASkipVerify: Boolean = False): Boolean;

    /// <summary>Returns true iff the given vector is an eigenvector of the
    ///  square matrix.</summary>
    /// <exception cref="EMathException">Raised if the matrix isn't square or
    ///  if the given vector is of a different dimension.</exception>
    function IsEigenvector(const u: TRealVector;
      const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns the associated eigenvalue of a given eigenvector.
    ///  </summary>
    /// <exception cref="EMathException">Raised if the given vector isn't an
    ///  eigenvector of the matrix. In particular, this happens if the matrix
    ///  isn't square or if the vector is the zero vector of its dimension or
    //// of a different dimension compared to the matrix.</exception>
    function EigenvalueOf(const u: TRealVector;
      const Epsilon: Extended = 0): TASR;

    /// <summary>Tries to find the associated eigenvalue of a given vector and
    ///  returns whether the given vector was an eigenvector or not.</summary>
    /// <param name="u">The given vector that might or might not be an eigenvector
    ///  of the matrix.</param>
    /// <param name="AEigenvalue">Receives the associated eigenvalue, if the
    ///  given vector turned out to be an eigenvector.</param>
    /// <returns>True iff the given vector is an eigenvector of the matrix.
    ///  </returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square
    ///  or if the dimension of the given vector doesn't match the size of the
    ///  matrix.</exception>
    function TryEigenvalueOf(const u: TRealVector; out AEigenvalue: TASR;
      const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether (lambda, u) is an eigenpair of the matrix.
    ///  </summary>
    function IsEigenpair(const lambda: TASR; const u: TRealVector;
      const Epsilon: Extended = 0): Boolean;

    /// <summary>Computes an eigenvector associated with a given eigenvalue.
    ///  </summary>
    /// <param name="lambda">The eigenvalue to find an associated eigenvalue for.
    ///  </param>
    /// <returns>A normalized eigenvector associated with the given eigenvalue.
    ///  </returns>
    /// <exception cref="EMathException">Raised if <c>lambda</c> isn't an
    ///  eigenvalue of the matrix. In particular, this happens if the matrix
    ///  isn't square. Also raised if the algorithm fails to compute an
    ///  eigenvector, for any reason.</exception>
    function EigenvectorOf(const lambda: TASR): TRealVector;

    /// <summary>Returns the spectral radius of the matrix, that is, the largest
    ///  absolute value of the eigenvalues.</summary>
    /// <remarks>The current implementation only works if all eigenvalues are
    ///  real; otherwise, an exception is raised.</remarks>
    function SpectralRadius: TASR; inline;

    /// <summary>Returns the array of singular values of the matrix.</summary>
    /// <remarks><para>The singular values are the eigenvalues of <c>sqrt(A* A)</c>,
    ///  or, equivalently, the square roots of the eigenvalues of the Hermitian
    ///  square <c>A* A</c>.</para>
    ///  <para>The singular values are listed in decreasing order.</para></remarks>
    function SingularValues: TRealVector;

    /// <summary>Returns a copy of the matrix with each element replaced with
    ///  its absolute value.</summary>
    function Abs: TRealMatrix;

    /// <summary>Replaces all elements that are extremely close to zero with
    ///  zero.</summary>
    /// <returns>A reference to the matrix.</returns>
    /// <remarks>The matrix is altered by this method; the result only points
    ///  to the original matrix that has been modified.</remarks>
    function Defuzz(const Eps: Double = 1E-8): TRealMatrix;

    /// <summary>Creates a new matrix that is identical to this matrix.</summary>
    function Clone: TRealMatrix; // empty matrix aware

    /// <summary>Returns the elements of the matrix as a real vector. The elements
    ///  are obtained from the matrix in column-major order. In particular, this
    ///  works well if the matrix is either a single row or a single column.
    ///  </summary>
    /// <see cref="AsVector" />
    function Vectorization: TRealVector; inline;

    /// <summary>Returns the elements of the matrix as a real vector. The elements
    ///  are obtained from the matrix in row-major order. In particular, this
    ///  works well if the matrix is either a single row or a single column.
    ///  </summary>
    /// <see cref="Vectorization" />
    /// <see cref="Data" />
    function AsVector: TRealVector; inline;

    /// <summary>Augments the matrix with a new matrix.</summary>
    /// <param name="A">The real matrix to augment the current matrix with. This
    ///  matrix must have the same number of rows as the current matrix.</param>
    /// <returns>A real matrix with A augmented to the right of the current
    ///  matrix.</return>
    /// <remarks>The current matrix is not modified by the method. A vector
    ///  may be used as the argument, since it will be automatically promoted to
    ///  the corresponding single-column matrix.</remarks>
    /// <see cref="AddRow" />
    function Augment(const A: TRealMatrix): TRealMatrix; overload; // empty matrix aware

    /// <summary>Augments the matrix with a zero column.</summary>
    /// <returns>A real matrix with a zero column augmented to the right of the
    ///  current matrix.</return>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Augment: TRealMatrix; overload; inline;

    /// <summary>Replaces a specific row with a new set of values.</summary>
    procedure SetRow(ARowIndex: Integer; const ARow: array of TASR);

    /// <summary>Replaces a specific column with a new set of values.</summary>
    procedure SetCol(AColIndex: Integer; const ACol: array of TASR);

    /// <summary>Returns the (m-1)×(n-1) submatrix obtained by removing a single
    ///  row and a single column from the matrix.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Submatrix(ARowToRemove: Integer = 0;
      AColToRemove: Integer = 0; AAllowEmpty: Boolean = False): TRealMatrix; overload; // empty matrix aware

    /// <summary>Returns a submatrix of the current matrix using the specified
    ///  rows and columns.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Submatrix(const ARows: array of Integer;
      const ACols: array of Integer): TRealMatrix; overload;

    /// <summary>Returns a principal submatrix of the current matrix using the
    ///  specified rows and columns.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Submatrix(const ARows: array of Integer): TRealMatrix; overload;

    /// <summary>Returns a leading principal submatrix of the matrix.</summary>
    function LeadingPrincipalSubmatrix(const ASize: Integer): TRealMatrix;

    /// <summary>Returns the submatrix obtained by removing the last column.
    ///  </summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Lessened: TRealMatrix;

    /// <summary>Returns the determinant of the submatrix obtained by removing
    ///  one specific row and one specific column.</summary>
    function Minor(ARow, ACol: Integer): TASR; inline;

    /// <summary>Returns (-1)^(ARow + ACol) times the corresponding minor.
    ///  </summary>
    function Cofactor(ARow, ACol: Integer): TASR; inline;

    /// <summary>Returns the cofactor matrix of the matrix.</summary>
    function CofactorMatrix: TRealMatrix;

    /// <summary>Returns the adjugate matrix of the matrix, that is, the
    ///  transpose of the cofactor matrix.</summary>
    function AdjugateMatrix: TRealMatrix; inline;

    /// <summary>Returns a copy of the matrix with each element replaced by
    ///  the result of a real-valued function applied to it.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Apply(AFunction: TRealFunctionRef): TRealMatrix;

    /// <summary>Replaces each entry of the matrix satisfying the given predicate
    ///  with a given value.</summary>
    /// <param name="APredicate">A predicate for real numbers. Each matrix entry
    ///  satisfying this predicate will be replaced with <c>ANewValue</c>.</param>
    /// <param name="ANewValue">The value to replace entries with.</param>
    /// <returns>A copy of the current matrix but with all entries satisfying
    ///  <c>APredicate</c> replaced with <c>ANewValue</c>.</returns>
    /// <remarks>The current matrix is not modified by this method.</remarks>
    function Replace(APredicate: TPredicate<TASR>;
      const ANewValue: TASR): TRealMatrix; overload;

    /// <summary>Replaces each entry with a specified value with a new value.
    ///  </summary>
    /// <param name="AOldValue">The value to search for in the original matrix.
    ///  </param>
    /// <param name="ANewValue">The value to replace with.</param>
    /// <returns>A copy of the current matrix but with each entry originally
    ///  eaual to <c>AOldValue</c> replaced with <c>ANewValue</c>.</returns>
    /// <remarks>The current matrix is not modified by this method.</remarks>
    function Replace(const AOldValue, ANewValue: TASR;
      const Epsilon: Extended = 0): TRealMatrix; overload;

    /// <summary>Sets each entry of the matrix to the same, new, value.</summary>
    /// <param name="ANewValue">The new value to set each entry to.</param>
    /// <returns>A matrix of the same size as the original matrix but with each
    ///  entry equal to <c>ANewValue</c>.</returns>
    /// <remarks>The current matrix is not modified by this method.</remarks>
    function Replace(const ANewValue: TASR): TRealMatrix; overload; {copy}

    /// <summary>Returns either a single-line textual representation of the
    /// matrix of the form ((M00, M01, ...), (M10, M11, ...), ...) or a multi-
    /// line tabular representation.</summary>
    function str(const AOptions: TFormatOptions): string; // mainly for debugging

    /// <summary>Adds a row to the matrix.</summary>
    /// <remarks>This is more of a computing-related function than a mathematical
    ///  operation. Notice that it does suffer from the
    ///  <c>SetLength(X, Length(X) + 1)</c> problem.</remarks>
    /// <see cref="Augment" />
    procedure AddRow(const ARow: array of TASR);

    /// <summary>Sorts the elemenets of the matrix.</summary>
    /// <remarks>This operation sorts the elements of the original matrix,
    ///  and returns a reference to it.</remarks>
    /// <returns>A reference to the original matrix, which has been sorted.
    ///  </returns>
    function Sort(AComparer: IComparer<TASR>): TRealMatrix;

    /// <summary>Shuffles the elemenets of the matrix.</summary>
    /// <remarks>This operation shuffles the elements of the original matrix,
    ///  and returns a reference to it.</remarks>
    /// <returns>A reference to the original matrix, which has been shuffled.
    ///  </returns>
    function Shuffle: TRealMatrix;

    /// <summary>Reverses the order of the elements in the matrix (in the
    /// row-major sense).</summary>
    /// <remarks>This operation reverses the elements of the original matrix,
    ///  and returns a reference to it.</remarks>
    /// <returns>A reference to the original matrix, which has been reversed.
    ///  </returns>
    function Reverse: TRealMatrix;

    // ***
    // Decompositions
    // ***

    /// <summary>Produces a LU decomposition of the matrix, which has to be square,
    ///  that is, finds a permutation matrix P, a unit lower triangular matrix L,
    ///  and an upper triangular matrix U such that Self = P.Transpose * L * U.
    ///  </summary>
    /// <param name="P">Receives the permutation matrix P.</param>
    /// <param name="L">Receives the unit lower triangular matrix L.</param>
    /// <param name="U">Receives the upper triangular matrix U.</param>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    procedure LU(out P, L, U: TRealMatrix);

    /// <summary>Computes the Cholesky decomposition of the square matrix, if
    ///  possible.</summary>
    /// <param name="R">Receives the upper triangular Cholesky factor R such
    ///  that Self = R.Transpose * R.</param>
    /// <returns>True iff the decomposition was successfully computed; this
    ///  happens precisely when the matrix is positive definite.</returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function Cholesky(out R: TRealMatrix): Boolean;

    /// <summary>Produces an orthogonal matrix Q and an upper triangular matrix
    ///  R such that Self = Q * R.</summary>
    /// <remarks>QR decomposition is only implemented for square and tall matrices.
    ///  </remarks>
    /// <exception cref="EMathException">Raised if the matrix is wide, that is,
    ///  if <c>Size.Cols &gt; Size.Rows</c>.</exception>
    procedure QR(out Q, R: TRealMatrix);

    /// <summary>Finds an upper Hessenberg matrix A and an orthogonal matrix
    ///  U such that Self = U A U.Transpose.</summary>
    procedure Hessenberg(out A, U: TRealMatrix);

  end;

/// <summary>Computes and returns the Nth power of the square matrix A.</summary>
/// <remarks>The zeroth power is the identity matrix of the same size as A.
///  Negative powers are defined iff A is invertible, and if so the result
///  is the |N|th power of the inverse of A (which is equal to the inverse
///  of the |N|th power of A).</remarks>
function mpow(const A: TRealMatrix; const N: Integer): TRealMatrix; overload;

/// <summary>Returns the unique positive semidefinite square root of a given
///  positive semidefinite symmetric matrix.</summary>
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;

/// <summary>Returns the zero matrix of a given size.</summary>
function ZeroMatrix(const ASize: TMatrixSize): TRealMatrix; inline;

/// <summary>Returns the identity matrix of a given size.</summary>
function IdentityMatrix(ASize: Integer): TRealMatrix; // empty matrix aware

/// <summary>Returns the reversal matrix of a given size.</summary>
function ReversalMatrix(ASize: Integer): TRealMatrix;

/// <summary>Returns a random matrix of a given size. The elements are real
///  numbers in [0, 1).</summary>
function RandomMatrix(const ASize: TMatrixSize): TRealMatrix;

/// <summary>Returns a random matrix of a given size. The elements are integers
///  in [a, b).</summary>
function RandomIntMatrix(const ASize: TMatrixSize; a, b: Integer): TRealMatrix;

/// <summary>Returns a diagonal matrix constructed using values from
///  a given array of real numbers.</summary>
/// <see cref="DirectSum" />
function diag(const AElements: array of TASR): TRealMatrix; overload;

/// <summary>Returns a diagonal matrix constructed using values from
///  a given real vector.</summary>
/// <see cref="DirectSum" />
function diag(const AElements: TRealVector): TRealMatrix; overload; inline;

/// <summary>Returns the outer product between two vectors.</summary>
function OuterProduct(const u, v: TRealVector): TRealMatrix; overload; inline;

/// <summary>Returns the n×n circulant matrix whose first row is a given vector.
///  </summary>
function CirculantMatrix(const AElements: array of TASR): TRealMatrix; overload;

/// <summary>Returns the Toeplitz matrix with a given first row and a given first
///  column.</summary>
function ToeplitzMatrix(const AFirstRow,
  AFirstCol: array of TASR): TRealMatrix; overload;

/// <summary>Returns the Hankel matrix with a given first row and a given last
///  column.</summary>
function HankelMatrix(const AFirstRow,
  ALastCol: array of TASR): TRealMatrix; overload;

/// <summary>Returns the backward shift matrix of a given size.</summary>
function BackwardShiftMatrix(ASize: Integer): TRealMatrix;

/// <summary>Returns the forward shift matrix of a given size.</summary>
function ForwardShiftMatrix(ASize: Integer): TRealMatrix;

/// <summary>Returns a Vandermonde matrix produced by a given vector.</summary>
/// <param name="AElements">The vector used as the second column in the matrix.
///  </param>
/// <param name="ACols">The number of columns in the matrix. If omitted,
///  the matrix will be square.</param>
function VandermondeMatrix(const AElements: array of TASR;
  ACols: Integer = 0): TRealMatrix; overload;

/// <summary>Returns a Hilbert matrix of a specific size, which need not be
///  square. Notice that the argument may be an integer, in which case a square
///  matrix is produced.</summary>
function HilbertMatrix(const ASize: TMatrixSize): TRealMatrix;

/// <summary>Returns the matrix for a rotation in a given coordinate plane in
///  some dimension.</summary>
/// <param name="AAngle">The angle of the rotation (in radians).</param>
/// <param name="ADim">The dimension of the vector space; the size of the
///  matrix.</param>
/// <param name="AIndex1">The index of the first coordinate that varies in the
///  plane in which the rotation takes place.</param>
/// <param name="AIndex2">The index of the second coordinate that varies in the
///  plane in which the rotation takes place.</param>
/// <remarks>A rotation matrix like this is also called a Givens rotation.</remarks>
function RotationMatrix(AAngle: TASR; ADim: Integer = 2; AIndex1: Integer = 0;
  AIndex2: Integer = 1): TRealMatrix; overload;

function RotationMatrix(AAngle: TASR; AAxis: TRealVector): TRealMatrix; overload;

/// <summary>Returns the matrix for a reflection in a hyperplane in some dimension.
///  </summary>
/// <param name="u">A vector perpendicular to the hyperplane.</param>
/// <remarks>A reflection matrix like this is also called a Householder matrix.
///  </remarks>
function ReflectionMatrix(const u: TRealVector): TRealMatrix; overload;

/// <summary>Returns the matrix for a reflection in a hyperplane in some dimension.
///  </summary>
/// <param name="u">A unit vector perpendicular to the hyperplane.</param>
/// <remarks>A reflection matrix like this is also called a Householder matrix.
///  The vector u MUST be a unit vector.</remarks>
function QuickReflectionMatrix(const u: TRealVector): TRealMatrix; overload; inline;

/// <summary>Returns the Hadamard product between two matrices of the same
///  size, that is, the matrix obtained by taking element-wise products.</summary>
function HadamardProduct(const A, B: TRealMatrix): TRealMatrix; overload;

/// <summary>Returns the direct sum of two matrices.</summary>
/// <see cref="diag" />
function DirectSum(const A, B: TRealMatrix): TRealMatrix; overload; inline;

/// <summary>Returns the direct sum of an open array of matrices.</summary>
/// <see cref="diag" />
function DirectSum(const Blocks: array of TRealMatrix): TRealMatrix; overload;

/// <summary>Determines whether two matrices commute or not.</summary>
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;

/// <summary>Returns the sum of all elements of the matrix.</summary>
function sum(const A: TRealMatrix): TASR; overload; inline;

/// <summary>Returns the arithmetic mean of all elements of the matrix.</summary>
function ArithmeticMean(const A: TRealMatrix): TASR; overload; inline;

/// <summary>Returns the geometric mean of all elements of the matrix.</summary>
function GeometricMean(const A: TRealMatrix): TASR; overload; inline;

/// <summary>Returns the harmonic mean of all elements of the matrix.</summary>
function HarmonicMean(const A: TRealMatrix): TASR; overload; inline;

/// <summary>Returns the product of all elements of the matrix.</summary>
function product(const A: TRealMatrix): TASR; overload; inline;

/// <summary>Returns the largest value in the matrix.</summary>
function max(const A: TRealMatrix): TASR; overload; inline;

/// <summary>Returns the smallest value in the matrix.</summary>
function min(const A: TRealMatrix): TASR; overload; inline;

/// <summary>Returns true iff there exists an entry of the matrix satisfying
///  the given predicate.</summary>
function exists(const A: TRealMatrix; APredicate: TPredicate<TASR>): Boolean; overload;

/// <summary>Returns the number of matrix entries that satisfy the given
///  predicate.</summary>
function count(const A: TRealMatrix; APredicate: TPredicate<TASR>): Integer; overload;

/// <summary>Returns the number of instances of a given value in the matrix.
///  </summary>
function count(const A: TRealMatrix; const AValue: TASR): Integer; overload; inline;

/// <summary>Returns true iff all entries of the matrix satisfy
///  the given predicate.</summary>
function ForAll(const A: TRealMatrix; APredicate: TPredicate<TASR>): Boolean; overload;

/// <summary>Returns true iff a specific value is found in the matrix.</summary>
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;

/// <summary>Solves the matrix equation AX = Y where A is a lower triangular
///  n×n matrix and X and Y are n-dimensional vectors using forward substitution.
///  </summary>
/// <param name="A">A square matrix of size n, which MUST be lower triangular.</param>
/// <param name="Y">An n-dimensional vector.</param>
/// <param name="IsUnitDiagonal">If true, the routine assumes that the main
///  diagonal of <c>A</c> consists of only 1's, without caring about the actual
///  elements on the main diagonal.</param>
/// <returns>The unique solution X of the equation AX = Y, if A is non-singular.</returns>
/// <exception cref="EMathException">Raised if A is not square and non-singular.
///  Also raised if the size of A differs from the dimension of Y.</exception>
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;

/// <summary>Solves the matrix equation AX = Y where A is an upper triangular
///  n×n matrix and X and Y are n-dimensional vectors using backward
///  substitution.</summary>
/// <param name="A">A square matrix of size n, which MUST be upper triangular.</param>
/// <param name="Y">An n-dimensional vector.</param>
/// <returns>The unique solution X of the equation AX = Y, if A is non-singular.</returns>
/// <exception cref="EMathException">Raised if A is not square and non-singular.
///  Also raised if the size of A differs from the dimension of Y.</exception>
function BackSubstitution(const A: TRealMatrix; const Y: TRealVector): TRealVector; overload;

/// <summary>Solves the matrix equation AX = Y.</summary>
/// <param name="AAugmented">The augmented matrix [A|Y].</param>
/// <returns>The unique solution X if a unique solution exists.</returns>
/// <exception cref="EMathException">Raised if a unique solution does not
///  exist.</exception>
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;

/// <summary>Solves the matrix equation AX = Y.</summary>
/// <param name="A">The matrix of coefficients, A.</param>
/// <param name="Y">The right-hand side of the equation, Y.</param>
/// <returns>The unique solution X if a unique solution exists.</returns>
/// <exception cref="EMathException">Raised if a unique solution does not
///  exist.</exception>
function SysSolve(const A: TRealMatrix; const Y: TRealVector): TRealVector; overload;
function SysSolve(const A: TRealMatrix; const Y: TRealMatrix): TRealMatrix; overload;

/// <summary>Given a set { (x_0, y_0), (x_1, y_1), ..., (x_(k-1), y_(k-1)) } of
///  points, this function determines the polynomial of degree at most n
///  whose graph best connects to the points in the least-squares sense.</summary>
/// <param name="X">The k-dimensional vector of X coordinates.</param>
/// <param name="Y">The k-dimensional vector of Y coordinates.</param>
/// <param name="ADegree">The maximum degree n >= 0 of the polynomial.</param>
/// <returns>A real vector consisting of the n + 1 coefficients c_0, c_1, ..., c_n
///  of the optimal polynomial c_0 + c_1 x + c_2 x^2 + ... + c_n x^n.</returns>
function LeastSquaresPolynomialFit(const X, Y: TRealVector;
  ADegree: Integer): TRealVector; overload;

/// <summary>A quick but unsafe procedure to copy data from one vector
///  to another.</summary>
/// <param name="ASource">The source vector.</param>
/// <param name="AFrom">The first component of the source vector to be copied.
///  </param>
/// <param name="ATo">The last component of the source vector to be copied.
///  </param>
/// <param name="ATarget">The vector to copy to.</param>
/// <param name="ATargetFrom">The index of the first component of the target
///  vector to be overwritten by the copied data.</param>
/// <remarks>The specified components in the source vector MUST exist, and the
///  target vector MUST be large enough to contain all data.</remarks>
procedure VectMove(const ASource: TRealVector; const AFrom, ATo: Integer;
  var ATarget: TRealVector; const ATargetFrom: Integer = 0); overload; inline;

/// <summary>A quick but unsafe procedure to copy data from a vector to a part
///  of a column of a matrix.</summary>
/// <param name="ASource">The source vector.</param>
/// <param name="AFrom">The first component of the source vector to be copied.
///  </param>
/// <param name="ATo">The last component of the source vector to be copied.
///  </param>
/// <param name="ATarget">The matrix to write data to.</param>
/// <param name="ATargetCol">The index of the column of the matrix to write to.
///  </param>
/// <param name="ATargetFirstRow">The index of the first row in the matrix
///  column to overwrite.</param>
/// <remarks>The specified components in the source vector MUST exist, and the
///  target matrix MUST be large enough to contain all data.</remarks>
procedure VectMoveToMatCol(const ASource: TRealVector; const AFrom, ATo: Integer;
  var ATarget: TRealMatrix; const ATargetCol: Integer;
  const ATargetFirstRow: Integer = 0); overload;

/// <summary>A quick but unsafe procedure to copy data from a vector to a part
///  of a row of a matrix.</summary>
/// <param name="ASource">The source vector.</param>
/// <param name="AFrom">The first component of the source vector to be copied.
///  </param>
/// <param name="ATo">The last component of the source vector to be copied.
///  </param>
/// <param name="ATarget">The matrix to write data to.</param>
/// <param name="ATargetRow">The index of the row of the matrix to write to.
///  </param>
/// <param name="ATargetFirstCol">The index of the first column in the matrix
///  row to overwrite.</param>
/// <remarks>The specified components in the source vector MUST exist, and the
///  target matrix MUST be large enough to contain all data.</remarks>
procedure VectMoveToMatRow(const ASource: TRealVector; const AFrom, ATo: Integer;
  var ATarget: TRealMatrix; const ATargetRow: Integer;
  const ATargetFirstCol: Integer = 0); overload; inline;

/// <summary>A quick but unsafe procedure to copy data from one matrix to
///  another.</summary>
/// <param name="ASource">The matrix to copy from.</param>
/// <param name="ARect">A rectangle describing the submatrix to copy; both
///  the top-left and the bottom-right elements are included in the submatrix.
///  </param>
/// <param name="ATarget">The matrix to copy to.</param>
/// <param name="ATargetTopLeft">The top-left element in the destination matrix
///  to be overwritten.</param>
/// <remarks>The specified elements in the source matrix MUST exist, and the
///  target matrix MUST be large enough to contain all data.</remarks>
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;

/// <summary>A quick but unsafe procedure to copy data from a column in a matrix
///  to a vector.</summary>
/// <param name="ASource">The matrix to copy data from.</param>
/// <param name="AColumn">The column of the matrix to copy data from.</param>
/// <param name="AFrom">The row index of the first element in the column to
///  copy.</param>
/// <param name="ATo">The row index of the last element in the column to
///  copy.</param>
/// <param name="ATarget">The vector to copy data into.</param>
/// <param name="ATargetFrom">The index of the first component of the vector
///  to overwrite.</param>
/// <remarks>The specified elements in the source matrix MUST exist, and the
///  target vector MUST be large enough to contain all data.</remarks>
procedure MatMoveColToVect(const ASource: TRealMatrix; const AColumn: Integer;
  const AFrom, ATo: Integer; var ATarget: TRealVector;
  const ATargetFrom: Integer = 0); overload;

/// <summary>A quick but unsafe procedure to copy data from a row in a matrix
///  to a vector.</summary>
/// <param name="ASource">The matrix to copy data from.</param>
/// <param name="ARow">The row of the matrix to copy data from.</param>
/// <param name="AFrom">The column index of the first element in the row to
///  copy.</param>
/// <param name="ATo">The column index of the last element in the row to
///  copy.</param>
/// <param name="ATarget">The vector to copy data into.</param>
/// <param name="ATargetFrom">The index of the first component of the vector
///  to overwrite.</param>
/// <remarks>The specified elements in the source matrix MUST exist, and the
///  target vector MUST be large enough to contain all data.</remarks>
procedure MatMoveRowToVect(const ASource: TRealMatrix; const ARow: Integer;
  const AFrom, ATo: Integer; var ATarget: TRealVector;
  const ATargetFrom: Integer = 0); overload; inline;

/// <summary>A quick but unsafe procedure to perform an inplace multiplication
///  of a submatrix with a square matrix from the left.</summary>
/// <param name="ATarget">The matrix to modify.</param>
/// <param name="ARect">A rectangle that identifies the submatrix to multiply.
///  Both the upper-left and the bottom-right corners belong to the submatrix.
///  The submatrix has format n × m. This parameter MUST specify a valid
///  submatrix.</param>
/// <param name="AFactor">The matrix to multiply the submatrix by, from the left.
///  This MUST be a matrix of size n × n.</param>
procedure MatMulBlockInplaceL(var ATarget: TRealMatrix; const ARect: TRect;
  const AFactor: TRealMatrix); overload;

/// <summary>A quick but unsafe procedure to perform an inplace multiplication
///  of a square submatrix with a square matrix from the left.</summary>
/// <param name="ATarget">The matrix to modify.</param>
/// <param name="APoint">The index of the upper-left element of the submatrix.</param>
/// <param name="AFactor">The matrix to multiply the submatrix by, from the left.
///  This MUST be a matrix of size n × n.</param>
/// <remarks>The implied n × n submatrix MUST be valid.</remarks>
procedure MatMulBlockInplaceL(var ATarget: TRealMatrix; const ATopLeft: TPoint;
  const AFactor: TRealMatrix); overload; inline;

/// <summary>
///  <para>Performs the operation</para>
///  <para><c>ATarget := ATarget * DirectSum(IdentityMatrix(k), AFactor)</c></para>
///  <para>where <c>k = n - m</c>, <c>n</c> is the size of the square matrix
///   <c>ATarget</c>, and <c>m <= n</c> is the size of the square matrix
///   <c>AFactor</c>.</para>
/// </summary>
/// <param name="ATarget">An n×n matrix. This MUST be square.</param>
/// <param name="AFactor">An m×m matrix which MUST be square and where m
///  MUST be <= n.</param>
procedure MatMulBottomRight(var ATarget: TRealMatrix; const AFactor: TRealMatrix); overload;

/// <summary>A quick but unsafe procedure to fill a submatrix with a given
///  real number.</summary>
/// <param name="ATarget">The matrix to modify.</param>
/// <param name="ARect">A rectangle describing the submatrix to fill; both
///  the top-left and the bottom-right elements of the rectangle belong to the
///  submatrix.</param>
/// <param name="Value">The real number to fill the submatrix with.</param>
/// <remarks>The specified submatrix MUST belong to the matrix.</remarks>
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;

  /// <summary>The AlgoSim complex matrix type.</summary>
  TComplexMatrix = record
  strict private
    const
      InitialRowOpSeqSize = 16; // strategy: then double
    var
      FSize: TMatrixSize;
      FElements: TASCArray;
      FRowOpSeq: TComplexRowOpSequence;
      FRowOpCount: Integer;
      FRowOpFactor: TASC;
      _FCollectRowOpData: Boolean; // never set directly -- use BeginCollectRowOpData and EndCollectRowOpData

    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;

    /// <summary>Clears any existing row op data, initializes the row op buffers
    /// and sets the collect row op data flag to true, so that row op data will
    /// be collected after this call.</summary>
    procedure BeginCollectRowOpData;

    /// <summary>Clears any existing row op data and sets the collect row op data
    ///  flag to false.</summary>
//    procedure EndCollectRowOpData;

    /// <summary>Adds a row operation record to the buffer.</summary>
    /// <remarks>Warning! This method MUST NOT be called if FCollectRowOpData
    ///  is set to false. Doing so will result in undefined behaviour.</remarks>
    procedure AddRowOpRecord(AType: TRowOperationType; ARow1, ARow2: Integer;
      AFactor: TASC);

    procedure GetHouseholderMap(const AVect: TComplexVector;
      out tau, gamma: TASC; out u: TComplexVector);

    /// <summary>Returns the two complex eigenvalues of the second leading
    ///  principal submatrix of the matrix.</summary>
    /// <remarks>The matrix MUST has at least two columns and it MUST have
    ///  at least two rows.</remarks>
    function Eigenvalues2x2: TComplexVector;

  private
    /// <summary>Creates a new complex matrix record but does not allocate data
    ///  for the matrix entries. The <c>Data</c> property is set to <c>nil</c>.
    ///  </summary>
//    constructor CreateWithoutAllocation(const ASize: TMatrixSize);

    /// <summary>Returns true iff the matrix is an empty matrix, that is, if
    ///  it has exactly zero elements.</summary>
    function IsEmpty: Boolean; inline;

    /// <summary>Performs a fast LU decomposition of the matrix, which has to be
    ///  an n×n square.</summary>
    /// <param name="A">An n×n matrix that receives the elements of L and U.
    ///  The upper triangular part, including the main diagonal, receives the
    ///  elments of U, and the lower triangular part, excluding the diagonal,
    ///  receives the elements of L that are below the main diagonal. The main
    ///  diagonal elements of L are all equal to 1, so A contains all
    ///  information about L and U.</param>
    /// <param name="P">An n-dimensional array of integers that contain all
    ///  information about the row exchanges made during the LU decomposition.
    ///  If the kth element of the array is m, then at the kth step, rows k and
    ///  m (>= k) were exchanged. If m = k, no row exchange was made at this
    ///  step. If m = -1, no row exchange was made at this step, and, in
    ///  addition, it was found that the matrix was singular at this step.
    ///  </param>
    /// <remarks>This method is as fast as possible. The public method
    ///  LU uses this method to construct the actual matrixes L, U, and P in
    ///  the LU decomposition Self = P.Transpose * L * U.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    procedure DoQuickLU(out A: TComplexMatrix; out P: TIndexArray);

    procedure RequireNonEmpty; inline;

  public

    // ***
    // Constructors
    // ***

    /// <summary>Allocates memory for a new complex matrix of a given size.
    ///  Equivalently, it creates a new complex matrix of a given size with
    ///  elements with undefined initial values.</summary>
    /// <remarks>Creating an uninitialized matrix of a given size is faster
    ///  compared to creating an initialized matrix of the same size.</remarks>
    constructor CreateUninitialized(const ASize: TMatrixSize);

    /// <summary>Copy constructor for <c>TComplexMatrix</c>.</summary>
    constructor Create(const AMatrix: TComplexMatrix); overload;

    /// <summary>Creates a new complex matrix from an array of complex-number
    ///  arrays. Each complex-number array will form a row in the matrix.</summary>
    /// <param name="Elements">An open array of complex-number arrays that will
    ///  be used to create the matrix.</param>
    /// <remarks>The complex-number arrays have to have the same length; this
    ///  common length becomes the number of columns of the matrix.</remarks>
    /// <see cref="CreateFromRows" />
    /// <see cref="CreateFromColumns" />
    constructor Create(const Elements: array of TASCArray); overload;

    /// <summary>Creates a new complex matrix from a single array of complex
    ///  numbers and a given matrix format; the complex numbers will populate
    ///  the matrix top to bottom and left to right. The complex number with
    ///  index i + 1 will be to the right of the element with index i unless the
    ///  latter is the right-most element of a row, in which case element i + 1
    ///  will become the first element of the next row.</summary>
    /// <param name="Elements">An array of complex numbers.</param>
    /// <param name="Cols">The number of columns in the matrix.</param>
    /// <remarks>The number of rows in the matrix will be the length of
    ///  <c>Elements</c> divided by <c>Cols</c>, which has to be an integer. The
    ///  default value of <c>Cols</c> is 1, which will create a single-column
    ///  matrix (a column vector) with the length of <c>Elements</c>.</remarks>
    constructor Create(const Elements: array of TASC; Cols: Integer = 1); overload;

    /// <summary>Creates a new complex matrix with a given format and fills it
    ///  with a given complex number.</summary>
    /// <param name="ASize">The format of the matrix.</param>
    /// <param name="AVal">The complex number to fill the matrix with.</param>
    /// <remarks>The default value of <c>AVal</c> is zero, which will create the
    ///  zero matrix of the given size.</remarks>
    constructor Create(const ASize: TMatrixSize; const AVal: TASC); overload;

    /// <summary>Creates a new complex matrix from a sequence of complex vectors,
    ///  which will form the rows of the matrix.</summary>
    /// <param name="Rows">An open array of complex vectors; each complex vector will
    ///  form a row in the matrix. The number of complex vectors will become the
    ///  number of rows in the matrix.</param>
    /// <remarks>The vectors need to have the same length, which will become the
    ///  number of columns of the matrix.</remarks>
    constructor CreateFromRows(const Rows: array of TComplexVector);

    /// <summary>Creates a new complex matrix from a sequence of complex vectors,
    ///  which will form the columns of the matrix.</summary>
    /// <param name="Columns">An open array of complex vectors; each complex vector will
    ///  form a column in the matrix. The number of complex vectors will become the
    ///  number of columns in the matrix.</param>
    /// <remarks>The vectors need to have the same length, which will become the
    ///  number of rows of the matrix.</remarks>
    constructor CreateFromColumns(const Columns: array of TComplexVector);

    /// <summary>Creates a new complex matrix A from two complex vectors u and v;
    ///  the matrix is the outer product of the vectors, that is, A_ij = u_i * v_j*.
    ///  </summary>
    /// <param name="u">The first complex vector.</param>
    /// <param name="v">The second complex vector.</param>
    /// <remarks>The vectors need not be of the same dimension.</remarks>
    constructor Create(const u, v: TComplexVector); overload;

    /// <summary>Creates a new complex matrix using given format and a
    ///  complex-valued function of valid matrix indices.</summary>
    /// <param name="ASize">The format of the matrix.</param>
    /// <param name="AFunction">A reference to a complex-valued function of
    ///  matrix indices (integers). This function is used to populate the matrix.
    ///  The element (i, j) in the matrix is the value AFunction(i, j). Here,
    ///  row and column indices are 1-based.</param>
    constructor Create(const ASize: TMatrixSize;
      AFunction: TMatrixIndexFunction<TASC>); overload;

    /// <summary>Creates a new diagonal complex matrix using a given array of
    ///  complex numbers.</summary>
    /// <param name="Elements">An open array of complex numbers that will be
    ///  used as the diagonal entries of the matrix.</param>
    constructor CreateDiagonal(const Elements: array of TASC); overload;

    /// <summary>Creates a new diagonal complex matrix using a given complex
    ///  vector.</summary>
    /// <param name="Elements">A complex vector that will be used to populate
    ///  the diagonal of the matrix.</param>
    constructor CreateDiagonal(const Elements: TComplexVector); overload;

    /// <summary>Creates a new diagonal complex matrix of a given size and a
    ///  (fixed) given value on the diagonal.</summary>
    /// <param name="ASize">The format of the matrix.</param>
    /// <param name="AVal">The value on the main diagonal.</param>
    constructor CreateDiagonal(ASize: Integer; AVal: TASC); overload;

    /// <summary>Creates a new complex matrix using pre-defined blocks.</summary>
    /// <param name="Blocks">An open array of the blocks of the matrix in row-
    ///  major order.</param>
    /// <param name="Cols">The number of blocks per row of the block matrix;
    ///  this value must be at least 1.</param>
    /// <remarks>Each block must contain at least a single element. In each row
    ///  of the block matrix, the number of rows in each block must be the same.
    ///  In each column of the block matrix, the number of columns in each block
    ///  must be the same. If the number of elements of <c>Blocks</c> is not
    ///  divisible by <c>Cols</c>, it will be padded with zero matrices until
    ///  the last row is filled.</remarks>
    /// <exception cref="EMathException">Raised if any of the constraints
    ///  specified above are violated.</exception>
    constructor Create(const Blocks: array of TComplexMatrix; Cols: Integer = 2); overload;

    // ***
    // Properties
    // ***

    /// <summary>The size of the matrix.</summary>
    property Size: TMatrixSize read FSize write SetSize;

    /// <summary>Gives access to the matrix elements using row and column indices.
    ///  </summary>
    /// <remarks>Row and column indices are zero-based.</remarks>
    property Elements[Row, Col: Integer]: TASC read GetElement write SetElement; default;

    /// <summary>Gives access to the rows of the matrix.</summary>
    /// <remarks>Row indices are zero-based.</remarks>
    property Rows[Index: Integer]: TComplexVector read GetRow write _DoSetRow;

    /// <summary>Gives access to the columns of the matrix.</summary>
    /// <remarks>Column indices are zero-based.</remarks>
    property Cols[Index: Integer]: TComplexVector read GetCol write _DoSetCol;

    /// <summary>Gives access to the first column.</summary>
    property FirstColumn: TComplexVector read GetFirstCol write SetFirstCol;

    /// <summary>Gives access to the last column.</summary>
    property LastColumn: TComplexVector read GetLastCol write SetLastCol;

    /// <summary>Gives access to the main diagonal of the matrix.</summary>
    /// <remarks>Diagonal indices are zero-based.</remarks>
    property MainDiagonal: TComplexVector read GetMainDiagonal write SetMainDiagonal;

    /// <summary>Gives access to the superdiagonal of the matrix.</summary>
    /// <remarks>Superdiagonal indices are zero-based.</remarks>
    property SuperDiagonal: TComplexVector read GetSuperDiagonal write SetSuperDiagonal;

    /// <summary>Gives access to the subdiagonal of the matrix.</summary>
    /// <remarks>Subdiagonal indices are zero-based.</remarks>
    property SubDiagonal: TComplexVector read GetSubDiagonal write SetSubDiagonal;

    /// <summary>Gives access to the antidiagonal of the matrix.</summary>
    /// <remarks>Antidiagonal indices are zero-based.</remarks>
    property AntiDiagonal: TComplexVector read GetAntiDiagonal write SetAntiDiagonal;

    /// <summary>Gives access to the underlying data as an array of complex numbers
    ///  in row-major order.</summary>
    /// <remarks>This may be cast to a <c>TComplexVector</c>. Notice that all
    ///  changes to the data will affect the original matrix. The <c>AsVector</c>
    ///  member function instead makes a copy of the data.</remarks>
    /// <see cref="AsVector" />
    property Data: TASCArray read FElements write FElements;

    /// <summary>Returns a pointer to the first element of a specific row.
    ///  The row index MUST be valid.</summary>
    property RowData[Index: Integer]: PASC read GetRowData;

    /// <summary>Returns a pointer to the first element of a specific row.
    ///  The function validates that the row index is valid, and raises an
    ///  exception if it is not.</summary>
    property SafeRowData[Index: Integer]: PASC read SafeGetRowData;

    /// <summary>The total size in bytes of this matrix.</summary>
    property MemorySize: Int64 read GetMemorySize;

    // ***
    // Operators
    // ***

    /// <summary>Creates a single-column matrix from a vector.</summary>
    class operator Implicit(const u: TComplexVector): TComplexMatrix;

    /// <summary>Promotes a real matrix to a complex matrix</summary>
    class operator Implicit(const A: TRealMatrix): TComplexMatrix;

    /// <summary>Treats a single-column matrix as a vector.</summary>
    class operator Explicit(const A: TComplexMatrix): TComplexVector; inline;

    /// <summary>Creates a single-entry matrix from a complex number.</summary>
    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;

    // ***
    // Tests
    // ***

    /// <summary>Returns true iff the matrix consists of a single row.</summary>
    function IsRow: Boolean; inline;

    /// <summary>Returns true iff the matrix consists of a single column.</summary>
    function IsColumn: Boolean; inline;

    /// <summary>Returns true iff the matrix is square.</summary>
    function IsSquare: Boolean; inline;

    /// <summary>Returns true iff the matrix is the N-dimensional identity
    ///  matrix.</summary>
    function IsIdentity(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is the m×n zero matrix (that is, if
    ///  it consists of only zeros).</summary>
    function IsZeroMatrix(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff all non-zero entries are located on the main
    ///  diagonal.</summary>
    function IsDiagonal(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is square and all non-zero entries
    ///  are located on the antidiagonal.</summary>
    function IsAntiDiagonal(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is the N-dimensional reversal
    ///  matrix.</summary>
    function IsReversal(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff all non-zero entries are located on or above
    ///  the main diagonal.</summary>
    /// <remarks>The matrix need not be square in order to be upper diagonal.
    ///  </remarks>
    function IsUpperTriangular(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff all non-zero entries are located on or below
    ///  the main diagonal.</summary>
    /// <remarks>The matrix need not be square in order to be lower diagonal.
    ///  </remarks>
    function IsLowerTriangular(const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns true iff the matrix is upper triangular or lower
    ///  triangular.</summary>
    function IsTriangular(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns the column index of the first non-zero element on a
    ///  given row.</summary>
    /// <returns>The column index of the first non-zero element on the given
    ///  row or -1 if the row consists of only zeros.</returns>
    function PivotPos(ARow: Integer; const Epsilon: Extended = 0): Integer;

    /// <summary>Returns whether or not a given row consists of only zeros.
    ///  </summary>
    function IsZeroRow(ARow: Integer; const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns whether or not a given row consists of only zeros except
    ///  for the right-most entry on the row.</summary>
    function IsEssentiallyZeroRow(ARow: Integer;
      const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is in row echelon form.</summary>
    function IsRowEchelonForm(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is in reduced row echelon form.
    ///  A matrix that is in reduced row echelon form is also in row echelon
    ///  form.</summary>
    function IsReducedRowEchelonForm(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determins whether the matrix is scalar, that is, if it is a
    ///  multiple of the identity matrix.</summary>
    function IsScalar(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is symmetric.</summary>
    function IsSymmetric(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is skew-symmetric.</summary>
    function IsSkewSymmetric(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is Hermitian.</summary>
    function IsHermitian(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is skew-Hermitian.</summary>
    function IsSkewHermitian(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is orthogonal.</summary>
    function IsOrthogonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is unitary.</summary>
    function IsUnitary(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is normal.</summary>
    function IsNormal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix consists only of ones and zeros.
    ///  </summary>
    function IsBinary(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is a permutation matrix.</summary>
    function IsPermutation(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is circulant.</summary>
    /// <remarks>A circulant matrix need not be square.</remarks>
    function IsCirculant(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is a Toeplitz matrix.</summary>
    /// <remarks>A Toeplitz matrix need not be square.</remarks>
    function IsToeplitz(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is a Hankel matrix.</summary>
    /// <remarks>A Hankel matrix need not be square.</remarks>
    function IsHankel(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is upper Hessenberg.</summary>
    /// <remarks>An upper Hessenberg matrix is always square.</remarks>
    function IsUpperHessenberg(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is lower Hessenberg.</summary>
    /// <remarks>A lower Hessenberg matrix is always square.</remarks>
    function IsLowerHessenberg(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix is tridiagonal.</summary>
    /// <remarks>A tridiagonal matrix is always square.</remarks>
    function IsTridiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is upper bidiagonal.</summary>
    /// <remarks>An upper bidiagonal matrix is always square.</remarks>
    function IsUpperBidiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is lower bidiagonal.</summary>
    /// <remarks>A lower bidiagonal matrix is always square.</remarks>
    function IsLowerBidiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is bidiagonal.</summary>
    /// <remarks>A bidiagonal matrix is always square.</remarks>
    function IsBidiagonal(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is centrosymmetric.</summary>
    function IsCentrosymmetric(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is a Vandermonde matrix.</summary>
    /// <remarks>A Vandermonde matrix need not be square.</remarks>
    function IsVandermonde(const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether the matrix commutes with a given matrix.
    ///  </summary>
    function CommutesWith(const A: TComplexMatrix;
      const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is idempotent, that is, if A^2 =
    ///  A.</summary>
    /// <remarks>An idempotent matrix is always square.</remarks>
    function IsIdempotent(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is an involution, that is, if A^2 =
    ///  I.</summary>
    /// <remarks>An involutory matrix is always square.</remarks>
    function IsInvolution(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is positive definite.</summary>
    /// <remarks>A positive definite matrix is always square and Hermitian.
    ///  </remarks>
    function IsPositiveDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is positive semidefinite.</summary>
    /// <remarks>A positive semidefinite matrix is always square and Hermitian.
    ///  </remarks>
    function IsPositiveSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is negative definite.</summary>
    /// <remarks>A negative definite matrix is always square and Hermitian.
    ///  </remarks>
    function IsNegativeDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is negative semidefinite.</summary>
    /// <remarks>A negative semidefinite matrix is always square and Hermitian.
    ///  </remarks>
    function IsNegativeSemiDefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the matrix is indefinite.</summary>
    /// <remarks>An indefinite matrix is always square and Hermitian.</remarks>
    function IsIndefinite(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Determines whether the square matrix is nilpotent.</summary>
    /// <exception cref="EMathException">Raised if the matrix isn't square or if
    ///  the method failed to determine nilpotency.</exception>
    function IsNilpotent(const Epsilon: Extended = 0): Boolean; inline;

    /// <summary>Returns the nilpotency index of the square matrix</summary>
    /// <returns>The nilpotency index if the matrix is nilpotent or -1 otherwise.
    ///  </returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square or if
    ///  the method failed to determine nilpotency.</exception>
    function NilpotencyIndex(const Epsilon: Extended = 0): Integer;

    /// <summary>Returns true iff the matrix is diagonally dominant, that is,
    ///  iff for every row the absolute value of the diagonal entry is greater
    ///  than or equal to the deleted absolute row sum.</summary>
    /// <remarks>Only defined for square matrices.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function IsDiagonallyDominant: Boolean;

    /// <summary>Returns true iff the matrix is strictly diagonally dominant,
    ///  that is, iff for every row the absolute value of the diagonal entry is
    ///  greater than the deleted absolute row sum.</summary>
    /// <remarks>Only defined for square matrices.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function IsStrictlyDiagonallyDominant: Boolean;

    // ***
    // Operations and computations
    // ***

    /// <summary>Returns the real part of the matrix.</summary>
    function RealPart: TRealMatrix;

    /// <summary>Returns the imaginary part of the matrix.</summary>
    function ImaginaryPart: TRealMatrix;

    /// <summary>Sets all elements to zero above the main diagonal.</summary>
    procedure MakeLowerTriangular;

    /// <summary>Sets all elements to zero below the main diagonal.</summary>
    procedure MakeUpperTriangular;

    /// <summary>Sets all elements to zero below the main subdiagonal.</summary>
    procedure MakeUpperHessenberg;

    /// <summary>Returns the square of the matrix.</summary>
    function Sqr: TComplexMatrix; inline;

    /// <summary>Returns the transpose of the matrix.</summary>
    function Transpose: TComplexMatrix;

    /// <summary>Returns the Hermitian conjugate (conjugate transpose, adjoint)
    ///  of the matrix.</summary>
    function Adjoint: TComplexMatrix;

    /// <summary>Returns the Hermitian square of the matrix.</summary>
    function HermitianSquare: TComplexMatrix; inline;

    /// <summary>Returns the modulus of the matrix.</summary>
    /// <remarks>The modulus is the unique positive-semidefinite square root of
    ///  the Hermitian square.</remarks>
    function Modulus: TComplexMatrix; inline;

    /// <summary>Returns the determinant of the matrix.</summary>
    function Determinant: TASC;

    /// <summary>Returns the trace of the matrix.</summary>
    function Trace: TASC;

    /// <summary>Returns the inverse of the matrix, if it exists.</summary>
    /// <exception cref="EMathException">Thrown if the matrix doesn't have an
    ///  inverse. This happens if the matrix isn't square or if its determinant
    ///  is zero (iff it doesn't have full rank).</exception>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    /// <see cref="TryInvert" />
    function Inverse: TComplexMatrix;

    /// <summary>Attempts to invert the matrix, which has to be square.</summary>
    /// <param name="AInverse">Receives the inverse of the matrix, if it can
    ///  be computed.</param>
    /// <returns>True iff the matrix can be inverted, that is, if and only if
    ///  it is non-singular.</returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    /// <see cref="Inverse" />
    function TryInvert(out AInverse: TComplexMatrix): Boolean;

    /// <summary>Returns the Moore-Penrose pseudoinverse of the matrix.</summary>
//    function Pseudoinverse: TComplexMatrix; experimental;

    /// <summary>Returns the rank of the matrix.</summary>
    function Rank: Integer; inline;

    /// <summary>Returns the nullity of the matrix.</summary>
    function Nullity: Integer; inline;

    /// <summary>Returns the condition number of the matrix using a specified
    ///  matrix norm.</summary>
    /// <exception cref="EMathException">Raised if the matrix isn't square and
    ///  invertible.</exception>
    function ConditionNumber(p: Integer = 2): TASR;

    /// <summary>Determines whether the matrix, which has to be square, is
    ///  singular.</summary>
    /// <returns>True iff the matrix is singular.</returns>
    /// <exception cref="EMathException">Thrown if the matrix isn't square.
    ///  </exception>
    function IsSingular: Boolean;

    /// <summary>Returns the Frobenius norm of the matrix, that is, the square
    ///  root of the sum of the squared moduli of all its elements.</summary>
    /// <remarks>This norm is submultiplicative.</remarks>
    function Norm: TASR; inline;

    /// <summary>Returns the square of the Frobenius norm of the matrix.</summary>
    function NormSqr: TASR; inline;

    /// <summary>Returns the (vector) p-norm of the matrix.</summary>
    function pNorm(const p: TASR): TASR; inline;

    /// <summary>Returns the max norm of the matrix, that is, the largest
    ///  modulus of its elements.</summary>
    ///  <remarks>This norm is NOT submultiplicative.</remarks>
    function MaxNorm: TASR; inline;

    /// <summary>Returns the sum norm of the matrix, that is, the sum of all
    ///  moduli of the entries.</summary>
    /// <remarks>This norm is submultiplicative.</remarks>
    function SumNorm: TASR; inline;

    /// <summary>Returns the k-norm of the matrix.</summary>
    function kNorm(const k: Integer): TASR; inline;

    /// <summary>Returns the maximum column sum norm of the matrix.</summary>
    /// <remarks>This norm is submultiplicative.</remarks>
    function MaxColSumNorm: TASR;

    /// <summary>Returns the maximum row sum norm of the matrix.</summary>
    /// <remarks>This norm is submultiplicative.</remarks>
    function MaxRowSumNorm: TASR;

    /// <summary>Returns the spectral norm of the matrix, that is, the largest
    ///  singular value.</summary>
    /// <remarks>This norm is submultiplicative.</remarks>
    function SpectralNorm: TASR; inline;

    /// <summary>Returns a deleted absolute row sum, that is, the sum of all
    ///  absolute-value entries on a given row except the row's diagonal entry
    ///  (if it has a diagonal entry).</summary>
    function DeletedAbsoluteRowSum(ARow: Integer): TASR;

    /// <summary>Swaps two rows in a matrix (in place).</summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowSwap(ARow1, ARow2: Integer): TComplexMatrix;

    /// <summary>Multiplies a single row in a matrix with a constant complex
    ///  number (in place).</summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowScale(ARow: Integer; AFactor: TASC): TComplexMatrix;

    /// <summary>Adds a multiple of one row to a second row in a matrix (in place).
    ///  </summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowAddMul(ATarget, ASource: Integer; AFactor: TASC;
      ADefuzz: Boolean = False; AFirstCol: Integer = 0): TComplexMatrix;

    /// <summary>Performs a given row operation to the matrix. The row operation
    ///  is given by a TRowOperationRecord.</summary>
    /// <returns>A reference to the original matrix.</returns>
    function RowOp(const ARowOp: TComplexRowOperationRecord): TComplexMatrix;

    /// <summary>Uses row operations to find a row echelon form of the matrix.
    ///  </summary>
    /// <param name="CollectRowOps">If set to true, the row operations performed
    ///  will be saved to the row op buffer.</param>
    /// <returns>Returns a row echelon form of the matrix.</returns>
    /// <remarks>The original matrix is not altered by this method; all row
    ///  operations are performed on a copy.</remarks>
    function RowEchelonForm(CollectRowOps: Boolean = False): TComplexMatrix;

    /// <summary>Uses row operations to find the reduced row echelon form of the
    ///  matrix.</summary>
    /// <param name="CollectRowOps">If set to true, the row operations performed
    ///  will be saved to the row op buffer.</param>
    /// <returns>Returns the reduced row echelon form of the matrix.</returns>
    /// <remarks>The original matrix is not altered by this method; all row
    ///  operations are performed on a copy.</remarks>
    function ReducedRowEchelonForm(CollectRowOps: Boolean = False): TComplexMatrix;

    /// <summary>Returns the number of zero rows in the matrix.</summary>
    function NumZeroRows(const AEpsilon: TASR = 0): Integer;

    /// <summary>Returns the number of zero rows at the bottom of the matrix.
    ///  </summary>
    function NumTrailingZeroRows(const AEpsilon: TASR = 0): Integer;

    /// <summary>Performs the Gram-Schmidt orthonormalization procedure on the
    ///  columns of the matrix; the result is returned.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function GramSchmidt: TComplexMatrix;

    procedure InplaceGramSchmidt(FirstCol, LastCol: Integer);

    /// <summary>Returns the matrix with zero or more columns removed so that
    ///  the remaining columns form a basis for the column space (range) of the
    ///  matrix.</summary>
    function ColumnSpaceBasis: TComplexMatrix;

    function ColumnSpaceProjection(const AVector: TComplexVector): TComplexVector;

    /// <summary>Returns the Euclidean distance between a vector and the column
    ///  space of the matrix.</summary>
    function DistanceFromColumnSpace(const AVector: TComplexVector): TASR;

    /// <summary>Computes and returns an upper Hessenberg matrix that is similar
    ///  to the original matrix, which must be square.</summary>
    /// <param name="A2x2Bulge">If true, the original matrix MUST be a PROPER
    ///  upper Hessenberg matrix except for a 2×2 bulge at the upper-left
    ///  position. In this case, a different algorithm suitable for this
    ///  particular case might be used, which is faster than the general-case
    ///  algorithm.</param>
    /// <remarks>If the original matrix is Hermitian, the result is tridiagonal.
    ///  </remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function SimilarHessenberg(A2x2Bulge: Boolean = False): TComplexMatrix;

    /// <summary>Returns the array of eigenvalues of the matrix.</summary>
    /// <returns>A complex vector of the n eigenvalues.</returns>
    /// <exception cref="EMathException">Raised if the eigenvalues couldn't be
    ///  computed, for any reason.</exception>
    function UnsortedEigenvalues: TComplexVector;

    function eigenvalues: TComplexVector;

    /// <summary>An alias for <c>eigenvalues</c>.</summary>
    function spectrum: TComplexVector; inline; // alias

    /// <summary>Computes the eigenvalues and eigenvectors of the matrix.
    ///  </summary>
    /// <param name="AEigenvalues">Receives the eigenvalues.</param>
    /// <param name="AEigenvectors">Receives a complex n×n matrix, the columns
    ///  of which are linearly independent eigenvectors of the matrix in the
    ///  same order as the eigenvalues.</param>
    /// <param name="ASkipVerify">If false (the default), the method verifies
    ///  that the obtained vectors really are linearly independent eigenvectors
    ///  of the matrix, in the same order as the eigenvalues. If true, this
    ///  verification is skipped. In any case, the eigenvalues are verified.
    ///  </param>
    /// <returns>True iff the eigenvalues and eigenvectors were successfully
    ///  computed and all checks complete successfully.</returns>
    /// <remarks>Currently only implemented for Hermitian matrices.</remarks>
    /// <exception cref="EMathException">Raised if the matrix isn't square
    ///  and Hermitian.</exception>
    function eigenvectors(out AEigenvalues: TComplexVector;
      out AEigenvectors: TComplexMatrix; ASkipVerify: Boolean = False): Boolean;

    /// <summary>Returns true iff the given vector is an eigenvector of the
    ///  square matrix.</summary>
    /// <exception cref="EMathException">Raised if the matrix isn't square or
    ///  if the given vector is of a different dimension.</exception>
    function IsEigenvector(const u: TComplexVector;
      const Epsilon: Extended = 0): Boolean;

    /// <summary>Returns the associated eigenvalue of a given eigenvector.
    ///  </summary>
    /// <exception cref="EMathException">Raised if the given vector isn't an
    ///  eigenvector of the matrix. In particular, this happens if the matrix
    ///  isn't square or if the vector is the zero vector of its dimension or
    //// of a different dimension compared to the matrix.</exception>
    function EigenvalueOf(const u: TComplexVector;
      const Epsilon: Extended = 0): TASC;

    /// <summary>Tries to find the associated eigenvalue of a given vector and
    ///  returns whether the given vector was an eigenvector or not.</summary>
    /// <param name="u">The given vector that might or might not be an eigenvector
    ///  of the matrix.</param>
    /// <param name="AEigenvalue">Receives the associated eigenvalue, if the
    ///  given vector turned out to be an eigenvector.</param>
    /// <returns>True iff the given vector is an eigenvector of the matrix.
    ///  </returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square
    ///  or if the dimension of the given vector doesn't match the size of the
    ///  matrix.</exception>
    function TryEigenvalueOf(const u: TComplexVector; out AEigenvalue: TASC;
      const Epsilon: Extended = 0): Boolean;

    /// <summary>Determines whether (lambda, u) is an eigenpair of the matrix.
    ///  </summary>
    function IsEigenpair(const lambda: TASC; const u: TComplexVector;
      const Epsilon: Extended = 0): Boolean;

    /// <summary>Computes an eigenvector associated with a given eigenvalue.
    ///  </summary>
    /// <param name="lambda">The eigenvalue to find an associated eigenvalue for.
    ///  </param>
    /// <returns>A normalized eigenvector associated with the given eigenvalue.
    ///  </returns>
    /// <exception cref="EMathException">Raised if <c>lambda</c> isn't an
    ///  eigenvalue of the matrix. In particular, this happens if the matrix
    ///  isn't square. Also raised if the algorithm fails to compute an
    ///  eigenvector, for any reason.</exception>
    function EigenvectorOf(const lambda: TASC): TComplexVector;

    /// <summary>Returns the spectral radius of the matrix, that is, the largest
    ///  absolute value of the eigenvalues.</summary>
    /// <remarks>The current implementation only works if all eigenvalues are
    ///  real; otherwise, an exception is raised.</remarks>
    function SpectralRadius: TASR; inline;

    /// <summary>Returns the array of singular values of the matrix.</summary>
    /// <remarks><para>The singular values are the eigenvalues of <c>sqrt(A* A)</c>,
    ///  or, equivalently, the square roots of the eigenvalues of the Hermitian
    ///  square <c>A* A</c>.</para>
    ///  <para>The singular values are listed in decreasing order.</para></remarks>
    function SingularValues: TRealVector;

    /// <summary>Returns a copy of the matrix with each element replaced with
    ///  its modulus.</summary>
    function Abs: TRealMatrix;

    /// <summary>Replaces all elements that are extremely close to zero with
    ///  zero.</summary>
    /// <returns>A reference to the matrix.</returns>
    /// <remarks>The matrix is altered by this method; the result only points
    ///  to the original matrix that has been modified.</remarks>
    function Defuzz(const Eps: Double = 1E-8): TComplexMatrix;

    /// <summary>Creates a new matrix that is identical to this matrix.</summary>
    function Clone: TComplexMatrix; // empty matrix aware

    /// <summary>Returns the elements of the matrix as a complex vector. The
    ///  elements are obtained from the matrix in column-major order. In
    ///  particular, this works well if the matrix is either a single row or a
    ///  single column.</summary>
    /// <see cref="AsVector" />
    function Vectorization: TComplexVector; inline;

    /// <summary>Returns the elements of the matrix as a complex vector. The
    ///  elements are obtained from the matrix in row-major order. In particular,
    ///  this works well if the matrix is either a single row or a single column.
    ///  </summary>
    /// <see cref="Vectorization" />
    /// <see cref="Data" />
    function AsVector: TComplexVector; inline;

    /// <summary>Augments the matrix with a new matrix.</summary>
    /// <param name="A">The complex matrix to augment the current matrix with.
    ///  This matrix must have the same number of rows as the current matrix.
    ///  </param>
    /// <returns>A complex matrix with A augmented to the right of the current
    ///  matrix.</return>
    /// <remarks>The current matrix is not modified by the method. A vector
    ///  may be used as the argument, since it will be automatically promoted to
    ///  the corresponding single-column matrix.</remarks>
    /// <see cref="AddRow" />
    function Augment(const A: TComplexMatrix): TComplexMatrix; overload; // empty matrix aware

    /// <summary>Augments the matrix with a zero column.</summary>
    /// <returns>A complex matrix with a zero column augmented to the right of the
    ///  current matrix.</return>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Augment: TComplexMatrix; overload; inline;

    /// <summary>Replaces a specific row with a new set of values.</summary>
    procedure SetRow(ARowIndex: Integer; const ARow: array of TASC);

    /// <summary>Replaces a specific column with a new set of values.</summary>
    procedure SetCol(AColIndex: Integer; const ACol: array of TASC);

    /// <summary>Returns the (m-1)×(n-1) submatrix obtained by removing a single
    ///  row and a single column from the matrix.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Submatrix(ARowToRemove: Integer = 0;
      AColToRemove: Integer = 0; AAllowEmpty: Boolean = False): TComplexMatrix; overload; // empty matrix aware

    /// <summary>Returns a submatrix of the current matrix using the specified
    ///  rows and columns.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Submatrix(const ARows: array of Integer;
      const ACols: array of Integer): TComplexMatrix; overload;

    /// <summary>Returns a principal submatrix of the current matrix using the
    ///  specified rows and columns.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Submatrix(const ARows: array of Integer): TComplexMatrix; overload;

    /// <summary>Returns a leading principal submatrix of the matrix.</summary>
    function LeadingPrincipalSubmatrix(const ASize: Integer): TComplexMatrix;

    /// <summary>Returns the submatrix obtained by removing the last column.
    ///  </summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Lessened: TComplexMatrix;

    /// <summary>Returns the determinant of the submatrix obtained by removing
    ///  one specific row and one specific column.</summary>
    function Minor(ARow, ACol: Integer): TASC; inline;

    /// <summary>Returns (-1)^(ARow + ACol) times the corresponding minor.
    ///  </summary>
    function Cofactor(ARow, ACol: Integer): TASC; inline;

    /// <summary>Returns the cofactor matrix of the matrix.</summary>
    function CofactorMatrix: TComplexMatrix;

    /// <summary>Returns the adjugate matrix of the matrix, that is, the
    ///  transpose of the cofactor matrix.</summary>
    function AdjugateMatrix: TComplexMatrix; inline;

    /// <summary>Returns a copy of the matrix with each element replaced by
    ///  the result of a complex-valued function applied to it.</summary>
    /// <remarks>The current matrix is not modified by the method.</remarks>
    function Apply(AFunction: TComplexFunctionRef): TComplexMatrix;

    /// <summary>Replaces each entry of the matrix satisfying the given predicate
    ///  with a given value.</summary>
    /// <param name="APredicate">A predicate for complex numbers. Each matrix entry
    ///  satisfying this predicate will be replaced with <c>ANewValue</c>.</param>
    /// <param name="ANewValue">The value to replace entries with.</param>
    /// <returns>A copy of the current matrix but with all entries satisfying
    ///  <c>APredicate</c> replaced with <c>ANewValue</c>.</returns>
    /// <remarks>The current matrix is not modified by this method.</remarks>
    function Replace(APredicate: TPredicate<TASC>;
      const ANewValue: TASC): TComplexMatrix; overload;

    /// <summary>Replaces each entry with a specified value with a new value.
    ///  </summary>
    /// <param name="AOldValue">The value to search for in the original matrix.
    ///  </param>
    /// <param name="ANewValue">The value to replace with.</param>
    /// <returns>A copy of the current matrix but with each entry originally
    ///  eaual to <c>AOldValue</c> replaced with <c>ANewValue</c>.</returns>
    /// <remarks>The current matrix is not modified by this method.</remarks>
    function Replace(const AOldValue, ANewValue: TASC;
      const Epsilon: Extended = 0): TComplexMatrix; overload;

    /// <summary>Sets each entry of the matrix to the same, new, value.</summary>
    /// <param name="ANewValue">The new value to set each entry to.</param>
    /// <returns>A matrix of the same size as the original matrix but with each
    ///  entry equal to <c>ANewValue</c>.</returns>
    /// <remarks>The current matrix is not modified by this method.</remarks>
    function Replace(const ANewValue: TASC): TComplexMatrix; overload;

    /// <summary>Returns either a single-line textual representation of the
    /// matrix of the form ((M00, M01, ...), (M10, M11, ...), ...) or a multi-
    /// line tabular representation.</summary>
    function str(const AOptions: TFormatOptions): string; // mainly for debugging

    /// <summary>Adds a row to the matrix.</summary>
    /// <remarks>This is more of a computing-related function than a mathematical
    ///  operation. Notice that it does suffer from the
    ///  <c>SetLength(X, Length(X) + 1)</c> problem.</remarks>
    /// <see cref="Augment" />
    procedure AddRow(const ARow: array of TASC);

    /// <summary>Sorts the elemenets of the matrix.</summary>
    /// <remarks>This operation sorts the elements of the original matrix,
    ///  and returns a reference to it.</remarks>
    /// <returns>A reference to the original matrix, which has been sorted.
    ///  </returns>
    function Sort(AComparer: IComparer<TASC>): TComplexMatrix;

    /// <summary>Shuffles the elemenets of the matrix.</summary>
    /// <remarks>This operation shuffles the elements of the original matrix,
    ///  and returns a reference to it.</remarks>
    /// <returns>A reference to the original matrix, which has been shuffled.
    ///  </returns>
    function Shuffle: TComplexMatrix;

    /// <summary>Reverses the order of the elements in the matrix (in the
    /// row-major sense).</summary>
    /// <remarks>This operation reverses the elements of the original matrix,
    ///  and returns a reference to it.</remarks>
    /// <returns>A reference to the original matrix, which has been reversed.
    ///  </returns>
    function Reverse: TComplexMatrix;

    // ***
    // Decompositions
    // ***

    /// <summary>Produces a LU decomposition of the matrix, which has to be square,
    ///  that is, finds a permutation matrix P, a unit lower triangular matrix L,
    ///  and an upper triangular matrix U such that Self = P.Adjoint * L * U.
    ///  </summary>
    /// <param name="P">Receives the permutation matrix P.</param>
    /// <param name="L">Receives the unit lower triangular matrix L.</param>
    /// <param name="U">Receives the upper triangular matrix U.</param>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    procedure LU(out P, L, U: TComplexMatrix);

    /// <summary>Computes the Cholesky decomposition of the square matrix, if
    ///  possible.</summary>
    /// <param name="R">Receives the upper triangular Cholesky factor R such
    ///  that Self = R.Adjoint * R.</param>
    /// <returns>True iff the decomposition was successfully computed; this
    ///  happens precisely when the matrix is positive definite.</returns>
    /// <exception cref="EMathException">Raised if the matrix isn't square.
    ///  </exception>
    function Cholesky(out R: TComplexMatrix): Boolean;

    /// <summary>Produces a unitary matrix Q and an upper triangular matrix
    ///  R such that Self = Q * R.</summary>
    /// <remarks>QR decomposition is only implemented for square and tall matrices.
    ///  </remarks>
    /// <exception cref="EMathException">Raised if the matrix is wide, that is,
    ///  if <c>Size.Cols > Size.Rows</c>.</exception>
    procedure QR(out Q, R: TComplexMatrix);

    /// <summary>Finds an upper Hessenberg matrix A and an orthogonal matrix
    ///  U such that Self = U A U.Transpose.</summary>
    procedure Hessenberg(out A, U: TComplexMatrix);

  end;

/// <summary>Computes and returns the Nth power of the square matrix A.</summary>
/// <remarks>The zeroth power is the identity matrix of the same size as A.
///  Negative powers are defined iff A is invertible, and if so the result
///  is the |N|th power of the inverse of A (which is equal to the inverse
///  of the |N|th power of A).</remarks>
function mpow(const A: TComplexMatrix; const N: Integer): TComplexMatrix; overload;

/// <summary>Returns the unique positive semidefinite square root of a given
///  positive semidefinite symmetric matrix.</summary>
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;

/// <summary>Returns the complex zero matrix of a given size.</summary>
function ComplexZeroMatrix(const ASize: TMatrixSize): TComplexMatrix; inline;

/// <summary>Returns the complex identity matrix of a given size.</summary>
function ComplexIdentityMatrix(ASize: Integer): TComplexMatrix; // empty matrix aware

/// <summary>Returns the complex reversal matrix of a given size.</summary>
function ComplexReversalMatrix(ASize: Integer): TComplexMatrix;

/// <summary>Returns a diagonal matrix constructed using values from
///  a given array of real numbers.</summary>
/// <see cref="DirectSum" />
function diag(const AElements: array of TASC): TComplexMatrix; overload;

/// <summary>Returns a diagonal matrix constructed using values from
///  a given complex vector.</summary>
/// <see cref="DirectSum" />
function diag(const AElements: TComplexVector): TComplexMatrix; inline; overload;

/// <summary>Returns the outer product between two vectors.</summary>
function OuterProduct(const u, v: TComplexVector): TComplexMatrix; overload; inline;

/// <summary>Returns the n×n circulant matrix whose first row is a given vector.
///  </summary>
function CirculantMatrix(const AElements: array of TASC): TComplexMatrix; overload;

/// <summary>Returns the Toeplitz matrix with a given first row and a given first
///  column.</summary>
function ToeplitzMatrix(const AFirstRow,
  AFirstCol: array of TASC): TComplexMatrix; overload;

/// <summary>Returns the Hankel matrix with a given first row and a given last
///  column.</summary>
function HankelMatrix(const AFirstRow,
  ALastCol: array of TASC): TComplexMatrix; overload;

/// <summary>Returns a Vandermonde matrix produced by a given vector.</summary>
/// <param name="AElements">The vector used as the second column in the matrix.
///  </param>
/// <param name="ACols">The number of columns in the matrix. If omitted,
///  the matrix will be square.</param>
function VandermondeMatrix(const AElements: array of TASC;
  ACols: Integer = 0): TComplexMatrix; overload;

/// <summary>Returns the matrix for a reflection in a hyperplane in some dimension.
///  </summary>
/// <param name="u">A vector perpendicular to the hyperplane.</param>
/// <remarks>A reflection matrix like this is also called a Householder matrix.
///  </remarks>
function ReflectionMatrix(const u: TComplexVector): TComplexMatrix; overload;

/// <summary>Returns the matrix for a reflection in a hyperplane in some dimension.
///  </summary>
/// <param name="u">A unit vector perpendicular to the hyperplane.</param>
/// <remarks>A reflection matrix like this is also called a Householder matrix.
///  The vector u MUST be a unit vector.</remarks>
function QuickReflectionMatrix(const u: TComplexVector): TComplexMatrix; overload; inline;

/// <summary>Returns the Hadamard product between two matrices of the same
///  size, that is, the matrix obtained by taking element-wise products.</summary>
function HadamardProduct(const A, B: TComplexMatrix): TComplexMatrix; overload;

/// <summary>Returns the direct sum of two matrices.</summary>
/// <see cref="diag" />
function DirectSum(const A, B: TComplexMatrix): TComplexMatrix; overload; inline;

/// <summary>Returns the direct sum of an open array of matrices.</summary>
/// <see cref="diag" />
function DirectSum(const Blocks: array of TComplexMatrix): TComplexMatrix; overload;

/// <summary>Determines whether two matrices commute or not.</summary>
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;

/// <summary>Returns the sum of all elements of the matrix.</summary>
function sum(const A: TComplexMatrix): TASC; overload; inline;

/// <summary>Returns the arithmetic mean of all elements of the matrix.</summary>
function ArithmeticMean(const A: TComplexMatrix): TASC; overload; inline;

/// <summary>Returns the geometric mean of all elements of the matrix.</summary>
function GeometricMean(const A: TComplexMatrix): TASC; overload; inline;

/// <summary>Returns the harmonic mean of all elements of the matrix.</summary>
function HarmonicMean(const A: TComplexMatrix): TASC; overload; inline;

/// <summary>Returns the product of all elements of the matrix.</summary>
function product(const A: TComplexMatrix): TASC; overload; inline;

/// <summary>Returns true iff there exists an entry of the matrix satisfying
///  the given predicate.</summary>
function exists(const A: TComplexMatrix; APredicate: TPredicate<TASC>): Boolean; overload;

/// <summary>Returns the number of matrix entries that satisfy the given
///  predicate.</summary>
function count(const A: TComplexMatrix; APredicate: TPredicate<TASC>): Integer; overload;

/// <summary>Returns the number of instances of a given value in the matrix.
///  </summary>
function count(const A: TComplexMatrix; const AValue: TASC): Integer; overload; inline;

/// <summary>Returns true iff all entries of the matrix satisfy
///  the given predicate.</summary>
function ForAll(const A: TComplexMatrix; APredicate: TPredicate<TASC>): Boolean; overload;

/// <summary>Returns true iff a specific value is found in the matrix.</summary>
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;

/// <summary>Solves the matrix equation AX = Y where A is a lower triangular
///  n×n matrix and X and Y are n-dimensional vectors using forward substitution.
///  </summary>
/// <param name="A">A square matrix of size n, which MUST be lower triangular.</param>
/// <param name="Y">An n-dimensional vector.</param>
/// <param name="IsUnitDiagonal">If true, the routine assumes that the main
///  diagonal of <c>A</c> consists of only 1's, without caring about the actual
///  elements on the main diagonal.</param>
/// <returns>The unique solution X of the equation AX = Y, if A is non-singular.</returns>
/// <exception cref="EMathException">Raised if A is not square and non-singular.
///  Also raised if the size of A differs from the dimension of Y.</exception>
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;

/// <summary>Solves the matrix equation AX = Y where A is an upper triangular
///  n×n matrix and X and Y are n-dimensional vectors using backward
///  substitution.</summary>
/// <param name="A">A square matrix of size n, which MUST be upper triangular.</param>
/// <param name="Y">An n-dimensional vector.</param>
/// <returns>The unique solution X of the equation AX = Y, if A is non-singular.</returns>
/// <exception cref="EMathException">Raised if A is not square and non-singular.
///  Also raised if the size of A differs from the dimension of Y.</exception>
function BackSubstitution(const A: TComplexMatrix; const Y: TComplexVector): TComplexVector; overload;

/// <summary>Solves the matrix equation AX = Y.</summary>
/// <param name="AAugmented">The augmented matrix [A|Y].</param>
/// <returns>The unique solution X if a unique solution exists.</returns>
/// <exception cref="EMathException">Raised if a unique solution does not
///  exist.</exception>
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;

/// <summary>Solves the matrix equation AX = Y.</summary>
/// <param name="A">The matrix of coefficients, A.</param>
/// <param name="Y">The right-hand side of the equation, Y.</param>
/// <returns>The unique solution X if a unique solution exists.</returns>
/// <exception cref="EMathException">Raised if a unique solution does not
///  exist.</exception>
function SysSolve(const A: TComplexMatrix; const Y: TComplexVector): TComplexVector; overload;
function SysSolve(const A: TComplexMatrix; const Y: TComplexMatrix): TComplexMatrix; overload;

/// <summary>Given a set { (x_0, y_0), (x_1, y_1), ..., (x_(k-1), y_(k-1)) } of
///  points, this function determines the polynomial of degree at most n
///  whose graph best connects to the points in the least-squares sense.</summary>
/// <param name="X">The k-dimensional vector of X coordinates.</param>
/// <param name="Y">The k-dimensional vector of Y coordinates.</param>
/// <param name="ADegree">The maximum degree n >= 0 of the polynomial.</param>
/// <returns>A complex vector consisting of the n + 1 coefficients c_0, c_1, ...,
///  c_n of the optimal polynomial c_0 + c_1 x + c_2 x^2 + ... + c_n x^n.</returns>
function LeastSquaresPolynomialFit(const X, Y: TComplexVector;
  ADegree: Integer): TComplexVector; overload;

/// <summary>A quick but unsafe procedure to copy data from one vector
///  to another.</summary>
/// <param name="ASource">The source vector.</param>
/// <param name="AFrom">The first component of the source vector to be copied.
///  </param>
/// <param name="ATo">The last component of the source vector to be copied.
///  </param>
/// <param name="ATarget">The vector to copy to.</param>
/// <param name="ATargetFrom">The index of the first component of the target
///  vector to be overwritten by the copied data.</param>
/// <remarks>The specified components in the source vector MUST exist, and the
///  target vector MUST be large enough to contain all data.</remarks>
procedure VectMove(const ASource: TComplexVector; const AFrom, ATo: Integer;
  var ATarget: TComplexVector; const ATargetFrom: Integer = 0); overload; inline;

/// <summary>A quick but unsafe procedure to copy data from a vector to a part
///  of a column of a matrix.</summary>
/// <param name="ASource">The source vector.</param>
/// <param name="AFrom">The first component of the source vector to be copied.
///  </param>
/// <param name="ATo">The last component of the source vector to be copied.
///  </param>
/// <param name="ATarget">The matrix to write data to.</param>
/// <param name="ATargetCol">The index of the column of the matrix to write to.
///  </param>
/// <param name="ATargetFirstRow">The index of the first row in the matrix
///  column to overwrite.</param>
/// <remarks>The specified components in the source vector MUST exist, and the
///  target matrix MUST be large enough to contain all data.</remarks>
procedure VectMoveToMatCol(const ASource: TComplexVector; const AFrom, ATo: Integer;
  var ATarget: TComplexMatrix; const ATargetCol: Integer;
  const ATargetFirstRow: Integer = 0); overload;

/// <summary>A quick but unsafe procedure to copy data from a vector to a part
///  of a row of a matrix.</summary>
/// <param name="ASource">The source vector.</param>
/// <param name="AFrom">The first component of the source vector to be copied.
///  </param>
/// <param name="ATo">The last component of the source vector to be copied.
///  </param>
/// <param name="ATarget">The matrix to write data to.</param>
/// <param name="ATargetRow">The index of the row of the matrix to write to.
///  </param>
/// <param name="ATargetFirstCol">The index of the first column in the matrix
///  row to overwrite.</param>
/// <remarks>The specified components in the source vector MUST exist, and the
///  target matrix MUST be large enough to contain all data.</remarks>
procedure VectMoveToMatRow(const ASource: TComplexVector; const AFrom, ATo: Integer;
  var ATarget: TComplexMatrix; const ATargetRow: Integer;
  const ATargetFirstCol: Integer = 0); overload; inline;

/// <summary>A quick but unsafe procedure to copy data from one matrix to
///  another.</summary>
/// <param name="ASource">The matrix to copy from.</param>
/// <param name="ARect">A rectangle describing the submatrix to copy; both
///  the top-left and the bottom-right elements are included in the submatrix.
///  </param>
/// <param name="ATarget">The matrix to copy to.</param>
/// <param name="ATargetTopLeft">The top-left element in the destination matrix
///  to be overwritten.</param>
/// <remarks>The specified elements in the source matrix MUST exist, and the
///  target matrix MUST be large enough to contain all data.</remarks>
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;

/// <summary>A quick but unsafe procedure to copy data from a column in a matrix
///  to a vector.</summary>
/// <param name="ASource">The matrix to copy data from.</param>
/// <param name="AColumn">The column of the matrix to copy data from.</param>
/// <param name="AFrom">The row index of the first element in the column to
///  copy.</param>
/// <param name="ATo">The row index of the last element in the column to
///  copy.</param>
/// <param name="ATarget">The vector to copy data into.</param>
/// <param name="ATargetFrom">The index of the first component of the vector
///  to overwrite.</param>
/// <remarks>The specified elements in the source matrix MUST exist, and the
///  target vector MUST be large enough to contain all data.</remarks>
procedure MatMoveColToVect(const ASource: TComplexMatrix; const AColumn: Integer;
  const AFrom, ATo: Integer; var ATarget: TComplexVector;
  const ATargetFrom: Integer = 0); overload;

/// <summary>A quick but unsafe procedure to copy data from a row in a matrix
///  to a vector.</summary>
/// <param name="ASource">The matrix to copy data from.</param>
/// <param name="ARow">The row of the matrix to copy data from.</param>
/// <param name="AFrom">The column index of the first element in the row to
///  copy.</param>
/// <param name="ATo">The column index of the last element in the row to
///  copy.</param>
/// <param name="ATarget">The vector to copy data into.</param>
/// <param name="ATargetFrom">The index of the first component of the vector
///  to overwrite.</param>
/// <remarks>The specified elements in the source matrix MUST exist, and the
///  target vector MUST be large enough to contain all data.</remarks>
procedure MatMoveRowToVect(const ASource: TComplexMatrix; const ARow: Integer;
  const AFrom, ATo: Integer; var ATarget: TComplexVector;
  const ATargetFrom: Integer = 0); overload; inline;

/// <summary>A quick but unsafe procedure to perform an inplace multiplication
///  of a submatrix with a square matrix from the left.</summary>
/// <param name="ATarget">The matrix to modify.</param>
/// <param name="ARect">A rectangle that identifies the submatrix to multiply.
///  Both the upper-left and the bottom-right corners belong to the submatrix.
///  The submatrix has format n × m. This parameter MUST specify a valid
///  submatrix.</param>
/// <param name="AFactor">The matrix to multiply the submatrix by, from the left.
///  This MUST be a matrix of size n × n.</param>
procedure MatMulBlockInplaceL(var ATarget: TComplexMatrix; const ARect: TRect;
  const AFactor: TComplexMatrix); overload;

/// <summary>A quick but unsafe procedure to perform an inplace multiplication
///  of a square submatrix with a square matrix from the left.</summary>
/// <param name="ATarget">The matrix to modify.</param>
/// <param name="APoint">The index of the upper-left element of the submatrix.</param>
/// <param name="AFactor">The matrix to multiply the submatrix by, from the left.
///  This MUST be a matrix of size n × n.</param>
/// <remarks>The implied n × n submatrix MUST be valid.</remarks>
procedure MatMulBlockInplaceL(var ATarget: TComplexMatrix; const ATopLeft: TPoint;
  const AFactor: TComplexMatrix); overload; inline;

/// <summary>
///  <para>Performs the operation</para>
///  <para><c>ATarget := ATarget * DirectSum(IdentityMatrix(k), AFactor)</c></para>
///  <para>where <c>k = n - m</c>, <c>n</c> is the size of the square matrix
///   <c>ATarget</c>, and <c>m <= n</c> is the size of the square matrix
///   <c>AFactor</c>.</para>
/// </summary>
/// <param name="ATarget">An n×n matrix. This MUST be square.</param>
/// <param name="AFactor">An m×m matrix which MUST be square and where m
///  MUST be <= n.</param>
procedure MatMulBottomRight(var ATarget: TComplexMatrix; const AFactor: TComplexMatrix); overload;

/// <summary>A quick but unsafe procedure to fill a submatrix with a given
///  complex number.</summary>
/// <param name="ATarget">The matrix to modify.</param>
/// <param name="ARect">A rectangle describing the submatrix to fill; both
///  the top-left and the bottom-right elements of the rectangle belong to the
///  submatrix.</param>
/// <param name="Value">The complex number to fill the submatrix with.</param>
/// <remarks>The specified submatrix MUST belong to the matrix.</remarks>
procedure MatBlockFill(var ATarget: TComplexMatrix; const ARect: TRect;
  const Value: TASC); overload;


// Various square roots
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;


// Common technical declarations and definitions
const
  /// <summary>Indicates whether the ASNum.pas unit was compiled with the
  ///  <c>REALCHECK</c> compiler directive. This compiler directive enables the
  ///  "real check" feature of the complex elementary functions. When turned on,
  ///  these functions first test to see if the given argument happens to be a real
  ///  number (zero imaginary part). If so, the complex elementary function calls
  ///  the standard real elementary function instead of using the general, complex
  ///  algorithm or formula to compute the value. In general, this improves
  ///  performance and accuracy if a real value is passed to a complex function.
  ///  </summary>
  /// <remarks>Normally, this functionality should be turned on. It may be turned
  ///  off for purposes of testing, however.</remarks>
  RealCheck = {$IFDEF REALCHECK}True{$ELSE}False{$ENDIF};

var
  /// <summary>The maximum number of threads to use for computations. The actual
  ///  number of threads used in a computation is never larger than this value,
  ///  but it may be smaller if the computer doesn't have enough "threads" (cores
  ///  or HT). The actual maximum number of threads is the smallest value of this
  ///  variable and the number of logical processors as reported by the operating
  ///  system.</summary>
  MaxThreadCount: Integer = 8;

function GetTotalCacheSize: Integer;
procedure ClearCaches;

const
  InvLn10: TASC = (Re: 0.434294481903251827651128918916605082294397005803666566114; Im: 0);

implementation

uses
  Classes, StrUtils, Windows, Character, WideStrUtils;

var
  FS: TFormatSettings;

procedure DoYield; inline;
begin
  if Assigned(GTYieldProc) then
    GTYieldProc;
end;


{ **************************************************************************** }
{ Utilities                                                                    }
{ **************************************************************************** }

// TODO: The TInt64Guard record can very much be improved.

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;

/// <summary>Returns (-1)^N.</summary>
/// <returns>(-1)^N</returns>
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;


{ **************************************************************************** }
{ TFormatStyleHelper                                                           }
{ **************************************************************************** }

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;


{ **************************************************************************** }
{ TRationalNumber                                                              }
{ **************************************************************************** }

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));

  // Get digit array
  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)));

        // Compute result length
        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); // decimal separator
        if (NumFracDigits > 0) and (AFormatOptions.Numbers.FracGrouping <> 0) then
          Inc(ResLen, Ceil(NumFracDigits / AFormatOptions.Numbers.FracGrouping) - 1);

        // Build result
        SetLength(Result, ResLen);

        c := 0; // index of last char written (valid at the beginning of each new part)

        // - Minus sign
        if (Numerator < 0) and (FirstSigDigitPos <> FirstSigDigitPos.MaxValue) then
        begin
          Inc(c);
          Result[c] := AFormatOptions.Numbers.MinusSign;
        end;

        // - Integer part
        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

          // - Decimal separator
          Inc(c);
          Result[c] := AFormatOptions.Numbers.DecimalSeparator;

          // - Fractional part
          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

        // Compute result length
        ResLen := NumDigits;
        if Numerator < 0 then
          Inc(ResLen, 1);
        if NumFracDigits >= 1 then
          Inc(ResLen, 1); // decimal separator
        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 {always a sign on the exp, even if >= 0});
        if (NumFracDigits > 0) and (AFormatOptions.Numbers.FracGrouping <> 0) then
          Inc(ResLen, Ceil(NumFracDigits / AFormatOptions.Numbers.FracGrouping) - 1);

        // Build result
        SetLength(Result, ResLen);

        c := 0; // index of last char written (valid at the beginning of each new part)

        // - Minus sign
        if Numerator < 0 then
        begin
          Inc(c);
          Result[c] := AFormatOptions.Numbers.MinusSign;
        end;

        // - MSD
        Inc(c);
        Result[c] := __DigChr(D[0]);

        // - Decimal separator
        if NumFracDigits > 0 then
        begin
          Inc(c);
          Result[c] := AFormatOptions.Numbers.DecimalSeparator;
        end;

        // - Fractional digits
        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;

        // - Exponent
        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;
  // INTEGER
  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
  // NOT INTEGER
  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); // will reduce
    tmp2 := RationalNumber(Y.Numerator, X.Denominator); // will reduce
    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;


{ **************************************************************************** }
{ TSimpleSymbolicForm                                                          }
{ **************************************************************************** }

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;


{ **************************************************************************** }
{ TASR: The AlgoSim real number type                                           }
{ **************************************************************************** }

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)); { "<=" is correct: sign bit and minus sign }
      {$WARN COMPARISON_TRUE OFF}
      Assert(AOptions.Numbers.MinLength.MaxValue < Length(NumDigits)); { "<" is correct: minus sign not included in min len }
      {$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;

  // Minus sign
  if (Result.Length > 0) and (Result[1] = '-') then
    Result[1] := AOptions.Numbers.MinusSign;

  // Digit grouping
  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 {sufficient} do
    if Result[i] = '-' then
      Result[i] := AOptions.Numbers.MinusSign;

  if AOptions.Numbers.PrettyExp then
    for i := 1 to Result.Length - 1 {required!} 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;

    // Digit grouping on the integral part
    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;

    // Digit grouping on the fractional part
    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;


{ **************************************************************************** }
{ TASRComparer                                                                 }
{ **************************************************************************** }

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;


{ **************************************************************************** }
{ Elementary functions for real numbers                                        }
{ **************************************************************************** }

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.'); // missed by Math.Power (which gladly returns 1...)
  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;


{ **************************************************************************** }
{ TASC: The AlgoSim complex number type                                        }
{ **************************************************************************** }

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;


{ **************************************************************************** }
{ TASCComparer                                                                 }
{ **************************************************************************** }

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;


{ **************************************************************************** }
{ Elementary functions for complex numbers                                     }
{ **************************************************************************** }

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

  // cpow(z, w) is defined for all z, w \in \C with z \ne 0 or w.Re > 0.

  // Warning! It is NOT true that z, w \in \R => cpow(z, w) = pow(z.Re, w.Re)
  // for all z and w for which cpow(z, w) is defined; a counterexample is
  // z = -2.3, w = 2.1.

  // Handle case z = 0
  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.');

  // Handle case z \ne 0: Result must "equal" cexp(w*cln(z));
  if IsZero(w.Im) and IsZero(frac(w.Re)) and (abs(w.Re) < 13) then
  begin
    // improve z^n where n \in [-12, 12] \cap \Z
    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 {0^0}
      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
    // improve real case
    Exit(pow(z.Re, w.Re))

  else if (z.Re = 0) and (w.Im = 0) and (frac(w.Re) = 0) then // notice that z.Im \ne 0
    // improve (bi)^n (notice that b \ne 0)
    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));   // Check
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));  // Check
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));                       // Check
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))                      // Check
  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)); // Check
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)); // Check
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));                                           // Check
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));                                // Check
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;                                       // Check
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;                                     // Check
  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));                                // Check
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));                                           // Check
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;


{ **************************************************************************** }
{ Integer functions                                                            }
{ **************************************************************************** }

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); // all digits after the decimal separator are 0 in an integer
  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

  // Preparation

  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 // initial zeros 0.00...
        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 // initial zeros 0.00...
        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;

  // Divmod loop

  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;

  // Rounding

  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;

  // Removal of trailing zeros

  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;
  // Now IndexHigh - PRECISION_DIGITS + 1 <= Index <= IndexHigh
  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]);
  // Unlike the GCD, the LCM might be greater than the arguments, so we need to check for overflow.
  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);

  // In general, we need primes up to Sqrt(N). But very often you are lucky,
  // and can tell it isn't prime by just having the first few primes available.
  // So don't spend a minute computing primes up to 1000000000 just in order to
  // find out that 1000000000000000000 is composite (the very first prime will do).
  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.');

  Result := _IsPrime[N];

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); // now N is the zero-based index

  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
//  Result := (IntPower(GoldenRatio, N) - IntPower(GoldenRatioConj, N)) / sqrt5; // less precise for large N and naive summation not much slower
  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; // to make compiler happy
  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; // to make sloppy compiler happy
  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 - floor(x / y) * y;
  Result := x mod y;
  if Result < 0 then
    Inc(Result, y);
end;

function imod(const x, y: TASI): TASI; inline;
begin
//  Result := x - floor(x / y) * y;
  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{empty product = 1};

  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

  N := Abs(N);

  if N = 0 then
    Exit(False);

  if Assigned(_IsPrime) and (N <= _IsPrime.Size - 1) and IsPrime(N) then
    Exit(True);

  a := N;

  ExpandPrimeCache(10000);

  for i := 0 to _Primes.Count - 1 do
  begin

    if _Primes[i] > Sqrt(a) then
      Exit(True);

    if a mod _Primes[i] = 0 then
    begin
      a := a div _Primes[i];
      if a = 1 then
        Exit(True);
      if a mod _Primes[i] = 0 then
        Exit(False);
    end

  end;

  if ExpandPrimeCache(Round(Sqrt(N))) then
    Result := IsSquareFree(N)
  else
    raise EMathException.Create('Couldn''t factor integer.');

end;

function factorize(const N: TASI): TArray<TASI>;
var
  neg: Boolean;
begin

  if (N = -1) or (N = 0) or (N = 1) then
  begin
    Result := [N];
    Exit;
  end;

  neg := N < 0;
  Result := PrimeFactors(Abs(N));
  if neg then
    Result[0] := -Result[0]; // exists since |N| <> 1

end;

function GetFactorizedString(const N: TASI; AMinusSign: Char = '-';
  AMultiplicationSign: Char = '*'; ASpace: Boolean = False): string;
begin
  Result := string.Join(
      IfThen(ASpace,