博客园  :: 首页  :: 新随笔  :: 联系 :: 订阅 订阅  :: 管理

高精度运算

Posted on 2010-10-20 11:30  桃子在路上  阅读(507)  评论(0编辑  收藏  举报

由于待处理的数据超过了任何一种数据类型所能容纳的范围,因此必须采用数串形式输入,并将其转化为数组。该数组的每一个元素对应一个十进制数,由其下标顺序指明位序号。由于高精度运算可能使得数据长度发生变化,因此除要用整数数组存储数据外,还需要一个整数变量纪录整数数组的元素个数,即数据的实际长度。

type
  numtype=array[1..255] of byte;
var
  a:numtype;
  la:byte; 
  s:string;
begin
  readln(s);
  la:=length(s);
  for i:=1 to la do
    a[la-i+1]:=ord(s[i])-ord('0');
end.

 

 高精度加法运算

    首先,确定a和b中的最大位数x,然后依照由低位至高位的顺序进行加法运算。在每一次运算中,a当前位加b当前位的和除以10,其整商即为进位,其余数即为和的当前位。在进行了x位的加法后,若最高位有进位(a[x+1]<>0),则a的长度为x+1。
以下只列出关键程序:

type
  numtype=array[1..255] of word;
var
  a,b,s:numtype;
  la,lb,ls:word; 
procedure plus(var a:numtype;var la:word;b:numtype;lb:word);  {利用过程实现}
  var
    i,x:word;
  begin
    if la>=lb
       then x:=la
       else x:=lb;
    for i:=1 to x do
      begin
        a[i]:=a[i]+b[i];
        a[i+1]:=a[i+1]+a[i] div 10;
        a[i]:=a[i] mod 10;
      end;
   while a[x+1]<>0 do
      x:=x+1;   
   la:=x;             {最高位若有进位,则长度增加}
  end;

 

 高精度减法运算(a>b)

    依照由低位至高位的顺序进行减法运算。在每一次位运算中,若出现不够减的情况,则向高位借位。在进行了la位的减法后,若最高位为零,则a的长度减1。
以下只列出关键程序:

type
  numtype=array[1..255] of longint;
var
  a,b:numtype;
  la,lb:word; 
procedure minus(var a:numtype;var la:word;b:numtype;);
  var
    i:word;
  begin
    for i:=1 to la do
      begin
        if a[i]<b[i]
           then begin
                  dec(a[i+1]);
                  a[i]:=a[i]+10;
                end;
        a[i]:=a[i]-b[i];
      end;
    while (a[la]=0) and (la>1) do
       dec(la);
  end;

 

   高精度乘法运算

    按照乘法规则,从a的第1位开始逐位与c(c为字节型)相乘。在第i位乘法运算中,a的i位与c的乘积必须加上i-1位的进位,然后规整积的i-1位。
以下只列出关键程序:

type
  numtype=array[1..1000] of word;
var
  a:numtype;
  la:word; 
procedure multiply(var a:numtype;var la:word;c:word);
  var
    i:word;
  begin
    a[1]:=a[1]*c;
    for i:=2 to la do
      begin
        a[i]:=a[i]*c;
        a[i]:=a[i]+a[i-1] div 10;
        a[i-1]:=a[i-1] mod 10;
      end;
    while a[la]>=10 do
      begin
       inc(la);
       a[la]:=a[la-1] div 10;
       a[la-1]:=a[la-1] mod 10;
      end;
  end;

 

  扩大进制数改善高精度运算效率

    用整数数组每一个元素表示一个十进制整数的方法存在的缺点是:如果十进制的位数很多,则对应的数组的长度会很长,并增加了高精度计算的时间。
    如果用一个元素记录2位数字、3位数字或更多位数字,则数组的长度可以明显缩小,但是还要考虑数的取值范围问题,必须保证程序运行过程中不越界。在权衡了两方面的情况后得出:用一个longint记录4位数字是最佳的方案。那么这个数组就相当于一个10000进制的数,其中每一个元素都是10000进制下的一位数。

一、数据类型定义:

type
  numtype=array[1..10000] of longint;  {可以存储40000位十进制数}
var
  a,n:numtype;
  la,ln:byte; 
  s:ansistring;       {任意长度的字符串类型}

二、整数数组的建立和输出

readln(s);
k:=length(s);
for i:=1 to k do
begin
j:=(k-i+4) div 4;
n[j]:=n[j]*10+ord(s[i])-48;
end;
ln:=(k+3) div 4;

    当得出最后结果n后,必须按照次高位到最低位的顺序,将每一位元素由10000进制数转换成十进制数,即必须保证每个元素对应4位十进制数。例如n[i]=0015(0<=i<=ln-2),对应的十进制数不能为15,否则会导致错误结果。可以按照如下方法输出n对应的十进制数:
write(n[ln]);
for i:=ln-1 downto 1 do
  write(n[i] div 1000,(n[i] div 100) mod 10,(n[i] div 10) mod 10,n[i] mod 10);

三、基本运算
    两个10000进制整数的加法和减法与前面的十进制运算方法类似,只是进制变成了10000进制。
    1、整数数组减1(n:=n-1,n为整数数组)
    从n[1]出发寻找第一个非零的元素,由于接受了低位的借位,因此减1,其后缀全为9999。如果最高位为0,则n的长度减1。

    j:=1;
    while (n[j]=0) do inc(j);             {寻找第一个非零的元素}
    dec(n[j]);                            {该位接受低位的借位,因此减1}
    for i:=1 to j-1 do
       n[i]:=9999;                        {其后缀全为9999}
    if (j=ln) and (n[j]=0)                {如果最高位为0,则n的长度减1}
       then dec(ln);

    2、整数数组除以整数(a:=a div i,a为整数数组,i为整数)
    按照从高位到低位的顺序,逐位相除,把余数乘进制后加到下一位继续相除。如果最高位为0,则a的长度减1。

    l:=0;
    for j:=la downto 1 do
      begin
        inc(a[j],l*10000);
        l:=a[j] mod i;
        a[j]:=a[j] div i;
      end;
    while a[la]=0 do dec(la);

    3、两个整数数组相乘(a:=a*n,a和n为整数数组)
    按照从高位到低位的顺序,将数组a的每一个元素与n相乘。

procedure multiply(a,b:numtype;la,lb:longint;var s:numtype;var ls:longint);
  var
    i,j:longint;
  begin
    for i:=1 to la do
      for j:=1 to lb do
          s[i+j-1]:=s[i+j-1]+a[i]*b[j];
    for i:=1 to la+lb-1 do
      begin
        s[i+1]:=s[i+1]+s[i] div 10000;
        s[i]:=s[i] mod 10000;
      end;
    if s[la+lb]=0
       then ls:=la+lb-1
       else ls:=la+lb;
  end;

 

习题与练习

一、用高精度计算出s=1!+2!+3!+...+100!。
    参考答案

program jiecheng;
type
  numtype=array[1..255] of longint;
var
  s,t:numtype;
  ls,lt,i:longint;
procedure plus(var a:numtype;var la:longint;b:numtype;lb:longint);
  var
    i,x:byte;
  begin
    if la>=lb
       then x:=la
       else x:=lb;
    for i:=1 to x do
      begin
        a[i]:=a[i]+b[i];
        a[i+1]:=a[i+1]+a[i] div 10000;
        a[i]:=a[i] mod 10000;
      end;
   while a[x+1]<>0 do
      x:=x+1;
   la:=x;
  end;
procedure multiply(var a:numtype;var la:longint;c:longint);
  var
    i:longint;
  begin
    a[1]:=a[1]*c;
    for i:=2 to la do
      begin
        a[i]:=a[i]*c;
        a[i]:=a[i]+a[i-1] div 10000;
        a[i-1]:=a[i-1] mod 10000;
      end;
    while a[la]>=10000 do
      begin
       inc(la);
       a[la]:=a[la-1] div 10000;
       a[la-1]:=a[la-1] mod 10000;
      end;
  end;
begin
  lt:=1;
  t[1]:=1;
  ls:=0;
  for i:=1 to 100 do
    begin
     multiply(t,lt,i);
     plus(s,ls,t,lt);
    end;
write(s[ls]);
for i:=ls-1 downto 1 do
  write(s[i] div 1000,(s[i] div 100) mod 10,(s[i] div 10) mod 10,s[i] mod 10);
writeln;
writeln;
end. 

二、输入任意两个整数a,b(a,b均在长整型范围内),计算a/b的结果,保留100位有效数字,最后一位要求四舍五入。

三、两个高精度数相乘。
    参考答案

program multiplys;
type
  numtype=array[1..1000] of int64;
var
  a,b,s:numtype;
  la,lb,ls,i:longint;
procedure inputa(var n:numtype;var ln:longint);
 var
   s:ansistring;
   k,j:longint;
 begin
   readln(s);
   k:=length(s);
   for i:=1 to k do
    begin
     j:=(k-i+4) div 4;
     n[j]:=n[j]*10+ord(s[i])-48;
    end;
   ln:=(k+3) div 4;
 end;{input}
procedure outputa(n:numtype;ln:longint);
  begin
   write(n[ln]);
   for i:=ln-1 downto 1 do
     write(n[i] div 1000,(n[i] div 100) mod 10,(n[i] div 10) mod 10,n[i] mod 10);
   writeln;
  end;
procedure multiply(a,b:numtype;la,lb:longint;var s:numtype;var ls:longint);
  var
    i,j:longint;
  begin
    for i:=1 to la do
      for j:=1 to lb do
          s[i+j-1]:=s[i+j-1]+a[i]*b[j];
    for i:=1 to la+lb-1 do
      begin
        s[i+1]:=s[i+1]+s[i] div 10000;
        s[i]:=s[i] mod 10000;
      end;
    if s[la+lb]=0
       then ls:=la+lb-1
       else ls:=la+lb;
  end;
begin
  fillchar(a,sizeof(a),0);
  fillchar(b,sizeof(b),0);
  fillchar(s,sizeof(s),0);
  inputa(a,la);
  inputa(b,lb);
  multiply(a,b,la,lb,s,ls);
  writeln(ls);
  outputa(s,ls);
  readln;
end.

四、2k进制数(digital.pas)(NIOP2006第四题)
【问题描述】设r是个2k 进制数,并满足以下条件:
   (1)r至少是个2位的2k 进制数。
   (2)作为2k 进制数,除最后一位外,r的每一位严格小于它右边相邻的那一位。
   (3)将r转换为2进制数q后,则q的总位数不超过w。
    在这里,正整数k(1≤k≤9)和w(k<W≤30000)是事先给定的。
    问:满足上述条件的不同的r共有多少个?
    我们再从另一角度作些解释:设S是长度为w 的01字符串(即字符串S由w个“0”或“1”组成),S对应于上述条件(3)中的q。将S从右起划分为若干个长度为k 的段,每段对应一位2k进制的数,如果S至少可分成2段,则S所对应的二进制数又可以转换为上述的2k 进制数r。
    例:设k=3,w=7。则r是个八进制数(23=8)。由于w=7,长度为7的01字符串按3位一段分,可分为3段(即1,3,3,左边第一段只有一个二进制位),则满足条件的八进制数有:
    2位数:高位为1:6个(即12,13,14,15,16,17),高位为2:5个,…,高位为6:1个(即67)。共6+5+…+1=21个。
    3位数:高位只能是1,第2位为2:5个(即123,124,125,126,127),第2位为3:4个,…,第2位为6:1个(即167)。共5+4+…+1=15个。
    所以,满足要求的r共有36个。
【输入文件】
    输入文件digital.in只有1行,为两个正整数,用一个空格隔开:
 k W
【输出文件】
    输出文件digital.out为1行,是一个正整数,为所求的计算结果,即满足条件的不同的r的个数(用十进制数表示),要求最高位不得为0,各数字之间不得插入数字以外的其他字符(例如空格、换行符、逗号等)。
(提示:作为结果的正整数可能很大,但不会超过200位)
【输入样例】
 3 7
【输出样例】
 36

【题目分析】
    考虑一个首位为r、位数为n且除最后一位每一位严格小于右边的2k进制数,它的种数等于从2k-(r+1)个自然数中选出n-1个升序排列。
    而题目所求答案等于首位为1..2w mod k-1、位数为w/k+1且符合要求的2k进制数种数,加上首位任意、位数不超过w/k不低于2的且符合要求的2k进制数种数。(当w mod k=0时直接考虑后者即可)
    因为k,W是指定的,则可以确定2k进制数r的最长位数是w div k +1,设d=w div k,则d位2k进制数组成的2位以上d位以下的升序排列数应该是:

这里还要比较d与的2k-1大小,必须保证d<=2k-1
    最后考虑首位的情况,显然首位的最大二进制位数是w mod k,则最大值m=2w mod k  - 1,则首位不大于m,总位数是d+1位的升序排列数应是:
       
这里也要比较d与的2k-1大小,必须保证d<=2k-1-m
         两者相加即所求,因为最终的运算结果可能会很大,所以必须使用高精度运算。
参考程序】

program digital1;
type
  numtype=array[1..1000] of longint;
const
  p:array[1..9] of longint=(2,4,8,16,32,64,128,256,512);
var
 k,w,i,j,m,d,ls,lc:longint;
 c,s:numtype;
procedure plus(var a:numtype;var la:longint;b:numtype;lb:longint);
  var
    i,x:longint;
  begin
    if la>=lb
       then x:=la
       else x:=lb;
    for i:=1 to x do
      begin
        a[i]:=a[i]+b[i];
        a[i+1]:=a[i+1]+a[i] div 10000;
        a[i]:=a[i] mod 10000;
      end;
   while a[x+1]<>0 do
      x:=x+1;
   la:=x;
  end;
procedure multiply(var a:numtype;var la:longint;c:longint);
  var
    i:longint;
  begin
    a[1]:=a[1]*c;
    for i:=2 to la do
      begin
        a[i]:=a[i]*c;
        a[i]:=a[i]+a[i-1] div 10000;
        a[i-1]:=a[i-1] mod 10000;
      end;
    while a[la]>=10000 do
      begin
       inc(la);
       a[la]:=a[la-1] div 10000;
       a[la-1]:=a[la-1] mod 10000;
      end;
  end;
procedure divice(var a:numtype;var la:longint;c:longint);
var
  l,j:longint;
begin
    l:=0;
    for j:=la downto 1 do
      begin
        inc(a[j],l*10000);
        l:=a[j] mod c;
        a[j]:=a[j] div c;
      end;
    while a[la]=0 do dec(la);
end;
procedure zuhe(m,n:longint;var c:numtype;var lc:longint);
 var
  i:longint;
 begin
  fillchar(c, sizeof(c),0);
  c[1]:=1;
  lc:=1;
  for i:=n downto n-m+1 do
    multiply(c,lc,i);
  for i:=m downto 2 do
    divice(c,lc,i);
 end;
begin
  assign(input,'digital.in');
  reset(input);
  assign(output,'digital1.out');
  rewrite(output);
  readln(k,w);
  fillchar(s,sizeof(s),0);
  ls:=0;
  d:=w div k;
  if d>p[k]-1 then  d:=p[k]-1;
  for i:= 2 to d do
    begin
     zuhe(i,p[k]-1,c,lc);
     plus(s,ls,c,lc);
    end;
  if w mod k<>0 then
   begin
     m:=p[w mod k]-1;
     d:=w div k;
     if m>p[k]-1-d then m:=p[k]-1-d;
     for i:=1 to m do
      begin
       zuhe(d,p[k]-1-i,c,lc);
       plus(s,ls,c,lc);
      end;
   end;
write(s[ls]);
for i:=ls-1 downto 1 do
  write(s[i] div 1000,(s[i] div 100) mod 10,(s[i] div 10) mod 10,s[i] mod 10);
  close(input);
  close(output);
end.


 【NIOP满分程序】

const
inf='digital.in';
ouf='digital.out';
maxn=512;
pw:array[1..9]of longint=(2,4,8,16,32,64,128,256,512);
mol=1000000;
type
 bignum=array[0..50]of longint;
var
 k,w,n,i,p,j:longint;
 f:array[boolean,0..maxn]of bignum;
 ans:bignum;
procedure init;
begin
   assign(input,inf);
   reset(input);
   readln(k,w);
   n:=1;
   for i:=1 to k do
     n:=n*2;
   close(input);
end;
procedure hiadd(var a,b:bignum);
var
  i,j,n,tmp:longint;
  ta:bignum;
begin
   if a[0]<b[0] then a[0]:=b[0];
   for i:=1 to a[0] do
     a[i]:=a[i]+b[i];
   for i:=1 to a[0] do
     begin
       tmp:=a[i] div mol; inc(a[i+1],tmp);
       a[i]:=a[i] mod mol;
     end;
   if a[a[0]+1]>0 then inc(a[0]);
end;
procedure cal;
var
  now,last:boolean;
  change:boolean;
  max:longint;
begin
   p:=k;
   now:=true; last:=false;
   fillchar(f,sizeof(f),0);
   for i:=0 to n do
     begin
       f[now,i][0]:=1; f[now,i][1]:=1;
     end;
   dec(n);
   if w>k*n then w:=k*n;
   ans[0]:=1; ans[1]:=0;
   repeat
     now:=not now; last:=not last;
     fillchar(f[now],sizeof(f[now]),0);
     change:=false;
     inc(p,k);
     if p>w then
      begin
        dec(p,k); max:=w-p; p:=w;
      end
     else max:=k;
     hiadd(f[now,n-1],f[last,n]);
     for i:=n-2 downto 0 do
       begin
         f[now,i]:=f[now,i+1];
         hiadd(f[now,i],f[last,i+1]);
       end;
     if p>2*k then hiadd(ans,f[now,0]);
   until p=w;
   for i:=1 to pw[max]-1 do
     hiadd(ans,f[now,i]);
   assign(output,ouf);
   rewrite(output);
   n:=ans[0];
   write(ans[n]);
   for i:=n-1 downto 1 do
     begin
       max:=mol div 10;
       while max<>0 do
        begin
          write(ans[i] div max mod 10);
          max:=max div 10;
        end;
     end;
   writeln;
   close(output);
end;
begin
  init;
  cal;
end.