(*
   The two units used should come after this message. Uncomment several write-
commands to get a "fully" operational program rather than this benchmark
version. You then also can skip the Timer unit and the two commands from that
unit (TimerOn and TimerOff) to make the program much smaller (no float math
linked into the program).
*)

program PiCalc;  { The fastest PI calculator you'll ever find... :) }

{ From bits and pieces picked up mainly from the FidoNet PASCAL echo }
{ Collected, optimized, unitized, etc. by Bjorn Felten @ 2:203:208 }
{ Public Domain  --  Nov 1994 }

{ Units needed are at the end !! }

uses HugeUtil, Timer; { use Crt if you want fast printout on screen }
                      { don't if you want to be able to redirekt o/p }

var
    words, number   : longint;
    nin, link, pii, a239    : HugePtr;

procedure ArcCoTan(n : integer; var angle : Huge);
var n2, del, remain : integer;
    positive : boolean;

begin                               { corresp. integer operations }
  ZeroHuge(angle,words);            { angle := 0 }
  ZeroHuge(nin^,words);             { nin   := 0 }
  ZeroHuge(link^,words);            { link  := 0 }
  angle.dat[angle.len] := 1;        { angle := 1 }
  DivHuge(angle,n,angle,remain);    { angle := angle div n }
  n2 := n*n;                        { n2    := n * n }
  del := 1;                         { del   := 1 }
  positive := true;
  CopyHuge(angle,nin^);             { nin   := angle }
  repeat
    DivHuge(nin^,n2,nin^,remain);   { nin   := nin div n2 }
    inc(del, 2);                    { del   := del + 2 }
    positive := not positive;
    DivHuge(nin^,del,link^,remain); { link  := nin div del }
    if positive then
      AddHuge(angle,link^)          { angle := angle + link }
    else
      SubHuge(angle,link^);         { angle := angle - link }
{    write(#13,word(del)) } { uncomment to see that program is not dead }
  until (link^.len <= 1) and (link^.dat[1] = 0);
{  writeln}                 { ... and this too }
end; { ArcCoTan }

begin
{  writeln('Program to get Pi (',pi:1:17,'...) with large precision.'); }
  write('Digits(max 40.000): '); readln(number);
  words := round(number / 4.7) + 3; { appr. 4.7 digits in one word }
  write(number:6,#9);
  TimerOn;
  GetHuge(pii,  words+2);
  GetHuge(a239, words+2);
  GetHuge(link, words+2);
  GetHuge(nin,  words+2);
  ArcCoTan(5,   pii^);        { ATan(1/5)  }
  AddHuge(pii^, pii^);
  AddHuge(pii^, pii^);        { * 4        }
  ArcCoTan(239, a239^);       { ATan(1/239)}
  SubHuge(pii^, a239^);
  AddHuge(pii^, pii^);
  AddHuge(pii^, pii^);        { * 4        }
  TimerOff;
{  WriteHuge(pii^, number)}     { uncomment if you want printout }
end.

unit HugeUtil;

interface

const HugeMax = $8000-16;

type  Huge = record
              len : word;
              dat : array[1..HugeMax] of word;
            end;
      HugePtr = ^Huge;

procedure AddHuge  (var Answer, Add : Huge);
procedure MulHuge  (var A : Huge; Mul : integer; var Answer : Huge);
procedure DivHuge  (var A : Huge; Del : integer; var Answer : Huge;
                    var Remainder : integer);
procedure SubHuge  (var Answer, Sub : Huge);
procedure ZeroHuge (var L : Huge; Size : word);
procedure CopyHuge (var Fra,Til : Huge);
procedure GetHuge  (var P : HugePtr; Size : word);
procedure WriteHuge(var L : Huge; Size: word);

implementation

procedure AddHuge; assembler; asm
  cld
  push  ds
  lds   di,Answer
  les   si,Add
  seges lodsw
  mov   cx,ax
  clc
@l1:
  seges lodsw
  adc   [si-2],ax
  loop  @l1
  jnb   @done
@l2:
  add   word [si],1
  inc   si
  inc   si
  jc    @l2
@done:
  mov   si,di
  lodsw
  shl   ax,1
  add   si,ax
  lodsw
  or    ax,ax
  je    @d2
  inc   word [di]
@d2:
  pop   ds
end;

procedure MulHuge; assembler; asm
  cld
  push  ds
  lds   si,A
  mov   bx,Mul
  les   di,Answer
  mov   cx,[si]
  mov   dx,si
  inc   di
  inc   di
  clc
@l1:
  mov   ax,[di]
  pushf
  mul   bx
  popf
  adc   ax,si
  stosw
  mov   si,dx
  loop  @l1
  adc   si,0
  mov   es:[di],si
  lds   di,A
  mov   di,[di]
  mov   ax,[di+2]
  or    ax,ax
  je    @l2
  inc   di
  inc   di
@l2:
  lds   si,Answer
  mov   [si],di
  pop   ds
end;

procedure DivHuge; assembler; asm
  std
  push  ds
  lds   si,A
  mov   bx,Del
  les   di,Answer
  mov   cx,[si]
  mov   di,cx
  add   di,cx
  xor   dx,dx
@l1:
  mov   ax,[di]
  div   bx
  stosw
  loop  @l1
  lds   si,Remainder
  mov   [si],dx
  lds   si,A
  mov   ax,[si]
  lds   di,Answer
  mov   [di],ax
  mov   si,[di]
  shl   si,1
@d3:
  lodsw
  or    ax,ax
  jne   @d2
  dec   word [di]
  jne   @d3
  inc   word [di]
@d2:
  pop   ds
end;

procedure SubHuge; assembler; asm
  cld
  push  ds
  lds   di,Answer
  les   si,Sub
  seges lodsw
  mov   cx,ax
  clc
@l1:
  seges lodsw
  sbb   [si-2],ax
  loop  @l1
  jnb   @done
@l2:
  sub   word [si],1
  inc   si
  inc   si
  jc    @l2
@done:
  mov   si,[di]
  shl   si,1
  std
@d3:
  lodsw
  or    ax,ax
  jne   @d2
  dec   word [di]
  jne   @d3
  inc   word [di]
@d2:
  pop   ds
end;


procedure WriteHuge;
var L1, L2, I, R, R1, X : integer;
begin
  with L do begin
    L1 := Len;
    L2 := L1 - 1;
    I := 1;
    write(dat[L1],'.');
    X := 0;
    for I := 1 to Size div 4 do begin
      Dat[L1] := 0;
      Len := L2;
      MulHuge(L,10000,L);
      R := dat[L1];
      R1 := R div 100;
      R  := R mod 100;
      write(chr(R1 div 10+48), chr(R1 mod 10+48),
            chr(R  div 10+48), chr(R  mod 10+48));
      inc(X);
      write(' ');
      if X > 14 then begin
        writeln; write('  ');
        X := 0
      end
    end
  end;
  writeln
end;                            { WriteHuge }

procedure ZeroHuge;
begin
  fillchar(L.Dat, Size * 2, #0);
  L.Len := Size
end;

procedure CopyHuge;
begin
  move(Fra, Til, Fra.Len * 2 + 2)
end;

procedure GetHuge;
var D : ^byte;
    Tries,
    Bytes : word;
begin
  Bytes := 2 * (Size + 1);
  Tries:=0;
  repeat
    getmem(P,Bytes);

{ To make it possible to use maximally large arrays, and to increase
  the speed of the computations, all records of type Huge MUST start
  at a segment boundary! }

    if ofs(P^) = 0 then begin
      ZeroHuge(P^,Size);
      exit
    end;
    inc(Tries);
    freemem(P,Bytes);
    new(D)
  until Tries>10;   { if not done yet, it's not likely we ever will be }
  writeln('Couldn''t get memory for array');
  halt(1)
end;                                   { GetHuge }

end.

unit Timer;

interface

procedure TimerOn;
procedure TimerOff;

implementation

var
  Time      : Longint absolute $0040:$006C;
  WaitTime,
  Temp      : Longint;

procedure TimerOn;
begin
  WaitTime:=Time
end;

procedure TimerOff;
begin
  Temp:=Time;
  writeln('Done! It took ',(Temp-WaitTime)/18.2:6:2,'s.')
end;

end.