unit ESBMaths;

{:
	ESBMaths 3.2.1 - contains useful Mathematical routines for Delphi 4, 5 & 6.

	Copyright ©1997-2001 ESB Consultancy<p>

	These routines are used by ESB Consultancy within the
	development of their Customised Applications, and have been
	under Development since the early Turbo Pascal days. Many of the routines were developed
	for specific needs.<p>

	ESB Consultancy retains full copyright.<p>

	ESB Consultancy grants users of this code royalty free rights
	to do with this code as they wish.<p>

	ESB Consultancy makes no guarantees nor excepts any liabilities
	due to the use of these routines<p>

	We does ask that if this code helps you in you development
	that you send as an email mailto:glenn@esbconsult.com.au or even
	a local postcard. It would also be nice if you gave us a
	mention in your About Box or Help File.<p>

	ESB Consultancy Home Page: http://www.esbconsult.com.au<p>

	Mail Address: PO Box 2259, Boulder, WA 6449 AUSTRALIA<p>

	Check out our new ESB Professional Computation Suite with
	3000+ Routines and 80+ Components for Delphi 4, 5 & 6.<p>
	http://www.esbconsult.com.au/esbpcs.html<p>

	Also check out Marcel Martin's HIT at:<p>
	http://www.esbconsult.com.au/esbpcs-hit.html<p>
	Marcel has been helping out to optimise and improve routines.<p>

	Rory Daulton has generously donated and helped with many optimised
	routines. Our thanks to him as well.<p>

	Marcel van Brakel has also been very helpful has includes ESBMaths
	into the Jedi Collection. http://www.delphi-jedi.org/

	Any mistakes made are mine rather than Rory's or the Marcels'.<p>

	v3.2.1 14 June 2001<p>
		- PiE constant renamed to PiToE
		- Fixed D3 compatibility

	v3.2 18 May 2001<p>
		- Added Delphi 6 Support

	v3.1 17 April 2001<p>
		- Added Rory Daulton's LongWord Routines<p>
		- Added Rory Daulton's IGreatestPowerOf2 Routine<p>
		- Added Gamma Function Calculation<p>
		- Added Inverse Gamma Function Calculation<p>
		- Added LnGamma (Natural Log of Gamma) Function Calculation<p>
		- Added Beta Function Calculation<p>
		- Added Incomplete Beta Function Calculation<p>

	v3.0 8 October 2000<p>
		- Many improvements to Optimisation.<p>
		- Improved and Added Constants.<p>
		- Michael Schnell developed BitsHighest that makes BitsNeeded much
			more efficient.<p>
		- Set Default ESBTolerance much smaller.<p>
		- Added FloatIsNegative & FloatIsPositive.<p>
		- Swaps changed back to Pascal since for D4 and above it optimises to
			better code.<p>
		- Added Swaps for the Floats and Int64.<p>
		- Improved ISqrt from Rory Daulton.<p>
		- Added ILog2 from Rory Daulton.<p>
		- Factorial, Permutations & Binomial Coefficients now correctly handle
			up to N = 1754.<p>
		- Added ESBArcSec and ESBArcCosec.<p>
		- Power & Log Routines completely overhauled - many thanks to
			Rory Daulton.<p>.

	v2.3 13 November 1999<p>
		- very nice GCD Routine donated by Marcel Martin.<p>
		- LCM and Relative Prime Routines to go with above added.<p>
		- InPlace Addition & Subtraction of Matrices added.<p>
		- Multiplication of Matrices added.<p>
		- Transposing Matrices added.<p>

	V2.2 18 September 1999<p>
		- Now over 130 Routines<p>
		- Added Delphi 5 Support<p>
		- All Optional Dependence on Math Unit removed by our own routines<p>
		- Added a collection of Constants for faster calculations<p>
		- Added several Trigonometric Routines<p>
		- Added Hyperbolic Routines<p>
		- Added Various Logarithm and Power Routines<p>

	v2.1 13 July 1999 <p>
		- Now over 110 Routines<p>
		- Added ESBDigits and ESBBitsNeeded for Cardinals to ESBMaths<p>
		- Added ESBMagnitude for Reals to ESBMaths<p>
		- Added GrandMean to ESBMath2<p>
		- Added Matrix Operations to ESBMath2, just the start<p>

	v2.0 16 Apr 1999<p>
		- Delphi 4 specific routines moved to ESBMaths2<p>
		- Added Vector operations to ESBMath2<p>

	v1.1 22 Jan 1999 - Improved Delphi 4 support<p>
		- Added: SumLWArray, SumSqEArray, SumSqDiffEArray, SumXYEArray<p>
		- Added: IsPositiveEArray<p>
		- Added: Geometric and Harmonic Means<p>
		- Added ESBMean, Median and Mode<p>
		- Added Routines for Variances and Means for Populations & Samples<p>
		- Added Quartile Calculations<p>
		- Fixed a couple of minor bugs<p>
		- Fixed bug in XY2Polar reported by Wolfgang Werner<p>
		- Added Help File<p>

	v1.01 17 July 1998 - Added some improvements suggested by Rory Daulton<p>
		- Added Factorials using Extendeds<p>
		- Added Permutations using Extendeds<p>
		- Added Combinations using Extendeds<p>
		- Added SwapI32 to Swap Integers<p>
		- Added SwapC to Swap Cardinals<p>
		- Added Factorials computed using Extendeds<p>
		- Added Permutations computed using Extendeds<p>
		- Added Binomial Coefficients (Combinations) computed using Extendeds<p>

	v1.00 17 April 1998 - first public/freeware release<p>
}

interface
{$IFDEF VER120}
{$DEFINE D4andAbove}
{$ENDIF}

{$IFDEF VER125}
{$DEFINE D4andAbove}
{$ENDIF}

{$IFDEF VER130}
{$DEFINE D4andAbove}
{$ENDIF}

{$IFDEF VER140}
{$DEFINE D4andAbove}
{$ENDIF}

{$J-} // Constants from here are not assignable
const
	//: Smallest Magnitude Single Available.
	MinSingle: Single = 1.5e-45;
	//: Largest Magnitude Single Available.
	MaxSingle: Single = 3.4e+38;
	//: Smallest Magnitude Double Available.
	MinDouble: Double = 5.0e-324;
	//: Largest Magnitude Double Available.
	MaxDouble: Double = 1.7e+308;
	//: Smallest Magnitude Extended Available.
	MinExtended: Extended = 3.6e-4951;
	//: Largest Magnitude Extended Available.
	MaxExtended: Extended = 1.1e+4932;
	// Smallest Delphi Currency Value.
	MinCurrency: Currency = -922337203685477.5807;
	// Largest Delphi Currency Value.
	MaxCurrency: Currency = 922337203685477.5807;

{--- Mathematical Constants ---}
{--- Taken from Abramowitz & Stegun: "Handbook of Mathematical Functions" }

const
	//: Square Root of 2.
	Sqrt2: Extended	= 1.4142135623730950488;

	//: Square Root of 3.
	Sqrt3: Extended	= 1.7320508075688772935;

	//: Square Root of 5.
	Sqrt5: Extended	= 2.2360679774997896964;

	//: Square Root of 10.
	Sqrt10: Extended = 3.1622776601683793320;

	//: Square Root of Pi.
	SqrtPi: Extended = 1.77245385090551602729;

	//: Cube Root of 2.
	Cbrt2: Extended  = 1.2599210498948731648;

	//: Cube Root of 3.
	Cbrt3: Extended  = 1.4422495703074083823;

	//: Cube Root of 10.
	Cbrt10: Extended = 2.1544346900318837219;

	//: Cube Root of 100.
	Cbrt100: Extended = 4.6415888336127788924;

	//: Cube Root of Pi.
	CbrtPi: Extended = 1.4645918875615232630;

	//: Inverse of Square Root of 2.
	InvSqrt2: Extended = 0.70710678118654752440;

	//: Inverse of Square Root of 3.
	InvSqrt3: Extended = 0.57735026918962576451;

	//: Inverse of Square Root of 5.
	InvSqrt5: Extended = 0.44721359549995793928;

	//: Inverse of Square Root of Pi.
	InvSqrtPi: Extended = 0.56418958354775628695;

	//: Inverse of Cube Root of Pi.
	InvCbrtPi: Extended	= 0.68278406325529568147;

	//: Natural Constant.
	ESBe: Extended	 = 2.7182818284590452354;

	//: Square of Natural Constant.
	ESBe2: Extended = 7.3890560989306502272;

	//: Natural Constant raised to Pi.
	ESBePi: Extended = 23.140692632779269006;

	//: Natural Constant raised to Pi/2.
	ESBePiOn2: Extended = 4.8104773809653516555;

	//: Natural Constant raised to Pi/4.
	ESBePiOn4: Extended = 2.1932800507380154566;

	//: Natural Log of 2.
	Ln2: Extended	= 0.69314718055994530942;

	//: Natural Log of 10.
	Ln10: Extended = 2.30258509299404568402;

	//: Natural Log of Pi.
	LnPi: Extended = 1.14472988584940017414;

	//: Log to Base 2 of 10.
	Log10Base2 = 3.3219280948873623478;

	//: Log to Base 10 of 2.
	Log2Base10: Extended = 0.30102999566398119521;

	//: Log to Base 10 of 3.
	Log3Base10: Extended = 0.47712125471966243730;

	//: Log to Base 10 of Pi.
	LogPiBase10: Extended	= 0.4971498726941339;

	//: Log to Base 10 of Natural Constant.
	LogEBase10: Extended = 0.43429448190325182765;

	//: Accurate Pi Constant.
	ESBPi: Extended = 3.1415926535897932385;

	//: Inverse of Pi.
	InvPi: Extended = 3.1830988618379067154e-1;

	//: Two * Pi.
	TwoPi: Extended = 6.2831853071795864769;

	//: Three * Pi.
	ThreePi: Extended = 9.4247779607693797153;

	//: Square of Pi.
	Pi2: Extended	= 9.8696044010893586188;

	//: Pi raised to the Natural Constant.
	PiToE: Extended	= 22.459157718361045473;

	//: Half of Pi.
	PiOn2: Extended = 1.5707963267948966192;

	//: Third of Pi.
	PiOn3: Extended = 1.0471975511965977462;

	//: Quarter of Pi.
	PiOn4: Extended = 0.7853981633974483096;

	//: Three Halves of Pi.
	ThreePiOn2: Extended = 4.7123889803846898577;

	//: Four Thirds of Pi.
	FourPiOn3: Extended =  4.1887902047863909846;

	//: 2^63.
	TwoToPower63: Extended = 9223372036854775808.0;

	//: One Radian in Degrees.
	OneRadian: Extended	= 57.295779513082320877;

	//: One Degree in Radians.
	OneDegree: Extended	= 1.7453292519943295769E-2;

	//: One Minute in Radians.
	OneMinute: Extended	= 2.9088820866572159615E-4;

	//: One Second in Radians.
	OneSecond: Extended	= 4.8481368110953599359E-6;

	//: Gamma Constant.
	ESBGamma: Extended	= 0.57721566490153286061;

	//: Natural Log of the Square Root of (2 * Pi)
	LnRt2Pi: Extended = 9.189385332046727E-1;

