[Delphi] Модуль Длинной Арифметики. , Готовые решения - выбераем лучшее!

[Navi1982]

В этой теме хотелось бы собрать готовые реализации и интересные идеи тех или иных алгоритмов, так или иначе связанных с вычислениями очень, очень-очень больших чисел!... Т.е. фактически было бы здорово собрать воедино комплекс решений в качестве одного модуля на языке Delphi (7).

Хотелось бы попросить участников данной темы придерживатся некоторых правил:
1. Свои реализации и их улучшения, по возможности, писать, изменять и дополнять в единном посте. Исключениями могут быть кардинально иные принцыпы реализации алгоритмов.
2. Было бы гораздо интересней, если перед каждым алгоритмом будет дано теоретическое объяснение, поясняющие принцып его работы.
3. Просьба тему не засорять коментариями, а лучше авторам постов создать отдельную тему где будет обсуждатся его идея и коды, а также история всех улучшений и модификаций. А затем выносить сюда последний наиболее шустрый вариант (или просто ссылку).

Ниже будет приведено оглавление лучших идей с именами их авторов и с ссылкой на их пост.

ОГЛАВЛЕНИЕ

1. THugeInt - конвертирование, сравнение, сложение, вычетание, умножение, деление
От Navi1982. Статус - незавершен, код - не документирован, скорость - частично оптимизирован.
Соавторы: sCreator, Bug Hunter.

Сообщение отредактировано: Navi1982 - 12.06.10, 19:54

[Navi1982]

Сдесь текущий статус моего модуля Длинной Арифметики. Конечно же мне помогали и подробную историю можно посмотреть сдесь или последнее сообщение.
Вот мои некоторые наработки:
1. Числа хранятся в виде динамических массивов в таком вот виде:

type
  THugeInt = record //HugeInt <=> HI
    flag: longint; //number (data) length and his sign
    data: array of longword; //every elem is a digit in range 0..4294967295
    end;

2. Перед реализацией основных функций, я задумался об удобном способе проверять результат. Следовательно пришлось начать с функций преобразования текстовой строки компонента TEdit в заданный тип и обратно.

{Форматирование числа в строке.
Сначала проверяет знак - если отсувствует, добавляет "+" по умолчанию.
Все остальные знаки считает цифрами, т.е принимает за "0".
Так-же удаляет все незначимые ноли в начале числа.}
Function FormatStr(s:string):String;
var
 num:String;
 i:Integer;
 rem0:boolean; //remove "0" while it's true
begin
 case s[1] of
 '-': begin num:='-'; i:=2; end;
 '+': begin num:='+'; i:=2; end;
 else begin num:='+'; i:=1; end;
 end;
 rem0:=true; //to remove not significant zeroes before number
 while i <= length(s) do
 begin
   case s[i] of
    '0': if rem0=false then num:=num+s[i];//to put or not to put? - that is question! :)
    '1'..'9': begin rem0:=false; num:=num+s[i]; end;
   else
    if rem0=false then num:=num+'0';//put if significant!
   end;
   i:=i+1;
 end;
 if length(num)=1 then num:=num+'0';
FormatStr:=num;
end;

