To FootThe term "Borland's Pascal" here includes "Turbo Pascal" and "Borland Pascal"; this page is largely for Borland's Pascal and similar.
In exact arithmetic, division by a number is exactly equivalent to multiplication by its reciprocal. In computer arithmetic, where it may not be possible to store the reciprocal of an unteger exactly, that is not always applicable.
In present-day computers, the multiplication and division operators are both about as quick as other operators, and there is only a little to be gained by preferring the former. In the past, that was not the case; division was significantly slower; therefore, compilers did, and still may, replace division by a constant with multiplication by the (possibly inexact) reciprocal.
In particular, therefore, for example :-
See Delphi RoundTo.
Prof. Timo Salmi's Pascal FAQ contains a section (#13) on Powers, and gives a general routine for calculating XY. This undoubtedly gives the correct results, but I beg to differ on a point of principle; I think that his routine is more general than normal usage requires.
Where complex arithmetic (a+jb) is NOT involved, I believe that, in real-world XY problems, one will always know that either ONE or BOTH of "Y is integer" OR "X is not negative" must hold. If Y is integer N, then multiplication (and maybe one division) will suffice. If X is positive, then the exp-log approach is valid. I believe that multiplication is preferable, if efficiently done, except perhaps when Y is large.
I would therefore use whichever one of the following functions fits the need - UNDERTESTED :-
function ExpFixed(const X : Xtype ; const N : longint) : Xtype { X^N } ;
function ExpNum(const A : Xtype ; const B : longint) : Xtype ;
begin if B=0 then ExpNum := 1.0 else
if Odd(B) then ExpNum := ExpNum(A, B-1) * A
else ExpNum := Sqr(ExpNum(A, B div 2)) ;
end {ExpNum} ;
begin
if N<0 then ExpFixed := ExpNum(1.0/X, -N)
else ExpFixed := ExpNum( X, +N) ;
end {ExpFixed} ;
function ExpFloat(const X, Y : Float) : Float { X^Y } ;
begin ExpFloat := Exp(Y*Ln(X)) end {ExpFloat} ;
Xtype and Float should be appropriate types; over-range checks may be added.
Using a well-designed, physically-reasonable program,
the above ExpFloat may well suffice.
For completeness, consider the special cases (AFAIR) :-
* 0.00=1.0,
* 0.0Y=0.0 (Ln fails),
* which makes sense for Y>0;
* but for Y<0, what?;
0.0Y then looks like 1.0/0.0, an attempt at the infinite.
function ExpFloat(const X, Y : Float) : Float { X^Y } ;
begin if X>0.0 then ExpFloat := Exp(Y*Ln(X))
else if X<0.0 then RunError(234) { or whatever }
else if Y>0.0 then ExpFloat := 0.0
else if Y=0.0 then ExpFloat := 1.0
else RunError(200) { (1/X)^Y } ;
end {ExpFloat} (* UNDERTESTED *) ;
Test program source: POWERS.PAS - under revision, Dec '97.
Here's another fast routine (not tested by me) :-
function RFpower(i, n : integer) : integer { i ** n } ;
{ Return i to the power of n. i and n are both of type integer.
The result is also of type integer. n must be non-negative. }
var y : integer ;
begin y := 1 ;
while n>0 do begin
while not odd(n) do begin n := n shr 1 ; i := sqr(i) end ;
n := pred(n) ; y := y*i end ;
RFpower:= y end {Bob Ferguson} ;
For X<0 and Y real, XY is not defined (AFAICS).
For X<0 and Y rational (Y=M/N; M, N integer),
XY is defined, and has N complex values.
After ensuring M, N are co-prime,
if M is even, the values certainly include a positive real;
if M is odd, they seemingly include a negative real.
The news:sci.math FAQ, I gather, refers.
There is normally no doubt as to the results of Div and Mod when both inputs are positive. If either may be negative: (1) think carefully what results you need; (2) check carefully what results you have.
Pascal Mod and Div should give the results of a single basic operation. With Numerator, Denominator, Quotient, and Remainder, one should certainly always have N=D*Q+R. This is not sufficient to define the results.
The following, at least, are possible (* |R| is Abs(R) *) :-
* Minimise |R| - |R| is no bigger than |D|/2
* Minimise |R|, subject to R≥0
* Minimise |R|, R has sign of N
* Minimise |R|, R has sign of D
Additionally, one may require :-
* Q≥0
* Q is zero or has same sign as N
* Q is zero or has same sign as D
I prefer that the Mod result should be invariant on translation by the Divisor, even across a Numerator of zero - that (N+kD) mod D should be independent of the signed integer k.
| Results for Div and Mod | |||||
|---|---|---|---|---|---|
| Inputs | Preferred | Actual | |||
| Numerator | Denominator | Quotient | Remainder | Quotient | Remainder |
| N | D | N div D | N mod D | N div D | N mod D |
| +33 | +5 | +6 | +3 | +6 | +3 |
| +33 | -5 | -7 | -2 | -6 | +3 |
| -33 | +5 | -7 | +2 | -6 | -3 |
| -33 | -5 | +6 | -3 | +6 | -3 |
This is consistent with the requirement that N=Q*D+R : it seems that Object Pascal is consistent with this, but not with my preferences (as in longcalc.pas), which give results as in the Table.
The real points, though, are that the operations are uncertainly defined unless both operands are positive; that the results needed depend on the circumstances; and that care must be taken. Similar considerations can apply to rounding, which is a sort of "Mod 1.0" operation, and may need to be towards -infinity, 0, or +infinity.
In news:borland.public.delphi.objectpascal, Henrick Hellström wrote: I would strongly prefer that (X mod Y) > 0 if Y > 0. The expression (((X mod Y) + Y) mod Y) one has to write from time to time just looks awful!
See also Float to Integer.
Johan Claeys has written "An introduction to complex numbers", part of a linked maths site.
There is no support for complex arithmetic in BP/TP. One can define a Record type TComplex as a pair of Floats of desired resolution, with maybe a length byte for passing it off as a string for function return as in Borland Pascal Normal Usage (and a pad-byte for word-alignment).
The standard operations can then be built up from available primitives. Example, compiled but not run :-
type Float = single {or double or extended} ;
{$IFNDEF CPXFN}
TComplex = record A, B : Float end ;
procedure AddCpx(var X : TComplex {:=}; const Y {+}, Z : TComplex) ;
begin X.A := Y.A + Z.A ; X.B := Y.B + Z.B end {AddCpx} ;
{$ELSE}
const S2 = 2*SizeOf(Float) ; type SComplex = string [S2] ;
TComplex = record
case byte of 0: (L : byte; A, B : Float) ; 1: (S : SComplex) end ;
function FAddCpx(const Y {+}, Z : SComplex) : SComplex ;
var TC : TComplex ;
begin
if (Length(Y)<>S2) or (Length(Z)<>S2) then begin
Writeln('FAddCpx para length error!') ; RunError(233) end ;
with TC do begin L := S2 ;
A := TComplex(Y).A + TComplex(Z).A ;
B := TComplex(Y).B + TComplex(Z).B end ;
FAddCpx := SComplex(TC) end {FAddCpx} ;
{$ENDIF}
Subtraction will be similar; multiplication and division just more complicated; general exponentiation will be worse, probably requiring polar conversion.
EFG has a Delphi unit for complex arithmetic; try also ESB.
This ties in with Principal Values. A non-negative number has two real square, fourth, (etc.) roots, one of each sign. Sometimes it is necessary to be quite careful to choose the correct result.
If complex numbers are allowed, the nth root of any number has n values, as will a number to the power of m/n - which suggests that a number to an irrational power may have an infinitude of values ... I don't know enough here.
It is clearly vital that, whenever an algorithm uses a multi-valued function, adequate care is taken to choose the correct value or values.
The arguments of such as sin and cos, and the result of such as arctan, are in radians; 2 π radians are 360 degrees of arc.
Some expressions have more than one inverse; some equations have more than one solution (X2 = 1, for example). In trigonometry, the Principal Value of an inverse function is taken as that satisfying :-
| -π/2 < PV ≤ +π/2 | --- | (sin, tan, cosec, cot) |
| 0.0 ≤ PV < π | --- | (cos, sec) |
|
(π is pi)
| ||
It is important to be sure to use the correct inverse value (or in some cases to use some or all), which is not necessarily the Principal Value. For example, a number X has TWO inverse tangents, ArcTan(X) and Pi+ArcTan(X). For ArcCos, one inverse is minus the other; for ArcSin, one inverse is Pi minus the other.
Also, Y = 2*N*Pi+ArcTrig(X) is a solution to X = Trig(Y), for any integer N, positive, zero, or negative.
ATAN2(Y, X), implemented in the FPU, does the angle part of cartesian-to-polar conversion; it returns values in -Pi < ATAN2 ≤ Pi (X=Y=0.0 is not an error; the value returned is 0.0). It is a better primitive than ArcTan(X) {which is implemented as ATAN2(1, X), anyway}.
ArcTan(X) cannot even cover the full range of Principal Value; it cannot return +Pi/2.
If necessary, ATAN2(Y, X) can be implemented as (untested translation from VBScript) :-
function AT2(const Y, X : extended) : extended ;
begin
if Y<>0.0
then AT2 := 2.0 * ArcTan( Y / ( X + Sqrt(X*X+Y*Y) ) )
else if X<0 then AT2 := Pi else AT2 := 0.0 ;
end {AT2} ;
Pascal's ArcTan function always returns an angle within -90 .. +90 degrees, in which range Cos is always positive. Therefore, ArcCos(x) := ArcTan(...) cannot be generally correct.
The expression ArcCos(x) = ArcTan (sqrt (1-sqr (x)) /x) given for ArcCos in TP7/BP7 on-line Help for ArcTan (and in Delphi Help up to at least D5) is trigonometrically correct, but it gives wrong results (for negative arguments) in Pascal and Delphi. ArcTan gives the principal value and Sqrt gives the positive root; but it is sometimes necessary to use the other value of one or the other of these. See INVCOS.PAS, which works in Pascal and Delphi, for details; it gives and tests various effective ways of calculating ArcTan, ArcSin and ArcCos. These programs
are in the programs directory (inter alia); INVCOS, FPATAN2 & ARCCOT run with "DCC32 -cc", and FSINCOS might possibly be made to do so.
See also #109 of Timo Salmi's Turbo Pascal FAQ, via my Pascal FAQ links.
function ArcCos2(x : float) : float ;
begin ArcCos2 := 0.5*Pi-ArcSin2(x) {OK} end ;
The above is satisfactory.
I read that ArcCos(x) = 2*ArcTan(Sqrt((1-x)/(1+x))) - needs thought.
The ArcSec and ArcCosec formulae presumably only need to be ArcCos(1/x) and ArcSin(1/x); ArcCot is *slightly* harder.
ArcCot := ArcTan(1/x) ; will fail at x=0;
it seems that ArcCot := HalfPi-ArcTan(x); is right & best.
Note that these two do not give the same value for x<0. They give equally valid values, but the first form is discontinuous, while the second form continuously covers the whole range.
Similar considerations may apply for inverse hyperbolic functions.
Soren wrote that:-
One can use FSIN and FCOS opcodes but since the internal assembler only supports 287 opcodes one needs to code them as
DB $D9 FE {SIN} ; DB $D9 FF {COS}The argument has to be (absolute) less than 2^63.
See also the FSINCOS link above.
Be very careful when averaging angles obtained by using inverse Trig functions; an unimportant input change can flip the angle-number from one end of the Principal range to the other. The appropriate average of 1° and 359° may be either 0° or 180°.
Remember that Pi is a predefined function in the System unit of Pascal and Delphi. Interestingly, I can find in Pascal no <digit string> that makes const Pi = 3.<digit string> ; give exactly the same value as the function System.Pi (which is right, I believe) does; I needed :-
const Pi = 3.1415926535897932385 + 6.5E-19 {!} ;
{ CRC #54: 3.14159265358979323846 26433 83279½ }
I have read that digits after the 19th are ignored. Pi as a const was useful in TP5 (and I think in TP6) for setting other const values such as HalfPi and TwoPi, where a function could not be used. The function can be used for this in BP7.
For Delphi 3, change 6.5E-19 to 2.5E-19 above; but the full CRC number above (omitting ½) sets the new constant equal to System.Pi.
Martin Harvey informs me that the coprocessor floating point instructions include FLDPI, more precise than one can (normally) obtain by trig operations, and setting the guard bits (81,82), used by the hardware to prevent cumulative rounding errors. But I have not seen evidence that BP7 uses this for System.Pi; Delphi, I read, does not do it.
As an Extended, I believe that Pi = $4000c90fdaa22168c235, which pattern can be entered as
const RR : array [0..9] of char = #$35#$c2#$68#$21#$a2#$da#$0f#$c9#$00#$40 ;
In any case, Pi, being transcendental, cannot be represented exactly by a floating-point number. Therefore, sin(Pi), sin(2*Pi), etc., are not necessarily (and are in fact not) computed as exactly zero. In some applications, particularly those using large angles (i.e. many Pi), it can be better to use rational arithmetic for the angles, reducing them to within the principal range or thereabouts before use as arguments to trigonometrical functions.
For Delphi Maths and Trig, see also in my Delphi Programming.
BP & TP provide only the fundamental base-e functions Ln(x) and Exp(x), from which all else can be derived.
| Loge(x) = Ln(x) | Log10(x) = Ln(x) / Ln(10) | Logy(x) = Ln(x) / Ln(y) |
| ex = Exp(x) | 10x = Exp(x*Ln(10)) | yx = Exp(x*Ln(y)) |
| To be re-checked | ||
See also Exponentiation.
"Constant folding" means that expressions involving true untyped constants - literals and those generated by such as const J = 3 ;, not by such as const K : integer = 4 ;, which are typed constants and can be varied - are evaluated, as far as is practicable, at compile time. With the above, Q := 7 ; and Q := J+4 ; generate identical code.
The rule for the run-time evaluation of A op B - that A and B are enlarged to the smallest type accommodating both their ranges - has no jurisdiction for compile-time arithmetic; but it may apply.
Whenever there is any doubt, check.
On Sun, 7 Jun 1998 in comp.lang.pascal.borland, Robert AH Prins wrote :-
... Borland made the RTL source of TP6/7 available for a fee, and gave it as a bonus with BP7. Norbert Juffa completely rewrote the TP6 RTL, and ported most of his changes to the BP7 RTL. Programs compiled with either are usually smaller and always faster than those using the Borland RTL. He also made "real" arithmetic as IEEE compliant as possible.
Since then, he has written his own version, 32-bit code only.
BPL70V20.ZIP contains optimized Run Time Libraries for Borland Pascal 7. The code is based on Norbert Juffa's BPL70N16, but incorporates changes Borland made in BP 7.01. Unlike the libraries in Norbert's BPL70N16, the ones include here only run on 386+ CPU's due to the extensive use of 32-bit instructions. Many of the original RTL routines have also been split up further to increase the smart-link granularity.
In addtion to the System unit, the RTLs also include optimized Dos, CRT (removal of RTE 200) and Overlay units.
459244 Oct 10 2004 ftp://garbo.uwasa.fi/pc/turbopa7/bpl70v20.zip
bpl70v20.zip Borland-Pascal 7.01 RT-Library Update Rel 2.0,
R.Prins
Wirth's Pascal was no doubt excellent for its purpose, of teaching; but Borland's Pascal and Delphi are sold as programming tools. Many of us have spent a modicum of time each implementing standard mathematical functions not in Wirth's minimal set; many of us have had to prove and debug these; some of us may still unknowingly need to debug them.
To my mind, it is a pity that Borland did not include with Pascal a standard unit, maybe MATHS.PAS/TPU/TPP/TPW, containing the rest of the standard functions, well implemented. We would then all have had access to the same small set of initial bugs, which by now would all have been killed off. Granted that others have written such, but only Borland could have packaged a truly standard one. What's in Delphi? - I think the Maths stuff is only in the more expensive or later editions.
See
ftp://garbo.uwasa.fi/pc/turbopas/ All sorts of good things, including ...
See
50898 Sep 6 07:38 ftp://garbo.uwasa.fi/pc/turbopas/tpmath1.zip tpmath1.zip Mathematical TP freeware library, units, Jean Debord 51028 Sep 6 07:48 ftp://garbo.uwasa.fi/pc/turbopas/tpmath2.zip tpmath2.zip Mathematical TP freeware library, demos, Jean Debord 19490 Apr 29 18:36 ftp://garbo.uwasa.fi/pc/turbopas/tpmath3.zip tpmath3.zip Mathematical TP freeware lib, interfacing, Jean Debord
See
1417 Sep 22 19:15 ftp://garbo.uwasa.fi/pc/turbopas/mathtp10.zip math10.zip Turbo Pascal 4+ Math unit uploaded, Michael Husted
The largest true integer type in TP/BP is longint, signed 32 bits. Some operations on unsigned 32 bits can be done with relevant checks off, and/or XORing with $80000000. 64-bit operations can with care be decomposed into 32-bit ones, in the same manner as one learnt at school to handle numbers in 10..99; this can be extended.
The comp type provides some operations on 64-bit integers.
Later versions of Delphi have proper unsigned 32-bit and signed 64-bit quantities.
I believe that there may be libraries on the Net and elsewhere to handle these.
My own code is used in LONGCALC.PAS, which is compatible with (at least) TP7, BP7, D3 "DCC32 -cc". This does integer arithmetic, bases 2..16, range over 65000 digits, and is scriptable in RPN. LONGCALC can be used as a Unit, by an experienced programmer.
The Prime Numbers are 2, 3, 5, 7, 11, 13, 17, 19, ... . The number 1 is not a prime.
I have terse Pascal implementations of the Sieve of Eratosthenes as eratost*.pas in my programs directory. They differ in that *1 marks in an array of booleans representing numbers, *2 in an array of booleans representing odd numbers, and *3 in an array of sets (effectively a bit array) representing odd numbers. The start at n2 optimisation is not yet included.
There is also primenos.pas, demonstrating a simple test on an individual (long) integer.
I have read that one should do a floating point division (4195835/3145727) to detect whether the Pentium (early versions) FDIV bug is present. The right result is 1.333820449136241; buggy CPUs give 1.333739068902037589.
I hear that the Merced chip "has no built-in FDIV at all".