Unit ABSDecimalRounding;

(* *****************************************************************************
Author: John Herbster 
DecimalRounding_JH1

 ROUTINES FOR ROUNDING IEEE-754 FLOATS TO SPECIFIED NUMBER OF DECIMAL FRACTIONS

  These routines round input values to fit as closely as possible to an
  output number with desired number of decimal fraction digits.

  Because, in general, numbers with decimal fractions cannot be exactly
  represented in IEEE-754 floating binary point variables, error limits
  are used to determine if the input numbers are intended to represent an
  exact decimal fraction rather than a nearby value.   Thus an error limit
  will be taken into account when deciding that a number input such as
  1.295, which internally is represented 1.29499 99999 …, really should
  be considered exactly 1.295 and that 0.29999 99999 ... really should
  be interpreted as 0.3 when applying the rounding rules.

  Some terminology:
    "NDFD" is used for Number of Decimal Fraction Digits.  If NDFD is
        negative, then the inputs will be rounded so that there are zeros
        on the left of the decimal point.  I.E. if NDFD = -3, then the
        output will be rounded to an integral multiple of a thousand.
    "MaxRelError" designates the maximum relative error to be allowed in the
        input values when deciding they are supposed to represent exact
        decimal fractions (as mentioned above).
    "Ctrl" determines the type of rounding to be done.  Nine kinds of
        rounding (plus no rounding) are defined.  They include almost
        every kind of rounding known.  See the definition of
        tDecimalRoundingCtrl below for the specific types.

 Rev. 2002.01.14 by John Herbster to add overloaded IsNaN() functions to
     solve problem with invalid operation when loading values for test.
 Rev. 2001.11.24 by John Herbster to add KnownErrorLimit, SafetyFactor,
     and MaxRelErrXxx constants. Also corrected DecimalRoundingCtrlStrs
     and comments per feedback from Marjan Venema.
 Pgm. 2001.09.19 by John Herbster.  Please send feedback and suggestions
     to author at herb-sci@swbell.net.

 Dedicated to the participants in the borland.public.delphi.objectpascal
     forum.  With special thanks to Brad White for his proofing and
     suggestions.

***************************************************************************** *)

Interface

uses Classes;

Type tDecimalRoundingCtrl =    {The defined rounding methods}
   (drNone,    {No rounding.}
    drHalfEven,{Round to nearest or to even whole number. (a.k.a Bankers) }
    drHalfPos, {Round to nearest or toward positive.}
    drHalfNeg, {Round to nearest or toward negative.}
    drHalfDown,{Round to nearest or toward zero.}
    drHalfUp,  {Round to nearest or away from zero.}
    drRndNeg,  {Round toward negative.                    (a.k.a. Floor) }
    drRndPos,  {Round toward positive.                    (a.k.a. Ceil ) }
    drRndDown, {Round toward zero.                        (a.k.a. Trunc) }
    drRndUp);  {Round away from zero.}

{ The following DecimalRound function is for doing the best possible job of
  rounding floating binary point numbers to the specified NDFD.  MaxRelError
  is the maximum relative error that will be allowed when determining the
  cut points for applying the rounding rules. }

Function DecimalRound
   (Value: extended;     {Input value to be rounded.}
    NDFD: integer;       {Number decimal fraction digits to figure in result.}
    MaxRelErr: double;   {Maximum relative error to assume in input value.}
    Ctrl: tDecimalRoundingCtrl = drHalfEven)  {Optional rounding rule}
    : extended;

{ The following functions have a two times "epsilon" error built in for the
  single, double, and extended argument respectively: }

Function DecimalRoundSgl (Value: single; NDFD: integer;
                          Ctrl: tDecimalRoundingCtrl = drHalfEven): extended;

Function DecimalRoundDbl (Value: double; NDFD: integer;
                          Ctrl: tDecimalRoundingCtrl = drHalfEven): extended;

Function DecimalRoundExt (Value: extended; NDFD: integer;
                          Ctrl: tDecimalRoundingCtrl = drHalfEven): extended;

{----- The following are utility functions and constants ----------------------}

{ The following "epsilon" values are representative of the resolution of the
  floating point numbers divided by the number being represented.
  These constants are supplied to the rounding routines to determine how much
  correction should be allowed for the natural errors in representing
  decimal fractions.  Using 2 times or higher multiples of these values may
  be advisable if the data have been massaged through arithmetic
  calculations. }
Const
  SglEps = 1.1920928955e-07;
  DblEps = 2.2204460493e-16;
  ExtEps = 1.0842021725e-19;
{ These epsilon values are smallest for which "1 + epsilon <> 1".
      For "1 - epsilon <> 1", divide the following values by 2. }