//Создана BugHunter'ом по мотивам sCreator'а
{Функция преобразования форматированой строки-10сс в массив-2^32сс}
{ Функция ожидает строку, содержащую корректную запись целого числа,  }
{ форматированую в соответствии со спецификацией: строка должна       }
{ содержать только цифры и возможно знак '+', '-' или один лидирующий }
{ пробел вместо знака '+'. Для входящей строки, не соответствующей    }
{ спецификации, результат работы функции непредсказуем.               }
Function StrToHugeInt(const s: string): THugeInt;
  Var
    i, rdlen, j, toend, k: integer;
    a: int64;
    alw: array[0..1] of longword absolute a;
    aa, inmind: longword;
  begin
    result.data:=nil;
    result.flag:=0;
 
    if (s[1] in ['+','-',' ']) then i:=2 else i:=1;
    toend:=length(s)-i+1;
    rdlen:=0;
 
    while (toend>0) do begin
      {йцукен фыва ячсми}
      if (toend>=9) then begin
        aa:=1000000000;
        k:=9;
        end
      else begin
        aa:=10; for k:=2 to toend do aa:=aa*10;
        k:=toend;
        end;
      {умножаем result.data на aa}
      inmind:=0;
      for j:=0 to rdlen-1 do begin
        a:=int64(result.data[j])*aa+inmind;
        result.data[j]:=alw[0];
        inmind:=alw[1];
        end;
      if (inmind>0) then begin
        setlength(result.data, rdlen+1); inc(rdlen);
        result.data[rdlen-1]:=inmind;
        end;
      {йцукен фыва ячсми}
      aa:=strtoint(copy(s, i, k));
      {складываем result.data и aa}
      inmind:=aa;
      for j:=0 to rdlen-1 do begin
        a:=int64(result.data[j])+inmind;
        result.data[j]:=alw[0];
        inmind:=alw[1];
        if (inmind=0) then break;
        end;
      if (inmind>0) then begin
        setlength(result.data, rdlen+1); inc(rdlen);
        result.data[rdlen-1]:=inmind;
        end;
      {йцукен фыва ячсми}
      inc(i, k); dec(toend, k);
      if (toend=0) then break;
      end;
 
    if (s[1]='-') then
      result.flag:=-rdlen
    else
      result.flag:=rdlen;
 
    //я тут добавил для обработки нолей.
    if (result.data=nil)and(result.flag=0) then
    begin
     SetLength(result.data,1);
     result.data[0]:=0;
     result.flag:=1;
    end;
  end;

{Функция преобразования числа HugeInt -> HexStr}
Function HugeIntToHexStr(Const op:THugeInt):String;
Var i:integer;
    t:String;
Begin
 if op.flag<0 then t:='-' else t:='+';
 if (op.flag<>0)or(op.data<>nil) then
 for i:=abs(op.flag)-1 downto 0 do t:=t+IntToHex(op.data[i],8);
result:=t;
End;

3. Далее идут вспомогательные функции и процедуры для работы основных.

{Процедура коррекции числа:
 убирает незначемые нулевые элементы с конца .data}
Procedure HICorrection(Var op:THugeInt);
Var i,sg:integer;
Begin
 if op.data <> nil then //or op.flag<>0
 begin
   if op.flag<0 then sg:=-1 else sg:=1;
   i:=length(op.data)-1;
   while (i>0)and(op.data[i]=0) do Dec(i); //пока не первый и равен нулю
   SetLength(op.data,i+1);
   op.flag:=length(op.data)*sg;
 end;
End;

{Функция сдвига разрядов в старшую сторону}
Function HIShl(const op:THugeInt; n:integer):THugeInt;
Var i,sg:integer;
begin
  if n<0 then n:=0;
  if op.flag<0 then sg:=-1 else sg:=1;
  result.flag:=abs(op.flag)+n*sg;
  SetLength(result.data,abs(result.flag));
  for i:=0 to abs(op.flag)-1 do result.data[i+n]:=op.data[i];
end;

4. Функции сравнения - так-же служат вспомогательными.

{функция сравнения по модулю, нужна для некоторых операций}
{возвращает: 1 if |op1|>|op2|; 0 if |op1|=|op2|; -1 if |op1|<|op2|}
{!!!доработать если операнды=nil}
Function HIModCompare(Const op1, op2:THugeInt):integer;
Var i,r:integer;
Begin
  if abs(op1.flag)>abs(op2.flag) then r:=1;
  if abs(op1.flag)<abs(op2.flag) then r:=-1;
  if abs(op1.flag)=abs(op2.flag) then
  begin
    i:=abs(op1.flag)-1;
    while (i>=0)and(op1.data[i]=op2.data[i]) do Dec(i);
    if i<0 then r:=0 else
    begin
     if op1.data[i]>op2.data[i] then r:=1;
     if op1.data[i]<op2.data[i] then r:=-1;
    end;
  end;
result:=r;
End;

