Когда команды нет: деление и корень из сдвигов и вычитаний

Если ты разделишь 100 на 7 столбиком в двоичной системе, найдёшь целый корень по одному биту и объяснишь, почему квадрат кандидата может переполниться, когда само число помещается, - листай до «Теперь сам».

Как думаешь, как делит процессор без команды деления

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

Прежде чем читать дальше, вспомни деление столбиком и ответь себе: какая часть этого способа самая трудная и что с ней станет, если цифры бывают только 0 и 1? Трудная часть - угадать очередную цифру частного. В двоичной системе угадывать нечего.

На что это похоже

Деление столбиком: сносишь очередную цифру делимого, смотришь, сколько раз делитель помещается в то, что получилось, пишешь эту цифру в частное, вычитаешь. «Сколько раз помещается» - от 0 до 9, и это приходится подбирать.

В двоичной записи вопрос «сколько раз» превращается в «помещается или нет». Поместился - пишешь единицу и вычитаешь делитель один раз. Не поместился - пишешь ноль и ничего не вычитаешь. Ни подбора, ни умножения делителя на цифру: одно сравнение и, может быть, одно вычитание на каждый бит.

Механизм: три алгоритма, один шаблон

Умножение. Множитель - сумма степеней двойки, по одной на единичный бит. Для каждого единичного бита i к ответу прибавляется множимое, сдвинутое на i. Это умножение на константу из прошлого урока, только константа неизвестна заранее и её биты проверяются в цикле.

Деление. Остаток начинается с нуля. На каждом круге в него справа приписывается очередной бит делимого, со старшего. Если остаток стал не меньше делителя - делитель вычитается, и в частное пишется единица.

100 / 7:    100 = 1100100₂

бит делимого    остаток    ≥ 7?         бит частного
     1             1        нет              0
     1             3        нет              0
     0             6        нет              0
     0            12        да, → 5          1
     1            11        да, → 4          1
     0             8        да, → 1          1
     0             2        нет              0

частное 0001110₂ = 14,  остаток 2

Корень. Целый корень из n - наибольшее r, для которого r² ≤ n. Биты r определяются со старшего: включаем очередной бит и смотрим, не перескочили ли. Если (r + b)² ≤ n - бит остаётся.

Все три устроены одинаково: на каждом круге - один бит, сравнение и одно сложение или вычитание.

Правило

Умножение: для каждого единичного бита i множителя прибавить a << i. Деление: для битов делимого со старшего - r = r << 1 | бит; если r ≥ d, то r = r − d и бит частного равен единице. Корень: для b от старшей степени двойки вниз - если (r + b)² ≤ n, то r = r + b. Умножение делает круг на каждый бит множителя, деление - на каждый бит делимого, корень - на каждый бит ответа. Умножению и делению хватает сдвига, сравнения, сложения и вычитания; корню - тоже, если вести остаток вместо квадрата.

Разбор: умножение и деление

def umnozhit(a, b):
    p = 0
    while b > 0:
        if b & 1:
            p = p + a
        a = a << 1
        b = b >> 1
    return p

def razdelit(n, d):
    q = 0
    r = 0
    for i in range(n.bit_length() - 1, -1, -1):
        r = r << 1 | n >> i & 1
        if r >= d:
            r = r - d
            q = q | 1 << i
    return [q, r]

print(umnozhit(13, 11), razdelit(100, 7), razdelit(1000000, 3), razdelit(5, 9))

143 [14, 2] [333333, 1] [0, 5]

umnozhit - та же программа, что быстрая степень из раздела про числа, только вместо умножения стоит сложение, а вместо возведения в квадрат - удвоение. Это не совпадение: степень - повторённое умножение, произведение - повторённое сложение, и обе раскладываются по битам второго числа.

В razdelit строка r << 1 | n >> i & 1 приписывает к остатку бит номер i: сдвиг освобождает младший разряд, | кладёт туда бит. Остаток всегда меньше делителя, поэтому после сдвига он меньше 2d, и вычесть делитель больше одного раза нельзя - цифра частного действительно только 0 или 1.

Деление 4 000 000 000 на 7 этим способом стоит 132 шага - 32 круга по четыре строки. Встроенное // тренажёр считает за один шаг - и здесь видно, что стоит за этим шагом. Схема деления в настоящем процессоре определяет за такт не один бит частного, а несколько, но всё равно работает в несколько раз дольше умножителя.

Разбор: корень по битам, и почему квадрат лучше не считать

def koren_prosto(n):
    r = 0
    b = 1 << 26
    while b > 0:
        if (r + b) * (r + b) <= n:
            r = r + b
        b = b >> 1
    return r

print(koren_prosto(180), koren_prosto(1000000), koren_prosto(999999))

13 1000 999

Начальный бит 2²⁶ - старший бит корня из числа меньше 2⁵³. Это работает, пока не встретится большое число. koren_prosto(4503599627370496) - корень из 2⁵² - останавливает программу в тренажёре: после того как бит 2²⁶ принят, следующий кандидат равен 2²⁶ + 2²⁵, и его квадрат больше 2⁵³, предела точного счёта. Само n помещается, ответ помещается, а промежуточный квадрат - нет. Настоящий Python это переживёт - его целые растут; тренажёр здесь ведёт себя как язык со словом фиксированной ширины. В C с 64-битными знаковыми числами то же случается у верхней границы типа, только без сообщения об ошибке.