{$IFNDEF D4andAbove}
type
	LongWord = Cardinal;
{$ENDIF}

type
	{: Used for a Bit List of 16 bits from 15 -> 0 }
	TBitList = Word;

var
	{: Tolerance used to decide when float close enough to zero.
		For MAXIMUM Tolerance use  MinExtended (3.6e-4951). }
	ESBTolerance: Extended = 5.0e-324; // MinDouble

{--- Integer Operations ---}

{$IFDEF D4andAbove}
{: Multiply two unsigned 32-bit integers, checking for overflow.
	EXCEPTIONS: EIntOverflow  if the product will not fit into a LongWord.<p>
	NOTES: 1. This is needed in Delphi 5.01, since an attempt to multiply two
		LongWords results in signed multiplication.<p>
	Developed By Rory Daulton - used with permission.
}
function UMul (const Num1, Num2: LongWord): LongWord;

{: Multiply two unsigned 32-bit integers then divide by 2**32 (ie
	return the upper DWORD of the product).
	EXCEPTIONS: None, since the answer always fits into 32 bits.<p>
	NOTES: This is useful in randomizing routines, to scale one integer by
		the other relative to 32 bits (reduce one integer according to the
		size of a second integer).<p>
	Developed By Rory Daulton - used with permission.
}
function UMulDiv2p32 (const Num1, Num2: LongWord): LongWord;


{: Multiply two unsigned 32-bit integers then divide by another unsigned 32-bit integer.
	EXCEPTIONS: EDivByZero  if the divisor is zero (see Note 1). <p>
	NOTES: 1. The product followed by division will cause an unflagged integer
		overflow if the answer will not fit into 32 bits. In this case, the
		returned value is the lower 32-bits of the true answer.  This will
		never be a problem if the two multipliers are below the divisor.<p>
		2. The result is truncated toward zero (as in most integer divides).<p>
		3. This is useful in randomizing routines, to scale one integer by
		the other relative to a third (reduce one integer according to
		the relative size of a second integer compared to a third integer).<p>
		4. The Win32 API has a similar function MulDiv(), but it works on
		signed integers and rounds the quotient.<p>
	Developed By Rory Daulton - used with permission.
}
function UMulDiv (const Num1, Num2, Divisor: LongWord): LongWord;

{: Multiply two unsigned 32-bit integers then take the modulus by another
	unsigned 32-bit integer.
	EXCEPTIONS: EDivByZero   if the divisor is zero (See Note 1).<p>
	NOTES: 1. There is never an overflow error, since the answer always fits
		into 32 bits.<p>
		2. This is useful in linear congruential generators.<p>
	Developed By Rory Daulton - used with permission.
}
function UMulMod (const Num1, Num2, Modulus: LongWord): LongWord;

{$ENDIF}

{: Returns True if X1 and X2 are within ESBTolerance of each other }
function SameFloat (const X1, X2: Extended): Boolean;

{: Returns True if X is within ESBTolerance of 0 }
function FloatIsZero (const X: Extended): Boolean;

{: Returns True if X is Positive, ie X > ESBTolerance.
}
function FloatIsPositive (const X: Extended): Boolean;

{: Returns True if X is Negative, ie X < -ESBTolerance.
}
function FloatIsNegative (const X: Extended): Boolean;

{: Increments a Byte up to Limit. If B >= Limit no increment occurs. }
procedure IncLim (var B: Byte; const Limit: Byte);

{: Increments a ShortInt up to Limit. If B >= Limit no increment occurs. }
procedure IncLimSI (var B: ShortInt; const Limit: ShortInt);

{: Increments a Word up to Limit. If B >= Limit no increment occurs. }
procedure IncLimW (var B: Word; const Limit: Word);

{: Increments an Integer up to Limit. If B >= Limit no increment occurs. }
procedure IncLimI (var B: Integer; const Limit: Integer);

{: Increments a LongInt up to Limit. If B >= Limit no increment occurs. }
procedure IncLimL (var B: LongInt; const Limit: LongInt);

{: Decrements a Byte down to Limit. If B <= Limit no increment occurs. BASM }
procedure DecLim (var B: Byte; const Limit: Byte);

{: Decrements a ShortInt down to Limit. If B <= Limit no increment occurs. }
procedure DecLimSI (var B: ShortInt; const Limit: ShortInt);

{: Decrements a Word down to Limit. If B <= Limit no increment occurs. }
procedure DecLimW (var B: Word; const Limit: Word);

{: Decrements an Integer down to Limit. If B <= Limit no increment occurs. }
procedure DecLimI (var B: Integer; const Limit: Integer);

{: Decrements a LongInt down to Limit. If B <= Limit no increment occurs. }
procedure DecLimL (var B: LongInt; const Limit: LongInt);

{: Returns the maximum value between two Bytes. BASM }
function MaxB (const B1, B2: Byte): Byte;

{: Returns the minimum value between two Bytes. BASM }
function MinB (const B1, B2: Byte): Byte;

{: Returns the maximum value between two ShortInts. }
function MaxSI (const B1, B2: ShortInt): ShortInt;

{: Returns the minimum value between two ShortInts. }
function MinSI (const B1, B2: ShortInt): ShortInt;

{: Returns the maximum value between two Words. BASM }
function MaxW (const B1, B2: Word): Word;

{: Returns the minimum value between two Words. BASM }
function MinW (const B1, B2: Word): Word;

{: Returns the maximum value between two Integers. }
function MaxI (const B1, B2: Integer): Integer;

{: Returns the minimum value between two Integers. }
function MinI (const B1, B2: Integer): Integer;

{: Returns the maximum value between two LongInts. }
function MaxL (const B1, B2: LongInt): LongInt;

{: Returns the minimum value between two LongInts. }
function MinL (const B1, B2: LongInt): LongInt;

{: Swap Two Bytes. }
procedure SwapB (var B1, B2: Byte);

{: Swap Two ShortInts. }
procedure SwapSI (var B1, B2: ShortInt); 

{: Swap Two Words. }
procedure SwapW (var B1, B2: Word);

{: Swap Two Integers. }
procedure SwapI (var B1, B2: SmallInt);

{: Swap Two LongInts. }
procedure SwapL (var B1, B2: LongInt);

{: Swap Two Integers (32-bit). }
procedure SwapI32 (var B1, B2: Integer);

{: Swap Two LongWords. }
procedure SwapC (var B1, B2: LongWord);