Const
  KnownErrorLimit = 1.234375;
  SafetyFactor = 2;
  MaxRelErrSgl = SglEps * KnownErrorLimit * SafetyFactor;
  MaxRelErrDbl = DblEps * KnownErrorLimit * SafetyFactor;
  MaxRelErrExt = ExtEps * KnownErrorLimit * SafetyFactor;
{ If MaxRelErr < XxxEps * KnownErrorLimit then errors can occur. }

Procedure LoadDecimalRoundingCtrlAbbrs(const Strings: tStrings);
{ This procedure loads the tDecimalRoundingCtrl descriptions and ordinals
  into the string list for such use as using a TCombobox to make rounding
  type selection. }

Procedure CalcEpsValues (var SglEps, DblEps, ExtEps: double);
{ This procedure was used to compute the values SglEps, DblEps, and ExtEps: }

Function IsNAN (const sgl: single)  : boolean; overload;
Function IsNAN (const dbl: double)  : boolean; overload;
Function IsNAN (const ext: extended): boolean; overload;
{ Each returns true if the value passed is not-a-number. }

Function GetX87CW: word;
{ Returns the FPU control word (which indicates interrupt masks and
      precision and rounding modes). }

Function IsFpuCwOkForRounding: boolean;
{ Returns true if floating point processor (FPU) is correctly set
    (1) to allow conversion from extended to double and double to single
        without creating the the loss-of-precision interrupt or exception,
    (2) to do arithmetic internal to FPU in extended precision, and
    (3) to internally use halves-to-even (a.k.a. bankers) rounding. }

Const
  DecimalRoundingCtrlStrs: array [tDecimalRoundingCtrl] of
      record Abbr: AnsiString; Dscr: AnsiString; end = (
    (Abbr: 'None'    ; Dscr: 'No rounding.'),
    (Abbr: 'HalfEven'; Dscr: 'Round to nearest or to even whole number'+
                             ' (a.k.a Bankers) '),
    (Abbr: 'HalfPos' ; Dscr: 'Round to nearest or toward positive'),
    (Abbr: 'HalfNeg' ; Dscr: 'Round to nearest or toward negative'),
    (Abbr: 'HalfDown'; Dscr: 'Round to nearest or toward zero'),
    (Abbr: 'HalfUp'  ; Dscr: 'Round to nearest or away from zero'),
    (Abbr: 'RndNeg'  ; Dscr: 'Round toward negative. (a.k.a. Floor) '),
    (Abbr: 'RndPos'  ; Dscr: 'Round toward positive. (a.k.a. Ceil ) '),
    (Abbr: 'RndDown' ; Dscr: 'Round toward zero. (a.k.a. Trunc) '),
    (Abbr: 'RndUp'   ; Dscr: 'Round away from zero.') );

Const
{ These FPU Control Word bit masks prevent interrupt when present: }
  IM = $0001; {Invalid op interrupt Mask}
  DM = $0002; {Denormalized op interrupt Mask}
  ZM = $0004; {Zero divide interrupt Mask}
  OM = $0008; {Overflow interrupt Mask}
  UM = $0010; {Underflow interrupt Mask}
  PM = $0020; {Loss of precision interrupt Mask}
{ The "pending interrupt" flags in status word have matching positions. }
  IntrM = IM or DM or ZM or OM or UM or PM;

{ These FPU Control Word bit fields change operation: }
  PC = $0300; {Precision Control mask}
  RC = $0C00; {Rounding Control mask}
  pcSingle  = $0000; pcDouble = $0200; pcExtended = $0300;
  rcBankers = $0000; rcFloor  = $0400; rcCeil     = $0800; rcChop = $0C00;

Implementation

Function DecimalRound(Value: extended; NDFD: integer; MaxRelErr: double;
                         Ctrl: tDecimalRoundingCtrl = drHalfEven): extended;
{ The DecimalRounding function is for doing the best possible job of rounding
  floating binary point numbers to the specified (NDFD) number of decimal
  fraction digits.  MaxRelErr is the maximum relative error that will allowed
  when determining when to apply the rounding rule.  }
var i64, j64: Int64; j: integer; m, ScaledVal, ScaledErr: extended;
begin

  If IsNaN(Value) or (Ctrl = drNone)
    then begin Result := Value; EXIT end;

{ Compute 10^NDFD and scale the Value and MaxError: }
  m := 1; For j := 1 to abs(NDFD) do m := m*10;
  If NDFD >= 0
    then begin
      ScaledVal := Value * m;
      ScaledErr := abs(MaxRelErr*Value) * m;
      end
    else begin
      ScaledVal := Value / m;
      ScaledErr := abs(MaxRelErr*Value) / m;
      end;

