Длинная арифметика, числа > longint Опции

[virt]

иногда нужно вычислять очень длинные числа ,например 14! = 87.178.291.200
в стандартные типы данных не лезет. Поэтому есть длинная арифметика.

Код

const _maxdig=1000;{максимальное количество 4х значных цифр}
     _osn=10000;{основание системы счисления}

type Tlong=array[0.._maxdig]of integer;{здесь храним само число}
    Plong=^Tlong;

  например число 1234567890 запишется так :
  | 0    | 1    | 2    | 3    | 4    | ...      -- индексы элементов массива
  ----------------------------------------------
  | 3    | 7890 | 3456 | 12   | 0    | ...      -- значения элементов
  ----------------------------------------------
     нулевой элемент массива -- количество 4х значных цифр.

чтение и запись длинных чисел :

procedure ReadLong(var f:text;a:Plong);
var ch:char;
    i:integer;
begin
   fillchar(a^,sizeof(a^),0);
   read(f,ch);
   while not (ch in ['0'..'9',#26]) do read(f,ch);
   while ch in ['0'..'9'] do
      begin
         for i:=a^[0] downto 1 do
            begin
               a^[i+1]:=a^[i+1]+(longint(a^[i])*10)div _osn;
               a^[i]:=(longint(a^[i])*10)mod _osn;
            end;
         a^[1]:=a^[1]+ord(ch)-ord('0');
         if a^[a^[0]+1]>0 then inc(a^[0]);
         read(f,ch);
      end;
end;


procedure WriteLong(var f:text;a:Plong);
var ls,s:string;
    i:integer;
begin
   str(_osn div 10,ls);
   write(f,a^[a^[0]]);
   for i:=a^[0]-1 downto 1 do
      begin
         str(a^[i],s);
         while length(s)<length(ls) do s:='0'+s;
         write(f,s);
      end;
   writeln(f);
end;

сложение и вычитание двух длинных чисел :

procedure SumLongTwo(a,b,c:Plong);
var i,k:integer;
begin
   fillchar(c^,sizeof(c^),0);
   if a^[0]>b^[0] then k:=a^[0] else k:=b^[0];
   for i:=1 to k do
      begin
         c^[i+1]:=(c^[i]+a^[i]+b^[i]) div _osn;
         c^[i]:=(c^[i]+a^[i]+b^[i]) mod _osn;
      end;
   if c^[k+1]=0 then c^[0]:=k else c^[0]:=k+1;
end;

procedure SubLongTwo(a,b:Plong;const sdvig:integer);
var i,j:integer;
begin
   for i:=1 to b^[0] do
      begin
         dec(a^[i+sdvig],b^[i]);
         j:=i;
         while (a^[j+sdvig]<0) and (j<=a^[0]) do
            begin
               inc(a^[j+sdvig],_osn);
               dec(a^[j+sdvig+1]);
               inc(j);
            end;
      end;
   i:=a^[0];
   while (i>1) and (a^[i]=0) do dec(i);
   a^[0]:=i;
end;
{a>=b*(_osn^sdvig); a<-- a-b*(_osn^sdvig) }

сравнение длинных чисел :

function EqLong(a,b:Plong):boolean;{a=b}
var i:integer;
begin
   EqLong:=false;
   if a^[0]=b^[0] then
      begin
         i:=1;
         while (i<=a^[0]) and (a^[i]=b^[i]) do inc(i);
         EqLong:=i=a^[0]+1;
      end;
end;



function MoreLong(a,b:Plong):boolean;{a>b}
var i:integer;
begin
   if a^[0]<b^[0] then MoreLong:=false else
      if a^[0]>b^[0] then MoreLong:=true else
         begin
            i:=a^[0];
            while (i>0) and (a^[i]=b^[i]) do dec(i);
            if i=0 then MoreLong:=false else
               if a^[i]>b^[i] then MoreLong:=true else MoreLong:=false;
         end;
end;


function LessLong(a,b:Plong):boolean;{a<b}
begin
   LessLong:=not(MoreLong(a,b) or EqLong(a,b));
end;


function More_EqLong(a,b:Plong):boolean;{a>=b}
begin
   More_EqLong:=MoreLong(a,b) or EqLong(a,b);
end;


function Less_EqLong(a,b:Plong):boolean;{a<=b}
begin
   Less_EqLong:=not(MoreLong(a,b));
end;


function MoreSdvigLong(a,b:Plong;const sdvig:integer):byte;{a>b*(_osn^sdvig) -- 0
                                                             a<b*(_osn^sdvig) -- 1
                                                             a=b*(_osn^sdvig) -- 2}
var i:integer; 
begin
   if a^[0]>(b^[0]+sdvig) then MoreSdvigLong:=0 else
      if a^[0]<(b^[0]+sdvig) then MoreSdvigLong:=1 else
         begin
            i:=a^[0];
            while (i>sdvig) and (a^[i]=b^[i-sdvig]) do dec(i);
            if i=sdvig then
               begin
                  MoreSdvigLong:=0;
                  for i:=1 to sdvig do
                     if a^[i]>0 then exit;
                  MoreSdvigLong:=2;
               end else
                  MoreSdvigLong:=byte(a^[i]<b^[i-sdvig]);{0 -- false(a>b);1 -- true(a<b)}
         end;
end;

умножение длинного числа на короткое :

procedure MulLongShort(a:Plong;const k:longint;c:Plong);
var i:integer;
begin
   fillchar(c^,sizeof(c^),0);
   if k=0 then inc(c^[0]) else
      begin
         for i:=1 to a^[0] do
            begin
               c^[i+1]:=(longint(a^[i])*k+c^[i]) div _osn;
               c^[i]:=(longint(a^[i])*k+c^[i]) mod _osn;
            end;
         if c^[a^[0]+1]>0 then c^[0]:=a^[0]+1 else c^[0]:=a^[0];
      end;
end;

(*
  при 0<=k<=_osn гарантированно работает правильно ,
  при _osn+1<=k<=67479 вроде тоже правильные ответы выдает (проверьте кто нить при всех ли)
  при k>=67480 во время исполнения проги возникает range check.
*)

умножение двух длинных чисел :

procedure MulLongTwo(a,b,c:Plong);
var i,j:integer;
    dv:longint;
begin
   fillchar(c^,sizeof(c^),0);
   for i:=1 to a^[0] do
      for j:=1 to b^[0] do
         begin
            dv:=longint(a^[i])*b^[j]+c^[i+j-1];
            c^[i+j]:=c^[i+j]+dv div _osn;
            c^[i+j-1]:=dv mod _osn;
         end;
   c^[0]:=a^[0]+b^[0];
   while (c^[0]>1) and (c^[c^[0]]=0) do dec(c^[0]);
end;

деление двух длинных чисел :

function FindBin(ost,b:Plong;const sp:integer):longint;
var up,down:word;
    c:Plong;
begin
   new©;
   down:=0;up:=_osn;
   while up-1>down do
      begin
         MulLongShort(b,(up+down) div 2,c);
         case MoreSdvigLong(ost,c,sp) of
          0:down:=(up+down) div 2;
          1:up:=(up+down) div 2;
          2:begin
               up:=(up+down) div 2;
               down:=up;
            end;
         end;
      end;
   MulLongShort(b,(up+down) div 2,c);
   if MoreSdvigLong(ost,c,0)=0 then SubLongTwo(ost,c,sp) else
      begin
         SubLongTwo(c,ost,sp);
         ost:=c;
      end;
   FindBin:=(up+down) div 2;
   dispose©;
end;


procedure MakeDel(a,b,res,ost:Plong);
var sp:integer;
begin
   ost^:=a^;
   sp:=a^[0]-b^[0];
   if MoreSdvigLong(a,b,sp)=1 then dec(sp);{!!!!!!!!!}
   res^[0]:=sp+1;
   while sp>=0 do
      begin
         res^[sp+1]:=FindBin(ost,b,sp);
         dec(sp);
      end;
end;


procedure DivLongTwo(a,b,res,ost:Plong);
begin
   fillchar(res^,sizeof(res^),0);res^[0]:=1;
   fillchar(ost^,sizeof(ost^),0);ost^[0]:=1;
   case MoreSdvigLong(a,b,0) of
    0:MakeDel(a,b,res,ost);
    1:ost^:=a^;
    2:res^[1]:=1;
   end;
end;

(*
  что бы разделить a на b вызовите : DivLongTwo(a,b,res,ost);
  res -- результат деления a на b; ost -- остаток от деления.
*)