{$IFDEF D4andAbove}
{: Swap Two Int64's }
procedure SwapInt64 (var X, Y: Int64);
{$ENDIF}

{: Returns: <p>
	-1  if B < 0 <p>
	 0  if B = 0 <p>
	 1  if B > 0  BASM }
function Sign (const B: LongInt): ShortInt;

{: Returns the Maximum of 4 Words - BASM }
function Max4Word (const X1, X2, X3, X4: Word): Word;

{: Returns the Minimum of 4 Words - BASM }
function Min4Word (const X1, X2, X3, X4: Word): Word;

{: Returns the Maximum of 3 Words - BASM }
function Max3Word (const X1, X2, X3: Word): Word;

{: Returns the Minimum of 3 Words - BASM }
function Min3Word (const X1, X2, X3: Word): Word;

{: Returns the Maximum of an array of Bytes }
function MaxBArray (const B: array of Byte): Byte;

{: Returns the Maximum of an array of Words }
function MaxWArray (const B: array of Word): Word;

{: Returns the Maximum of an array of ShortInts }
function MaxSIArray (const B: array of ShortInt): ShortInt;

{: Returns the Maximum of an array of Integers }
function MaxIArray (const B: array of Integer): Integer;

{: Returns the Maximum of an array of LongInts }
function MaxLArray (const B: array of LongInt): LongInt;

{: Returns the Minimum of an array of Bytes }
function MinBArray (const B: array of Byte): Byte;

{: Returns the Minimum of an array of Words }
function MinWArray (const B: array of Word): Word;

{: Returns the Minimum of an array of ShortInts }
function MinSIArray (const B: array of ShortInt): ShortInt;

{: Returns the Minimum of an array of Integers }
function MinIArray (const B: array of Integer): Integer;

{: Returns the Minimum of an array of LongInts }
function MinLArray (const B: array of LongInt): LongInt;

{: Returns the Sum of an array of Bytes. All Operation in Bytes }
function SumBArray (const B: array of Byte): Byte;

{: Returns the Sum of an array of Bytes. All Operation in Words }
function SumBArray2 (const B: array of Byte): Word;

{: Returns the Sum of an array of ShortInts. All Operation in ShortInts }
function SumSIArray (const B: array of ShortInt): ShortInt;

{: Returns the Sum of an array of ShortInts. All Operation in Integers }
function SumSIArray2 (const B: array of ShortInt): Integer;

{: Returns the Sum of an array of Words. All Operation in Words }
function SumWArray (const B: array of Word): Word;

{: Returns the Sum of an array of Words. All Operation in Longints }
function SumWArray2 (const B: array of Word): LongInt;

{: Returns the Sum of an array of Integers. All Operation in Integers }
function SumIArray (const B: array of Integer): Integer;

{: Returns the Sum of an array of Longints. All Operation in Longints }
function SumLArray (const B: array of LongInt): LongInt;

{: Returns the Sum of an array of LongWord. All Operation in LongWords }
function SumLWArray (const B: array of LongWord): LongWord;

{: Returns the number of digits (Magnitude) of a Positive Integer}
function ESBDigits (const X: LongWord): Byte;

{: Returns the highest Bit set within a LongWord. Donated by Michael Schnell }
function BitsHighest (const X: LongWord): Integer;

{: Returns the number of Bits needed to represent a given positive integer}
function ESBBitsNeeded (const X: LongWord): Integer;

{: Returns the Greatest Common (Positive) Divisor (GCD)of two Integers. Also
	Refered to as the Highest Common Factor (HCF). Uses Euclid's Algorithm.
     BASM Routine donated by Marcel Martin }
function GCD (const X, Y: LongWord): LongWord;

{$IFDEF D4andAbove}
{: Returns the Least Common Multiple of two Integers. }
function LCM (const X, Y : LongInt): Int64;
{$ELSE}
{: Returns the Least Common Multiple of two Integers. }
function LCM (const X, Y : LongInt): LongInt;
{$ENDIF}

{: If two Integers are Relative Prime to each other then GCD (X, Y) = 1.
	CoPrime is another term for Relative Prime. Some interpretive
     problems may arise when '0' and/or '1' are used. }
function RelativePrime (const X, Y: LongWord): Boolean;

{--- Floating Point Operations ---}

{: Returns the 80x87 Control Word <p>
  15-12 Reserved  <p>
	On 8087/80287 12 was Infinity Control <p>
	 0 Projective <p>
	 1 Affine <p>
  11-10 Rounding Control <p>
    00 Round to nearest even <p>
    01 Round Down <p>
    10 Round Up <p>
    11 Chop - Truncate towards Zero <p>
  9-8  Precision Control <p>
    00 24 bits Single Precision <p>
    01 Reserved <p>
    10 53 bits Double Precision <p>
    11 64 bits Extended Precision (Default) <p>
  7-6  Reserved <p>
	On 8087 7 was Interrupt Enable Mask <p>
  5  Precesion Exception Mask <p>
  4  Underflow Exception Mask <p>
  3  Overflow Exception Mask <p>
  2  Zero Divide Exception Mask <p>
  1  Denormalised Operand Exception Mask <p>
  0  Invalid Operation Exception Mask <p>
  BASM }
function Get87ControlWord: TBitList;

{: Sets the 80x87 Control Word <p>
  15-12 Reserved <p>
	On 8087/80287 12 was Infinity Control <p>
	 0 Projective <p>
	 1 Affine <p>
  11-10 Rounding Control <p>
    00 Round to nearest even <p>
    01 Round Down <p>
    10 Round Up <p>
    11 Chop - Truncate towards Zero <p>
  9-8  Precision Control <p>
    00 24 bits Single Precision <p>
    01 Reserved <p>
    10 53 bits Double Precision <p>
    11 64 bits Extended Precision (Default) <p>
  7-6  Reserved <p>
	On 8087 7 was Interrupt Enable Mask <p>
  5  Precesion Exception Mask <p>
  4  Underflow Exception Mask <p>
  3  Overflow Exception Mask <p>
  2  Zero Divide Exception Mask <p>
  1  Denormalised Operand Exception Mask <p>
  0  Invalid Operation Exception Mask <p>
  BASM }
procedure Set87ControlWord (const CWord: TBitList);

{: Swap two Extendeds. }
procedure SwapExt (var X, Y: Extended);

{: Swap two Doubles. }
procedure SwapDbl (var X, Y: Double);

{: Swap two Singles. }
procedure SwapSing (var X, Y: Single);

{: Returns  <p>
	-1 if X < 0 <p>
	 0 if X = 0 <p>
	 1 if X > 0 <p>
 }
function Sgn (const X: Extended): ShortInt;

{: Returns the straight line Distance between (X1, Y1) and (X2, Y2) }
function Distance (const X1, Y1, X2, Y2: Extended): Extended;

{: Performs Floating Point Modulus. ExtMod := X - Floor ( X / Y ) * Y }
function ExtMod (const X, Y: Extended): Extended;

{: Performs Floating Point Remainder. ExtRem := X - Int ( X / Y ) * Y }
function ExtRem (const X, Y: Extended): Extended;

{: Returns X mod Y for Comp Data Types }
function CompMOD (const X, Y: Comp): Comp;

{: Converts Polar Co-ordinates into Cartesion Co-ordinates }
procedure Polar2XY (const Rho, Theta: Extended; var X, Y: Extended);

{: Converts Cartesian Co-ordinates to Polar Co-ordinates }
procedure XY2Polar (const X, Y: Extended; var Rho, Theta: Extended);

{: Converts Degrees/Minutes/Seconds into an Extended Real }
function DMS2Extended (const Degs, Mins, Secs: Extended): Extended;

{: Converts an Extended Real into Degrees/Minutes/Seconds }
procedure Extended2DMS (const X: Extended; var Degs, Mins, Secs: Extended);

{: Returns the Maximum of two Extended Reals }
function MaxExt (const X, Y: Extended): Extended;

{: Returns the Minimum of two Extended Reals }
function MinExt (const X, Y: Extended): Extended;

{: Returns the Maximum of an array of Extended Reals }
function MaxEArray (const B: array of Extended): Extended;

{: Returns the Minimum of an array of Extended Reals }
function MinEArray (const B: array of Extended): Extended;

{: Returns the Maximum of an array of Single Reals }
function MaxSArray (const B: array of Single): Single;

{: Returns the Maximum of an array of Single Reals }
function MinSArray (const B: array of Single): Single;

{: Returns the Maximum of an array of Comp Reals }
function MaxCArray (const B: array of Comp): Comp;

{: Returns the Minimum of an array of Comp Reals }
function MinCArray (const B: array of Comp): Comp;

{: Returns the Sum of an Array of Single Reals }
function SumSArray (const B: array of Single): Single;

{: Returns the Sum of an Array of Extended Reals }
function SumEArray (const B: array of Extended): Extended;

{: Returns the Sum of the Square of an Array of Extended Reals }
function SumSqEArray (const B: array of Extended): Extended;

{: Returns the Sum of the Square of the difference of
	an Array of Extended Reals from a given Value }
function SumSqDiffEArray (const B: array of Extended; Diff: Extended): Extended;

{: Returns the Sum of the Pairwise Product of two
	Arrays of Extended Reals }
function SumXYEArray (const X, Y: array of Extended): Extended;

{: Returns the Sum of an Array of Comp Reals }
function SumCArray (const B: array of Comp): Comp;

{: Returns A! i.e Factorial of A - only values up to 1754 are handled
	returns 0 if larger }
function FactorialX (A: LongWord): Extended;

{: Returns nPr i.e Permutation of r objects from n.
	Only values of N up to 1754 are handled	returns 0 if larger
	If R > N  then 0 is returned }
function PermutationX (N, R: LongWord): Extended;

{: Returns nCr i.e Combination of r objects from n.
	These are also known as the Binomial Coefficients
	Only values of N up to 1754 are handled	returns 0 if larger
	If R > N  then 0 is returned }

function BinomialCoeff (N, R: LongWord): Extended;

{: Returns True if all elements of X > ESBTolerance }
function IsPositiveEArray (const X: array of Extended): Boolean;

{: Returns the Geometric Mean of the values }
function GeometricMean (const X: array of Extended): Extended;

{: Returns the Harmonic Mean of the values }
function HarmonicMean (const X: array of Extended): Extended;

{: Returns the Arithmetic Mean of the Values }
function ESBMean (const X: array of Extended): Extended;

{: Returns the Variance of the Values, assuming a Sample.
	Square root this value to get Standard Deviation }
function SampleVariance (const X: array of Extended): Extended;

{: Returns the Variance of the Values, assuming a Population.
	Square root this value to get Standard Deviation }
function PopulationVariance (const X: array of Extended): Extended;

{: Returns the Mean and Variance of the Values, assuming a Sample.
	Square root the Variance to get Standard Deviation }
procedure SampleVarianceAndMean (const X: array of Extended;
	var Variance, Mean: Extended);

{: Returns the Mean and Variance of the Values, assuming a Population.
	Square root the Variance to get Standard Deviation }
procedure PopulationVarianceAndMean (const X: array of Extended;
	var Variance, Mean: Extended);

{: Returns the Median (2nd Quartiles) of the Values. The array
	MUST be sorted before using this operation }
function GetMedian (const SortedX: array of Extended): Extended;

{: Returns the Mode (most frequent) of the Values. The array
	MUST be sorted before using this operation. Function is False
	if no Mode exists }
function GetMode (const SortedX: array of Extended; var Mode: Extended): Boolean;

{: Returns the 1st and 3rd Quartiles - Median is 2nd Quartiles - of the Values.
	The array	MUST be sorted before using this operation }
procedure GetQuartiles (const SortedX: array of Extended; var Q1, Q3: Extended);

{: Returns the Magnitued of a given Value }
function ESBMagnitude (const X: Extended): Integer;

{-- Trigonometric Functions ---}

//: Returns Tangent of Angle given in Radians.
function ESBTan (Angle : Extended): Extended;

//: Returns CoTangent of the Angle given in Radians.
function ESBCot (Angle : Extended): Extended;

//: Returns CoSecant of the Angle given in Radians.
function ESBCosec (const Angle : Extended): Extended;

//: Returns Secant of the Angle given in Radians.
function ESBSec (const Angle : Extended): Extended;

//: Returns the ArcTangent of Y / X - Result is in Radians
function ESBArcTan (X, Y: Extended): Extended;

//: Fast Computation of Sin and Cos, where Angle is in Radians
procedure ESBSinCos (Angle: Extended; var SinX, CosX: Extended);

//: Given a Value returns the Angle whose Cosine it is, in Radians
function ESBArcCos (const X: Extended): Extended;

//: Given a Value returns the Angle whose Sine it is, in Radians
function ESBArcSin (const X: Extended): Extended;

{: Given a Value returns the Angle whose Secant it is, in Radians.
}
function ESBArcSec (const X: Extended): Extended;

{: Given a Value returns the Angle whose Cosecant it is, in Radians.
}
function ESBArcCosec (const X: Extended): Extended;

{-- Logarithm & Power Functions ---}

//: Returns Logarithm of X to Base 10
function ESBLog10 (const X: Extended): Extended;

//: Returns Logarithm of X to Base 2
function ESBLog2 (const X: Extended): Extended;

//: Returns Logarithm of X to Given Base
function ESBLogBase (const X, Base: Extended): Extended;

{: Calculate 2 to the given floating point power. Developed by Rory Daulton
	and used with permission. December 1998.<p>
	The algorithm used is to scale the power of the fractional part
	of X, using FPU commands. Although the FPU (Floating Point Unit) is used,
	the answer is exact for integral X, since the FSCALE FPU command is.<p>
	EOverflow Exception when X > Log2(MaxExtended) = 11356.5234062941439494
	(if there was no other FPU error condition, such as underflow or denormal,
	before entry to this routine)<p>
	EInvalidOp Exception on some occasions when EOverflow would be expected,
	due to some other FPU error condition (such as underflow) before entry to
	this routine.
}
function Pow2 (const X: Extended): Extended;

{: Calculate any float to non-negative integer power. Developed by Rory Daulton
	and used with permission. Last modified December 1998.<p>
}
function IntPow (const Base: Extended; const Exponent: LongWord): Extended;

{: Raises Values to an Integer Power. Thanks to Rory Daulton for improvements.
}
function ESBIntPower (const X: Extended; const N: LongInt): Extended; 

//: Returns X^Y - handles all cases
function XtoY (const X, Y: Extended): Extended;

//: Returns 10^Y - handles all cases
function TenToY (const Y: Extended): Extended;

//: Returns 2^Y - handles all cases
function TwoToY (const Y: Extended): Extended;

//: Returns log X using base Y - handles all valid cases
function LogXtoBaseY (const X, Y: Extended): Extended;

{: ISqrt (I) computes INT (SQRT (I)), that is, the integral part of the
	square root of integer I.
	Code  originally developed by Marcel Martin, used with permission.
	Rory Daulton introduced a faster routine (based on Marcel's) for most
	occassions and this is now used with Permission.
}
function ISqrt (const I: LongWord): Longword;

{: Calculate the integer part of the logarithm base 2 of an integer.
	Developed by Rory Daulton and used with Permission.<p>
	An Exception is raised if I is Zero.
}
function ILog2 (const I: LongWord): LongWord;

{: Calculate the greatest integral power of two that is less than or equal to
	the parameter.
	EXCEPTIONS:  EInvalidArgument for argument zero.<P>
	Developed By Rory Daulton - used with permission.
}
function IGreatestPowerOf2 (const N: LongWord): LongWord;

{--- Hyperbolic Functions ---}

//: Returns the inverse hyperbolic cosine of X
function ESBArCosh (X : Extended) : Extended;

//: Returns the inverse hyperbolic sine of X
function ESBArSinh (X : Extended) : Extended;

//: Returns the inverse hyperbolic tangent of X
function ESBArTanh (X : Extended) : Extended;

//: Returns the hyperbolic cosine of X
function ESBCosh (X : Extended) : Extended;

//: Returns the hyperbolic sine of X
function ESBSinh (X : Extended) : Extended;

//: Returns the hyperbolic tangent of X
function ESBTanh (X : Extended) : Extended;

{--- Special Functions ---}

{: Returns 1/Gamma(X) using a Series Expansion as defined in Abramowitz &
	Stegun. Defined for all values of X.<p>
	Accuracy: Gives about 15 digits.
}
function InverseGamma (const X: Extended): Extended;

{: Returns Gamma(X) using a Series Expansion for 1/Gamma (X) as defined in
	Abramowitz & Stegun. Defined for all values of X except negative
	integers and 0.<p>
	Accuracy: Gives about 15 digits.
}
function Gamma (const X: Extended): Extended;

{: Logarithm to base e of the gamma function.
	Accurate to about 1.e-14.
	Programmer: Alan Miller - developed for Fortan 77, converted with permission.
}
function LnGamma (const X: Extended): Extended;

{: Returns Beta(X,Y) using a Series Expansion for 1/Gamma (X) as defined in
	Abramowitz & Stegun. Defined for all values of X, Y except negative
	integers and 0.
	Accuracy: Gives about 15 digits.
}
function Beta (const X, Y: Extended): Extended;

{: Returns the Incomplete Beta Ix(P, Q), where 0 <= X <= 1 and
	P and Q are positive.
	Accuracy: Gives about 17 digits.<p>
	Adapted From Collected Algorithms from CACM,
	Algorithm 179 - Incomplete Beta Function Ratios,
	Oliver G Ludwig.
}
function IncompleteBeta (X: Extended; P, Q: Extended): Extended;

implementation

uses
	SysUtils;

procedure IncLim (var B: Byte; const Limit: Byte);
begin
	if B < Limit then
		Inc (B);
end;

procedure IncLimSI (var B: ShortInt; const Limit: ShortInt);
begin
	if B < Limit then
		Inc (B);
end;

procedure IncLimW (var B: Word; const Limit: Word);
begin
	if B < Limit then
		Inc (B);
end;

procedure IncLimI (var B: Integer; const Limit: Integer);
begin
	if B < Limit then
		Inc (B);
end;

procedure IncLimL (var B: LongInt; const Limit: LongInt);
begin
	if B < Limit then
		Inc (B);
end;

procedure DecLim (var B: Byte; const Limit: Byte);
begin
	if B > Limit then
		Dec (B);
end;

procedure DecLimSI (var B: ShortInt; const Limit: ShortInt);
begin
	if B > Limit then
		Dec (B);
end;

procedure DecLimW (var B: Word; const Limit: Word);
begin
	if B > Limit then
		Dec (B);
end;

procedure DecLimI (var B: Integer; const Limit: Integer);
begin
	if B > Limit then
		Dec (B);
end;

procedure DecLimL (var B: LongInt; const Limit: LongInt);
begin
	if B > Limit then
		Dec (B);
end;

function MaxB (const B1, B2: Byte): Byte;
begin
	if B1 > B2 then
		Result := B1
	 else
		Result := B2;
end;

function MinB (const B1, B2: Byte): Byte;
begin
	if B1 < B2 then
		Result := B1
	else
		Result := B2;
end;

function MaxSI (const B1, B2: ShortInt): ShortInt;
begin
	if B1 > B2 then
		Result := B1
	else
		Result := B2;
end;

function MinSI (const B1, B2: ShortInt): ShortInt;
begin
	if B1 < B2 then
		Result := B1
	else
		Result := B2;
end;

function MaxW (const B1, B2: Word): Word;
begin
	if B1 > B2 then
		Result := B1
	else
		Result := B2;
end;

function MinW (const B1, B2: Word): Word;
begin
	if B1 < B2 then
		Result := B1
	else
		Result := B2;
end;

function MaxI (const B1, B2: Integer): Integer;
begin
	if B1 > B2 then
		Result := B1
	else
		Result := B2;
end;

function MinI (const B1, B2: Integer): Integer;
begin
	if B1 < B2 then
		Result := B1
	else
		Result := B2;
end;

function MaxL (const B1, B2: LongInt): LongInt;
begin
	if B1 > B2 then
		Result := B1
	else
		Result := B2;
end;

function MinL (const B1, B2: LongInt): LongInt;
begin
	if B1 < B2 then
		Result := B1
	else
		Result := B2;
end;

procedure SwapB (var B1, B2: Byte);
var
	Temp: Byte;
begin
	Temp := B1;
	B1 := B2;
	B2 := Temp;
end;

procedure SwapSI (var B1, B2: ShortInt);
var
	Temp: ShortInt;
begin
	Temp := B1;
	B1 := B2;
	B2 := Temp;
end;

procedure SwapW (var B1, B2: Word);
var
	Temp: Word;
begin
	Temp := B1;
	B1 := B2;
	B2 := Temp;
end;
procedure SwapI (var B1, B2: SmallInt);
var
	Temp: SmallInt;
begin
	Temp := B1;
	B1 := B2;
	B2 := Temp;
end;

procedure SwapL (var B1, B2: LongInt);
var
	Temp: LongInt;
begin
	Temp := B1;
	B1 := B2;
	B2 := Temp;
end;

procedure SwapI32 (var B1, B2: Integer);
var
	Temp: Integer;
begin
	Temp := B1;
	B1 := B2;
	B2 := Temp;
end;

procedure SwapC (var B1, B2: LongWord);
var
	Temp: LongWord;
begin
	Temp := B1;
	B1 := B2;
	B2 := Temp;
end;

function Sign (const B: LongInt): ShortInt;
begin
	if B < 0 then
		Result := -1
	else if B = 0 then
		Result := 0
	else
		Result := 1;
end;

function Max4Word (const X1, X2, X3, X4: Word): Word;
begin
	Result := X1;
	if X2 > Result then
		Result := X2;
	if X3 > Result then
		Result := X3;
	if X4 > Result then
		Result := X4;
end;

function Min4Word (const X1, X2, X3, X4: Word): Word;
begin
	Result := X1;
	if X2 < Result then
		Result := X2;
	if X3 < Result then
		Result := X3;
	if X4 < Result then
		Result := X4;
end;

function Max3Word (const X1, X2, X3: Word): Word; 
begin
	Result := X1;
	if X2 > Result then
		Result := X2;
	if X3 > Result then
		Result := X3;
end;

function Min3Word (const X1, X2, X3: Word): Word; 
begin
	Result := X1;
	if X2 < Result then
		Result := X2;
	if X3 < Result then
		Result := X3;
end;

function MaxBArray (const B: array of Byte): Byte;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		if B [I] > Result then
			Result := B [I];
end;

function MaxWArray (const B: array of Word): Word;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		if B [I] > Result then
		Result := B [I];
end;

function MaxIArray (const B: array of Integer): Integer;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		if B [I] > Result then
			   Result := B [I];
end;

function MaxSIArray (const B: array of ShortInt): ShortInt;
var
	 I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
	if B [I] > Result then
		Result := B [I];
end;

function MaxLArray (const B: array of LongInt): LongInt;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] > Result then
			Result := B [I];
end;

function MinBArray (const B: array of Byte): Byte;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] < Result then
		   Result := B [I];
