Длинное деление: проблема с остатком

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

Код (Text):

1. ...
2. Num = record
3.      x: Digit;
4.      l: Word;
5. end;
6. ...
7. function LongDiv(u, v: Num; var q: Num; var r: Num): byte;
8. label again;
9. const
10.   b = 4294967296;
11. var
12.   s: Cardinal;
13.   un, vn: Num;
14.   i, j: Word;
15.   qhat: int64;
16.   rhat: int64;
17.   p, t: int64;
18.   temp1, temp2, k: uint64;
19.   owf: Boolean;
20.  
21. begin
22.   if (v.l <= 0) or (v.x[v.l-1] = 0) then
23.   begin
24.     LongDiv := 1;
25.     exit;
26.   end;
27.   if More(v, u) then
28.   begin
29.     SetLngth(q, 1);
30.     q.x[0] := 0;
31.     Move(r, u);
32.     LongDiv := 0;
33.     exit;
34.   end;
35.   SetLngth(q, u.l);
36.   for i := 0 to u.l - 1 do q.x[i] := 0;
37.   SetLngth(r, v.l);
38.   SetTo0(r, v.l);
39.  
40.   if v.l = 1 then
41.   begin
42.       {$Q-}
43.       {$R-}
44.       k := 0;
45.       for j := u.l-1 downto 0 do
46.       begin
47.         q.x[j] := (k*b+u.x[j]) div v.x[0];
48.         k := (k*b+u.x[j]) - q.x[j]*v.x[0];
49.       end;
50.       SetLngth(r, 1);
51.       r.x[0] := k;
52.       {$Q+}
53.       {$R+}
54.       LongDiv := 0;
55.       exit;
56.   end;
57. //--------------------------
58.   {$R-}
59.   s := nlz(v.x[v.l-1]);
60.   {$R+}
61.   SetLngth(vn, v.l);
62.   if s <> 0 then
63.   begin
64.     for j := v.l-1 downto 1 do
65.       vn.x[j] := (v.x[j] shl s) or (v.x[j-1] shr (32-s));
66.     vn.x[0] := v.x[0] shl s;
67.     SetLngth(un, u.l+1);
68.     un.x[u.l] := u.x[u.l-1] shr (32-s);
69.     for j := u.l-1 downto 1 do
70.       un.x[j] := (u.x[j] shl s) or (u.x[j-1] shr (32-s));
71.    un.x[0] := u.x[0] shl s;
72.   end
73.   else
74.   begin
75.     for j := v.l-1 downto 0 do
76.       vn.x[j] := v.x[j];
77.     SetLngth(un, u.l+1);
78.     for j := u.l-1 downto 0 do
79.       un.x[j] := u.x[j];
80.   end;
81. //--------------------------
82.   {$Q-}
83.   {$R-}
84.   for j := u.l-v.l downto 0 do
85.   begin
86.     qhat := (un.x[j+v.l]*b + un.x[j+v.l-1]) div vn.x[v.l-1];
87.     rhat := (un.x[j+v.l]*b + un.x[j+v.l-1]) - qhat*vn.x[v.l-1];
88. again:
89.     temp1 := qhat*vn.x[v.l-2];
90.     temp2 := b*rhat + un.x[j+v.l-2];
91.     if ((qhat >= b) or (temp1 > temp2)) then
92.     begin
93.       dec(qhat);
94.       rhat := rhat + vn.x[v.l-1];
95.       if (rhat < b) then goto again;
96.     end;
97.     k := 0;
98.  
99.     for i := 0 to v.l - 1 do
100.     begin
101.       p := qhat*vn.x[i];
102.       t := un.x[i+j] - k - (p and $0FFFFFFFF);
103.       un.x[i+j] := t;
104.       k :=  (p shr 32) - (t shr 32) ;
105.     end;
106.  
107.     t := un.x[j+v.l] - k;
108.     un.x[j+v.l] := t;
109.     q.x[j] := qhat;
110.     if t < 0 then
111.     begin
112.       q.x[j] := q.x[j] - 1;
113.       k := 0;
114.       for i := 0 to v.l-1 do
115.       begin
116.         t := un.x[i+j] + vn.x[i] + k;
117.         un.x[i+j] := t;
118.         k := t shr 32;
119.       end;
120.       un.x[j+v.l] := un.x[j+v.l] + k;
121.     end;
122.   end;
123.   {$R+}
124.   {$Q+}
125.   for j := 0 to v.l - 1 do
126.   begin
127.       r.x[j] := (un.x[j] shr s) or (un.x[j+1] shl (32-s));
128.   end;
129.  
130.   SetLngth(q, GetLngth(q.x, q.l));
131.   SetLngth(r, GetLngth(r.x, r.l));
132.   LongDiv := 0;
133. 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-ей итерации главного цикла в блоке

