q[j] = (k*b + u[j] )/v[0] ;

к = (k*b + U[j] ) - q[j]*v[0] ;

if (r != NOLL) r[0] = k; return 0;

Нормализация путем сдвига v влево, такого, что

старший бит становится единичным, и сдвига и влево

на ту же величину. Нам может потребоваться добавить

старшую цифру к частному,- мы делаем это безусловно

S = nlz(v[n-l]) - 16; О <= S <= 16.

vn = (unsigned short *)alloca(2*n);

for (i = n - 1; i > 0; i--)

vn[i] = (v[i] « s) I (v[i-l] >> 16-s) ; vn[0] = v[0] << s;

un = (unsigned short *) alloca (2* (m + D); un[m] = u[m-l] >> 16-s; for (i = m - 1; i > 0; i--)

un[i] = (u[i] « s) I (u[i-l] » 16-s); un[0] = u[0] << s;

for (j = m - n; j >= 0; j--) Главный цикл

{ Вычисляем оценку q[j]

qhat = (un[j-fn]*b + un[j+n-l])/vn [n-1] ;

rhat = (un[j-fn]*b + unij-fn-l]) - qhat*vn [n-1] ; again:

if (qhat >= b qhat*vn[n-2] > b*rhat + un[j+n-2])

qhat = qhat - 1;

rhat = rhat + vn[n-l];

if (rhat < b) goto again;

Умножение и вычитание к = 0;

for (i = 0; i < п; i-n-) p = qhat*vn [i] ;

t = un[i-fj] - к - (p & OxFFFF) ; un[i-fj] = t;

к = (p >> 16) - (t >> 16) ;

t = un[j-bn] - k; un[j-)-n] = t;

q[j] = qhat; Сохранение цифры частного if (t < 0) Если мы вычли слишком много, { вернем назад

q[j] = q[j] - 1;

к = 0;

for (i = 0; i < п;

t = un[i-bj] + vn[i] + k; un[i-fj] = t;

к = t >> 16;

unCj-fn] = un[j-fn] + k; } for j

Если вызывающей функции нужно значение остатка, денормализуем и возвращаем его if (г != NULL)

for (i =0; i < n; i++)

r[i] = (un[i] >> s) I (un[i+l] « 16-s);

return 0;

Алгоритм обрабатывает входные и выходные данные по полслова. Конечно, бьшо бы предпочтительнее обрабатывать данные по слову, но такой алгоритм требует наличия команды деления 64+32 =* 32. Однако предполагается, что такой команды в компьютере нет (или она недоступна из языка высокого уровня). Хотя при этом предполагается наличие коммщы

деления 32+32 => 32, для решения нашей задачи достаточно команды 32+16 => 16.

Таким образом, в данной реализации алгоритма Кнута основание счисления b равно 65536. В [39] вы найдете подробное пояснение данного алгоритма.

Делимое и и делитель v используют прямой порядок, т.е. и [0] и v [0] представляют собой младшие цифры числа (впрочем, данный код корректно работает при использовании как прямого, так обратного порядка). Параметры тип представляют собой количество полуслов в и и V соответственно (Кнут определяет m как длину частного). Вызывающая функция должна обеспечить память для хранения частного q и (необязательно) для остатка г. Пространство для частного должно иметь размер, как минимум, m-n+l полуслов и п полуслов для остатка. Если значение остатка нас не интересует, вместо адреса для его размещения можно передать параметр NULL.

Алгоритм требует, чтобы старшая цифра делителя v[n-i] была ненулевой. Это упрощает нормализацию и помогает убедиться, что вызывающая функция вьщелила достаточное количество памяти для частного. Код проверяет, является ли значение v[n-l] ненулевым, а также выполнение условий п > 1 и m 2 п. Если любое из этих условий нарушено, возвращается код ошибки (значение 1).

После этих проверок код выполняет деление для частного случая, когда делитель имеет длину 1. Это выделение частного случая не для ускорения работы; просто оставшаяся часть алгоритма требует, чтобы длина делителя была не менее 2.

Если делитель имеет длину 2 или большую, алгоритм нормализует делитель, сдвигая его влево на величину, достаточную, чтобы старший бит делителя был равен 1. На ту же величину сдвигается и делимое, так что эти сдвиги никак не влияют на величину частного. Как поясняется у Кнута, эти шаги необходимы для того, чтобы облегчить оценку каждой цифры частного с хорошей точностью. Для определения величины сдвига используется функция количества ведущих нулевых битов nlz(je).