end;

function MinWArray (const B: array of Word): Word;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] < Result then
			  Result := B [I];
end;

function MinIArray (const B: array of Integer): Integer;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] < Result then
			 Result := B [I];
end;

function MinSIArray (const B: array of ShortInt): ShortInt;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] < Result then
			Result := B [I];
end;

function MinLArray (const B: array of LongInt): LongInt;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] < Result then
			Result := B [I];
end;

function ISqrt (const I: LongWord): LongWord;
const
	Estimates: array[0..31] of Word = (
	// for index i, constant := Trunc(Sqrt((Int64(2) shl i) - 1.0)), which is
	//   the largest possible ISqrt(n) for  2**i <= n < 2**(i+1)
	    1,     1,     2,     3,     5,     7,    11,    15,
	   22,    31,    45,    63,    90,   127,   181,   255,
	  362,   511,   724,  1023,  1448,  2047,  2896,  4095,
	 5792,  8191, 11585, 16383, 23170, 32767, 46340, 65535);
	// eax  // ebx  // ecx  // edx
asm  // entry:   // eax = I
  // calc the result quickly for zero or one (sqrt equals the argument)
	   cmp     eax, 1
	   jbe     @@end
  // save registers and the argument
	   push    ebx
	   mov     ebx, eax                // ebx = I
  // use the logarithm base 2 to load an initial estimate, which is greater
  //   than or equal to the actual value
	   bsr     eax, ebx        // eax = ILog2(I) (note upper WORD is now zero)
	   mov     ax, [word ptr Estimates + eax * 2]
						  // eax = X
  // repeat ...
