Announcement

Collapse
No announcement yet.

calculate with large numbers

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • calculate with large numbers

    Hello,

    i need some help to translate the following Pascal code to
    pb/dll-code. i need the fast calculating routines, because i want
    to calculate with numbers with up to 100 digits or more. the timer
    functions are not needed.
    Code:
    program PiCalc;
    
    { 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.
    ------------------

  • #2
    'k. I gather Bjorn Felten's HugeUtil unit is intended for free public
    distribution. We have the occasional request for this sort of functionality.
    I'll take a look at it.

    ------------------
    Tom Hanlin
    PowerBASIC Staff

    Comment


    • #3
      The DivHuge implementation appears to be outright broken, which doesn't
      inspire a lot of confidence in these routines. I'll keep digging, though.


      ------------------
      Tom Hanlin
      PowerBASIC Staff

      Comment

      Working...
      X