Big Integers in Turbo Pascal 3.01A (3)

Division

The easy way to implement integer division is to use repeated subtraction. I tried this – but even on a Raspberry Pi B under Free Pascal it’s (unsurprisingly) incredibly slow, especially if the denominator is very much smaller than the numerator,

...
      quotient := '0';
      remainder := num;
      while ge(remainder,den) do
      begin
        remainder := sub(remainder,den);
        quotient := add(quotient, '1')
      end;
...

There are many, many algorithms that improve on this – almost all of which need serious study to make them comprehensible. I intend to revisit these if I ever make a version 2 of this library as they should be faster than the method I eventually implemented, based on a recursive algorithm written in C by Bastian Molkenthin. (The Turbo Pascal 3.0 reference manual for CP/M warns against the use of recursion as it’s slow. Quite how slow will be revealed in part 4!).

Here’s the Pascal realisation of Bastian’s method.

var quott, remt: BigInt; { Used as var parameters so can't be declared  }
                         { locally in fastdiv - see TP reference manual }
procedure fastdiv(num: BigInt; den: BigInt; var quot: BigInt; var rem: BigInt);
var dent: BigInt;
begin
  if lt(num,den) = TRUE then   { numerator < denominator, therefore quotient is 0 }
  begin
    quot := '0';
    rem := num
  end
  else
  begin
    dent := multiply(den,'2');  { double the denominator and make recursive call }
    fastdiv(num,dent,quott,remt);
    if lt(remt,den) then        { remainder < denominator, so quotient is doubled }
    begin
      quot := multiply(quott,'2');
      rem := remt
    end
    else                           { remainder >= denominator, so quotient is double }
    begin                          { the remainder plus 1 and the remainder is }
      quot := multiply(quott,'2'); { subtracted from the denominator }
      quot := add(quot, '1');
      rem := sub(remt,den)
    end;
  end;
end;

Modulo

While there are a number of different definitions of modulo, this implements the same definition as Turbo Pascal’s native integer mod function. This is simply the remainder from an integer division, which is positive if the numerator is positive and negative if the numerator is negative, regardless of the sign of the denominator (which is always treated as being positive).

Square root

Straightforward now that integer division is implemented, while remembering to check for attempts to find the square roots of negative numbers, 0 or 1.

...
    if (num[1]='-') then                       { Square roots of negative numbers are complex ! }
      isqrt := 'NaN'
    else
    begin
      if (eq(num,'0') or eq(num,'1')) then     { Trap square root of 0 and 1 }
        isqrt := num
      else
      begin
        t0 := divide(num,'2');
        t1 := divide(num,t0);
        t1 := add(t0,t1);
        t1 := divide(t1,'2');
        while (lt(t1,t0)=TRUE) do              { Keep dividing until the square root is found }
        begin
          t0 := t1;
          t1 := divide(num,t0);
          t1 := add(t0,t1);
          t1 := divide(t1,'2');
        end;
        isqrt := t0
...

The mildly exciting news is that we now have enough of a library implemented to benchmark a Sieve of Eratosthenes prime number algorithm.