精确计算2的平方根

精确计算2的平方根(小数点后10000位) 计算2的平方根，用到下面的公式： 1/sqrt(1-x)=1+ 1/2 * x + (1*3)/(2*4) * x^2+ (1*3*5)/(2*4*6) * x^3 + ……… 理论上代入x=1/2 时就得到了 2的平方根。但是收敛会太慢。为了使得收敛的速度更快，这里， 程序选择了 x=1/57122 可以代入左边的式子看看它和sqrt(2)有什么关系，57122是怎么找出来的？呵呵！当然不是我 想出来的，而是用编程的找到的。类似的，如让 x= 1/50 , x= 1/1682 ,都远比 x=1/2 收敛速度要快。 计算出来的最后，我们只要乘回一个系数就达到目的了。 program sqrt_2; label ext; const      dn=2504; var     i,k:longint;     sum,a:array[1..dn] of integer;     ip:integer; procedure outp; var   i,m:integer; procedure testm; begin      if m mod 10=0 then write(' ');      if (m mod 50=0) and (m mod 1000<>0) then               writeln(':',m:8);      if m mod 1000<>0 then exit;      writeln(':',m:8,'  Press Enter..');      readln; end; procedure writep(num:integer); begin     write(num div 1000);   m:=m+1;   testm;     write(num div 100 mod 10); m:=m+1;  testm;     write(num div 10 mod 10);  m:=m+1;  testm;     write(num mod 10);  m:=m+1;  testm; end; begin    writeln('sqrt(2)=');    writeln(sum[1],'.');    m:=0;    for i:=2 to dn do         writep(sum[i]) ;    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 sm_div(k:longint); var    i:integer;    r1,c:longint; begin    c:=0;    for i:=1 to dn do       begin          r1:=c*10000+sum[i];          sum[i]:=r1 div k;          c:=r1 mod k;       end; end; procedure m_mul(k:longint); var    i:integer;    r1,c:longint; begin    c:=0;    for i:=dn downto ip do       begin          r1:=k*a[i]+c;          a[i]:=r1 mod 10000;          c:=r1 div 10000;       end;      if ip=1 then         begin            a[1]:=c;            exit;         end;      i:=ip-1;      while c>0 do       begin          if i=1 then             begin                a[1]:=c;                exit             end;          a[i]:=c mod 10000;          c:=c div 10000;          i:=i-1;       end; end; procedure sm_mul(k:longint); var    i:integer;    r1,c:longint; begin    c:=0;    for i:=dn downto 2 do       begin          r1:=k*sum[i]+c;          sum[i]:=r1 mod 10000;          c:=r1 div 10000;       end;       r1:=1;       r1:=r1*k*sum[1]+c;       sum[1]:=r1 ; end; procedure m_add; var    i:integer;    c:longint; begin    c:=0;    for i:=dn downto 1 do      begin         sum[i]:=sum[i]+a[i]+c;         if sum[i]>=10000 then           begin              c:=sum[i] div 10000;              sum[i]:=sum[i] mod 10000;           end         else c:=0;      end; end; begin     writeln('Please wait...');     for i:=1 to dn do a[i]:=0;     a[1]:=1;     sum:=a;     k:=1;     ip:=1;     repeat        i:=ip;        while (a[i]=0) do          begin             i:=i+1;             if i>dn then goto ext;          end;         ip:=i;         m_div(2*k);         m_mul(2*k-1);         m_div(57122);         m_add;         k:=k+1;     until false; ext:     sm_div(169);     sm_mul(239);     outp; end.