Using the Math Unit

Problem/Question/Abstract:

I've heard of the Math unit that's included with the Developer level version of Delphi. How do I use it?

Answer:

An Unsung Hero?

If you haven't heard of the Math unit, you're not alone. It's one of those units that's kind of buried in the myriad directories under the Delphi directory. Also, it's only included in the Delphi 2.0 Developer version and above. For those of you who have Delphi 2.0 Developer, you can find the source code file in the \\Borland\Delphi 2.0\Source\RTL\SYS directory.

This unit is one of those things I've heard a lot of people mention in the past but unfortunately, I haven't seen many examples of using it. I'm not sure whether it's because developers don't use it much or don't know about it. In any case, it's a shame because what the Math unit has to offer could be helpful to people needing to use mathematical functions in their code.

Not only are they helpful by their mere existence, but several of the member functions are written in Assembly language that is optimized for the Pentium FPU, so they're incredibly fast and efficient. For example, the procedure SinCos, which is one of the procedures written in Assembly, will produce the Sin and Cos simultaneously of an any angle faster than if you called the Sin and Cos functions individually. I've seen it at work and it's amazing.

The Math unit includes several categories of mathematical functions you can incorporate into your code. These include:

Trigonometric, Hyperbolic and Angle Conversion Functions
Logorithmic Functions
Exponential Functions
Some Miscellaneous Computing Functions
Several Statistical Functions
The Standard Set of Quattro Pro Financial Functions

Mind you, not all of the functions are coded in Assembly, but the mere fact that they've already been written means you don't have to, so that's a real time saver. Below are the function prototypes for the unit, so you can see what's in the file:

{-----------------------------------------------------------------------
Copyright (c) Borland International.

Most of the following trig and log routines map directly to Intel 80387 FPU
floating point machine instructions.  Input domains, output ranges, and
error handling are determined largely by the FPU hardware.
Routines coded in assembler favor the Pentium FPU pipeline architecture.
-----------------------------------------------------------------------}

{ Trigonometric functions }
function ArcCos(X: Extended): Extended;  { IN: |X| <= 1  OUT: [0..PI] radians }
function ArcSin(X: Extended): Extended;  { IN: |X| <= 1  OUT: [-PI/2..PI/2] radians }

{ ArcTan2 calculates ArcTan(Y/X), and returns an angle in the correct quadrant.
  IN: |Y| < 2^64, |X| < 2^64, X <> 0   OUT: [-PI..PI] radians }
function ArcTan2(Y, X: Extended): Extended;

{ SinCos is 2x faster than calling Sin and Cos separately for the same angle }
procedure SinCos(Theta: Extended; var Sin, Cos: Extended) register;
function Tan(X: Extended): Extended;
function Cotan(X: Extended): Extended;           { 1 / tan(X), X <> 0 }
function Hypot(X, Y: Extended): Extended;        { Sqrt(X**2 + Y**2) }

{ Angle unit conversion routines }
function DegToRad(Degrees: Extended): Extended;  { Radians := Degrees * PI / 180}
function RadToDeg(Radians: Extended): Extended;  { Degrees := Radians * 180 / PI }
function GradToRad(Grads: Extended): Extended;   { Radians := Grads * PI / 200 }
function RadToGrad(Radians: Extended): Extended; { Grads := Radians * 200 / PI }
function CycleToRad(Cycles: Extended): Extended; { Radians := Cycles * 2PI }
function RadToCycle(Radians: Extended): Extended;{ Cycles := Radians / 2PI }

{ Hyperbolic functions and inverses }
function Cosh(X: Extended): Extended;
function Sinh(X: Extended): Extended;
function Tanh(X: Extended): Extended;
function ArcCosh(X: Extended): Extended;   { IN: X >= 1 }
function ArcSinh(X: Extended): Extended;
function ArcTanh(X: Extended): Extended;   { IN: |X| <= 1 }

{ Logorithmic functions }
function LnXP1(X: Extended): Extended;   { Ln(X + 1), accurate for X near zero }
function Log10(X: Extended): Extended;                     { Log base 10 of X}
function Log2(X: Extended): Extended;                      { Log base 2 of X }
function LogN(Base, X: Extended): Extended;                { Log base N of X }

{ Exponential functions }

{ IntPower: Raise base to an integral power.  Fast. }
function IntPower(Base: Extended; Exponent: Integer): Extended register;