@@repeat:
  // --  save the last estimate [ X ]
	   mov     ecx, eax                        // ecx = X
  // --  calc the new estimate [ (I/X + X) / 2 ; the Newton-Raphson formula ]
	   xor     edx, edx                                // edx = 0
	   mov     eax, ebx        // eax = I (so edx:eax = I)
	   div     ecx             // eax = I/X            // edx = I mod X
	   add     eax, ecx        // eax = I/X + X
	   shr     eax, 1          // eax = XNew = (I/X+X)/2
  // until the new estimate >= the last estimate
  //   [which can never happen in exact floating-point arithmetic, and can
  //   happen due to truncation only if the last estimate <= Sqrt(I) ]
	   cmp     eax, ecx
	   jb      @@repeat
  // use the next-to-last estimate as the result
	   mov     eax, ecx        // eax = X
  // restore registers
	   pop     ebx
@@end:              //exit:     // eax = Result
end {ISqrt};

function ILog2 (const I: LongWord): LongWord;
	procedure BadILog2;
	begin
		raise EMathError.Create('Division By Zero');
	end {BadILog2};
asm
    bsr     eax,eax
    jz      BadILog2
end {ILog2};

function SumBArray (const B: array of Byte): Byte;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumBArray2 (const B: array of Byte): Word;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumSIArray (const B: array of ShortInt): ShortInt;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumSIArray2 (const B: array of ShortInt): Integer;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumWArray (const B: array of Word): Word;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= 0 to High (B) do
		Result := Result + B [I];
end;

function SumWArray2 (const B: array of Word): LongInt;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumIArray (const B: array of Integer): Integer;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumLArray (const B: array of LongInt): LongInt;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumLWArray (const B: array of LongWord): LongWord;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I:= Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function Get87ControlWord: TBitList; 
var
	Temp: Word;
asm
	fstcw [Temp]  		{ Get '87 Control Word }
	mov  ax, [Temp] 	{ Leave in AX for function }
end;

procedure Set87ControlWord (const CWord: TBitList);
var
	Temp: Word;
asm
	mov   [Temp], ax
	fldcw [Temp] 		{ Load '87 Control Word }
end;

procedure Polar2XY (const Rho, Theta: Extended; var X, Y: Extended);
begin
	ESBSinCos (Theta, Y, X); { quite fast }
	X := Rho * X;
	Y := Rho * Y;
end;

procedure XY2Polar (const X, Y: Extended; var Rho, Theta: Extended);
begin
	Rho := Sqrt (Sqr (X) + Sqr (Y));
	if Abs (X) > ESBTolerance then
	begin
		Theta := ArcTan (abs (Y) / abs (X));
		if Sgn (X) = 1 then
		begin
			if Sgn (Y) = -1 then
				Theta := TwoPi - Theta
		end
		else
		begin
			if Sgn (Y) = 1 then
				Theta := ESBPi - Theta
			else
				Theta := ESBPi + Theta
		end;
	end
	else
		Theta := Sgn (Y) * Pion2
end;

function Sgn (const X: Extended): ShortInt;
begin
	if X < 0.0 then
		Result := -1
	else if X = 0.0 then
		Result := 0
	else
		Result := 1
end;

function DMS2Extended (const Degs, Mins, Secs: Extended): Extended;
begin
	Result := Degs + Mins / 60.0 + Secs / 3600.0
end;

procedure Extended2DMS (const X: Extended; var Degs, Mins, Secs: Extended);
var
	Y: Extended;
begin
	Degs := Int (X);
	Y := Frac (X) * 60;
	Mins := Int (Y);
	Secs := Frac (Y) * 60;
end;

function Distance (const X1, Y1, X2, Y2: Extended): Extended;
{ Rory Daulton suggested this more tolerant routine }
var
	X, Y: Extended;
begin
	X := Abs (X1 - X2);
	Y := Abs (Y1 - Y2);
	if X > Y then
		Result := X * Sqrt (1 + Sqr (Y / X))
	else if Y <> 0 then
		Result := Y * Sqrt (1 + Sqr (X / Y))
	else
		Result := 0
end;

function ExtMod (const X, Y: Extended): Extended;
var
	Z: Extended;
begin
	 Result := X / Y;
	 Z := Int (Result);
	 if Result < 0 then
		Z := Z - 1.0;
	 { Z now has Floor (X / Y) }
	 Result := X - Z * Y
end;

function ExtRem (const X, Y: Extended): Extended;
begin
	Result := X - Int (X / Y) * Y
end;

function MaxExt (const X, Y: Extended): Extended;
begin
	if X > Y then
		Result := X
	else
		Result := Y
end;

function MinExt (const X, Y: Extended): Extended;
begin
	if X < Y then
		Result := X
	else
		Result := Y
end;

function CompMOD (const X, Y: Comp): Comp;
begin
	Result := X - Y * Int (X / Y)
end;

function MaxEArray (const B: array of Extended): Extended;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] > Result then
			Result := B [I];
end;

function MinEArray (const B: array of Extended): Extended;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
	if B [I] < Result then
		Result := B [I];
end;

function MaxSArray (const B: array of Single): Single;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] > Result then
			Result := B [I];
end;

function MinSArray (const B: array of Single): Single;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] < Result then
			  Result := B [I];
end;

function MaxCArray (const B: array of Comp): Comp;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] > Result then
			Result := B [I];
end;

function MinCArray (const B: array of Comp): Comp;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		if B [I] < Result then
			Result := B [I];
end;

function SumSArray (const B: array of Single): Single;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumEArray (const B: array of Extended): Extended;
var
	I: Integer;
begin
	Result := B [Low (B)];
	for I := Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SumSqEArray (const B: array of Extended): Extended;
var
	I: Integer;
begin
	Result := Sqr (B [Low (B)]);
	for I := Low (B) + 1 to High (B) do
		Result := Result + Sqr (B [I]);
end;

function SumSqDiffEArray (const B: array of Extended; Diff: Extended): Extended;
var
	I: Integer;
begin
	Result := Sqr (B [Low (B)] - Diff);
	for I := Low (B) + 1 to High (B) do
		Result := Result + Sqr (B [I] - Diff);
end;

function SumXYEArray (const X, Y: array of Extended): Extended;
var
	I: Integer;
	M, N: Integer;
begin
	M := MaxL (Low (X), Low (Y));
	N := MinL (High (X), High (Y));
	Result := X [M] * Y [M];
	for I := M + 1 to N do
		Result := Result + X [I] * Y [I];
end;

function SumCArray (const B: array of Comp): Comp;
var
	I: Integer;
begin
	Result := Low (B);
	for I := Low (B) + 1 to High (B) do
		Result := Result + B [I];
end;

function SameFloat (const X1, X2: Extended): Boolean;
begin
	Result := abs (X1 - X2) < ESBTolerance
end;

function FloatIsZero (const X: Extended): Boolean;
begin
	Result := abs (X) < ESBTolerance
end;

function FloatIsPositive (const X: Extended): Boolean;
begin
	Result := (X >= ESBTolerance);
end;

function FloatIsNegative (const X: Extended): Boolean;
begin
	Result := (X <= -ESBTolerance);
end;

function FactorialX (A: LongWord): Extended;
var
	I: Integer;
begin
	if A  > 1754 then
	begin
		Result := 0.0;
		Exit;
	end;
	Result := 1.0;
	for I := 2 to A do
		Result := Result * I;
end;

function PermutationX (N, R: LongWord): Extended;
var
	I : Integer;
begin
	if (N = 0) or (R > N) or (N > 1754) then
	begin
		Result := 0.0;
		Exit;
	end;
	Result := 1.0;
	if (R = 0) then
		Exit;
	try
		for I := N downto N - R + 1 do
			Result := Result * I;
		Result := Int (Result + 0.5);
	except
		Result := -1.0
	end;
end;

function BinomialCoeff (N, R: LongWord): Extended;
var
	I: Integer;
	K: LongWord;
begin
	if (N = 0) or (R > N) or (N > 1754) then
	begin
		Result := 0.0;
		Exit;
	end;
	Result := 1.0;
	if (R = 0) or (R = N) then
		Exit;
	if R > N div 2 then
		R := N - R;
	K := 2;
	try
		for I := N - R + 1 to N do
		begin
			Result := Result * I;
			if K <= R then
			begin
				Result := Result / K;
				Inc (K);
			end;
		end;
		Result := Int (Result + 0.5);
	except
		Result := -1.0
	end;
end;

function IsPositiveEArray (const X: array of Extended): Boolean;
var
	I: Integer;
begin
	Result := False;
	for I := 0 to High (X) do
		if X [I] <= ESBTolerance then
			Exit;
	Result := True;
end;

function GeometricMean (const X: array of Extended): Extended;
var
	I: Integer;
begin
	if High (X) < 0 then
		raise Exception.Create ('Array is Empty!')
	else if not IsPositiveEArray (X) then
		raise Exception.Create ('Array contains values <= 0!')
	else
	begin
		Result := 1;
		for I := 0 to High (X) do
			Result := Result * X [I];
		Result := XtoY (Result, 1 / (High (X) + 1));
	end;
end;

function HarmonicMean (const X: array of Extended): Extended;
var
	I: Integer;
begin
	if High (X) < 0 then
		raise Exception.Create ('Array is Empty!')
	else if not IsPositiveEArray (X) then
		raise Exception.Create ('Array contains values <= 0!')
	else
	begin
		Result := 0;
		for I := 0 to High (X) do
			Result := Result + 1 / X [I];
		Result := Result / (High (X) + 1);
		Result := 1 / Result;
	end;