{функция сравнения с учетом знака}
{возвращает: 1 if op1>op2; 0 if op1=op2; -1 if op1<op2}
{!!!доработать если операнды=nil}
Function HICompare(Const op1, op2:THugeInt):integer;
Var i,r:integer;
Begin
  if op1.flag>op2.flag then r:=1;
  if op1.flag<op2.flag then r:=-1;
  if op1.flag=op2.flag then
  begin
    i:=abs(op1.flag)-1;
    while (i>=0)and(op1.data[i]=op2.data[i]) do Dec(i);
    if i<0 then r:=0 else
    begin
     if op1.data[i]>op2.data[i] then r:=1;
     if op1.data[i]<op2.data[i] then r:=-1;
     if op1.flag<0 then r:=r*-1;
    end;
  end;
result:=r;
End;

5. И основные арифметические функции:
Сложение и вычитание:

{Функция сложения по модулю}
Function HIModAdd(const op1,op2:THugeInt):THugeInt;
Var
  i,len,lenmax:integer;
  t:Int64;
  cr:array [0..1] of longword absolute t;
begin
  //возьмем число за основу
  if abs(op1.flag)<abs(op2.flag) then
  begin
    len:=abs(op1.flag); lenmax:=abs(op2.flag);
    SetLength(result.data,lenmax); //определим размер результата
    //начинаем сложение с младших элементов до минимальной длины
    t:=0; //<=> cr[0..1]:=0;
    for i:=0 to len-1 do begin
      t:=int64(op1.data[i])+int64(op2.data[i])+int64(cr[1]);
      result.data[i]:=cr[0];
    end;
    //далее прибавляем 0 с переносами до конца
    for i:=len to lenmax-1 do begin
      t:=int64(op2.data[i])+int64(cr[1]);
      result.data[i]:=cr[0];
    end;
  end else //abs(op1.flag)>=abs(op2.flag)
  begin
    len:=abs(op2.flag); lenmax:=abs(op1.flag);
    SetLength(result.data,lenmax); //определим размер результата
    //начинаем сложение с младших элементов до минимальной длины
    t:=0; //<=> cr[0..1]:=0;
    for i:=0 to len-1 do begin
      t:=int64(op1.data[i])+int64(op2.data[i])+int64(cr[1]);
      result.data[i]:=cr[0];
    end;
    //далее прибавляем 0 с переносами до конца
    for i:=len to lenmax-1 do begin
      t:=int64(op1.data[i])+int64(cr[1]);
      result.data[i]:=cr[0];
    end;
  end;
  //учтем значение последнего переноса
  if t>$FFFFFFFF then begin
    lenmax:=lenmax+1;
    SetLength(result.data,lenmax);
    result.data[lenmax-1]:=cr[1];
  end;
  result.flag:=lenmax; //длина результата
end;

{Функция вычетания по модулю}
{result = op1 - op2; only if |op1| >= |op2|}
Function HIModSub(const op1,op2:THugeInt):THugeInt;
Var
  len,rif:integer;
  P2:THugeInt;
begin
  result.data:=nil;
  result.flag:=0;
  //проверка, если |op1| >= |op2|
  rif:=HIModCompare(op1,op2);
  if rif>=0 then
  begin
   len:=length(op1.data);
   P2.data:=Copy(op2.data);
   SetLength(P2.data,len); P2.flag:=len;
//начинаем вычитане
asm
        PUSH    EAX
        PUSH    EDI
        PUSH    ESI
        PUSH    ECX
        MOV     EDI,P2.data
        MOV     ESI,op1 // ESI <- указатель на op1 (thx to leo!)
        MOV     ESI,[ESI+THugeInt.data] // ESI <- указатель на op1.data
        MOV     ECX,len
        CLD
        CLC
@@1:    LODSD
        SBB     EAX,[EDI]
        STOSD
        LOOP    @@1
        POP     ECX
        POP     ESI
        POP     EDI
        POP     EAX
end;
HICorrection(P2);
result.data:=Copy(P2.data);
result.flag:=P2.flag;
end;
end;

Теперь с учетом положительных и отрицательных знаков:

