精确计算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
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;
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;
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;
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;
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.。