В процессе обновления своей библиотеки для RSA-шифрования столкнулся с ошибкой в распространённом алгоритме деления больших чисел. Причём эта ошибка была обнаружена мной совершенно случайно, ибо возникает она лишь в очень редких случаях и нигде не упоминается.
Алгоритм деления больших чисел на основании догадки qGuess в том виде, в котором он ходит по интернету:
Обозначим делитель B=(b0, b1, …, bn-1 ), делимое A=(a0, a1, …, an+m-1), числа имеют длины n и n+m соответственно.
Общий цикл остается почти таким же: алгоритм состоит из последовательно выполняемых шагов, где шаг представляет собой деление (n+1)-разрядной части A на n-разрядный делитель B(исключение может составлять первый шаг, но это укладывается в общую схему, которая будет обсуждаться ниже). При этом i-й шаг состоит из двух действий:
1. Угадать i-ю цифру частного q[i]
2. Не создавая временных чисел, вычесть “сдвинутое” произведения B*q[i] из A.
При этом сдвиг B относительно A на каждом шаге уменьшается.Как заставить компьютер генерировать правильное частное q ? Или, хотя бы, достаточно близкое к нему число ?
Довольно интересный способ состоит в высказывании догадки qGuess по первым цифрам делителя и делимого. Понятно, что этих нескольких цифр недостаточно для гарантированно правильного результата, однако неплохое приближение все же получится.
Пусть очередной шаг представляет собой деление некоторого U = (u0, u1, …, un) на B=(b0, b1, …, bn- 1). Если bn-1 >= BASE/2, то можно доказать следующие факты:
1. Если положить qGuess = (un*BASE + un-1) / bn-1 , то qGuess-2 >= q >= qGuess. Иначе говоря, вычисленная таким способом “догадка” будет не меньше искомого частного, но может быть больше на 1 или 2.
2. Если же дополнительно выполняется неравенство qGuess*bn-2 > BASE*r + un-2 , где r – остаток при нахождении qGuess и qGuess != BASE, то qGuess-1 >= q >= qGuess, причем вероятность события qGuess = q+1 приблизительно равна 2/BASE.Таким образом, если bn-1 > BASE/2, то можно вычислить qGuess = (un*BASE + un-1) / bn-1 и уменьшать на единицу до тех пор, пока не станут выполняться условия (2). Получившееся
значение будет либо правильным частным q, либо, с вероятностью 2/BASE, на единицу большим числом.Что делать, если bn-1 слишком мало, чтобы пользоваться таким способом ? Например, можно домножить делитель и делимое на одно и то же число scale = BASE / ( bn-1 +1 ). При этом несколько изменится способ вычисления остатка, а частное останется прежним. Такое домножение иногда называют нормализацией числа.
На тот случай, если qGuess получилось все же на единицу большим q, будем использовать в пункте (b) вычитание, аналогичное функции Sub, которое вместо отрицательного числа даст дополнение до следующей степени основания. Если такое произошло, то последний перенос будет равен borrow = -1. Это сигнал, что необходимо прибавить одно B назад. Заметим, что в конце сложения будет лишний перенос carry=1, о котором нужно забыть (он компенсирует borrow=-1).
Однако в этом алгоритме есть как минимум одна ошибка и один подводный камень для реализации.
Начну с ошибки. Пусть мы работает с основанием BASE=10. Возьмём два числа A и B, состоящие из 4-х разрядов и равные 9999. Перемножим их, R = 99980001. А теперь разделим R на A (в результате должно получиться B и остаток равный 0). Попробуйте выполнить это деление в столбик (советую обратить внимание что qGuess может быть равен основанию или даже на единицу больше его). На первом же шаге деления вы столкнётесь с тем, что алгоритм показывает нам такое: 10 >= q >= 11 (это при том, что верное q равно 9). А вот этого “Таким образом, если bn-1 > BASE/2, то можно вычислить qGuess = (un*BASE + un-1) / bn-1 и уменьшать на единицу до тех пор, пока не станут выполняться условия (2).” вообще никогда нельзя делать. Чтобы понять почему, просто обратите внимание на последнее деление в моём примере.
Как мы видим выше qGuess может запросто быть больше основания и с этим и сопряжен подводный камень. Дело в том что если операция умножения имеет вид RAX * [R/M64] = RDX:RAX (EAX * [R/M32] = EDX: EAX у 32-х биных систем), то операция деления выглядит как RDX:RAX / [R/M64] = RAX (остаток от деления помещается в RDX). Как вы понимаете поместить число, равное или превосходящее BASE (в нашем случае 2^64), просто некуда. Деление с результатом равным или превосходящим BASE попросту заканчивается исключением. А так как такой случай крайне редок и нигде специально не заостряют внимание на этом факте (хотя бы потому что во всех примерах используются очень маленькие основания), то в этом месте можно запросто сделать ошибку, которую потом будет трудно обнаружить.