Unit RRMath

*********************************************************************** Author: Robert Rossmair Wagelwitz, Rodaer Str. 6 04688 Mutzschen Germany e-mail: Robert.Rossmair@t-online.de Module: RRMath Version: 02-OCT-1997 Description: 32-bit (not at all!) math routines Copyright (C) 1997 Robert Rossmair This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program (file COPYING); if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

Classes

Functions

Average - uses ExpData; (* const { atanRadix[k] := Round($2000000 * (1 shl k) * ArcTan(1/(1 shl k))); } atanRadix: array[0.
BitCount -
ClipInteger - binaere Quersumme
ConvCartesian -
ConvPolar -
DivComplex -
exp2 -
Factorial - 2^X

function exp2scale(P: byte; X: Smallint): Integer; assembler; { computes 2^(2^P * X/32768) * 32768 } asm push ebx movzx ebx, X mov al, bh mov cl, 7 sub cl, P sar al, cl sub cl, 3 shr bx, cl mov cl, al and bx, 2047 shl bx, 1 mov ax, WORD PTR ExpTable[ebx] dec cl js @2 mov dx, 1 { MSBit substituieren } mov ch, 0 jcxz @exit @1: shl ax, 1 rcl dx, 1 loop @1 jmp @exit @2: not cl stc rcr ax, 1 shr ax, cl adc ax, 0 { Runden } mov dx, 0 @exit: pop ebx end;
FPU_RoundMethod - X!

berechnet X!
I0 - clear pending exceptions
IMax -
IMin -
IMulDiv - Round(Min(Max(X, LongMin), LongMax)
Interpolate - berechnet modifizierte Besselfunktion 1.
Interpolate_F -
Long - ======================================== aus c't 2/92, berechnet Normalverteilung ========================================
Max -
MaxD -
MaxS - Interpolate_F := y1 + (y2 - y1) * (x - x1) / (x2 - x1);
Min -
MinD -
MinS -
MulComplex -
NextMultipleTo -
NormDistribution - idiv Y
Sign - Die Matisse von single ist bloss 23 bit lang, so dass Werte >= $7FFFFFC0 auf $80000000 gerundet werden.
SqrFixed -
tan - Result in edx:eax
UMax -
UMin -
UMulDiv -

Types

Coordinates
IntFraction
NaturalNumber
PComplex
PComplexArray
PComplexL
PComplexLArray
PFixed
PFixedArray
T3DVector
TComplex
TComplexArray
TComplexL
TComplexLArray
TFixedArray
WordFraction

Constants

divSqrt2Pi
exDenorm
exErrorSum
exInvOp
exOverflow
exPrecision
exStackFault
exUnderflow
exZeroDiv
ld_10
ld_e
lg_2
ln_2
M_1_PI
M_1_SQRTPI
M_2_PI
M_2_SQRTPI
M_E
M_LN10
M_LN2
M_LOG10E
M_LOG2E
M_PI
M_PI_2
M_PI_4
M_SQRT2
M_SQRT_2
prDouble
prExtended
prSingle
rcChop
rcDown
rcNearOrEven
rcUp
RC_Mask
stBusy
stC0
stC1
stC2
stC3
stTOP
TwoExp31
TwoExp32

Variables


Functions


function Average(X, Y: Integer): Integer;

uses ExpData; (* const { atanRadix[k] := Round($2000000 * (1 shl k) * ArcTan(1/(1 shl k))); } atanRadix: array[0..23] of longint = ( $01921FB5, $01DAC670, $01F5B760, $01FD5BAA, $01FF55BB, $01FFD55C, $01FFF556, $01FFFD55, $01FFFF55, $01FFFFD5, $01FFFFF5, $01FFFFFD, $01FFFFFF, $02000000, $02000000, $02000000, $02000000, $02000000, $02000000, $02000000, $02000000, $02000000, $02000000, $02000000);

function BitCount(X: Integer): byte;


function ClipInteger(X: Integer): Smallint;

binaere Quersumme

procedure ConvCartesian(var Z: TComplex);


procedure ConvPolar(var Z: TComplex);


procedure DivComplex(Z1, Z2: TComplex; var Result: TComplex);


function exp2(X: Single): Single;


function Factorial(X: word): extended;

2^X

function exp2scale(P: byte; X: Smallint): Integer; assembler; { computes 2^(2^P * X/32768) * 32768 } asm push ebx movzx ebx, X mov al, bh mov cl, 7 sub cl, P sar al, cl sub cl, 3 shr bx, cl mov cl, al and bx, 2047 shl bx, 1 mov ax, WORD PTR ExpTable[ebx] dec cl js @2 mov dx, 1 { MSBit substituieren } mov ch, 0 jcxz @exit @1: shl ax, 1 rcl dx, 1 loop @1 jmp @exit @2: not cl stc rcr ax, 1 shr ax, cl adc ax, 0 { Runden } mov dx, 0 @exit: pop ebx end;


procedure FPU_RoundMethod(Method: word);

X!

berechnet X!


function I0(X: single): single;

clear pending exceptions

function IMax(Arg1, Arg2: Integer): Integer;


function IMin(Arg1, Arg2: Integer): Integer;


function IMulDiv(Factor1, Factor2, Divisor: Integer): Integer;

Round(Min(Max(X, LongMin), LongMax)

function Interpolate(X, X1, Y1, X2, Y2: Integer): Integer;

berechnet modifizierte Besselfunktion 1. Gattung, 0. Ordnung } // Interpolate?? berechnen die Funktion // // y := y1 + (y2 - y1) * (x - x1) / (x2 - x1

berechnet modifizierte Besselfunktion 1. Gattung, 0. Ordnung. Quelle: Oppenheim / Schafer: Zeitdiskrete Signalverarbeitung, 7.4.3 (S. 519), Oldenburg 1992


function Interpolate_F(X, X1, Y1, X2, Y2: Single): Single;


function Long(X: single): longint;

======================================== aus c't 2/92, berechnet Normalverteilung ========================================

function Max(Arg1, Arg2: Extended): Extended;


function MaxD(Arg1, Arg2: Double): Double;


function MaxS(Arg1, Arg2: Single): Single;

Interpolate_F := y1 + (y2 - y1) * (x - x1) / (x2 - x1);

function Min(Arg1, Arg2: Extended): Extended;


function MinD(Arg1, Arg2: Double): Double;


function MinS(Arg1, Arg2: Single): Single;


procedure MulComplex(Z1, Z2: TComplex; var Result: TComplex);


function NextMultipleTo(X, Y: Integer): Integer;


function NormDistribution(rx, ru, rs: single): single;

idiv Y

function Sign(X: Integer): Integer;

Die Matisse von single ist bloss 23 bit lang, so dass Werte >= $7FFFFFC0 auf $80000000 gerundet werden. Der groesste Wert, der nicht zum Runtime Error 207 "Invalid floating point operation" fuehrt, wird auf $7FFFFF80 gerundet. } { >= $7FFFFFBF ? [longint(single = 2^31) = $4F000000]

function SqrFixed(X: TFixed): TFixed;


function tan(x: single): single;

Result in edx:eax

function UMax(Arg1, Arg2: UInt): UInt;


function UMin(Arg1, Arg2: UInt): UInt;


function UMulDiv(Factor1, Factor2, Divisor: UInt): UInt;


Types


Coordinates = (Cartesian, Polar);
Nenner
IntFraction = record
Numerator : Smallint;
Denominator : Smallint;
end;
Bruch
NaturalNumber = 1..High(Integer)

PComplex = ^TComplex

PComplexArray = ^TComplexArray

PComplexL = ^TComplexL

PComplexLArray = ^TComplexLArray

PFixed = ^TFixed

PFixedArray = ^TFixedArray

T3DVector = array[0..2] of Double;

TComplex = record
end;

TComplexArray = array[Word] of TComplex;

TComplexL = record
end;

TComplexLArray = array[Word] of TComplexL;

TFixedArray = array[Word] of TFixed;

WordFraction = record
Numerator : word;
Denominator : word;
end;
Nenner

Constants

divSqrt2Pi = 3.98942280401433E-0001

lg(2) = 1/ld(10)

exDenorm = $0002

exErrorSum = $0080

exInvOp = $0001

Konstanten fuer FPU-Exceptions (Status- und Kontrollregister)

exOverflow = $0008

exPrecision = $0020

exStackFault = $0040

exUnderflow = $0010

exZeroDiv = $0004

ld_10 = 3.32192809488736E+0000

ld(e) = 1/ln(2)

ld_e = 1.44269504088896E+0000

lg_2 = 3.01029995663981E-0001

ln_2 = 6.93147180559945E-0001

ld(10) = 1/lg(2)

M_1_PI = 0.318309886183790671538

M_1_SQRTPI = 0.564189583547756286948

M_2_PI = 0.636619772367581343076

M_2_SQRTPI = 1.12837916709551257390

M_E = 2.71828182845904523536

Constants rounded for 21 decimals.

M_LN10 = 2.30258509299404568402

M_LN2 = 0.693147180559945309417

M_LOG10E = 0.434294481903251827651

M_LOG2E = 1.44269504088896340736

M_PI = 3.14159265358979323846

M_PI_2 = 1.57079632679489661923

M_PI_4 = 0.785398163397448309116

M_SQRT2 = 1.41421356237309504880

M_SQRT_2 = 0.707106781186547524401

prDouble = $0002

prExtended = $0003

prSingle = $0000

FPU precision control

rcChop = $0C00

rcDown = $0400

rcNearOrEven = $0000

Konstanten fuer RC-Feld des FPU - Kontrollregisters:

rcUp = $0800

RC_Mask = $0C00

stBusy = $8000

stC0 = $0100

FPU - Statusregister

stC1 = $0200

stC2 = $0400

stC3 = $4000

stTOP = $3800

TwoExp31 = 2.147483648E+0009

1/Sqrt(2*Pi)

TwoExp32 = 4.294967296E+0009


Variables