{Функция сложения со знаком}
Function HIAdd(const op1,op2:THugeInt):THugeInt;
Var cmp:integer;
begin
 cmp:=HIModCompare(op1,op2);
 if cmp<0 then
 begin
  if (op1.flag>0)and(op2.flag>0) then
   result:=HIModAdd(op2,op1);
  if (op1.flag>0)and(op2.flag<0) then
  begin
   result:=HIModSub(op2,op1);
   result.flag:=result.flag*(-1);
  end;
  if (op1.flag<0)and(op2.flag>0) then
   result:=HIModSub(op2,op1);
  if (op1.flag<0)and(op2.flag<0) then
  begin
   result:=HIModAdd(op2,op1);
   result.flag:=result.flag*(-1);
  end;
 end else
 begin
  if (op1.flag>0)and(op2.flag>0) then
   result:=HIModAdd(op1,op2);
  if (op1.flag>0)and(op2.flag<0) then
   result:=HIModSub(op1,op2);
  if (op1.flag<0)and(op2.flag>0) then
  begin
   result:=HIModSub(op1,op2);
   result.flag:=result.flag*(-1);
  end;
  if (op1.flag<0)and(op2.flag<0) then
  begin
   result:=HIModAdd(op1,op2);
   result.flag:=result.flag*(-1);
  end;
 end;
end;

{Функция вычитания со знаком}
Function HISub(const op1,op2:THugeInt):THugeInt;
Var cmp:integer;
begin
 cmp:=HIModCompare(op1,op2);
 if cmp<0 then
 begin
  if (op1.flag>0)and(op2.flag>0) then
  begin
   result:=HIModSub(op2,op1);
   result.flag:=result.flag*(-1);
  end;
  if (op1.flag>0)and(op2.flag<0) then
   result:=HIModAdd(op2,op1);
  if (op1.flag<0)and(op2.flag>0) then
  begin
   result:=HIModAdd(op2,op1);
   result.flag:=result.flag*(-1);
  end;
  if (op1.flag<0)and(op2.flag<0) then
   result:=HIModSub(op2,op1);
 end else
 begin
  if (op1.flag>0)and(op2.flag>0) then
   result:=HIModSub(op1,op2);
  if (op1.flag>0)and(op2.flag<0) then
   result:=HIModAdd(op1,op2);
  if (op1.flag<0)and(op2.flag>0) then
  begin
   result:=HIModAdd(op1,op2);
   result.flag:=result.flag*(-1);
  end;
  if (op1.flag<0)and(op2.flag<0) then
  begin
   result:=HIModSub(op1,op2);
   result.flag:=result.flag*(-1);
  end;
 end;
end;

Умножение и деление:

{Функция умножения по модулю}
Function HIModMul(const op1,op2:THugeInt):THugeInt;
Var
  i,j,len1,len2:integer;
  t:Int64;
  cr:array [0..1] of longword absolute t;
  rp,r:THugeInt; //rp - is precedent result of multiply
begin
  len1:=abs(op1.flag);
  len2:=abs(op2.flag);
  //begin rp:=0;
  rp.flag:=1; SetLength(rp.data,rp.flag);
  rp.data[rp.flag-1]:=0; //end rp:=0;
  j:=0;
  while j<len2 do
  begin
    i:=0; t:=0;
    //begin r:=0;
    SetLength(r.data,0);
    SetLength(r.data,len1+j+1); //size of op1.data + shift + 1 elem of .data
    //end r:=0;
    r.flag:=0+j;
    while i<len1 do
    begin
      t:=int64(op2.data[j])*int64(op1.data[i])+t;
      inc(r.flag); //SetLength(r.data,r.flag);
      r.data[r.flag-1]:=cr[0];
      t:=int64(cr[1]);
      inc(i);
    end;
      inc(r.flag);
      r.data[r.flag-1]:=t;
    //Add r with rp
    r:=HIModAdd(r,rp);
    //begin rp:=r;
    rp.flag:=r.flag; SetLength(rp.data,rp.flag);
    rp.data:=Copy(r.data); //end rp:=r;
    inc(j);
  end;
  HICorrection(rp); // убирает незначемые нулевые элементы из массива числа
  result:=rp;
end;

{Функция деления по модулю без остатка op1 DIV op2}
{разряд результата угадывается по методу дихотомии}
Function HIModDiv(const op1,op2:THugeInt):THugeInt;
Var
  i,len1:integer;
  r,t,y,rHI:THugeInt;
  rk:longword;
  base32,r64,x1,x2:Int64;