При выполнении нормализации для делимого и делителя выделяется новая память. Это делается в связи с тем, что с точки зрения вызывающей функции крайне нежелательно изменение входных данных, которые вообще могуг оказаться константами, расположенными в памяти( доступной только для чтения. Кроме того, из-за сдвига делимому может понадобиться дополнительное полуслово для хранения старшей цифры. Для выделения этой памяти в С идеально подходит функция а Носа (), которая обычно выделяет память с помощью двух-трех команд и не требует явного ее освобождения. Эта память выделяется в программном стеке, так что она автоматически освобождается при возврате из функции.

В основном цикле происходит быстрое вычисление цифр частного, по одной цифре за итерацию; делимое прн этом уменьшается и в конечном счете становится равным остатку. Оценка значения очередной цифры частного qhat после уточнения в цикле, помеченном меткой again, всегда оказывается либо точной, либо больше точного значения на 1.

Следующие шаги состоят в умножении qhat на делитель и вычитании произведения из полученного остатка, как и в методе деления столбиком. Если остаток оказывается отрицательным, необходимо уменьшить цифру частного на 1 и либо заново вычислить произведение и вычесть его из остатка, либо, что гораздо проще, прибавить к отрицательному значению остатка делитель. Такой шаг делается не более одного раза, поскольку цифра частного либо точна, либо больше точного значения на 1.

Последнее действие состоит в возврате остатка вызывающей функции, если переданный адрес памяти для ?фанения остатка не NULL. Остаток при этом должен быть сдвинут вправо на ту же величину, на которую в целях нормализации сдвигались влево делитель и делимое.

Шаги, связанные с неверной оценкой цифры частного, выполняются достаточно редко. Чтобы убедиться в этом, заметим, что первое вычисление оценки каждой цифры частного qhat выполняется путем деления двух старших цифр текущего остатка на старшую цифру делителя. Шаги цкла again уточняют оценку путем деления трех старших цифр текущего остатка на две старшие цифры делителя (доказательство опущено; убедиться в корректности этих действий можно, проведя несколько тестовых вычислений с основанием системы счисления b = ю). Заметим, что из-за выполненной нормализации делитель имеет значение, не меньшее ь/2, а делимое не более чем в b раз превышает делитель (поскольку любой остаток меньше делителя).

Насколько точна оценка частного при использовании только трех цифр делимого р двух - делителя? Можно показать, что в связи с выполненной нормализацией такая оценю довольно точна. Чтобы увидеть это в какой-то мере интуитивно (без формального доказатель сгва), рассмотрим оценку u/v в арифметике с основанием счисления 10. Можно показать, чтс такая оценка всегда несколько завышена (или точна). Следовательно, наихудшим случаем яв ляется тот, при котором усечение делителя до двух цифр максимально уменьшает его, а усе чение делимого до трех цифр не уменьшает его величину, которая должна быть максимальго возможной. Это происходит в случае 49900.. .0/5099.. .9, что приводит к оценке 499/50=9.98 Точный результат составляет примерно 499/51=9.7843. Разница в 0.1957 показывает, чт( оценка цифры частного отличается от истинного значения не более чем на 1, причем это про исходит примерно в 20% случаев (в предположении, что цифры частного распределены рав номерно). Это, в свою очередь, означает, что дополнительные действия, связанные с неточно) оценкой цифры частного, будуг выполняться примерно в 20% случаев.

Проведение этого (нестрогого) анализа для общего случая системы счисления b при водит к выводу о том, что оценочное и истинное значение цифры частного отличается н более чем на 2/Ь. В случае /?=65536 мы получим, что оценочное и истинное значени

цифры частного отличаются не более чем на 1 и происходит это с вероятностью пример но 2/65536=«0.00003. Следовательно, дополнительные действия выполняются тольк примерно для 0.003% цифр частного.

В качестве конкретного примера, когда требуются дополнительные корректн рующие действия, в десятичной системе счисления можно привести делени 4500/501. Аналогичный пример для системы счисления по основанию 65536-Ox7FFF 8 ООО 0000 0000/0x8000 0000 ООО 1.

Не будем пытаться оценить время работы программ в целом, ограничимся лишь з: мечанием, что для больших тип время выполнения определяется циклом умнож* ния/вычитания. Хороший компилятор в состоянии использовать для этого цикла всего 1 базовых RISC-команд, одна из которых - умножение. Цикл for j выполняется п раз, цикл умножения/вычитания - т-п+1 раз, что дает время выполнения этой части прс граммы, равное (15 +mul)п(т-п + \) тактов, где mul- время умножения двух К