正文

整理我以前的PASCAL源程序-高精度计算(1)计算圆周率2010-08-19 22:05:00

【评论】 【打印】 【字体: 】 本文链接:http://blog.pfan.cn/jtchang/51679.html

分享到:

序言   看着这些程序,使我想起了以前DOS时代。这些PASCAL程序是我多年前编的。有些在论坛贴过,有些没贴出来。最近在整理,希望对以后学PASCAL的人有用。    高精度计算(1)计算圆周率   以前就一个16位的编译器TP7,数组开大一点还不行。当时是想开一个数组,把一些常数计算得位数更从一些。没有用别人的算法,没有用别人的代码,就自己编,不求太高的效率。   最先计算的是圆周率,失败了几次,最后程序成功地算出一千位来。当时心情真的很兴奋。以前网络没有普及,拿着一本书,在那里对着数位,直到书上数位没有了,程序仍然显示着后面的数位……     后来改了几次。下面的程序,用Machin公式,在TP7下可以算出一万位圆周率来。 program spi; { pi=16acrtg(1/5)-4arctg(1/239)   } { pi/4=1/1*(20/25-239/57121)-1/3(20/25^2-239/57121^2)+... } { Programmed by j.t.Chang.} label ext,ext2; const       dn=2504; var    i,ip,c:integer;    k:longint;    a,b,sum:array[1..dn] of integer; procedure oupt; var i:integer;     k:longint; procedure testk; var    ch:char; begin      if k mod 10=0 then write(' ');      if (k mod 50=0) and (k mod 1000<>0) then               writeln(':',k:8);      if k mod 1000<>0 then exit;      writeln(':',k:8,'  Press Enter..');      readln; end; begin    k:=0;    writeln('PI=');    writeln(sum[1],'.');    for i:=2 to dn do      begin        write(sum[i] div 1000);  k:=k+1; testk;        write(sum[i] div 100 mod 10); k:=k+1; testk;        write(sum[i] div 10 mod 10 ); k:=k+1; testk;        write(sum[i] mod 10);  k:=k+1; testk;      end;    writeln; writeln;    writeln('Programmed by j.t.chang.'); end; procedure m_div(k:longint); var     i:integer;     r1,c:longint; begin    c:=0;    for i:=ip to dn do      begin         r1:=c*10000+a[i];         a[i]:=r1 div k;         c:=r1 mod k;      end; end; procedure m_div2(k:longint); var     i:integer;     r1,c:longint; begin    c:=0;    for i:=ip to dn do      begin         r1:=c*10000+a[i];         b[i]:=r1 div k;         c:=r1 mod k;      end; end; procedure sum_add; var    i:integer;    r1,c:longint; begin    c:=0;    for i:= dn downto ip do     begin         r1:=c+b[i]+sum[i];         if r1>9999 then           begin               sum[i]:=r1 - 10000;               c:=1;           end          else            begin              sum[i]:=r1;              c:=0;            end;     end;    if c=0 then exit;    i:=ip-1;    while c>0 do       begin          if i=1 then             begin                sum[1]:=c+sum[1];                exit;             end          else             begin               r1:=c+sum[i];               if r1>9999 then                 begin                     sum[i]:=r1 - 10000;                     c:=1;                  end                else                   begin                     sum[i]:=r1;                     c:=0;                   end;             end;           i:=i-1;       end; end; procedure sum_sub; var    i:integer;    r1,c:longint; begin    c:=0;    for i:= dn downto ip do     begin         r1:=c-b[i]+sum[i];         if r1<0 then           begin               sum[i]:=r1+10000;               c:=-1;           end          else            begin              sum[i]:=r1;              c:=0;            end;     end;    if c=0 then exit;    i:=ip-1;    while c<0 do       begin          if i=1 then             begin                sum[1]:=c+sum[1];                exit;             end          else             begin                r1:=c+sum[i];                if r1<0 then                  begin                     sum[i]:=r1+10000;                     c:=-1;                  end               else                   begin                      sum[i]:=r1;                      c:=0;                   end;             end;           i:=i-1;       end; end; procedure proc; var   r1,c:longint;   i:integer; begin     c:=0;     for i:=dn downto 1 do       begin         r1:=1;         r1:=c+r1*sum[i]*4;         if r1>9999 then           begin              sum[i]:=r1 mod 10000;              c:=r1 div 10000;           end         else            begin              sum[i]:=r1;              c:=0;            end;       end; end; begin    writeln('Please wait...');    for i:=1 to dn do sum[i]:=0;    a:=sum;    a[1]:=20;    k:=1;    ip:=1;    c:=1;    repeat       i:=ip;       while  a[i]=0 do         begin             i:=i+1;             if i>dn then goto ext;         end;       ip:=i;       m_div(25);       m_div2(k);       if c=1 then  sum_add        else sum_sub;       k:=k+2;       c:=-c;    until false; ext:    for i:=1 to dn do a[i]:=0;    a[1]:=239;    b:=a;    k:=1;    ip:=1;    c:=-1;    repeat       i:=ip;       while  a[i]=0 do         begin             i:=i+1;             if i>dn then goto ext2;         end;       ip:=i;       m_div(57121);       m_div2(k);       if c=1 then  sum_add        else sum_sub;       k:=k+2;       c:=-c;    until false; ext2:       proc;       oupt; end.

阅读(2039) | 评论(0)


版权声明:编程爱好者网站为此博客服务提供商,如本文牵涉到版权问题,编程爱好者网站不承担相关责任,如有版权问题请直接与本文作者联系解决。谢谢!

评论

暂无评论
您需要登录后才能评论,请 登录 或者 注册