home *** CD-ROM | disk | FTP | other *** search
/ PC-X 1997 October / pcx14_9710.iso / swag / math.swg / 0099_Calculation PI.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1995-09-04  |  6.9 KB  |  325 lines
  1. (*
  2.    The two units used should come after this message. Uncomment several write-
  3. commands to get a "fully" operational program rather than this benchmark
  4. version. You then also can skip the Timer unit and the two commands from that
  5. unit (TimerOn and TimerOff) to make the program much smaller (no float math
  6. linked into the program).
  7. *)
  8.  
  9. program PiCalc;  { The fastest PI calculator you'll ever find... :) }
  10.  
  11. { From bits and pieces picked up mainly from the FidoNet PASCAL echo }
  12. { Collected, optimized, unitized, etc. by Bjorn Felten @ 2:203:208 }
  13. { Public Domain  --  Nov 1994 }
  14.  
  15. { Units needed are at the end !! }
  16.  
  17. uses HugeUtil, Timer; { use Crt if you want fast printout on screen }
  18.                       { don't if you want to be able to redirekt o/p }
  19.  
  20. var
  21.     words, number   : longint;
  22.     nin, link, pii, a239    : HugePtr;
  23.  
  24. procedure ArcCoTan(n : integer; var angle : Huge);
  25. var n2, del, remain : integer;
  26.     positive : boolean;
  27.  
  28. begin                               { corresp. integer operations }
  29.   ZeroHuge(angle,words);            { angle := 0 }
  30.   ZeroHuge(nin^,words);             { nin   := 0 }
  31.   ZeroHuge(link^,words);            { link  := 0 }
  32.   angle.dat[angle.len] := 1;        { angle := 1 }
  33.   DivHuge(angle,n,angle,remain);    { angle := angle div n }
  34.   n2 := n*n;                        { n2    := n * n }
  35.   del := 1;                         { del   := 1 }
  36.   positive := true;
  37.   CopyHuge(angle,nin^);             { nin   := angle }
  38.   repeat
  39.     DivHuge(nin^,n2,nin^,remain);   { nin   := nin div n2 }
  40.     inc(del, 2);                    { del   := del + 2 }
  41.     positive := not positive;
  42.     DivHuge(nin^,del,link^,remain); { link  := nin div del }
  43.     if positive then
  44.       AddHuge(angle,link^)          { angle := angle + link }
  45.     else
  46.       SubHuge(angle,link^);         { angle := angle - link }
  47. {    write(#13,word(del)) } { uncomment to see that program is not dead }
  48.   until (link^.len <= 1) and (link^.dat[1] = 0);
  49. {  writeln}                 { ... and this too }
  50. end; { ArcCoTan }
  51.  
  52. begin
  53. {  writeln('Program to get Pi (',pi:1:17,'...) with large precision.'); }
  54.   write('Digits(max 40.000): '); readln(number);
  55.   words := round(number / 4.7) + 3; { appr. 4.7 digits in one word }
  56.   write(number:6,#9);
  57.   TimerOn;
  58.   GetHuge(pii,  words+2);
  59.   GetHuge(a239, words+2);
  60.   GetHuge(link, words+2);
  61.   GetHuge(nin,  words+2);
  62.   ArcCoTan(5,   pii^);        { ATan(1/5)  }
  63.   AddHuge(pii^, pii^);
  64.   AddHuge(pii^, pii^);        { * 4        }
  65.   ArcCoTan(239, a239^);       { ATan(1/239)}
  66.   SubHuge(pii^, a239^);
  67.   AddHuge(pii^, pii^);
  68.   AddHuge(pii^, pii^);        { * 4        }
  69.   TimerOff;
  70. {  WriteHuge(pii^, number)}     { uncomment if you want printout }
  71. end.
  72.  
  73. unit HugeUtil;
  74.  
  75. interface
  76.  
  77. const HugeMax = $8000-16;
  78.  
  79. type  Huge = record
  80.               len : word;
  81.               dat : array[1..HugeMax] of word;
  82.             end;
  83.       HugePtr = ^Huge;
  84.  
  85. procedure AddHuge  (var Answer, Add : Huge);
  86. procedure MulHuge  (var A : Huge; Mul : integer; var Answer : Huge);
  87. procedure DivHuge  (var A : Huge; Del : integer; var Answer : Huge;
  88.                     var Remainder : integer);
  89. procedure SubHuge  (var Answer, Sub : Huge);
  90. procedure ZeroHuge (var L : Huge; Size : word);
  91. procedure CopyHuge (var Fra,Til : Huge);
  92. procedure GetHuge  (var P : HugePtr; Size : word);
  93. procedure WriteHuge(var L : Huge; Size: word);
  94.  
  95. implementation
  96.  
  97. procedure AddHuge; assembler; asm
  98.   cld
  99.   push  ds
  100.   lds   di,Answer
  101.   les   si,Add
  102.   seges lodsw
  103.   mov   cx,ax
  104.   clc
  105. @l1:
  106.   seges lodsw
  107.   adc   [si-2],ax
  108.   loop  @l1
  109.   jnb   @done
  110. @l2:
  111.   add   word [si],1
  112.   inc   si
  113.   inc   si
  114.   jc    @l2
  115. @done:
  116.   mov   si,di
  117.   lodsw
  118.   shl   ax,1
  119.   add   si,ax
  120.   lodsw
  121.   or    ax,ax
  122.   je    @d2
  123.   inc   word [di]
  124. @d2:
  125.   pop   ds
  126. end;
  127.  
  128. procedure MulHuge; assembler; asm
  129.   cld
  130.   push  ds
  131.   lds   si,A
  132.   mov   bx,Mul
  133.   les   di,Answer
  134.   mov   cx,[si]
  135.   mov   dx,si
  136.   inc   di
  137.   inc   di
  138.   clc
  139. @l1:
  140.   mov   ax,[di]
  141.   pushf
  142.   mul   bx
  143.   popf
  144.   adc   ax,si
  145.   stosw
  146.   mov   si,dx
  147.   loop  @l1
  148.   adc   si,0
  149.   mov   es:[di],si
  150.   lds   di,A
  151.   mov   di,[di]
  152.   mov   ax,[di+2]
  153.   or    ax,ax
  154.   je    @l2
  155.   inc   di
  156.   inc   di
  157. @l2:
  158.   lds   si,Answer
  159.   mov   [si],di
  160.   pop   ds
  161. end;
  162.  
  163. procedure DivHuge; assembler; asm
  164.   std
  165.   push  ds
  166.   lds   si,A
  167.   mov   bx,Del
  168.   les   di,Answer
  169.   mov   cx,[si]
  170.   mov   di,cx
  171.   add   di,cx
  172.   xor   dx,dx
  173. @l1:
  174.   mov   ax,[di]
  175.   div   bx
  176.   stosw
  177.   loop  @l1
  178.   lds   si,Remainder
  179.   mov   [si],dx
  180.   lds   si,A
  181.   mov   ax,[si]
  182.   lds   di,Answer
  183.   mov   [di],ax
  184.   mov   si,[di]
  185.   shl   si,1
  186. @d3:
  187.   lodsw
  188.   or    ax,ax
  189.   jne   @d2
  190.   dec   word [di]
  191.   jne   @d3
  192.   inc   word [di]
  193. @d2:
  194.   pop   ds
  195. end;
  196.  
  197. procedure SubHuge; assembler; asm
  198.   cld
  199.   push  ds
  200.   lds   di,Answer
  201.   les   si,Sub
  202.   seges lodsw
  203.   mov   cx,ax
  204.   clc
  205. @l1:
  206.   seges lodsw
  207.   sbb   [si-2],ax
  208.   loop  @l1
  209.   jnb   @done
  210. @l2:
  211.   sub   word [si],1
  212.   inc   si
  213.   inc   si
  214.   jc    @l2
  215. @done:
  216.   mov   si,[di]
  217.   shl   si,1
  218.   std
  219. @d3:
  220.   lodsw
  221.   or    ax,ax
  222.   jne   @d2
  223.   dec   word [di]
  224.   jne   @d3
  225.   inc   word [di]
  226. @d2:
  227.   pop   ds
  228. end;
  229.  
  230.  
  231. procedure WriteHuge;
  232. var L1, L2, I, R, R1, X : integer;
  233. begin
  234.   with L do begin
  235.     L1 := Len;
  236.     L2 := L1 - 1;
  237.     I := 1;
  238.     write(dat[L1],'.');
  239.     X := 0;
  240.     for I := 1 to Size div 4 do begin
  241.       Dat[L1] := 0;
  242.       Len := L2;
  243.       MulHuge(L,10000,L);
  244.       R := dat[L1];
  245.       R1 := R div 100;
  246.       R  := R mod 100;
  247.       write(chr(R1 div 10+48), chr(R1 mod 10+48),
  248.             chr(R  div 10+48), chr(R  mod 10+48));
  249.       inc(X);
  250.       write(' ');
  251.       if X > 14 then begin
  252.         writeln; write('  ');
  253.         X := 0
  254.       end
  255.     end
  256.   end;
  257.   writeln
  258. end;                            { WriteHuge }
  259.  
  260. procedure ZeroHuge;
  261. begin
  262.   fillchar(L.Dat, Size * 2, #0);
  263.   L.Len := Size
  264. end;
  265.  
  266. procedure CopyHuge;
  267. begin
  268.   move(Fra, Til, Fra.Len * 2 + 2)
  269. end;
  270.  
  271. procedure GetHuge;
  272. var D : ^byte;
  273.     Tries,
  274.     Bytes : word;
  275. begin
  276.   Bytes := 2 * (Size + 1);
  277.   Tries:=0;
  278.   repeat
  279.     getmem(P,Bytes);
  280.  
  281. { To make it possible to use maximally large arrays, and to increase
  282.   the speed of the computations, all records of type Huge MUST start
  283.   at a segment boundary! }
  284.  
  285.     if ofs(P^) = 0 then begin
  286.       ZeroHuge(P^,Size);
  287.       exit
  288.     end;
  289.     inc(Tries);
  290.     freemem(P,Bytes);
  291.     new(D)
  292.   until Tries>10;   { if not done yet, it's not likely we ever will be }
  293.   writeln('Couldn''t get memory for array');
  294.   halt(1)
  295. end;                                   { GetHuge }
  296.  
  297. end.
  298.  
  299. unit Timer;
  300.  
  301. interface
  302.  
  303. procedure TimerOn;
  304. procedure TimerOff;
  305.  
  306. implementation
  307.  
  308. var
  309.   Time      : Longint absolute $0040:$006C;
  310.   WaitTime,
  311.   Temp      : Longint;
  312.  
  313. procedure TimerOn;
  314. begin
  315.   WaitTime:=Time
  316. end;
  317.  
  318. procedure TimerOff;
  319. begin
  320.   Temp:=Time;
  321.   writeln('Done! It took ',(Temp-WaitTime)/18.2:6:2,'s.')
  322. end;
  323.  
  324. end.
  325.