{ Power: Raise base to any power.
  For fractional exponents, or exponents > MaxInt, base must be > 0. }
function Power(Base, Exponent: Extended): Extended;


{ Miscellaneous Routines }

{ Frexp:  Separates the mantissa and exponent of X. }
procedure Frexp(X: Extended; var Mantissa: Extended; var Exponent: Integer) register;

{ Ldexp: returns X*2**P }
function Ldexp(X: Extended; P: Integer): Extended register;

{ Ceil: Smallest integer >= X, |X| < MaxInt }
function Ceil(X: Extended):Integer;

{ Floor: Largest integer <= X,  |X| < MaxInt }
function Floor(X: Extended): Integer;

{ Poly: Evaluates a uniform polynomial of one variable at value X.
    The coefficients are ordered in increasing powers of X:
    Coefficients[0] + Coefficients[1]*X + ... + Coefficients[N]*(X**N) }
function Poly(X: Extended; const Coefficients: array of Double): Extended;

{-----------------------------------------------------------------------
Statistical functions.

Common commercial spreadsheet macro names for these statistical and
financial functions are given in the comments preceding each function.
-----------------------------------------------------------------------}

{ Mean:  Arithmetic average of values.  (AVG):  SUM / N }
function Mean(const Data: array of Double): Extended;

{ Sum: Sum of values.  (SUM) }
function Sum(const Data: array of Double): Extended register;
function SumOfSquares(const Data: array of Double): Extended;
procedure SumsAndSquares(const Data: array of Double;
  var Sum, SumOfSquares: Extended) register;

{ MinValue: Returns the smallest signed value in the data array (MIN) }
function MinValue(const Data: array of Double): Double;

{ MaxValue: Returns the largest signed value in the data array (MAX) }
function MaxValue(const Data: array of Double): Double;

{ Standard Deviation (STD): Sqrt(Variance). aka Sample Standard Deviation }
function StdDev(const Data: array of Double): Extended;

{ MeanAndStdDev calculates Mean and StdDev in one pass, which is 2x faster than
  calculating them separately.  Less accurate when the mean is very large
  (> 10e7) or the variance is very small. }
procedure MeanAndStdDev(const Data: array of Double; var Mean, StdDev: Extended);

{ Population Standard Deviation (STDP): Sqrt(PopnVariance).
  Used in some business and financial calculations. }
function PopnStdDev(const Data: array of Double): Extended;

{ Variance (VARS): TotalVariance / (N-1). aka Sample Variance }
function Variance(const Data: array of Double): Extended;

{ Population Variance (VAR or VARP): TotalVariance/ N }
function PopnVariance(const Data: array of Double): Extended;

{ Total Variance: SUM(i=1,N)[(X(i) - Mean)**2] }
function TotalVariance(const Data: array of Double): Extended;

{ Norm:  The Euclidean L2-norm.  Sqrt(SumOfSquares) }
function Norm(const Data: array of Double): Extended;

{ MomentSkewKurtosis: Calculates the core factors of statistical analysis:
  the first four moments plus the coefficients of skewness and kurtosis.
  M1 is the Mean.  M2 is the Variance.
  Skew reflects symmetry of distribution: M3 / (M2**(3/2))
  Kurtosis reflects flatness of distribution: M4 / Sqr(M2) }
procedure MomentSkewKurtosis(const Data: array of Double;
  var M1, M2, M3, M4, Skew, Kurtosis: Extended);

{ RandG produces random numbers with Gaussian distribution about the mean.
  Useful for simulating data with sampling errors. }
function RandG(Mean, StdDev: Extended): Extended;

{-----------------------------------------------------------------------
Financial functions.  Standard set from Quattro Pro.

Parameter conventions:

From the point of view of A, amounts received by A are positive and
amounts disbursed by A are negative (e.g. a borrower's loan repayments
are regarded by the borrower as negative).

Interest rates are per payment period.  11% annual percentage rate on a
loan with 12 payments per year would be (11 / 100) / 12 = 0.00916667

-----------------------------------------------------------------------}

type
  TPaymentTime = (ptEndOfPeriod, ptStartOfPeriod);

{ Double Declining Balance (DDB) }
function DoubleDecliningBalance(Cost, Salvage: Extended;
  Life, Period: Integer): Extended;

{ Future Value (FVAL) }
function FutureValue(Rate: Extended; NPeriods: Integer; Payment, PresentValue:
  Extended; PaymentTime: TPaymentTime): Extended;