end;

function ESBMean (const X: array of Extended): Extended;
begin
	Result := SumEArray (X) / (High (X) - Low (X) + 1)
end;

function SampleVariance (const X: array of Extended): Extended;
var
	I: Integer;
	SumSq: Extended;
	Mean: Extended;
begin
	Mean := ESBMean (X);
	SumSq := 0.0;
	for I := Low (X) to High (X) do
		SumSq := SumSq + Sqr (X [I] - Mean);
	Result := SumSq / (High (X) - Low (X))
end;

function PopulationVariance (const X: array of Extended): Extended;
var
	I: Integer;
	SumSq: Extended;
	Mean: Extended;
begin
	Mean := ESBMean (X);
	SumSq := 0.0;
	for I := Low (X) to High (X) do
		SumSq := SumSq + Sqr (X [I] - Mean);
	Result := SumSq / (High (X) - Low (X) + 1)
end;

procedure SampleVarianceAndMean (const X: array of Extended;
	var Variance, Mean: Extended);
var
	I: Integer;
	SumSq: Extended;
begin
	Mean := ESBMean (X);
	SumSq := 0.0;
	for I := Low (X) to High (X) do
		SumSq := SumSq + Sqr (X [I] - Mean);
	if High (X) > Low (X) then
		Variance := SumSq / (High (X) - Low (X))
	else
		Variance := 0;
end;

procedure PopulationVarianceAndMean (const X: array of Extended;
	var Variance, Mean: Extended);
var
	I: Integer;
	SumSq: Extended;
begin
	Mean := ESBMean (X);
	SumSq := 0.0;
	for I := Low (X) to High (X) do
		SumSq := SumSq + Sqr (X [I] - Mean);
	Variance := SumSq / (High (X) - Low (X) + 1)
end;

function GetMedian (const SortedX: array of Extended): Extended;
var
	N: Integer;
begin
	N := High (SortedX) + 1;
	if N <= 0 then
		raise Exception.Create ('Array is Empty!')
	else if N = 1 then
		Result := SortedX [0]
	else if Odd (N) then
		Result := SortedX [N div 2]
	else
		Result := (SortedX [N div 2 - 1] + SortedX [N div 2]) / 2;
end;

function GetMode (const SortedX: array of Extended; var Mode: Extended): Boolean;
var
	I, Freq, HiFreq: Integer;
	Matched: Boolean;
begin
	if High (SortedX) < 0 then
	begin
		raise Exception.Create ('Array is Empty!')
	end
	else if High (SortedX) = 0 then
	begin
		Mode := SortedX [0];
		Result := True;
	end
	else
	begin
		Mode := -MaxExtended;
		Freq := 1;
		HiFreq := 0;
		Matched := False;
		for I := 1 to High (SortedX) do
		begin
			if SameFloat (SortedX [I - 1], SortedX [I]) then
				Inc (Freq)
			else
			begin
				if Freq <> 1 then
				begin
					if Freq = HiFreq then
						Matched := True
					else if Freq > HiFreq then
					begin
						Mode := SortedX [I - 1];
						HiFreq := Freq;
						Matched := False;
					end;
					Freq := 1;
				end;
			end;
		end;
		if HiFreq > 0 then
		begin
			if Freq = HiFreq then
				Matched := True
			else if Freq > HiFreq then
			begin
				Mode := SortedX [High (SortedX)];
				Matched := False;
			end;
		end
		else if Freq > 1 then
		begin
			HiFreq := Freq;
			Mode := SortedX [High (SortedX)];
			Matched := False;
		end;
		Result := (HiFreq > 0) and not Matched;
	end;
end;

procedure GetQuartiles (const SortedX: array of Extended; var Q1, Q3: Extended);
var
	N: Single;
	I: Integer;
begin
	if High (SortedX) < 0 then
		raise Exception.Create ('Array is Empty!')
	else if High (SortedX) = 0 then
	begin
		Q1 := SortedX [0];
		Q3 := SortedX [0];
	end
	else
	begin
		N := (High (SortedX) + 1) / 4 + 0.5;
		I := Trunc (N);
		N := Frac (N);
		if I - 1 < High (SortedX) then
			Q1 := SortedX [I - 1] + (SortedX [I] - SortedX [I - 1]) * N
		else
			Q1 := SortedX [I - 1];

		N := 3 * (High (SortedX) + 1) / 4 + 0.5;
		I := Trunc (N);
		N := Frac (N);
		if I - 1 < High (SortedX) then
			Q3 := SortedX [I - 1] + (SortedX [I] - SortedX [I - 1]) * N
		else
			Q3 := SortedX [I - 1];
	end;
end;

function ESBDigits (const X: LongWord): Byte;
// Returns the number of digits in a given positive integer
begin
	if X = 0 then
		Result := 0
	else
		Result := Trunc (ESBLog10 (X)) + 1;
end;

function ESBMagnitude (const X: Extended): Integer;
// Returns the number of digits in a given positive integer
begin
	if X = 0 then
		Result := 0
	else
	begin
		if X >= 1 then
			Result := Trunc (ESBLog10 (Abs (X))) + 1
		else
			Result := Trunc (ESBLog10 (Abs (X)) - 0.00000000001) - 1
	end;
end;

function BitsHighest (const X: LongWord): Integer;
asm
	   MOV     ECX, EAX
	   MOV     EAX, -1
	   BSR     EAX, ECX
end;

function ESBBitsNeeded (const X: LongWord): Integer;
begin
	Result := BitsHighest (X) + 1;
end;


{--- Trigonometric Functions ---}

function ESBTan (Angle : Extended): Extended;
	function FTan (Angle : Extended): Extended;
	asm
		fld		[Angle]	// St(0) <- Angle
		ffree     st(7)	// Ensure st(7) is free
		fptan		     // St(1) <- Tan (Angle), St(0) <- 1
		fstp		st(0)	// Dispose of 1
		fwait
	end;
begin
	if abs (Angle) >= TwoToPower63 then // must be less then 2^63
		raise EMathError.Create ('Angle Magnitude too large for ESBTan')
	else
		Result := FTan (Angle);
end;

function ESBCot (Angle : Extended): Extended;
	function FCot (Angle : Extended): Extended;
	asm
		fld		[Angle]	// St(0) <- Angle
		ffree     st(7)	// Ensure st(7) is free
		fptan		     // St(1) <- Tan (Angle), St(0) <- 1
		fdivrp			// St(0) <- St(0)/St(1) which is Cot
		fwait
	end;
begin
	if abs (Angle) >= TwoToPower63 then // must be less then 2^63
		raise EMathError.Create ('Angle Magnitude too large for ESBCot')
	else
		Result := FCot (Angle);
end;

function ESBArcTan (X, Y: Extended): Extended;
	function FArcTan (X, Y: Extended): Extended;
	asm
		fld		[Y]		// St(0) <- Y
		fld	     [X]		// St(0) <- X, St (1) <- Y
		fpatan			// St(0) <- ArcTan (Y/X)
		fwait
	end;
begin
	if X = 0 then
     begin
     	if Y = 0 then
          	Result := 0
          else if Y > 0 then
          	Result := PiOn2
          else
          	Result := - PiOn2
     end
     else if Y = 0 then
	begin
		if X > 0 then
			Result := 0
		else
			Result := ESBPi
	end
	else
		Result := FArcTan (X, Y);
end;

procedure ESBSinCos (Angle: Extended; var SinX, CosX: Extended);
	procedure FSinCos (Angle: Extended; var SinX, CosX: Extended);
	asm
		fld		[Angle]			// St(0) <- Angle
		fsincos
		fstp		tbyte ptr [edx]    // St(0) -> CosX
		fstp    	tbyte ptr [eax]    // St(0) -> SinX
		fwait
	end;
begin
	if abs (Angle) >= TwoToPower63 then // must be less then 2^63
		raise EMathError.Create ('Angle Magnitude too large for ESBSinCos')
	else
		FSinCos (Angle, SinX, CosX);
end;

function ESBArcCos (const X: Extended): Extended;
var
	Y: Extended;
begin
	if abs (X) > 1 then
		raise EMathError.Create ('ESBArcCos requires values with magnitude <= 1')
	else
	begin
		Y := Sqrt (1 - Sqr (X));
		if abs (X) <= Y then
			Result := PiOn2 - ESBArcTan (Y, X)
		else	if X < 0 then
			Result := ESBArcTan (X, Y) + ESBPi
		else
			Result := ESBArcTan (X, Y)
	end;
end;

function ESBArcSin (const X: Extended): Extended;
var
	Y: Extended;
begin
	if abs (X) > 1 then
		raise EMathError.Create ('ESBArcSin requires values with magnitude <= 1')
	else
	begin
		Y := Sqrt (1 - Sqr (X));
		if abs (X) <= Y then
			Result :=  ESBArcTan (Y, X)
		else
			Result := Sgn (X) * PiOn2 - ESBArcTan (X, Y)
	end;
end;

function ESBCosec (const Angle: Extended): Extended;
var
	Y: Extended;
begin
	Y := Sin (Angle);
	if Y = 0 then
		raise EMathError.Create ('ESBCosec is Undefined')
	else
		Result := 1 / Y;
end;

function ESBSec (const Angle: Extended): Extended;
var
	Y: Extended;
begin
	Y := Cos (Angle);
	if Y = 0 then
		raise EMathError.Create ('ESBSec is Undefined')
	else
		Result := 1 / Y;
end;

function ESBArcSec (const X: Extended): Extended;
begin
	if FloatIsZero (X) then
		raise EMathError.Create ('Divide By Zero');
	Result := ESBArcCos (1 / Abs (X));
	if X < 0 then
		Result := Result + ESBPI;
end;

function ESBArcCosec (const X: Extended): Extended;
begin
	if FloatIsZero (X) then
		raise EMathError.Create ('Divide By Zero');
	Result := ESBArcSin (1 / Abs (X));
	if X < 0 then
		Result := Result + ESBPI;
end;

{--- Logarithm & Power Functions ---}