Выход - не считать квадрат вовсе. Квадрат суммы: (r + b)² = r² + 2rb + b². Если хранить не n, а остаток n − r², то проверка (r + b)² ≤ n превращается в остаток ≥ 2rb + b². А 2rb и b² можно вести сдвигами: b² - степень четвёрки, она сдвигается на два разряда за круг. Когда бит принят, значение 2rb для следующего, вдвое меньшего веса - это 2(r + b) · b/2 = rb + b²: старое пополам плюс b². Когда не принят - просто старое пополам.

def koren(n):
    r = 0
    bit = 1 << 52
    while bit > n:
        bit = bit >> 2
    while bit > 0:
        if n >= r + bit:
            n = n - (r + bit)
            r = (r >> 1) + bit
        else:
            r = r >> 1
        bit = bit >> 2
    return r

print(koren(180), koren(4503599627370496), koren(9007199136250224), int(9007199136250224 ** 0.5))

13 67108864 94906264 94906265

Здесь bit - это b², остаток n − r² хранится в самой переменной n, а в переменной r во время работы лежит не сам корень, а 2rb - набранная часть корня, умноженная на удвоенный текущий вес. К концу, когда b доходит до единицы и ещё раз делится пополам, множитель исчезает, и в r остаётся корень. Для 180:

остаток    r      bit    остаток ≥ r + bit?
  180       0      64     да:  остаток 116, r = 64
  116      64      16     да:  остаток 36,  r = 48
   36      48       4     нет:              r = 24
   36      24       1     да:  остаток 11,  r = 13

Ответ 13, остаток 11: 13² + 11 = 180. Ни одно число по дороге не вышло за 180. В общем случае сумма r + bit бывает больше n, но не больше чем на четверть - до переполнения далеко.

Последнее число в выводе - граница другого рода. 9 007 199 136 250 224 - это 94 906 265² − 1, и настоящий корень из него - 94 906 264. Корень через дробные числа возвращает 94 906 265: точный результат меньше целого на 1 / (2 · 94 906 265), пять миллиардных, а соседние дробные числа такой величины стоят в 2⁻²⁶, пятнадцати миллиардных, друг от друга, и результат округляется до целого. Ошибка в одной единице, зато в самом неудобном месте - у полного квадрата.

Есть и способ быстрее побитового - метод Ньютона: x = (x + n // x) // 2, пока x убывает. На том же числе от начального значения 2²⁷ он сходится за 6 кругов вместо 27, потому что каждый круг удваивает число верных битов. Но в круге у него деление, и там, где деления нет, выбирают сдвиги.

Что ломается без этого

Первое: int(n ** 0.5) и math.sqrt. До 2⁵² результат верен, дальше - лотерея у полных квадратов. В Python для целого корня есть точная math.isqrt.

Второе: проверка r * r <= n в цикле. Привычная запись «пока квадрат не больше» переполняется у границы типа. Безопасная форма - r <= n // r: делением вместо умножения.

Третье: деление на ноль. В razdelit при d = 0 условие r >= d истинно всегда, и функция молча возвращает частное из одних единиц. Команда процессора в этом случае вызывает аварийное прерывание; самодельное деление обязано проверить делитель само.

Четвёртое: знаки. Все три алгоритма - для неотрицательных чисел. Знаковое деление делают так: запомнить знаки, разделить модули, поставить знак частному. И тут возвращается −2³¹, у которого нет модуля.

Зачем это в работе

Микроконтроллеры и встраиваемые системы: если команды деления нет, компилятор подставляет вместо / вызов библиотечной функции, внутри которой - цикл из разбора. Поэтому в коде прошивок деления избегают, а делители выбирают степенями двойки.

Длинная арифметика: деление числа из тысячи разрядов - тот же столбик, только «цифра» - это 32 или 64 бита, и её уже приходится подбирать. Деление в Python над длинными числами устроено так.

Целый корень нужен чаще, чем кажется: граница перебора делителей в проверке простоты, сторона квадратной сетки, в которую укладывается n элементов, проверка «полный ли это квадрат» в задачах по теории чисел. Везде, где важен точный ответ на границе, дробный корень не годится.

Теперь сам

1. Раздели 45 на 6 двоичным столбиком: выпиши остаток после каждого приписанного бита.
2. Сколько кругов сделает umnozhit(a, 1000) и сколько в них будет сложений?
3. Найди koren(50) по таблице: остаток, r, bit на каждом круге.
4. Почему в razdelit делитель никогда не вычитается дважды за круг?
5. Тип - 64 бита со знаком, предел около 9,2 · 10¹⁸. При n = 10¹⁸ кандидат r + b не поднимается выше 1,1 · 10⁹. Поместится ли его квадрат? А при n = 9 · 10¹⁸, когда кандидат равен 2³¹ + 2³⁰?