begin
 //учтем деление на ноль!
  t.flag:=1; SetLength(t.data,t.flag); // t=0;
  t.data[t.flag-1]:=0;
 if HICompare(t,op2)<>0 then //op2<>0 ? //иначе result=nil <=> бесконечность
 begin
  base32:=$100000000; //2^32
  len1:=abs(op1.flag); //"лишнее наследство" //вобщем: проход начнем со старших разрядов делимого,
  i:=len1;            //добавляя в промежуточное делимое по одному разряду справа...
  r.flag:=0; SetLength(r.data,r.flag); // r=nil обнуляем часное
  t.flag:=0; SetLength(t.data,t.flag); // t=nil обнуляем промежуточное делимое
  while i>0 do        //...пока не дойдем до последнего разряда из делимого!
  begin
    t:=HIShl(t,1);    //Т.е. в промежуточном делимом сдвигаем разряды влево (умножаем на базу 2^32),
    HICorrection(t);  //убеждаемся что после такой процедуры неокажется "левых" нолей,
    t.data[0]:=op1.data[i-1]; //и "прибавляем" разряд справа от делимого.
    //gues r[k] - начинаем угадывание разряда в часном (макс.подх.знач.) по дихотомии...
    x1:=0; x2:=base32;                            //...в отрезке 0..2^32
    r64:=(x1+x2) div 2;                           //делим отрезок на две половины
    rk:=Int64Rec(r64).Lo;       //rk и есть текущий разряд в часном и = только что вычесленной середине отрезка.
    rHI.flag:=1; SetLength(rHI.data,rHI.flag); //временная переменная THugeInt для rk
    while x1<>rk do   //будим разделять отрезки пока x1<>rk
    begin
     rHI.data[0]:=rk;
     y:=HIModMul(rHI,op2);//rk*делитель
                          //теперь проверим в каком отрезке находится
                          //максимально подходящее значение для rk (текущего разряда часного)
     if HICompare(t,y)>=0 then x1:=rk else x2:=rk; //If t>=y Then [значение в старшем отр.] else [младшем]
     r64:=(x1+x2) div 2;                           //и выбраный отрезок снова на пополам
     rk:=Int64Rec(r64).Lo;                                   //rk=середина отрезка.
    end;
    rHI.data[0]:=rk;
    y:=HIModMul(rHI,op2); // y=[текущий разряд часного]*[делитель]
    t:=HIModSub(t,y);     // [промежуточное делимое]=[промежуточное делимое]-y
    HICorrection(t);      // убираем "левые" ноли, если имеются.
    r:=HIShl(r,1);        // фиксируем часное умножая его на базу
    r.data[0]:=rk;        //  и "прибавим" только что вычеслиный разряд
    dec(i);               // опускаем следующий разряд делимого справа к промежуточному
  end;
  HICorrection(r);  // в результате могут оказаться ноли и мы их почистим кроме значемого
 end else //if op2=0 => r=nil ! or infinite (см. в начало)
 begin
  r.flag:=0; SetLength(r.data,r.flag); // r=nil
 end;
  result:=r;
end;

Теперь со знаком:

{Функция умножения со знаком}
Function HIMul(const op1,op2:THugeInt):THugeInt;
Var f1,f2,m:integer;
begin
  f1:=op1.flag;
  f2:=op2.flag;
  m:=0;
  if (f1>0)and(f2>0) then m:=1 else
   if (f1<0)and(f2<0) then m:=1 else
    if (f1>0)and(f2<0) then m:=-1 else
     if (f1<0)and(f2>0) then m:=-1;
  result:=HIModMul(op1,op2);
  result.flag:=result.flag*m;
end;

{Функция деления op1 DIV op2 со знаком}
Function HIDiv(const op1,op2:THugeInt):THugeInt;
Var f1,f2,m:integer;
begin
  f1:=op1.flag;
  f2:=op2.flag;
  m:=0;
  if (f1>0)and(f2>0) then m:=1 else
   if (f1<0)and(f2<0) then m:=1 else
    if (f1>0)and(f2<0) then m:=-1 else
     if (f1<0)and(f2>0) then m:=-1;
  result:=HIModDiv(op1,op2);
  result.flag:=result.flag*m;
end;

Сообщение отредактировано: Navi1982 - 07.07.10, 19:50