{------------------------------------------------------------------------------
PURPOSE:     Calculate 2 to the given floating point power.
AUTHOR:      Rory Daulton
DATE:        December 1998
PARAMETERS:  X: power to take of 2 [so the result is 2**X]
EXCEPTIONS:  EOverflow  when X > Log2(MaxExtended) = 11356.5234062941439494
				    (if there was no other FPU error condition, such as
				    underflow or denormal, before entry to this routine)
		   EInvalidOp on some occasions when EOverflow would be expected, due
				    to some other FPU error condition (such as underflow)
				    before entry to this routine
NOTES:       1. The algorithm used is to scale the power of the fractional part
			 of X, using FPU commands.
		   2. Although the FPU (Floating Point Unit) is used, the answer is
			 exact for integral X, since the FSCALE FPU command is.
------------------------------------------------------------------------------}
function Pow2 (const X: Extended): Extended;
asm
	   fld     X
// find Round(X)
	   fld     st
	   frndint
// find _Frac(X) [minimal fractional part of X, between -0.5 and 0.5]
	   fsub    st(1),st
	   fxch    st(1)
// Find 2**_Frac(X)
	   f2xm1
	   fld1
	   fadd
// Result := 2**_Frac(X) * 2**Round(X)
	   fscale
	   fstp    st(1)
	   fwait
end {Pow2};

function IntPow (const Base: Extended; const Exponent: LongWord): Extended;
  { Heart of Rory Daulton's IntPower: assumes valid parameters &
	non-negative exponent }
asm
		fld1                      { Result := 1 }
		cmp     eax, 0            { eax := Exponent }
		jz      @@3
		fld     Base
		jmp     @@2
  @@1:    fmul    ST, ST            { X := Base * Base }
  @@2:    shr     eax,1
		jnc     @@1
		fmul    ST(1),ST          { Result := Result * X }
		jnz     @@1
		fstp    st                { pop X from FPU stack }
  @@3:
		fwait
  end {P};


function ESBLog10 (const X: Extended): Extended;
	function FLog10 (X: Extended): Extended;
	asm
		fldlg2			// St(0) <- Log10 of 2
		fld		[X]		// St(0) <- X, St(1) <- Log10 of 2
		fyl2x			// St(0) <- log10 (2) * Log2 (X)
		fwait
	end;
	function AltFLog10 (X: Extended): Extended;
	asm
		fldlg2			// St(0) <- Log10 of 2
		fld		[X]		// St(0) <- X, St(1) <- Log10 of 2
		fyl2xp1			// St(0) <- log10 (2) * Log2 (X+1)
		fwait
	end;
begin
	if not FloatIsPositive (X) then // must be Positive
		raise EMathError.Create ('Value must be > 0')
	else if abs (X - 1) < 0.1 then
		Result := AltFLog10 (X - 1)
	else
		Result := FLog10 (X);
end;

function ESBLog2 (const X: Extended): Extended;
	function FLog2 (X: Extended): Extended;
	asm
		fld1				// St(0) <- 1
		fld		[X]		// St(0) <- X, St(1) <-1
		fyl2x			// St(0) <- 1 * Log2 (X)
		fwait
	end;
	function AltFLog2 (X: Extended): Extended;
	asm
		fld1				// St(0) <- 1
		fld		[X]		// St(0) <- X, St(1) <-1
		fyl2xp1			// St(0) <- 1 * Log2 (X+1)
		fwait
	end;
begin
	if not FloatIsPositive (X) then // must be Positive
		raise EMathError.Create ('Value must be > 0')
	else if abs (X - 1) < 0.1 then
		Result := AltFLog2 (X - 1)
	else
		Result := FLog2 (X);
end;

function ESBLogBase (const X, Base: Extended): Extended;
begin
	if not FloatIsPositive (X) then // must be Positive
		raise EMathError.Create ('Value must be > 0')
	else if not FloatIsPositive (Base) then // must be Positive
		raise EMathError.Create ('Value must be > 0')
	else
		Result := ESBLog2 (X) / ESBLog2 (Base);
end;

function LogXtoBaseY (const X, Y: Extended): Extended;
begin
	Result := ESBLog2 (X) / ESBLog2 (Y);
end;

function ESBIntPower (const X: Extended; const N: LongInt): Extended;
var
	P: LongWord;
begin
	if N = 0 then
		Result := 1
	else if (X = 0) then
	begin
		if N < 0 then
			raise EMathError.Create ('Zero cannot be raised to a Negative Power')
		else
			Result := 0
	end
	else if (X = 1) then
		Result := 1
	else if X = -1 then
	begin
		if Odd (N) then
			Result := -1
		else
			Result := 1
	end
	else if N > 0 then
		Result := IntPow (X, N)
	else
	begin
		if N <> Low (LongInt) then
			P := abs (N)
		else
			{$IFDEF D4andAbove}
			P := LongWord (High (LongInt)) + 1;
			{$Else}
			raise EMathError.Create ('Cannot be raised to Low (LongInt)');
			{$ENDIF}
		try
			Result := IntPow (X, P);
		except
			on EMathError do
			begin
				Result := IntPow (1 / X, P); { try again with another method, }
				Exit;                           {   perhaps less precise         }
			end {on};
		end {try};
		Result := 1 / Result;
	end;
end;

function XtoY (const X, Y: Extended): Extended;
	function PowerAbs: Extended; // Routine developed by Rory Daulton
	var
		ExponentPow2: Extended; // equivalent exponent to power 2
	begin
		try
			ExponentPow2 := ESBLog2 (Abs (X)) * Y;
		except
			on EMathError do
			   // allow underflow, when ExponentPow2 would have been negative
				if (Abs (X) > 1) <> (Y > 0) then
				begin
					Result := 0;
					Exit;
				end {if}
				else
					raise;
		end {try};
		Result := Pow2 (ExponentPow2);
    end;
begin
	if FloatIsZero (Y) then
		Result := 1
	else if FloatIsZero (X) then
	begin
		if Y < 0 then
			raise EMathError.Create ('Zero cannot be raised to a Negative Power')
		else
			Result := 0
	end
	else if FloatIsZero (Frac (Y)) then
	begin
		if (Y >= Low (LongInt)) and (Y <= High (LongInt)) then
			Result := ESBIntPower (X, LongInt (Round (Y)))
		else
		begin
			if (X > 0) or FloatIsZero (Frac (Y / 2.0)) then
				Result := PowerAbs
			else
				Result := -PowerAbs
		end;
	end
	else if X > 0 then
		Result := PowerAbs
	else
		raise EMathError.Create ('Invalid X^Y')
end;

function TenToY (const Y: Extended): Extended;
begin
	if FloatIsZero (Y) then
		Result := 1
	else if Y < -3.58e4931 then
		Result := 0
	else
		Result := Pow2 (Y * Log10Base2)
end;

function TwoToY (const Y: Extended): Extended;
begin
	if FloatIsZero (Y) then
		Result := 1
	else
		Result := Pow2 (Y)
end;

{--- Hyperbolic Functions ---}

function ESBArCosh (X: Extended): Extended;
begin
	Result := Ln (X + Sqrt (Sqr (X) - 1));
end;

function ESBArSinh (X: Extended): Extended;
begin
	Result := Ln (X + Sqrt (Sqr (X) + 1));
end;

function ESBArTanh (X: Extended): Extended;
var
	Y: Extended;
begin
	if X = 1 then
		raise EMathError.Create ('ESBArTanh is not Defined for X = 1')
	else
	begin
		Y :=(1 + X) / (1 - X);
		if Y <= 0 then
			raise EMathError.Create ('ESBArTanh is not Defined for that Value of X')
		else
			Result := Ln (Y) / 2;
	end;
end;

function ESBCosh (X: Extended): Extended;
begin
	Result := (Exp (X) + Exp (-X)) / 2.0;
end;

function ESBSinh (X: Extended): Extended;
begin
	Result := (Exp (X) - Exp (-X)) / 2.0;
end;

function ESBTanh (X: Extended): Extended;
var
	Y: Extended;
begin
	Y := ESBCosh (X);
	if Y = 0 then
		raise EMathError.Create ('X has a Cosh of 0, illegal for ESBTanh')
	else
		Result := ESBSinh (X) / Y;
end;

function GCD (const X, Y: LongWord): LongWord;
asm
     jmp   @01      // We start with EAX <- X, EDX <- Y, and check to see if Y = 0
@00: mov   ecx, edx	// ECX <- EDX prepare for division
     xor   edx, edx // clear EDX for Division
     div   ecx		// EAX <- EDX:EAX div ECX, EDX <- EDX:EAX mod ECX
     mov   eax, ecx	// EAX <- ECX, and repeat if EDX <> 0
@01: and   edx, edx // test to see if EDX is zero, without changing EDX
     jnz   @00		// when EDX is zero EAX has the result
end;

{$IFDEF D4andAbove}
function LCM (const X, Y : LongInt): Int64;
begin
	Result:= (X div LongInt (GCD (Abs (X), Abs (Y)))) * Int64 (Y);
end;
{$ELSE}
function LCM (const X, Y : LongInt): LongInt;
begin
	Result:= (X div LongInt (GCD (Abs (X), Abs (Y)))) * Y;
end;
{$ENDIF}

function RelativePrime (const X, Y: LongWord): Boolean;
begin
	Result := GCD (X, Y) = 1;
end;

{$IFDEF D4andAbove}
procedure SwapInt64 (var X, Y: Int64);
var
	Temp: Int64;
begin
	Temp := X;
	X := Y;
	Y := Temp;
end;
{$ENDIF}

procedure SwapExt (var X, Y: Extended);
var
	Temp: Extended;
begin
	Temp := X;
	X := Y;
	Y := Temp;
end;

procedure SwapDbl (var X, Y: Double);
var
	Temp: Double;
begin
	Temp := X;
	X := Y;
	Y := Temp;
end;

procedure SwapSing (var X, Y: Single);
var
	Temp: Single;
begin
	Temp := X;
	X := Y;
	Y := Temp;
end;

{$IFDEF D4andAbove}
{---------------------------------------------------------------------------
---
PURPOSE: Show an unsigned-integer overflow error with an exception
NOTES:   This is an internal routine, not declared in the interface section
----------------------------------------------------------------------------
--}
procedure FlagUIntOverflow;
begin
	raise EIntOverflow.Create('Unsigned integer overflow in RDMath');
end {FlagUIntOverflow};


