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

Ограничимся делением двухразрядных чисел без знака. Деление чисел большей разрядности можно обобщить, при необходимости обратившись к первоисточнику [1, п.4.3.1]. Описываемый алгоритм назовем «программный 128/128». Заметим, что во многих 64-битных компиляторах он уже реализован (GCC, Clang, Intel Compiler) и может быть использован напрямую без изобретения велосипеда.

Цель данной статьи — подробно объяснить детали алгоритма, чтобы снизить порог входа в общеизвестные труды Д. Кнута, в том числе объяснить почему деление в процессоре дает лишь одноразрядное частное (конкретно для 64-битных процессоров можно делить 128-битное число на 64-битное, получая лишь 64-битные частное и остаток). Назовем процессорный алгоритм деления как «аппаратный 128/64».

Ключевым моментом в понимании алгоритмов деления является процесс нормализации чисел, который и позволяет смело пользоваться встроенным в процессор делением 128/64, не опасаясь непредсказуемого поведения процессора.

Алгоритм деления двухразрядных чисел можно классифицировать по фактической разрядности делителя на два алгоритма:

1. Половинчатое программное деление, когда делитель одноразрядный,

2. Полное программное деление, когда делитель двухразрядный.

Заметим, что аппаратный алгоритм 128/64 является половинчатым; он будет использован в обеих ветках программного алгоритма. Несмотря на привязку к 64-битным процессорам, в описываемых алгоритмах будем использовать обозначения для -битного процессора.

Половинчатое программное деление

Задача - разделить два числа

Здесь делитель занимает один разряд (пусть величина разряда равна , а битовая ширина равна ), а делимое - два разряда. Перед делением делаем нормализацию делимого и делителя, которая не меняет частного ; остаток при этом придется масштабировать обратно, чтобы получить оригинальный остаток . Процесс нормализации заключается в одновременном умножении двух чисел (делимого и делителя) на некоторое число до тех пор, пока делитель не станет больше (или равным) половине разряда:

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

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

Это и позволяет автоматически воспользоваться аппаратным алгоритмом "128/64", что упрощает (и ускоряет) программный алгоритм деления.

Получившийся остаток, который в данном граничном случае равен (смотри пояснение ниже*), если к нему приставить младший разряд делимого, будет двухразрядным числом, которое также можно разделить на делитель аппаратным алгоритмом "128/64":

*Величину остатка можно получить, принимая частное за , и прямо вычитая результат умножения (частного на делитель) из делимого:

при этом получившийся здесь остаток следует разделить на два, потому что мы взяли делимое и делитель из середины формулы, где они умножены на два:

Два получившихся частных следует объединить в двухразрядное число, которое даст итоговое частное от деления, а остаток от последнего (аппаратного) деления следует нормализовать в обратную сторону, то есть привести к исходному масштабу (делением на два Q раз, если нормализация заключалась в умножении на два Q раз).

Отметим, что половинчатое деление (помимо нормализации) в итоге требует два аппаратных деления.

Алгоритм половинчатого деления:

1. Умножением на два нормализовать делимое и делитель , так, чтобы старший бит делителя был равен единице. После нормализации делимое обозначается как , а делитель - как .

2. Разделить двухразрядное делимое, образованное из старших разрядов трехразрядного делимого , на одноразрядный делитель, используя аппаратный алгоритм деления.

3. Скомпоновать из найденного в пункте 2 остатка и младшего разряда делимого новое двухразрядное делимое.

4. Разделить новое двухразрядное делимое на одноразрядный делитель, используя аппаратный алгоритм деления.

5. Из найденных на этапах 2 и 4 частных скомпоновать двухразрядное частное, которое будет итоговым частным.

6. Найденный на этапе 4 остаток разделить на два столько раз, сколько раз было умножение на два в процедуре нормализации. Получившийся остаток будет итоговым остатком.

Полное программное деление

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

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

Отметим, что полное программное деление (помимо нормализации) требует одно аппаратное деление и одно расширяющее умножение с одним-двумя вычитаниями.

Алгоритм полного деления:

1. Умножением на два нормализовать делимое и делитель , так, чтобы старший бит делителя был равен единице. После нормализации делимое обозначается как , а делитель - как .

2. Разделить двухразрядное делимое, образованное из двух старших разрядов трехразрядного делимого , на старший разряд делителя, используя аппаратный алгоритм деления.

3. Скомпоновать из найденного в пункте 2 остатка и младшего разряда делимого двухразрядный остаток .

4. Умножить младший разряд делителя на найденное в пункте 2 частное , формируя некоторое двухразрядное число .

5. Если , то сформировать остаток . Перейти к пункту 7.

6. Если , то сформировать остаток и декрементировать частное.

7. Найденный на этапах 5 или 6 остаток разделить на два столько раз, сколько раз было умножение на два в процедуре нормализации. Получившийся остаток будет итоговым остатком. Найденное частное будет итоговым частным.

import random

def generate_random_unsigned_int(bitwidth):
    """
    Генерирует случайное беззнаковое число с заданной битовой шириной.
    Генерируемое число N удовлетворяет неравенству 0 <= N < 2^bitwidth.
    """
    if bitwidth < 0:
        raise ValueError("Битовая ширина должна быть неотрицательной.")
    return random.getrandbits(bitwidth)