Код (Text):

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

А в чём именно, я понять не могу. Подскажите, кто знает.
Спасибо.

 

А почему не на С? Переходи, не пожалеешь.

 

Ой, совсем забыл про вспомогательные функции:

Код (Text):

1. procedure SetLngth(var a: Num; m: Word);
2. begin
3.   SetLength(a.x, m);
4.   a.l := m;
5. end;
6.  
7. procedure Move(var a: Num; b: Num);
8. var
9.   i: Word;
10. begin
11.   SetLngth(a, b.l);
12.   a.l := b.l;
13.   for i := 0 to a.l - 1 do
14.     a.x[i] := b.x[i];
15. end;
16.  
17. procedure SetTo0(var a: Num; m: Word); overload;
18. var
19.   i: Word;
20. begin
21.   for i := 0 to a.l - 1 do a.x[i] := 0;
22.   SetLngth(a, m);
23. end;
24.  
25. function Eq(a, b: Num): Boolean;
26. var i: Word;
27. begin
28.     Eq:=False;
29.     if a.l <> b.l then exit
30.     else
31.     begin
32.         i := 0;
33.         while (i <= a.l) and (a.x[i] = b.x[i]) do
34.           Inc(i);
35.         Eq := (i = a.l) ;
36.     end;
37. end;
38.  
39. function More(a, b: Num): Boolean;
40. var i:Integer;
41. begin
42.   if a.l < b.l then
43.     More := false
44.   else if a.l > b.l then
45.     More := true
46.   else
47.   begin
48.     i := a.l - 1 ;
49.     while (i > 1) and (a.x[i] = b.x[i]) do dec(i);
50.     if a.x[i] > b.x[i] then
51.         More := true
52.     else
53.         More := false;
54.   end;
55. end;

Вроде всё. Если чего не хватает, пишите.

 

Phuntik, да я изначально хотел на С, но пришлось на Delphi. Не по своей воле.

 

Ну что вы все молчите?

Нашёл одну ошибку: при вычислении qhat возникало переполнение. Надо писать

Код (Text):

1.  
2. qhat := un.x[j+v.l];
3. qhat := (qhat*b + un.x[j+v.l-1]) div vn.x[v.l-1];
4. rhat := un.x[j+v.l];
5. rhat := (rhat*b + un.x[j+v.l-1]) - qhat*vn.x[v.l-1];

А при вычитании un - qhat*vn так и не могу понять в чём дело, вроде переносы и заёмы правильно учитываются, а в старшем разряде всё равно единица.
Надеюсь на вашу помощь.

 

Ах да, забыл про GetLngth:

Код (Text):

1. ...
2. type
3. ...
4. Digit = array of Cardinal;
5. ...
6.  
7. function GetLngth(x : Digit; m: Word): Word;
8. var i : Integer;
9. begin
10.   i := m-1;
11.   while (i > 0) and (x[i] = 0) do
12.     i := i - 1;
13.   GetLngth := i+1;
14. end;

 

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