Текущий архив: 2010.02.21;
Скачать: CL | DM;

Вниз

Длинное деление: проблема с остатком   Найти похожие ветки 

← →
AssemblerIA64   (2009-03-25 23:02) [0]

Имеется код, выполняющий деление длинных чисел (взят из книги Уоррена):

...
Num = record
    x: Digit;
    l: Word;
end;
...
function LongDiv(u, v: Num; var q: Num; var r: Num): byte;
label again;
const
 b = 4294967296;
var
 s: Cardinal;
 un, vn: Num;
 i, j: Word;
 qhat: int64;
 rhat: int64;
 p, t: int64;
 temp1, temp2, k: uint64;
 owf: Boolean;

begin
 if (v.l <= 0) or (v.x[v.l-1] = 0) then
 begin
   LongDiv := 1;
   exit;
 end;
 if More(v, u) then
 begin
   SetLngth(q, 1);
   q.x[0] := 0;
   Move(r, u);
   LongDiv := 0;
   exit;
 end;
 SetLngth(q, u.l);
 for i := 0 to u.l - 1 do q.x[i] := 0;
 SetLngth(r, v.l);
 SetTo0(r, v.l);

 if v.l = 1 then
 begin
     {$Q-}
     {$R-}
     k := 0;
     for j := u.l-1 downto 0 do
     begin
       q.x[j] := (k*b+u.x[j]) div v.x[0];
       k := (k*b+u.x[j]) - q.x[j]*v.x[0];
     end;
     SetLngth(r, 1);
     r.x[0] := k;
     {$Q+}
     {$R+}
     LongDiv := 0;
     exit;
 end;
//--------------------------
 {$R-}
 s := nlz(v.x[v.l-1]);
 {$R+}
 SetLngth(vn, v.l);
 if s <> 0 then
 begin
   for j := v.l-1 downto 1 do
     vn.x[j] := (v.x[j] shl s) or (v.x[j-1] shr (32-s));
   vn.x[0] := v.x[0] shl s;
   SetLngth(un, u.l+1);
   un.x[u.l] := u.x[u.l-1] shr (32-s);
   for j := u.l-1 downto 1 do
     un.x[j] := (u.x[j] shl s) or (u.x[j-1] shr (32-s));
  un.x[0] := u.x[0] shl s;
 end
 else
 begin
   for j := v.l-1 downto 0 do
     vn.x[j] := v.x[j];
   SetLngth(un, u.l+1);
   for j := u.l-1 downto 0 do
     un.x[j] := u.x[j];
 end;
//--------------------------
 {$Q-}
 {$R-}
 for j := u.l-v.l downto 0 do
 begin
   qhat := (un.x[j+v.l]*b + un.x[j+v.l-1]) div vn.x[v.l-1];
   rhat := (un.x[j+v.l]*b + un.x[j+v.l-1]) - qhat*vn.x[v.l-1];
again:
   temp1 := qhat*vn.x[v.l-2];
   temp2 := b*rhat + un.x[j+v.l-2];
   if ((qhat >= b) or (temp1 > temp2)) then
   begin
     dec(qhat);
     rhat := rhat + vn.x[v.l-1];
     if (rhat < b) then goto again;
   end;
   k := 0;

   for i := 0 to v.l - 1 do
   begin
     p := qhat*vn.x[i];
     t := un.x[i+j] - k - (p and $0FFFFFFFF);
     un.x[i+j] := t;
     k :=  (p shr 32) - (t shr 32) ;
   end;

   t := un.x[j+v.l] - k;
   un.x[j+v.l] := t;
   q.x[j] := qhat;
   if t < 0 then
   begin
     q.x[j] := q.x[j] - 1;
     k := 0;
     for i := 0 to v.l-1 do
     begin
       t := un.x[i+j] + vn.x[i] + k;
       un.x[i+j] := t;
       k := t shr 32;
     end;
     un.x[j+v.l] := un.x[j+v.l] + k;
   end;
 end;
 {$R+}
 {$Q+}
 for j := 0 to v.l - 1 do
 begin
     r.x[j] := (un.x[j] shr s) or (un.x[j+1] shl (32-s));
 end;

 SetLngth(q, GetLngth(q.x, q.l));
 SetLngth(r, GetLngth(r.x, r.l));
 LongDiv := 0;