def _div_cpu(x: int, y: int, b: int):
    """
    Тестовая заглушка (mock) деления двухразрядного числа x = L1 + H1*b на 
    одноразрядное y = L2, таких, что частное и остаток - одноразрядные 
    числа. Данная функция должна быть заменена соответствующим интринсиком 
    (например, интринсик _div128 MSVC), который компилятор заменит 
    соответствующими инструкциями процессора для заданной целевой 
    архитектуры.

    Args
        x - делимое.
        y - делитель.
        b - основание позиционной системы счисления 
        (например, 2^64 при делении 128-битного числа на 64-битное).

    Returns
        Частное q и остаток r.
    """
    assert(x >= 0)
    assert(y > 0)
    assert(x < b*b)
    assert(y < b)
    q = x // y # CPU: (a, b) / (c) => Одноразрядные частное и остаток.
    r = x % y  # CPU
    assert(q < b)
    return (q, r)

def _divh(x: int, y: int, b: int):
    """
    Делит двухразрядное число x = L1 + H1*b на одноразрядное y = L2.

    Args
        x - делимое.
        y - делитель.
        b - основание позиционной системы счисления.

    Returns
        Частное q и остаток r.
    """
    c = 0
    # Нормализация
    while y < (b // 2):
        c += 1
        x <<= 1
        y <<= 1
    assert(y < b)     # y - некоторое 1-разрядное число.
    assert(x < b*b*b) # x - некоторое 3-разрядное число.
    #
    x_2d = x // b       # Двухразрядное число.
    low_x = x & (b - 1) # Младший разряд x.
    q, r = _div_cpu(x_2d, y, b)
    r *= b
    r += low_x
    q1, r = _div_cpu(r, y, b)
    q *= b
    q += q1
    r >>= c # Денормализация.
    return (q, r)

def div(x: int, y: int, w: int):
    """
    Делит двухразрядное число x = L1 + H1*b на двухразрядное y = L2 + H2*b.
    Здесь b - "лимб" в терминологии библиотеки GMP (здесь будет 
    упоминаться как "разряд"). Основан на алгоритме деления из тома 2 
    "Искусства программирования" Дональда Кнута.
    В итоге, используется аппаратное деление двухразрядного числа на 
    одноразрядное, которое дает одноразрядные частное и остаток 
    (например, интринсик _div128 MSVC).
    Некоторые компиляторы (GCC, Clang, Intel Compiler) включают готовый 
    алгоритм (тип unsigned __int128, поддержка проверяется условием 
    #if defined(__SIZEOF_INT128__)).

    Args
        x - делимое.
        y - делитель.
        w - битовая ширина (делимое и делитель не должны превышать 
        эту ширину).

    Returns
        Частное q и остаток r.
    """
    assert( y > 0 )
    assert( x >= 0 )
    assert( (w % 2) == 0 )
    half_w = w >> 1      # Битовая полуширина.
    high_y = y >> half_w # Старший разряд y.
    b = 1 << half_w
    if high_y == 0:
        return _divh(x, y, b)
    c = 0
    # Нормализация
    while high_y < (b // 2):
        c += 1
        x <<= 1
        y <<= 1
        high_y = y >> half_w
    #
    x_2d = x >> half_w # 2-разрядное число.
    assert(x_2d < b*b) # x_2d - некоторое 2-разрядное число.
    q, r = _div_cpu(x_2d, high_y, b)
    low_x = x & (b - 1) # Младший разряд x.
    low_y = y & (b - 1) # Младший разряд y.
    x_2d = low_x + (r << half_w) # Новое 2-разрядное число.
    assert(x_2d < b*b)
    T = q * low_y # Учет младшего разряда y.
    assert(T < b*b)
    r = x_2d - T
    if x_2d < T: # Коррекция частного и остатка.
        q -= 1
        r += y
    r >>= c # Денормализация.
    return (q, r)

bit_width = 12 # Задаваемая битовая ширина чисел.
for i in range(10_000_000): # Тест на случайных данных.
    x = generate_random_unsigned_int(bit_width)
    y = generate_random_unsigned_int(bit_width)
    if y == 0:
        continue
    q1, r1 = div(x, y, bit_width) # Проверяемый алгоритм.
    q2, r2 = x // y, x % y # Эталонный алгоритм.
    assert((q1, r1) == (q2, r2))

Кое что про обобщение деления на многоразрядные числа

Например, требуется разделить два трехразрядных числа (делитель имеет ненулевой третий разряд). После процедуры нормализации выходит так, что следует разделить четырехразрядное число на трехразрядное. Мы можем взять три старших разряда делимого и разделить на два старших разряда делителя, используя описанное выше полное программное деление (правда, без этапа нормализации, так как он уже сделан). И далее применить коррекцию неучтенного младшего разряда делителя аналогично пунктам 4-6 алгоритма полного программного деления. И так далее рекурсивно... В итоге, базовым при реализации арифметики произвольной точности будет алгоритм деления разрядного числа на нормализованное разрядное, что согласуется с [1].

Вывод

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

Ключевым моментом в алгоритмах деления является процедура нормализации, которая позволяет использовать аппаратные возможности процессора и упрощает алгоритм программного деления.

Источники

1. Дональд Э. Кнут. Искусство программирования, Т.2., Получисленные алгоритмы, 3-е издание, 2018.

UPD. Добавил код на Python.