[Back to MATH SWAG index] [Back to Main SWAG index] [Original]

    (*
       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.

[Back to MATH SWAG index] [Back to Main SWAG index] [Original]