end;

В нём при делении (2851389442 + 2361951239*2^32 + 2247011709*2^64 + 50788274*2^96)/(4294967295 + 536870911*2^32) остаток определяется неправильно. Последнее значение un (это ещё не денормализованный остаток) равняется (3411107704 + 671223006*2^32 + 1*2^64), а должно быть (3411107704 + 671223006*2^32). Ошибка происходит на 3-ей итерации главного цикла в блоке
for i := 0 to v.l - 1 do
begin
     p := qhat*vn.x[i];
     t := un.x[i+j] - k - (p and $0FFFFFFFF);
     un.x[i+j] := t;
     k :=  (p shr 32) - (t shr 32) ;
end;
А в чём именно, я понять не могу. Подскажите, кто знает.
Спасибо.

+ вспомогательные функции

...
type
...
Digit = array of Cardinal;
...

function GetLngth(x : Digit; m: Word): Word;
var i : Integer;
begin
 i := m-1;
 while (i > 0) and (x[i] = 0) do
   i := i - 1;
 GetLngth := i+1;
end;

procedure SetLngth(var a: Num; m: Word);
begin
 SetLength(a.x, m);
 a.l := m;
end;

procedure Move(var a: Num; b: Num);
var
 i: Word;
begin
 SetLngth(a, b.l);
 a.l := b.l;
 for i := 0 to a.l - 1 do
   a.x[i] := b.x[i];
end;

procedure SetTo0(var a: Num; m: Word); overload;
var
 i: Word;
begin
 for i := 0 to a.l - 1 do a.x[i] := 0;
 SetLngth(a, m);
end;

function Eq(a, b: Num): Boolean;
var i: Word;
begin
   Eq:=False;
   if a.l <> b.l then exit
   else
   begin
       i := 0;
       while (i <= a.l) and (a.x[i] = b.x[i]) do
         Inc(i);
       Eq := (i = a.l) ;
   end;
end;

function More(a, b: Num): Boolean;
var i:Integer;
begin
 if a.l < b.l then
   More := false
 else if a.l > b.l then
   More := true
 else
 begin
   i := a.l - 1 ;
   while (i > 1) and (a.x[i] = b.x[i]) do dec(i);
   if a.x[i] > b.x[i] then
       More := true
   else
       More := false;
 end;
end;

+ upd
Нашёл одну ошибку: при вычислении qhat возникало переполнение. Надо писать
qhat := un.x[j+v.l];
qhat := (qhat*b + un.x[j+v.l-1]) div vn.x[v.l-1];
rhat := un.x[j+v.l];
rhat := (rhat*b + un.x[j+v.l-1]) - qhat*vn.x[v.l-1];

---

← →
Slym ©   (2009-03-26 05:00) [1]

Чем готовый FGInt не угодил?

---

← →
AssemblerIA64   (2009-03-26 23:29) [2]

Всё, проблема исчерпала себя.
Ошибки, в основном, были в операторах присвоения, в которых в левой части стояла переменная типа int64, а в правой выражение, в котором были переменные типа Cardinal. Из-за этого не происходило приведения младшего типа данных к старшиму до вычисления.

---

← →
Германн ©   (2009-03-27 00:28) [3]

> а в правой выражение, в котором были переменные типа Cardinal

Уж сколько раз твердили миру!  Особенно АП.

---

Страницы: 1 вся ветка

Текущий архив: 2010.02.21;
Скачать: CL | DM;

Наверх

Память: 0.48 MB
Время: 0.014 c