{---------------------------------------------------------------------------
---
PURPOSE:     Multiply two unsigned 32-bit integers, checking for overflow
DATE:        March 2001
EXCEPTIONS:  EIntOverflow  if the product will not fit into a LongWord
NOTES:       1. This is needed in Delphi 5.01, since an attempt to multiply two
			 LongWords results in signed multiplication.
----------------------------------------------------------------------------
--}
function UMul (const Num1, Num2: LongWord): LongWord;
asm
	// Num1 is in eax, Num2 is in edx
	mul     edx               // product now in edx:eax
	jc      FlagUIntOverflow  // error on overflow
end {UMul};


{---------------------------------------------------------------------------
---
PURPOSE:     Multiply two unsigned 32-bit integers then divide by 2**32 (i.e.
		   return the upper DWORD of the product)
DATE:        Modified 13 Apr 2001 to use LongWords instead of Integers
EXCEPTIONS:  None, since the answer always fits into 32 bits
NOTES:       This is useful in randomizing routines, to scale one integer by
		   the other relative to 32 bits (reduce one integer according to the
		   size of a second integer).
----------------------------------------------------------------------------
--}
function UMulDiv2p32 (const Num1, Num2: LongWord): LongWord;
asm
	// Num1 is in eax, Num2 is in edx
	mul     edx           // result is in edx:eax
	mov     eax,edx       // return high dword
end {UMulDiv2p32};


{---------------------------------------------------------------------------
---
PURPOSE:     Multiply two unsigned 32-bit integers then divide by another
		   unsigned 32-bit integer.
DATE:        Created December 1998
		   Modified 13 Apr 2001 to use LongWords instead of Integers
EXCEPTIONS:  EDivByZero  if the divisor is zero
		   (see Note 1)
NOTES:       1. The product followed by division will cause an unflagged
			 integer overflow if the answer will not fit into 32 bits. In
			 this case, the returned value is the lower 32-bits of the true
			 answer.  This will never be a problem if the two multipliers
			 are below the divisor.
		   2. The result is truncated toward zero (as in most integer
			 divides).
		   3. This is useful in randomizing routines, to scale one integer by
			 the other relative to a third (reduce one integer according to
			 the relative size of a second integer compared to a third
			 integer).
		   4. The Win32 API has a similar function MulDiv(), but it works on
			 signed integers and rounds the quotient.
----------------------------------------------------------------------------
--}
function UMulDiv (const Num1, Num2, Divisor: LongWord): LongWord;
asm
	// Num1 is in eax, Num2 is in edx, Divisor is in ecx
	mul     edx           // product now in edx:eax
	div     ecx           // quotient now in eax, remainder in edx
end {UMulMod};


{---------------------------------------------------------------------------
---
PURPOSE:     Multiply two unsigned 32-bit integers then take the modulus by
		   another unsigned 32-bit integer.
DATE:        Created December 1998
		   Modified 13 Apr 2001 to use LongWords instead of Integers
EXCEPTIONS:  EDivByZero   if the divisor is zero
		   (See Note 1)
NOTES:       1. There is never an overflow error, since the answer always fits
			 into 32 bits.
		   2. This is useful in linear congruential generators.
----------------------------------------------------------------------------
--}
function UMulMod (const Num1, Num2, Modulus: LongWord): LongWord;
asm
	// Num1 is in eax, Num2 is in edx, Modulus is in ecx
	mul     edx           // product now in edx:eax
	div     ecx           // quotient now in eax, remainder in edx
	mov     eax,edx       // return the remainder
end {UMulMod};
{$ENDIF}

{--- Special Functions ---}

function InverseGamma (const X: Extended): Extended;
var
	C: array [1..26] of Extended;
	Z: Extended;
	XF: Extended;
	I: Integer;
begin
	C [1]  :=  1;
	C [2]  :=  0.5772156649015329;
	C [3]  := -0.6558780715202538;
	C [4]  := -0.0420026350340952;
	C [5]  :=  0.1665386113822915;
	C [6]  := -0.0421977345555443;
	C [7]  := -0.0096219715278770;
	C [8]  :=  0.0072189432466630;
	C [9]  := -0.0011651675918591;
	C [10] := -0.0002152416741149;
	C [11] :=  0.0001280502823882;
	C [12] := -0.0000201348547807;
	C [13] := -0.0000012504934821;
	C [14] :=  0.0000011330272320;
	C [15] := -0.0000002056338417;
	C [16] :=  0.0000000061160950;
	C [17] :=  0.0000000050020075;
	C [18] := -0.0000000011812746;
	C [19] :=  0.0000000001043427;
	C [20] :=  0.0000000000077823;
	C [21] := -0.0000000000036968;
	C [22] :=  0.0000000000005100;
	C [23] := -0.0000000000000206;
	C [24] := -0.0000000000000054;
	C [25] :=  0.0000000000000014;
	C [26] :=  0.0000000000000001;
	Result := 0;
	Z := 1;
	XF := Frac (X);
	if XF = 0 then
		XF := Sgn (X);
	for I := 1 to 26 do
	begin
		Z := Z * XF;
		Result := Result + C [I] * Z;
	end;
	if X > 0 then
	begin
		while XF < X do
		begin
			Result := Result / XF;
			XF := XF + 1;
		end;
	end
	else if X < 0 then
	begin
		while XF > X do
		begin
			XF := XF - 1;
			Result := XF * Result;
		end;
	end
end;

function Gamma (const X: Extended): Extended;
begin
	Result := InverseGamma (X);
	if abs (Result) < 1e-4000 then
		raise EMathError.Create ('Not Defined for given Value')
	else
		Result := 1 / Result;
end;

{ Logarithm to base e of the gamma function.

  Accurate to about 1.e-14.
  Programmer: Alan Miller

  Latest revision of Fortran 77 version - 28 February 1988
}
function LnGamma (const X: Extended): Extended;
const
	A1 = -4.166666666554424E-02;
	A2 = 2.430554511376954E-03;
	A3 = -7.685928044064347E-04;
	A4 = 5.660478426014386E-04;
var
	Temp, Arg, Product: Extended;
	Reflect: Boolean;
begin
//  lngamma is not defined if x = 0 or a negative integer.
	if FloatIsZero (X) or FloatIsNegative (X) and (Abs (X - Int (X)) < ESBTolerance) then
		raise EMathError.Create ('Invalid Value');

// If X < 0, use the reflection formula:
//        gamma(x) * gamma(1-x) = pi * cosec(pi.x)

	Reflect := X < 0.0;
	if Reflect then
		Arg := 1.0 - X
	else
		Arg := X;

// Increase the argument, if necessary, to make it > 10.

	Product := 1.0;
	while (Arg <= 10.0) do
	begin
		Product := Product * Arg;
		Arg := Arg + 1.0;
	end;

// Use a polynomial approximation to Stirling's formula.
// N.B. The real Stirling's formula is used here, not the simpler, but less
//     accurate formula given by De Moivre in a letter to Stirling, which
//     is the one usually quoted.

	Arg := Arg - 0.5;
	Temp := 1.0 / Sqr (Arg);
	Result := LnRt2Pi + Arg * (Ln (Arg) - 1.0 +
		(((A4 * Temp + A3) * Temp + A2) * Temp + A1) * Temp) - Ln (Product);

	if Reflect then
	begin
		Temp := Sin (ESBPi * X);
		Result := Ln (ESBPi / Temp) - Result;
	end;
end;

function Beta (const X, Y: Extended): Extended;
var
	R1, R2: Extended;
begin
	R1 := InverseGamma (X);
	R2 := InverseGamma (Y);
	if FloatIsZero (R1) or FloatIsZero (R2) then
		raise EMathError.Create ('Not Defined for given Value')
	else
		Result := InverseGamma (X + Y) / (R1 * R2);
end;

function IncompleteBeta (X: Extended; P, Q: Extended): Extended;
{ Adapted From Collected Algorithms from CACM
	Algorithm 179 - Incomplete Beta Function Ratios
	Oliver G Ludwig
}
const
	Epsilon: Extended = 0.5e-18;
	MaxIterations = 1000;
var
	FinSum, InfSum, Temp, Temp1, Term, Term1, QRecur, Index: Extended;
	I : Integer;
	Alter: Boolean;
begin
	if (P <= 0) or (Q <= 0) or (X < 0) or (X > 1) then
		raise EMathError.Create ('Not Defined for given Value')
	else
	begin
		if (X = 0) or (X = 1) then
			Result := X
		else
		begin
			// Interchange arguments if necessary to get better convergence
			if X <= 0.5 then
				Alter := False
			else
			begin
				Alter := True;
				SwapExt (P, Q);
				X := 1 - X;
			end;

			// Recurs on the (effective) Q until the Power Series doesn't alternate
			FinSum := 0;
			Term := 1;
			Temp := 1 - X;
			QRecur := Q;
			Index := Q;
			repeat
				Index := Index - 1;
				if Index <= 0 then
					Break;
				QRecur := Index;
				Term := Term * (QRecur + 1) / (Temp * (P + QRecur));
				FinSum := FinSum + Term;
			until False;

			// Sums a Power Series for non-integral effective Q and yields unity for integer Q
			InfSum := 1;
			Term := 1;
			for I := 1 to MaxIterations do
			begin
				if Term <= Epsilon then
					Break;
				Index := I;
				Term := Term * X * (Index - QRecur) * (P + Index - 1) /
					(Index * (P + Index));
				InfSum := InfSum + Term;
			end;

			// Evaluates Gammas
			Temp := Gamma (QRecur);
			Temp1 := Temp;
			Term := Gamma (QRecur + P);
			Term1 := Term;
			Index := QRecur;
			repeat
				Temp1 := Temp1 * Index;
				Term1 := Term1 * (Index + P);
				Index := Index + 1;
			until Index >= Q - 0.5;

			Temp := XtoY (X, P) * (InfSum * Term / (P * Temp) + FinSum * Term1
				* XtoY (1 - X, Q) / (Q * Temp1))/Gamma (P);

			if Alter then
				Result := 1 - Temp
			else
				Result := Temp
		end;
	end;
end;

{------------------------------------------------------------------------
PURPOSE:     Calculate the greatest integral power of two that is less
		   than or equal to the parameter
DATE:        April 2000
EXCEPTIONS:  EInvalidArgument for argument zero
------------------------------------------------------------------------}
function IGreatestPowerOf2 (const N: LongWord): LongWord;
begin
	Result := LongWord(1) shl ILog2 (N);
end {GreatestPowerOf2};



end.