{ Interest Payment (IPAYMT)  }
function InterestPayment(Rate: Extended; Period, NPeriods: Integer; PresentValue,
  FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

{ Interest Rate (IRATE) }
function InterestRate(NPeriods: Integer;
  Payment, PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

{ Internal Rate of Return. (IRR) Needs array of cash flows. }
function InternalRateOfReturn(Guess: Extended;
  const CashFlows: array of Double): Extended;

{ Number of Periods (NPER) }
function NumberOfPeriods(Rate, Payment, PresentValue, FutureValue: Extended;
  PaymentTime: TPaymentTime): Extended;

{ Net Present Value. (NPV) Needs array of cash flows. }
function NetPresentValue(Rate: Extended; const CashFlows: array of Double;
  PaymentTime: TPaymentTime): Extended;

{ Payment (PAYMT) }
function Payment(Rate: Extended; NPeriods: Integer;
  PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

{ Period Payment (PPAYMT) }
function PeriodPayment(Rate: Extended; Period, NPeriods: Integer;
  PresentValue, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

{ Present Value (PVAL) }
function PresentValue(Rate: Extended; NPeriods: Integer;
  Payment, FutureValue: Extended; PaymentTime: TPaymentTime): Extended;

{ Straight Line depreciation (SLN) }
function SLNDepreciation(Cost, Salvage: Extended; Life: Integer): Extended;

{ Sum-of-Years-Digits depreciation (SYD) }
function SYDDepreciation(Cost, Salvage: Extended; Life, Period: Integer): Extended;

Clearing Things Up: Usage

Whew! That's a lot of information to digest! I listed it here to impress upon you just how much there is. For those of you creating financial applications, the financial functions will come in handy (I sure wish I had these functions available when I was writing financial software).

Listing the functions is one thing - actually using them is another. As you can see, most of the input parameters require an Extended numeric type, which is a 10-byte number ranging from 3.4 * 10e-4932 to 1.1 * 10e4932 in scientific notation. In other words, you can have incredibly huge numbers as input values.

Take a moment to look at the statistical functions. Notice anything odd about almost all of the functions' input parameters? Most of them take a constant open array of double! This implies you can pass any size array as an input parameter, but it must be passed as a constant array, which means that you have to pass the array in the form of (1, 2, 3, 4, 5, 6 ..). That's not so difficult with small known sets of numbers; just hard code them in. But arrays in Pascal are typically of the variable type, where you define a finite number of elements, then fill in the element values. This poses a bit of a problem. Fortunately, there's a solution.

Buried in the System unit is a function called Slice: function Slice(var A: array; Count: Integer): array; Essentially, what Slice does is take a certain number of elements from an array, starting at the beginning, and passes the slice of the array as an open array parameter, fooling the compiler into thinking a constant array is being passed. This means that you can pass the entire array or smaller subset. In fact, Slice can only be used within the context of being passed as an open array parameter. Using it outside of this context will create a compiler error. What a convenient function! So, we can define a variable type array as we're used to in Delphi, fill it up, put it into Slice, which is then used in one of the functions. For example: MyExtendedNumber := Mean(Slice(MyArray, NumElementsPassed));

At this point, you're probably thinking this is pretty incredible stuff. But there's one thing that still bothers me about it: The place where the statistical functions would be most useful is on columnar calculations on tables. Unfortunately, you never know how many records are in a table until runtime. Granted, depending upon the amount of RAM in your system, you could make an incredibly large array of let's say 100K elements, fill it up to the record count of your table, then apply Slice to grab only those you need. However, that's pretty inefficient. Also, in my immediate experience, many of my tables have well over 100K records, which means I'd have to hard code an even greater upper limit. But there will also be tables that have far fewer records than 100K - more like 10K. So the idea then is to strike a balance. No thanks!

Doing the DynArray Thing

So where am I if I can't accept defining a huge array, or making some sort of size compromise? I guess I need to create a variable sized array whose size can be defined at runtime.

Wait a minute! You're not supposed to be able to do that in Delphi!

You can, but it takes some pointers to be able to pull it off. For an in-depth discussion of the technique, I'm going to point you to an article on the enquiry site entitled Runtime Resizeable Arrays, which will show you how to create an array that has an element count you don't know about until runtime. I highly suggest reading the article before continuing, if you're not familiar with the technique.

What gives you the ability to create a dynamic array is a function like the following (this in the article):

type
  TResizeArr = array[0..0] of string;
  PResizeArr = ^TResizeArr;
  ...

  {============================================================================
   Procedure which defines the dynarray. Note that the THandle and Pointer to
   the array are passed by reference. This is so that they may defined outside
   the scope of this proc.
   ============================================================================}

procedure DefineDynArray(var h: THandle; {Handle to mem pointer}
  NumElements: LongInt; {Number of items in array}
  var PArr: PResizeArr); {Pointer to array struct}
begin
  {Allocate Windows Global Heap memory}
  h := GlobalAlloc(GMEM_FIXED, NumElements * sizeof(TResizeArr));
  PArr := GlobalLock(h);
end;

Note that you can just as easily replace the array of String with an array of Double. In any case, the gist of what the procedure does is allocate memory for the number of elements that you want to have in your array (TResizeArr), then locks that area of the heap and assigns it to an array pointer (PResizeArr). To load values into the array, you de-reference the pointer as follows: MyPointerArray^[0] := 1234.1234;

To pass the entire array into the Mean function as above, all we have to do is de-reference the entire array as follows: MyExtendedNumber := Mean(Slice(MyPointerArray^, NumElementsPassed));

Putting It All Together

I mentioned above that the best place to employ the statistical functions is in performing statistics on columnar data in a table. The unit code below provides a simple example of loading a DynArray with columnar data, then performing the Mean function on the loaded array. Note the table that I used had about 80K records in it.

unit parrmain;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  DB, DBTables, ComCtrls, Math, StdCtrls, ds_stats;

type
  TResizeArr = array[0..0] of Double;
  PResizeArr = ^TResizeArr;
  TForm1 = class(TForm)
    StatusBar1: TStatusBar;
    Button1: TButton;
    Table1: TTable;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
    PA: PResizeArr;
    Hndl: THandle;
    procedure DefineDynArray(var H: THandle;
      NumElem: LongInt;
      var PArr: PResizeArr);
    procedure FillArray;
    procedure SaveScreen;
  public
    { Public declarations }

  end;

var
  Form1: TForm1;

implementation

{$R *.DFM}

procedure TForm1.DefineDynArray(var H: THandle;
  NumElem: LongInt;
  var PArr: PResizeArr);
begin
  H := GlobalAlloc(GMEM_FIXED, NumElem * SizeOf(TResizeArr));
  PArr := GlobalLock(H);
end;

procedure TForm1.LoadArray;
var
  tbl: TTable;
  recs: LongInt;
begin
  tbl := TTable.Create(Application); //Create a TTable instance
  with tbl do
  begin //and set properties
    Active := False;
    DatabaseName := 'Primary';
    TableName := 'Test';
    TableType := ttParadox;
    Open;

    recs := RecordCount; //Get number of records in table

    DefineDynArray(Hndl, recs, PA); //Now, define our dynarray
    recs := 0; //Reset recs for reuse

    StatusBar1.SimpleText := 'Now filling array';

    while not EOF do
    begin
      Application.ProcessMessages; //allow background processing
      try
        PA^[recs] := FieldByName('Charge').AsFloat;
        StatusBar1.SimpleText := 'Grabbed value of: ' + FloatToStr(PA^[recs]);
        StatusBar1.Invalidate;
        Inc(recs);
        Next;
      except
        GlobalUnlock(Hndl);
        GlobalFree(Hndl);
        Exit;
      end;
    end;

    //Pop up a message to show what was calculated.
    ShowMessage(FloatToStr(Mean(Slice(PA^, RecordCount))));
      //pass the array using Slice
    GlobalUnlock(Hndl); //Unlock and Free memory and TTable instance.
    GlobalFree(Hndl);
    tbl.Free;
  end;
end;

procedure TForm1.Button1Click(Sender: TObject);
begin
  LoadArray;
end;

end.

You could get pretty complex with this by creating a component that encapsulates the statistical functions and grabs data off a table. Using the principles of the code above, it shouldn't be too hard to do. Follow the code; better yet, try it out, supply your own values, and see what you come up with.

We covered a lot of ground here. I wasn't happy to tell you merely about the Math unit and all the wonderful routines it contains; I wanted to show you a way to employ a major portion of it in as flexible a way as possible.

In my opinion, it's not enough just to know about something in programming; you have to know how to use it. With the material I've presented, you should be able to employ the functions of the Math unit in very little time.