Умножая вычеты
В известном алгоритме замены деления умножением произведение имеет двойную длину. Существует ли вариант с "коротким" произведением?
Молчит наука
Мы рассмотрим алгоритм деления целых чисел на константу, основанный на вычислении "короткого" произведения делимого на число, обратное делителю. При этом частное и остаток вычисляются с использованием таблицы предварительно рассчитанных значений, соответствующих различным значениям старших битов произведения. Ниже приводится обоснование алгоритма, которое можно пропустить и сразу перейти к описанию алгоритма в последнем разделе.
Полезные факты
Используя операции сдвига, умножение и деление на четные числа в арифметике по модулю 2^N можно свести к умножению и делению на нечетные числа, для которых существует обратный по умножению элемент (modular multiplicative inverse).
Из того, что для взаимно обратных чисел M и W выполняется равенство M*W=W*M=1 (mod 2^N) следует, что операция квазиделения Y на M, т.е. умножения Y на W=M^-1, взаимно однозначна.
Очевидно, что если заранее известно, что Y делится на M (Y=X*M), то для нахождения частного X можно использовать квазиделение: Y*W=X*M*W=X.
В каждом умножении есть доля деления
Давайте посмотрим, как ведут себя при умножении на обратный элемент W, числа, дающие при делении на M одинаковый остаток R:
(0*M+R)*W=0*M*W+R*W=R*W+0=S(R)+0
(1*M+R)*W=1*M*W+R*W=R*W+1=S(R)+1
(2*M+R)*W=2*M*W+R*W=R*W+2=S(R)+2
Легко заметить, что произведения занимают все целые точки некоторого диапазона значений, начиная с некоторого числа S(R). Причем вследствие взаимной однозначности умножения диапазоны значений произведений, соответствующие разным значениям остатка, не пересекаются, т.е. весь дипазон [0..2^N-1] разбивается на M непересекающихся поддиапазонов примерно равной длины. Заметим также, что первый поддиапазон начинается с S(0)=0.
Назовем произведение остатка на обратный элемент следом остатка S(R)=R*W. Зная величину следа остатка S(R), мы могли бы найти частное вычитанием:
пусть Y=X*M+R, тогда Y*W=X*M*W+R*W=X+R*W=X+S(R), X=Y*W-S(R)
Гадая на остатках
Пусть длина делителя M равна B бит. Тогда, очевидно, длина следа ненулевого остатка должна превышать N-B, иначе в нулевом поддиапазоне окажется слишком мало чисел. То же самое справедливо для любого поддиапазона. Следовательно, старшие B бит любой пары следов различных остатков различны.
Это значит, что, заранее вычислив следы всех остатков, мы могли бы угадать след любого остатка по его старшим B битам. Но проблема в том, что старшие B бит результата квазиделения могут не совпадать со старшими битами следа, например, из-за переноса.
Элегантно избавиться от этой проблемы не удается, поэтому мы просто уменьшим в два раза максимальное допустимое значение делимого и будем угадывать след остатка по его старшим B+1 битам. С одной стороны, это ограничение может добавить в алгоритм деления лишнюю проверку для нечетного делителя. С другой стороны, т.к. теперь результаты квазиделения покрывают нижнюю половину каждого поддиапазона, то старшие B+1 бит могут совпасть только у результатов, соответствующих следам одинаковых остатков.
Алгоритм и примеры
Сначала давайте посмотрим примеры деления с остатком на 10 (проверка переполнения должна быть отключена):
function ShaDiv10(ADividend: cardinal): cardinal;
const
Shift= 1;
Inverse= $CCCCCCCD;
ToIndex= 28;
Stamp: array[0..15] of cardinal= ($00000000, $00000000, $00000001, $33333334,
$33333334, $00000001, $66666667, $66666667,
$66666667, $9999999A, $9999999A, $9999999A,
$CCCCCCCD, $CCCCCCCD, $CCCCCCCD, $00000001);
begin;
Result:=ADividend shr Shift * Inverse;
Result:=Result - Stamp[Result shr ToIndex];
end;
function ShaMod10(ADividend: cardinal): cardinal;
const
Shift= 1;
Inverse= $CCCCCCCD;
ToIndex= 28;
Rest: array[0..15] of cardinal= (0, 0, 10, 8, 8, 10, 6, 6, 6, 4, 4, 4, 2, 2, 2, 10);
begin;
Result:=Rest[ADividend shr Shift * Inverse shr ToIndex] + ADividend and (1 shl Shift - 1);
end;
function ShaDivMod10(ADividend: cardinal; var ARest: cardinal): cardinal;
const
Shift= 1;
Inverse= $CCCCCCCD;
ToIndex= 28;
Stamp: array[0..15] of cardinal= ($00000000, $00000000, $00000001, $33333334,
$33333334, $00000001, $66666667, $66666667,
$66666667, $9999999A, $9999999A, $9999999A,
$CCCCCCCD, $CCCCCCCD, $CCCCCCCD, $00000001);
Rest: array[0..15] of cardinal= (0, 0, 10, 8, 8, 10, 6, 6, 6, 4, 4, 4, 2, 2, 2, 10);
var
i: cardinal;
begin;
Result:=ADividend shr Shift * Inverse;
i:=Result shr ToIndex;
Result:=Result - Stamp[i];
ARest:=Rest[i] + ADividend and (1 shl Shift - 1);
end;и на 100:
function ShaDivMod100(ADividend: cardinal; var ARest: cardinal): cardinal;
const
Shift= 2;
Inverse= $C28F5C29;
ToIndex= 26;
Stamp: array[0..63] of cardinal= (
$00000000, $00000000, $0A3D70A4, $0A3D70A4, $00000001, $147AE148, $147AE148, $1EB851EC,
$1EB851EC, $00000001, $28F5C290, $28F5C290, $33333334, $33333334, $33333334, $3D70A3D8,
$3D70A3D8, $47AE147B, $47AE147B, $47AE147B, $51EB851F, $51EB851F, $00000001, $5C28F5C3,
$5C28F5C3, $66666667, $66666667, $00000001, $70A3D70B, $70A3D70B, $7AE147AF, $7AE147AF,
$7AE147AF, $851EB852, $851EB852, $8F5C28F6, $8F5C28F6, $8F5C28F6, $9999999A, $9999999A,
$A3D70A3E, $A3D70A3E, $A3D70A3E, $AE147AE2, $AE147AE2, $00000001, $B851EB86, $B851EB86,
$C28F5C29, $C28F5C29, $00000001, $CCCCCCCD, $CCCCCCCD, $D70A3D71, $D70A3D71, $D70A3D71,
$E147AE15, $E147AE15, $EB851EB9, $EB851EB9, $EB851EB9, $F5C28F5D, $F5C28F5D, $00000001);
Rest: array[0..63] of cardinal= (
0, 0, 16, 16,100, 32, 32, 48, 48,100, 64, 64, 80, 80, 80, 96,
96, 12, 12, 12, 28, 28,100, 44, 44, 60, 60,100, 76, 76, 92, 92,
92, 8, 8, 24, 24, 24, 40, 40, 56, 56, 56, 72, 72,100, 88, 88,
4, 4,100, 20, 20, 36, 36, 36, 52, 52, 68, 68, 68, 84, 84,100);
var
i: cardinal;
begin;
Result:=ADividend shr Shift * Inverse;
i:=Result shr ToIndex;
Result:=Result - Stamp[i];
ARest:=Rest[i] + ADividend and (1 shl Shift - 1);
end;Инициализация и деление 32-битных чисел в общем виде:
type
TSlotRecord= record
Stamp: cardinal;
Rest: cardinal;
end;
TDivisorRecord= record
Divisor: cardinal;
Shift: cardinal;
Mask: cardinal;
Inverse: cardinal;
Quot32: cardinal;
Num32: cardinal;
ToIndex: cardinal;
Slots: array of TSlotRecord;
end;
PDivisorRecord= ^TDivisorRecord;
procedure InitDivisor(p: PDivisorRecord; ADivisor: cardinal);
const
cMaxDividend= $7FFFFFFF;
cMaxDivisor= $000FFFFF;
var
i, j: integer;
c, m: cardinal;
begin;
if (ADivisor<=0) or (ADivisor>cMaxDivisor) then ADivisor:=ADivisor div (ADivisor-ADivisor);
with p^ do begin;
Divisor:=ADivisor;
Shift:=0; while Divisor and (1 shl Shift)=0 do inc(Shift);
Mask:=1 shl Shift - 1;
m:=Divisor shr Shift;
Inverse:=1; c:=m;
while c<>1 do begin;
Inverse:=Inverse*c;
c:=c*c;
end;
Quot32:=cMaxDividend div m + 1;
Num32:=Quot32 * m;
j:=1; while m>=1 shl j do inc(j);
ToIndex:=SizeOf(Inverse)*8-j-1;
SetLength(Slots, 1 shl (j+1));
for j:=0 to Length(Slots)-1 do with Slots[j] do begin;
Stamp:=1;
Rest:=Divisor;
end;
for i:=0 to m-1 do begin;
c:=Inverse*cardinal(i);
for j:=c shr ToIndex to (c+Quot32-1) shr ToIndex do with Slots[j] do begin;
Stamp:=c;
Rest:=i shl Shift;
end;
end;
if Divisor and 1=0 then begin;
Quot32:=0;
Num32:=0;
end;
end;
end;
function ShaDivMod(p: PDivisorRecord; ADividend: cardinal; var ARest: cardinal): cardinal;
var
c: cardinal;
begin;
with p^ do begin;
Result:=0;
//need for big divisor only
if Divisor shr (SizeOf(Divisor)*8-1)<>0 then begin;
ARest:=ADividend;
if ADividend>=Divisor then begin;
dec(ARest, Divisor);
inc(Result);
end;
exit;
end;
//need for odd divisor only
if (Divisor and 1<>0) and (ADividend>=Num32) then begin;
ADividend:=ADividend-Num32;
Result:=Quot32;
end;
c:=ADividend shr Shift * Inverse;
with Slots[c shr ToIndex] do begin;
Result:=Result + c - Stamp;
ARest:=Rest + ADividend and Mask;
end;
end;
end;Экспериментальная проверка корректности алгоритма была выполнена для типов byte, word, cardinal, uint64. В двух последних случаях проверка выполнялась только для делителей в диапазоне 1..2^20:
function Validate(TypeSize, MaxDivisor: integer): string;
var
b, e, i, j, m, mlast, len: integer;
n, d, x, umask, smask: int64;
a: array of integer;
label
err1, err2;
begin;
if (TypeSize<>1) and (TypeSize<>2) and (TypeSize<>4) and (TypeSize<>8) then begin;
Result:=Format('invalid type size = %d', [TypeSize]);
exit;
end;
TypeSize:=TypeSize*8;
if MaxDivisor=0 then MaxDivisor:=1 shl (TypeSize-1) - 1;
if (MaxDivisor<3) or (MaxDivisor>1024*1024) then begin;
Result:=Format('invalid max divisor = %d', [MaxDivisor]);
exit;
end;
e:=0;
umask:=int64(-1) shr (64-TypeSize);
smask:=umask shr 1;
mlast:=smask and MaxInt; if mlast>MaxDivisor then mlast:=MaxDivisor;
len:=1; while len<MaxDivisor do len:=len*2;
len:=len*2; SetLength(a, len);
m:=3;
repeat;
d:=smask div m * m;
n:=1; x:=m;
while x<>1 do begin;
n:=n*x and umask;
x:=x*x and umask;
end;
b:=1; while m>=1 shl b do inc(b);
inc(b);
for i:=0 to 1 shl b-1 do a[i]:=m;
for i:=0 to m-1 do begin;
j:=n*i and umask shr (TypeSize-b);
if a[j]<>m then goto err1;
a[j]:=i;
j:=n*(d+i) and umask shr (TypeSize-b);
if (a[j]<>i) and (a[j]<>m) then goto err2;
a[j]:=i;
end;
inc(m,2);
until m>=mlast;
Result:='done';
exit;
err2: inc(e);
err1: inc(e);
Result:=Format('error%d on divisor %d',[e, m]);
end;
procedure TForm1.bValidate8Click(Sender: TObject);
begin;
Memo1.Lines.Add(Validate(SizeOf(byte), 0));
end;
procedure TForm1.bValidate16Click(Sender: TObject);
begin;
Memo1.Lines.Add(Validate(SizeOf(word), 0));
end;
procedure TForm1.bValidate32Click(Sender: TObject);
begin;
Memo1.Lines.Add(Validate(SizeOf(cardinal), 64*1024)); //1024*1024
end;
procedure TForm1.bValidate64Click(Sender: TObject);
begin;
Memo1.Lines.Add(Validate(SizeOf(int64), 64*1024)); //1024*1024
end;Описанный алгоритм может применяться в функциях перевода числа из одной системы счисления в другую и им подобных. Например, в функции IntToStr 64-битное число, превышающее MaxInt, удобно разрезать на два фрагмента, длина младшего из которых составляет 10 десятичных цифр, а после в цикле применить квазиделение на 100.