{ Do the diferent basic types separately: }
  Case Ctrl of
    drHalfEven: begin
      i64 := round((ScaledVal - ScaledErr));
      j64 := round((ScaledVal + ScaledErr));
      if Odd(i64)
        then i64 := j64;
      end;
    drHalfDown:  {Round to nearest or toward zero.}
      i64 := round((abs(ScaledVal) - ScaledErr));
    drHalfUp:    {Round to nearest or away from zero.}
      i64 := round((abs(ScaledVal) + ScaledErr));
    drHalfPos:   {Round to nearest or toward positive.}
      i64 := round((ScaledVal + ScaledErr));
    drHalfNeg:   {Round to nearest or toward negative.}
      i64 := round((ScaledVal - ScaledErr));
    drRndNeg:    {Truncate toward negative. (a.k.a. Floor)}
      i64 := round((ScaledVal + (ScaledErr - 1/2)));
    drRndPos:    {Truncate toward positive. (a.k.a. Ceil)}
      i64 := round((ScaledVal - (ScaledErr - 1/2)));       {}
    drRndDown:   {Truncate toward zero (a.k.a. Trunc).}
      i64 := round((abs(ScaledVal) + (ScaledErr - 1/2)));
    drRndUp:     {Truncate away from zero.}
      i64 := round((abs(ScaledVal) - (ScaledErr - 1/2)));
    else i64 := round(ScaledVal);
    end{cases};

{ Finally convert back to the right order: }
  If NDFD >= 0
    then Result := i64 / m
    else Result := i64 * m;

  If (Ctrl in [drHalfDown,drHalfUp,drRndDown,drRndUp]) and
     (Value < 0)
    then Result := -Result;

end{DecimalRound};

Function DecimalRoundSgl(Value: single; NDFD: integer;
                         Ctrl: tDecimalRoundingCtrl = drHalfEven): extended;
begin
  Result := DecimalRound(Value,NDFD,MaxRelErrSgl,Ctrl);
end;

Function DecimalRoundDbl(Value: double; NDFD: integer;
                         Ctrl: tDecimalRoundingCtrl = drHalfEven): extended;
begin
  Result := DecimalRound(Value,NDFD,MaxRelErrDbl,Ctrl);
end;

Function DecimalRoundExt(Value: extended; NDFD: integer;
                         Ctrl: tDecimalRoundingCtrl = drHalfEven): extended;
begin
  Result := DecimalRound(Value,NDFD,MaxRelErrExt,Ctrl);
end;

Procedure CalcEpsValues (var SglEps, DblEps, ExtEps: double);
{ Compute smallest 1/(2^n) epsilon values for which "1 + epsilon <> 1".
  For "1 - epsilon <> 1", divide these computed values by 2. }
var s: single; d: double; e,f: extended;
begin
{ Compute for single, s: }
  f := 1;
  Repeat f := f/2; s := 1 + f/2; Until s = 1;
  SglEps := f;
{ Compute for double, d: }
  f := 1;
  Repeat f := f/2; d := 1 + f/2; Until d = 1;
  DblEps := f;
{ Compute with extended, e: }
  f := 1;
  Repeat f := f/2; e := 1 + f/2; Until e = 1;
  ExtEps := f;
end;

{$WRITEABLECONST OFF}
type  TExtPackedRec = packed record Man: Int64; Exp: word end;
const SglExpBits: LongInt = $7F800000;          { 8 bits}
      DblExpBits: Int64   = $7FF0000000000000;  {11 bits}
      ExtExpBits: word    = $7FFF;              {15 bits}

Function IsNAN (const sgl: single): boolean;
var InputX: LongInt absolute sgl;
begin
  Result := (InputX <> 0) and ((InputX and SglExpBits)=SglExpBits);
end;

Function IsNAN (const dbl: double): boolean;
var InputX: Int64 absolute dbl;
begin
  Result := (InputX <> 0) and ((InputX and DblExpBits)=DblExpBits);
end;

Function IsNAN (const ext: extended): boolean;
var InputX: TExtPackedRec absolute ext;
begin
  Result := (InputX.Man <> 0) and ((InputX.Exp and ExtExpBits)=ExtExpBits);
end;

Function GetX87CW: word;
{ Returns the FPU control word (which indicates interrupt masks and
      precision and rounding modes). }
asm
  FStCW [Result]
end;

Function IsFpuCwOkForRounding: boolean;
{ Checks to see that floating point processor (FPU) is correctly set to
    (1) allow conversion from extended to double and double to single
        without creating the the loss of precision interrupt or exception,
    (2) do arithmetic internal to FPU in extended precision, and
    (3) use round halves-to-even (a.k.a. bankers rounding) internally. }
var CW: word;
begin
  CW := GetX87CW;
  Result := ((CW and (PC or RC or PM)) = (rcBankers or pcExtended or PM));
end;

Procedure LoadDecimalRoundingCtrlAbbrs(const Strings: tStrings);
var dr: tDecimalRoundingCtrl;
begin
  Assert(Strings<>nil);
  Strings.Clear;
  For dr := low(dr) to high(dr)
      do Strings.AddObject(DecimalRoundingCtrlStrs[dr].Abbr,pointer(dr));
end;

End.
