高精度计算
高精度计算
高精度运算,是指参与运算的数(加数,减数,因子……)范围大大超出了标准数据类型(整型,实型)能表示的范围的运算。例如,求两个200位的数的和。这时,就要用到高精度算法了。
概述
所谓的高精度运算,是指参与运算的数(加数,减数,因子……)范围大大超出了标准数据类型(整型,实型)能表示的范围的运算。例如,求两个200位的数的和。这时,就要用到高精度算法了。在这里,我们先讨论高精度加法。
高精度加法
高精度运算主要解决以下三个问题:
一、加数、减数、运算结果的输入和存储
运算因子超出了整型、实型能表示的范围,肯定不能直接用一个数的形式来表示。在Pascal中,能表示多个数的数据类型有两种:数组和字符串。
数组:每个数组元素存储1位(在优化时,这里是一个重点!),有多少位就需要多少个数组元素;用数组表示数的优点:每一位都是数的形式,可以直接加减;运算时非常方便。用数组表示数的缺点:数组不能直接输入;输入时每两位数之间必须有分隔符,不符合数值的输入习惯;
字符串:String型字符串的最大长度是255,可以表示255位。Ansistring型字符串长度不受限制。用字符串表示数的优点:能直接输入输出,输入时,每两位数之间不必分隔符,符合数值的输入习惯;用字符串表示数的缺点:字符串中的每一位是一个字符,不能直接进行运算,必须先将它转化为数值再进行运算;运算时非常不方便;
综合以上所述,对上面两种数据结构取长补短:用字符串读入数据,用数组存储数据:
var st:string;
x,y:array[0..255]of integer;{定义两个数组,X和Y,用来储存数}
i,j,l1,l2:integer;
begin
readln(st);
l1:=length(st);{------length(x),该函数是获取字符串X的长度,返回为整型}
for i:=0 to 255 do x[i]:=0;{数组初始化,该句等价于‘fillchar(x,sizeof(x),o);’,即给一数组整体赋值,但运行速度快于用‘for’语句对数组中的每一个数赋值}
for i:=l1 downto 1 do
x[l1-i+1]:=ord(st[i])-ord(''0'');{------这里是重点,把字符串转换为数值,储存在数组中}
readln(st);
l2:=length(st);{------length(x),该函数是获取字符串X的长度,返回为整型}
for i:=0 to 255 do y[i]:=0;{数组初始化,该句等价于‘fillchar(y,sizeof(y),o);’}
for i:=l2 downto 1 do
y[l2-i+1]:=ord(st[i])-ord(''0'');{------这里是重点,把字符串转换为数值,储存在数组中}
对字符串转为数值原理补充:ord(x)-48,如果X=''1'',因为''1''的ASCLL码是49,所以减去48就等于1,间接地把字符转换为数值了,各位初手要好好体会.
二、运算过程
在往下看之前,大家先列竖式计算35+86。
注意的问题:
(1)运算顺序:两个数靠右对齐;从低位向高位运算;先计算低位再计算高位;
(2)运算规则:同一位的两个数相加再加上从低位来的进位,成为该位的和;这个和去掉向高位的进位就成为该位的值;如上例:3+8+1=12,向前一位进1,本位的值是2;可借助MOD、DIV运算完成这一步;
(3)最后一位的进位:如果完成两个数的相加后,进位位值不为0,则应添加一位;
(4)如果两个加数位数不一样多,则按位数多的一个进行计算;
if l1<l2 then l1:=l2;
for i:=1 to l1 do
begin
x[i]:=x[i]+y[i];
x[i+1]:=x[i+1]+x[i] div 10;
x[i]:=x[i] mod 10;
end;
三、结果的输出(这也是优化的一个重点)
按运算结果的实际位数输出
var st:string;
x,y:array[0..255]of integer;
i,j,l1,l2:integer;
begin
readln(st);
l1:=length(st);
for i:=0 to 255 do x[i]:=0;
for i:=l1 downto 1 do
x[l1-i+1]:=ord(st[i])-ord(''0'');
readln(st);
l2:=length(st);
for i:=0 to 255 do y[i]:=0;
for i:=l2 downto 1 do
y[l2-i+1]:=ord(st[i])-ord(''0'');
if l1<l2 then l1:=l2;
for i:=1 to l1 do
begin
x[i]:=x[i]+y[i];
x[i+1]:=x[i+1]+x[i] div 10;
x[i]:=x[i] mod 10;
end;
write(''x+y='');
j:=255;
while x[j]=0 do j:=j-1;
for i:=j downto 1 do write(x[i]);
readln;
end.
四、优化:
以上的方法的有明显的缺点:
(1)浪费空间:一个整型变量(-32768~32767)只存放一位(0~9);
(2)浪费时间:一次加减只处理一位;
针对以上问题,我们做如下优化:一个数组元素存放四位数;(integer的最大范围是32767,5位的话可能导致出界)将标准数组改为紧缩数组。第一步的具体方法:
l:=length(s1);
k1:=260;
repeat {————有关字符串的知识}
s:=copy(s1,l-3,4);
val(s,a[k1],code);
k1:=k1-1;
s1:=copy(s1,1,l-4);
l:=l-4;
until l<=0;
k1:=k1+1;
而因为这个改进,算法要相应改变:
(1)运算时:不再逢十进位,而是逢万进位(mod 10000; div 10000);
(2)输出时:最高位直接输出,其余各位,要判断是否足够4位,不足部分要补0;例如:1,23,2345这样三段的数,输出时,应该是100232345而不是1232345。
改进后的算法:
program add;
var s,s1,s2:string;
a,b,c:array[1..260]of integer; {260是随便取的}
x,y,k,i,l,k1,k2,code:longint;
begin
readln(s1);
readln(s2);
l:=length(s1);
k1:=260;
repeat
s:=copy(s1,l-3,4);
val(s,a[k1],code);
dec(k1);
s1:=copy(s1,1,l-4);
l:=l-4;
until l<=0;
k1:=k1+1;
l:=length(s2);
k2:=260;
repeat
s:=copy(s2,l-3,4);
val(s,b[k2],code);
dec(k2);
s2:=copy(s2,1,l-4);
dec(l,4);
until l<=0;
inc(k2);
if k1<k2 then k:=k1 else k:=k2;
y:=0;
for i:=260 downto k do
begin
x:=a[i]+b[i]+y;
c[i]:=x mod 10000;
y:=x div 10000;
end;
write(c[k]);
for i:=k+1 to 260 do
begin
if c[i]<1000 then write(''0'');
if c[i]<100 then write(''0'');
if c[i]<10 then write(''0'');
write(c[i]);
end;
end.
高精度减法
和高精度加法相比,减法在差为负数时处理的细节更多一点:当被减数小于减数时,差为负数,差的绝对值是减数减去被减数;在程序实现上用一个变量来存储符号位,用另一个数组存差的绝对值。
算法流程:(1).读入被减数S1,S2(字符串);
(2).置符号位:判断被减数是否大于减数:大则将符号位置为空;小则将符号位置为“- ”,交换减数与被减数;
(3).被减数与减数处理成数值,放在数组中;
(4).运算:A、取数;
B、判断是否需要借位;
C、减,将运算结果放到差数组相应位中;
D、判断是否运算完成:是,转5;不是,转A;
(5).打印结果:符号位,第1位,循环处理第2到最后一位;
细节:▲如何判断被减数与减数的大小?
如果位数一样,直接比较字符串大小;否则,位数多的大。
k1:=length(s1); k2:=length(s2);
if k1=k2 then
if s1<s2 then begin fh:=''-''; s:=s1;s1:=s2; s2:=s;end
else if k1<k2 then begin fh:=''-'';s:=s1;s1:=s2;s2:=s;end;{s1存被减数,fh存符号}
将字符串处理成数值:
l:=length(s1);{求出s1的长度,也即s1的位数;有关字符串的知识。}
k1:=260;
for i:=l downto 1 do
begin
a[k1]:=ord(s1[i])-48;{将字符转成数值}
k1:=k1-1;
end;
k1:=k1+1;
运算(减法跟加法比较,减法退位处理跟加法进位处理不一样):
处理退位: 跟加法一样,在for语句外面先将退位清零,用被减数再减去退位,{注意:由于每一个数位不一定都得向前一位借位,所以这里退位得清零。例如,234-25,个位需借位,而十位不用} 接着,再判断,当被减数某一位不够减时,则需加上前一位退位过来的数。注意:由于这里采用优化方法,所以退一位,就等于后一位加上10000。)最后,再拿一个数组来存储两个减数的差。
jw:=0;
for i:=260 downto k1 do
begin
a[i]:=a[i]-jw;{此处jw为从刚处理的那一位上从本一位上的借位}
jw:=0; {此处jw为I 位准备向高一位的借位}
if a[i]<b[i] then
begin
jw:=1;
a[i]:=a[i]+10000;
end;
c[i]:=a[i]-b[i]
end;
打印结果: 先找到差的第一个非零数,如果差的所有位数都为零,就直接输出零; 如果不是,就输出符号位和差的第一位。剩下部分,打印补足零;因为优化后的高精度减法,是把每四个数位分成一段,而每一段则必须有四个数,当有一段不足四个数时,就得用"0"补足.(如:第一位是''1'',第二位是''34'',第三位是''345'',第四位是''8'', 则应写为''1003403450008'').注意:第一位不用补零,(如:第一位为''3'',则写成''3'').
while (c[k]=0) and (k<=260) do k:=k+1;
if k>260 then write(''0'')
else begin
write(fh,c[k]);{k是差的第1位;}
for i:=k+1 to 260 do
begin
if c[i]<100 then write(''0'');
if c[i]<10 then write(''0'');
write(c[i]);
end;
end;
参考程序:
program ZDloveQC;
var s1,s2,s3,s4,s:string;
a,b,c:array[1..260]of integer;
i,k1,k2,l,code,jw:longint;
fh:string;
begin
readln(s1); readln(s2);
k1:=length(s1); k2:=length(s2); fh:='''';
if k1=k2 then
if s1<s2 then begin fh:=''-'';s:=s1; s1:=s2; s2:=s; end
else
if k1<k2 then begin fh:=''-'';s:=s1; s1:=s2; s2:=s; end;
k1:=260;
l:=length(s1);
repeat
s3:=copy(s1,l-3,4);
val(s3,a[k1],code);
dec(k1);
s1:=copy(s1,1,l-4);
l:=l-4;
until l<=0;
inc(k1);
l:=length(s2);
k2:=260;
repeat
s4:=copy(s2,l-3,4);
val(s4,b[k2],code);
dec(k2);
s2:=copy(s2,1,l-4);
l:=l-4;
until l<=0;
inc(k2);
jw:=0;
for i:=260 downto k1 do
begin
a[i]:=a[i]-jw;
jw:=0;
if a[i]<b[i] then
begin
jw:=1;
a[i]:=a[i]+10000;
end;
c[i]:=a[i]-b[i];
end;
while (c[k1]=0)and(k1<260) do inc(k1);
if k1>260 then writeln(''0'')
else begin
write(fh,c[k1]);
for i:=k1+1 to 260 do
begin
if c[i]<1000 then write(''0'');
if c[i]<100 then write(''0'');
if c[i]<10 then write(''0'');
write(c[i]);
end;
end;
end.
高精度乘法
高精度乘法基本思想和加法一样。其基本流程如下:
①读入被乘数s1,乘数s2
②把s1、s2分成4位一段,转成数值存在数组a,b中;记下a,b的长度k1,k2;
③i赋为b中的最低位;
④从b中取出第i位与a相乘,累加到另一数组c中;(注意:累加时错开的位数应是多少位?)
⑤i:=i-1;检测i值:小于k2则转⑥,否则转④
⑥打印结果
参考程序:
program chengfa;
const n=100;
type ar=array [1..n] of integer;
var a,b:ar; k1,k2,k:integer;
c:array [1..200] of integer;
s1,s2:string;
procedure fenge(s:string;var d:ar; var kk:integer); {将s分割成四位一组存放在d中,返回的kk值指向d的最高位}
var ss:string;
i,code:integer;
begin
i:=length(s);
kk:=n;
repeat
ss:=copy(s,i-3,4);
val(ss,d[kk],code);
kk:=kk-1;
s:=copy(s,1,i-4);
i:=i-4;
until i<0;
kk:=kk+1;
end;
procedure init;
var i:integer;
begin
for i:=1 to n do begin a:=0; b:=0; end;
for i:=1 to 2*n do c:=0;
write(''input 2 numbers:'');
readln(s1);
readln(s2);
fenge(s1,a,k1);
fenge(s2,b,k2);
end;
procedure jisuan;
var i,j,m:integer; x,y,z,jw:longint;
begin
i:=n; k:=2*n;
repeat
x:=b; z:=0; m:=k; jw:=0;
for j:=n downto k1 do
begin
y:=a[j];
z:=c[m];
x:=x*y+z+jw;
jw:=x div 10000;
c[m]:=x mod 10000;
m:=m-1;
x:=b;
end;
if jw<>0 then c[m]:=jw else m:=m+1;
i:=i-1;
k:=k-1;
until i<k2;
k:=m;
end;
procedure daying;
var i:integer;
begin
write(c[k]);
for i:=k+1 to 2*n do
begin
if c<1000 then write(''0'');
if c<100 then write(''0'');
if c<10 then write(''0'');
write(c);
end;
writeln;
end;
begin
init;
jisuan;
daying;
end.
教材“基础编”P87高精乘法参考程序:
program ex3_1;
var
a,b,c:array[0..1000] of word;
procedure init;
var
s:string;
ok,i,j:integer;
begin
readln(s);
a[0]:=length(s);
for i:=1 to a[0] do
val(s[a[0]-i+1],a,ok);
readln(s);
b[0]:=length(s);
b[0]:=length(s);
for i:=1 to b[0] do
val(s[b[0]-i+1],b,ok);
end;
procedure highmul;
var i,j,k:integer;
begin
c[0]:=a[0]+b[0];
for i:=1 to b[0] do
for j:=1 to a[0]+1 do
begin
inc(c[i+j-1],a[j]*b mod 10);
c[i+j]:=c[i+j]+(a[j]*b div 10)+(c[i+j-1] div 10);
c[i+j-1]:=c[i+j-1] mod 10;
end;
end;
procedure print;
var i:integer;
begin
while c[c[0]]=0 do dec(c[0]);
for i:=c[0] downto 1 do
write(c);
writeln;
end;
begin
init;
highmul;
print;
end.
高精度除法
高精度除法:
1).高精度除以整型数据(integer);
程序如下:
program HighPrecision3_Multiply1;
const
fn_inp=''hp5.inp'';
fn_out=''hp5.out'';
maxlen=100; { max length of the number }
type
hp=record
len:integer; { length of the number }
s:array[1..maxlen] of integer
{ s[1] is the lowest position
s[len] is the highest position }
end;
var
x,y:hp;
z,w:integer;
procedure PrintHP(const p:hp);
var i:integer;
begin
for i:=p.len downto 1 do write(p.s[i]);
end;
procedure init;
var
st:string;
i:integer;
begin
assign(input,fn_inp);
reset(input);
readln(st);
x.len:=length(st);
for i:=1 to x.len do { change string to HP }
x.s:=ord(st[x.len+1-i])-ord(''0'');
readln(z);
close(input);
end;
procedure Divide(a:hp;b:integer;var c:hp;var d:integer);
{ c:=a div b ; d:=a mod b }
var i,len:integer;
begin
fillchar(c,sizeof(c),0);
len:=a.len;
d:=0;
for i:=len downto 1 do { from high to low }
begin
d:=d*10+a.s[i];
c.s:=d div b;
d:=d mod b;
end;
while(len>1) and (c.s[len]=0) do dec(len);
c.len:=len;
end;
procedure main;
begin
Divide(x,z,y,w);
end;
procedure out_;
begin
assign(output,fn_out);
rewrite(output);
PrintHP(y);
writeln;
writeln(w);
close(output);
end;
begin
init;
main;
out_;
end.
2).高精度除以高精度
程序如下:
program HighPrecision4_Multiply2;
const
fn_inp=''hp6.inp'';
fn_out=''hp6.out'';
maxlen=100; { max length of the number }
type
hp=record
len:integer; { length of the number }
s:array[1..maxlen] of integer
{ s[1] is the lowest position
s[len] is the highest position }
end;
var
x:array[1..2] of hp;
y,w:hp; { x:input ; y:output }
procedure PrintHP(const p:hp);
var i:integer;
begin
for i:=p.len downto 1 do write(p.s);
end;
procedure init;
var
st:string;
j,i:integer;
begin
assign(input,fn_inp);
reset(input);
for j:=1 to 2 do
begin
readln(st);
x[j].len:=length(st);
for i:=1 to x[j].len do { change string to HP }
x[j].s:=ord(st[x[j].len+1-i])-ord(''0'');
end;
close(input);
end;
procedure Subtract(a,b:hp;var c:hp); { c:=a-b, suppose a>=b }
var i,len:integer;
begin
fillchar(c,sizeof(c),0);
if a.len>b.len then len:=a.len { get the bigger length of a,b }
else len:=b.len;
for i:=1 to len do { subtract from low to high }
begin
inc(c.s,a.s-b.s);
if c.s<0 then
begin
inc(c.s,10);
dec(c.s[i+1]); { add 1 to a higher position }
end;
end;
while(len>1) and (c.s[len]=0) do dec(len);
c.len:=len;
end;
function Compare(const a,b:hp):integer;
{
1 if a>b
0 if a=b
-1 if a < b
}
var len:integer;
begin
if a.len>b.len then len:=a.len { get the bigger length of a,b }
else len:=b.len;
while(len>0) and (a.s[len]=b.s[len]) do dec(len);
{ find a position which have a different digit }
if len=0 then compare:=0 { no difference }
else compare:=a.s[len]-b.s[len];
end;
procedure Multiply10(var a:hp); { a:=a*10 }
var i:Integer;
begin
for i:=a.len downto 1 do
a.s[i+1]:=a.s;
a.s[1]:=0;
inc(a.len);
while(a.len>1) and (a.s[a.len]=0) do dec(a.len);
end;
procedure Divide(a,b:hp;var c,d:hp); { c:=a div b ; d:=a mod b }
var i,j,len:integer;
begin
fillchar(c,sizeof(c),0);
len:=a.len;
fillchar(d,sizeof(d),0);
d.len:=1;
for i:=len downto 1 do
begin
Multiply10(d);
d.s[1]:=a.s; { d:=d*10+a.s}
{ c.s:=d div b ; d:=d mod b; }
{ while(d>=b) do begin d:=d-b;inc(c.s) end }
while(compare(d,b)>=0) do
begin
Subtract(d,b,d);
inc(c.s);
end;
end;
while(len>1)and(c.s[len]=0) do dec(len);
c.len:=len;
end;
procedure main;
begin
Divide(x[1],x[2],y,w);
end;
procedure out;
begin
assign(output,fn_out);
rewrite(output);
PrintHP(y);
writeln;
PrintHP(w);
writeln;
close(output);
end;
begin
init;
main;
out;
end.
高精度阶乘
作为一种高精度乘法的扩展算法,实质为高精度乘低精度,算法如下: var
a:array[1..10000] of longint;
i,j,k,l,p,o,q,x,y,w:integer;
begin
read(i);
a[1]:=1;
w:=1;
for j:=1 to i do
begin
y:=0; //到“For”前可省,但改为for k:=1 to 10000 do
x:=j;
while x>0 do
begin
y:=y+1;
x:=x div 10;
end;
o:=0;
for k:=w to l+y+1 do
begin
q:=a[k]*j+o;
o:=q div 10;
a[k]:=q mod 10;
end;
l:=10000;
while (a[l]=0) and (l>1) do l:=l-1;
w:=1;
while (a[w]=0) and (w<9999) do w:=w+1;
end;
for p:=l downto 1 do
write(a[p]);
writeln;
end.