Быстрое нахождениe остатка от деления больших чисел для делителей специального вида
В этой статье я расскажу об одном способе вычисления x mod p, для p вида (2 ** n — omega), причём omega значительно меньше 2 ** n. Напишу генератор констант на Python. Приведу пару игрушечных примеров на С++, для которых может быть выполнено исчерпывающее тестирование для всех возможных аргументов. А в качестве серьёзной проверки — вычислю 97! mod (2 ** 256 — 2 ** 32 — 977).
Зачем это нужно?
Если у вас под рукой нет подходящей библиотеки для длинной арифметики (или чудесного Python’а);
Если вам доступны только сложения, умножения и битовые сдвиги, а нужно реализовать вычисление остатка от деления;
Для оптимизации под конкретные делители, если скорость работы универсальных библиотек не устраивает.
Теория
Дано:
Делитель p вида (2 ** n — omega), причём omega значительно меньше (2 ** n);
Делимое x, размера m бит, причём m > n;
Размер битового слова s, s <= n, n и m нацело делятся на s.
Для начала, заметим что (2 ** n) mod p = (2 ** n) mod (2 ** n — omega) == omega.
Пусть x_low — младшие n бит числа x, а x_high — старшие (m — n) бит числа x, т.е. x можно представить в виде x_low + x_high * (2 ** n).
Тогда x mod p = (x_low + x_high * (2 ** n)) mod p == x_low + x_high * omega (mod p). Применение этой формулы не гарантирует, что результат окажется в диапазоне [0; p — 1], однако при достаточно большом количестве рекурсивных применений результат окажется в диапазоне [0; (2 ** n) — 1]. С вероятностью примерно omega / (2 ** n) потребуется вычесть p из результата, чтобы он оказался в диапазоне [0; p — 1].
Разобьем x на слова w[i] размера s бит: x = sum (i = 0 … (m / s — 1); w[i] * 2 ** (i * s)). К такому представлению делимого также можно применить формулу x mod p == x_low + x_high * omega (mod p); в результате показатель степени у старших слов уменьшится, а коэффициент при слове будет домножен на omega. Применив формулу рекурсивно достаточное количество раз можем добиться того, что коэффициент при каждом битовом слове будет меньше (2 ** n). То есть, будет получено представление x mod p в виде суммы битовых слов w[i] размера s бит, домноженных на коэффициенты размером не более n бит.
Генератор коэффициентов на Python
def build_reducer(input_bits, target_bits, limb_size_bits, omega):
num_input_limbs = input_bits // limb_size_bits
target_limit = 2 ** target_bits
def split_and_reduce(coeff):
low_part = coeff % target_limit
high_part = coeff // target_limit
return low_part + high_part * omega
coeffs = [2 ** (limb_size_bits * i) for i in range(num_input_limbs)]
while any([k >= target_limit for k in coeffs]):
coeffs = [split_and_reduce(k) for k in coeffs]
return coeffs
def print_reducer(coeffs, target_bits):
for k in coeffs:
print('%0*x' % (target_bits // 4, k))Аргументы:
input_bits — количество бит на входе (m);
target_bits — желаемое количество бит на выходе (n);
limb_size_bits — размер битового слова (s);
omega — одна из констант, задающих делитель (2 ** n — omega).
Принцип работы:
Определяем количество битовых слов размера s, необходимых для представления числа из m бит;
Находим степень двойки (2 ** n), по границе которой будем разбивать каждый коэффициент на «младшие» (1…n) и «старшие» (n+1…m) биты;
Находим степени двойки, являющиеся коэффициентами битовых слов в разложении делимого;
Пока хотя бы один из коэффициентов превышает (2 ** n) — дробим такие коэффициенты на «младшие» и «старшие» биты и пересобираем из них новый коэффициент: («младшие биты» + «старшие биты» * omega).
Примеры работы генератора коэффициентов:
Коэффициенты для сокращения 32-битного числа до 8-битного по модулю (2 ** 8 — 17) = 239; размер слова — 8 бит:
>>> r = build_reducer(32, 8, 8, 17)
>>> print_reducer(r, 8)
01
11
32
85Коэффициенты для сокращения 32-битного числа до 16-битного по модулю (2 ** 16 — 666) = 64870; размер слова — 8 бит:
>>> r = build_reducer(32, 16, 8, 666)
>>> print_reducer(r, 16)
0001
0100
029a
9f34Коэффициенты для сокращения 512-битного числа до 256-битного по модулю (2 ** 256 — 2 ** 32 — 977); размер слова — 32 бита; для удобства разряды сгруппированы по 32 бита:
>>> r = build_reducer(512, 256, 32, 2 ** 32 + 977)
>>> print_reducer(r, 256)
00000000_00000000_00000000_00000000_00000000_00000000_00000000_00000001
00000000_00000000_00000000_00000000_00000000_00000000_00000001_00000000
00000000_00000000_00000000_00000000_00000000_00000001_00000000_00000000
00000000_00000000_00000000_00000000_00000001_00000000_00000000_00000000
00000000_00000000_00000000_00000001_00000000_00000000_00000000_00000000
00000000_00000000_00000001_00000000_00000000_00000000_00000000_00000000
00000000_00000001_00000000_00000000_00000000_00000000_00000000_00000000
00000001_00000000_00000000_00000000_00000000_00000000_00000000_00000000
00000000_00000000_00000000_00000000_00000000_00000000_00000001_000003d1
00000000_00000000_00000000_00000000_00000000_00000001_000003d1_00000000
00000000_00000000_00000000_00000000_00000001_000003d1_00000000_00000000
00000000_00000000_00000000_00000001_000003d1_00000000_00000000_00000000
00000000_00000000_00000001_000003d1_00000000_00000000_00000000_00000000
00000000_00000001_000003d1_00000000_00000000_00000000_00000000_00000000
00000001_000003d1_00000000_00000000_00000000_00000000_00000000_00000000
000003d1_00000000_00000000_00000000_00000000_00000000_00000001_000003d1Коэффициенты для сокращения 512-битного числа до 256-битного по модулю (2 ** 256 — 2 ** 32 — 977); размер слова — 64 бита; для удобства разряды сгруппированы по 64 бита:
>>> r = build_reducer(512, 256, 64, 2 ** 32 + 977)
>>> print_reducer(r, 256)
0000000000000000_0000000000000000_0000000000000000_0000000000000001
0000000000000000_0000000000000000_0000000000000001_0000000000000000
0000000000000000_0000000000000001_0000000000000000_0000000000000000
0000000000000001_0000000000000000_0000000000000000_0000000000000000
0000000000000000_0000000000000000_0000000000000000_00000001000003d1
0000000000000000_0000000000000000_00000001000003d1_0000000000000000
0000000000000000_00000001000003d1_0000000000000000_0000000000000000
00000001000003d1_0000000000000000_0000000000000000_0000000000000000Коэффициенты для сокращения 512-битного числа до 256-битного по модулю (2 ** 256 — 432420386565659656852420866394968145599); размер слова — 32 бита; для удобства разряды сгруппированы по 32 бита:
>>> r = build_reducer(512, 256, 32, 432420386565659656852420866394968145599)
>>> print_reducer(r, 256)
00000000_00000000_00000000_00000000_00000000_00000000_00000000_00000001
00000000_00000000_00000000_00000000_00000000_00000000_00000001_00000000
00000000_00000000_00000000_00000000_00000000_00000001_00000000_00000000
00000000_00000000_00000000_00000000_00000001_00000000_00000000_00000000
00000000_00000000_00000000_00000001_00000000_00000000_00000000_00000000
00000000_00000000_00000001_00000000_00000000_00000000_00000000_00000000
00000000_00000001_00000000_00000000_00000000_00000000_00000000_00000000
00000001_00000000_00000000_00000000_00000000_00000000_00000000_00000000
00000000_00000000_00000000_00000001_45512319_50b75fc4_402da173_2fc9bebf
00000000_00000000_00000001_45512319_50b75fc4_402da173_2fc9bebf_00000000
00000000_00000001_45512319_50b75fc4_402da173_2fc9bebf_00000000_00000000
00000001_45512319_50b75fc4_402da173_2fc9bebf_00000000_00000000_00000000
45512319_50b75fc4_402da173_2fc9bec0_45512319_50b75fc4_402da173_2fc9bebf
50b75fc4_402da173_2fc9bec0_9d671cd5_1b343a1b_66926b57_d2a4c1c6_1536bda7
402da173_2fc9bec0_9d671cd5_81c69bc5_9509b0b0_74ec0aea_8f564d66_7ec7eb3c
2fc9bec0_9d671cd5_81c69bc5_e697f5e4_1f12c33a_0a7b6f4e_3302b92e_a029cecdКоэффициенты для сокращения 512-битного числа до 256-битного по модулю (2 ** 256 — 432420386565659656852420866394968145599); размер слова — 64 бита; для удобства разряды сгруппированы по 64 бита:
>>> r = build_reducer(512, 256, 64, 432420386565659656852420866394968145599)
>>> print_reducer(r, 256)
0000000000000000_0000000000000000_0000000000000000_0000000000000001
0000000000000000_0000000000000000_0000000000000001_0000000000000000
0000000000000000_0000000000000001_0000000000000000_0000000000000000
0000000000000001_0000000000000000_0000000000000000_0000000000000000
0000000000000000_0000000000000001_4551231950b75fc4_402da1732fc9bebf
0000000000000001_4551231950b75fc4_402da1732fc9bebf_0000000000000000
4551231950b75fc4_402da1732fc9bec0_4551231950b75fc4_402da1732fc9bebf
402da1732fc9bec0_9d671cd581c69bc5_9509b0b074ec0aea_8f564d667ec7eb3cp.s. (2 ** 256 — 2 ** 32 — 977) — параметр p эллиптической кривой SECP256K1; (2 ** 256 — 432420386565659656852420866394968145599) — параметр N эллиптической кривой SECP256K1. (См. https://en.bitcoin.it/wiki/Secp256k1)
Что можно сказать, глядя на эти коэффициенты:
Чем меньше единичных бит в представлении числа omega — тем проще выглядят коэффициенты;
Чем меньше размер битового слова — тем быстрее происходит перенос в младшие разряды и тем самым сокращается длина битового представления x mod p;
Коэффициенты при младших битовых словах содержат по одному единичному биту на позициях, соответствующим началам битовых слов. Т.е. первые n бит числа берутся как есть, и к ним прибавляются группы старших битов, помноженные на вычисленные коэффициенты.
Игрушечные примеры
Для первого примера коэффициентов:
const uint32_t modulus = (1 << 8) - 17;
uint32_t mod_exact(uint32_t x) {
return x % modulus;
}
uint32_t mod_manual(uint32_t x) {
while (x & 0xFFFFFF00) {
uint32_t buffer = 0;
buffer += 0x01 * (x & 0xFF); x >>= 8;
buffer += 0x11 * (x & 0xFF); x >>= 8;
buffer += 0x32 * (x & 0xFF); x >>= 8;
buffer += 0x85 * (x & 0xFF); x >>= 8;
x = buffer;
}
return (x < modulus) ? x : (x - modulus);
}Для второго примера коэффициентов:
const uint32_t modulus = (1 << 16) - 666;
uint32_t mod_exact(uint32_t x) {
return x % modulus;
}
uint32_t mod_manual(uint32_t x) {
while (x & 0xFFFF0000) {
uint32_t buffer = (x & 0xFFFF); x >>= 16;
buffer += 0x029a * (x & 0xFF); x >>= 8;
buffer += 0x9f34 * (x & 0xFF); x >>= 8;
x = buffer;
}
return (x < modulus) ? x : (x - modulus);
}Программа для исчерпывающего тестирования двух примеров выше:
#include
// put some implementation of mod_exact() and mod_manual() here
int main() {
uint32_t number = 0, fails = 0;
do {
uint32_t exact = mod_exact(number);
uint32_t manual = mod_manual(number);
fails += ((exact == manual) ? 0 : 1);
if ((number & 0x00FFFFFF) == 0) {
std::cout << number << "; fails=" << fails << std::endl;
}
number++;
} while (number != 0);
std::cout << "done; fails=" << fails << std::endl;
return 0;
} Серьёзный пример: 97! mod (2 ** 256 — 2 ** 32 — 977)
Возьмем третий пример работы генератора коэффициентов:
3. Коэффициенты для сокращения 512-битного числа до 256-битного по модулю (2 ** 256 - 2 ** 32 - 977); размер слова - 32 бита; для удобства разряды сгруппированы по 32 бита:
>>> r = build_reducer(512, 256, 32, 2 ** 32 + 977)
>>> print_reducer(r, 256)
00000000_00000000_00000000_00000000_00000000_00000000_00000000_00000001
00000000_00000000_00000000_00000000_00000000_00000000_00000001_00000000
00000000_00000000_00000000_00000000_00000000_00000001_00000000_00000000
00000000_00000000_00000000_00000000_00000001_00000000_00000000_00000000
00000000_00000000_00000000_00000001_00000000_00000000_00000000_00000000
00000000_00000000_00000001_00000000_00000000_00000000_00000000_00000000
00000000_00000001_00000000_00000000_00000000_00000000_00000000_00000000
00000001_00000000_00000000_00000000_00000000_00000000_00000000_00000000
00000000_00000000_00000000_00000000_00000000_00000000_00000001_000003d1
00000000_00000000_00000000_00000000_00000000_00000001_000003d1_00000000
00000000_00000000_00000000_00000000_00000001_000003d1_00000000_00000000
00000000_00000000_00000000_00000001_000003d1_00000000_00000000_00000000
00000000_00000000_00000001_000003d1_00000000_00000000_00000000_00000000
00000000_00000001_000003d1_00000000_00000000_00000000_00000000_00000000
00000001_000003d1_00000000_00000000_00000000_00000000_00000000_00000000
000003d1_00000000_00000000_00000000_00000000_00000000_00000001_000003d1Входное 512-битное число разобьём на 16 групп по 32 бита, обозначим их w[0…15]. Тогда в результате умножения этих групп на ненулевые кусочки коэффициентов получим такое представление x mod p:
(2 ** 0) * (w[0] + 977 * w[8] + 977 * w[15]) +
(2 ** 32) * (w[1] + w[8] + 977 * w[9] + w[15]) +
(2 ** 64) * (w[2] + w[9] + 977 * w[10]) +
(2 ** 96) * (w[3] + w[10] + 977 * w[11]) +
(2 ** 128) * (w[4] + w[11] + 977 * w[12]) +
(2 ** 160) * (w[5] + w[12] + 977 * w[13]) +
(2 ** 192) * (w[6] + w[13] + 977 * w[14]) +
(2 ** 224) * (w[7] + w[14] + 977 * w[15])При каждой степени двойки (соответствующей одному битовому слову результата) стоит сумма чисел, которая наверняка должна укладываться в 64 бита, с учётом переноса в старшие разряды. Степени двойки (2 ** 256) и выше отсутствуют, т.е. с учётом возможного переноса результат будет состоять максимум из 9 битовых слов по 32 бита, а остальные будут равны нулю. Таким образом, для последующих итераций можно использовать упрощенные выражения:
(2 ** 0) * (w[0] + 977 * w[8]) +
(2 ** 32) * (w[1] + w[8]) +
(2 ** 64) * w[2] +
(2 ** 96) * w[3] +
(2 ** 128) * w[4] +
(2 ** 160) * w[5] +
(2 ** 192) * w[6] +
(2 ** 224) * w[7]Используя эти результаты напишем функцию для сокращения 512-битного числа до 256-битного по модулю (2 ** 256 — 2 ** 32 — 977):
void reduce_512(const uint32_t x[16], uint32_t y[8]) {
// check if any of bits[257..512] is set
uint32_t mask = 0;
for (int i = 8; i < 16; mask |= x[i++]);
if (mask == 0) return;
uint64_t buffer = 0;
// stage 1: reduce 16 limbs into 9
// (2 ** 0) * (w[0] + 977 * w[8] + 977 * w[15]) +
buffer += (uint64_t)x[0] + 977 * (uint64_t)x[8] + 977 * (uint64_t)x[15];
y[0] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 32) * (w[1] + w[8] + 977 * w[9] + w[15]) +
buffer += (uint64_t)x[1] + (uint64_t)x[8] + 977 * (uint64_t)x[9] + (uint64_t)x[15];
y[1] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 64) * (w[2] + w[9] + 977 * w[10]) +
buffer += (uint64_t)x[2] + (uint64_t)x[9] + 977 * (uint64_t)x[10];
y[2] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 96) * (w[3] + w[10] + 977 * w[11]) +
buffer += (uint64_t)x[3] + (uint64_t)x[10] + 977 * (uint64_t)x[11];
y[3] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 128) * (w[4] + w[11] + 977 * w[12]) +
buffer += (uint64_t)x[4] + (uint64_t)x[11] + 977 * (uint64_t)x[12];
y[4] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 160) * (w[5] + w[12] + 977 * w[13]) +
buffer += (uint64_t)x[5] + (uint64_t)x[12] + 977 * (uint64_t)x[13];
y[5] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 192) * (w[6] + w[13] + 977 * w[14]) +
buffer += (uint64_t)x[6] + (uint64_t)x[13] + 977 * (uint64_t)x[14];
y[6] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 224) * (w[7] + w[14] + 977 * w[15])
buffer += (uint64_t)x[7] + (uint64_t)x[14] + 977 * (uint64_t)x[15];
y[7] = buffer & 0xFFFFFFFF; buffer >>= 32;
// stage 2: reduce 9 limbs into 8
while (buffer != 0) {
// save 9th limb (10 bits max) and flush buffer
const uint64_t overflow256 = buffer & 0xFFFFFFFF; buffer >>= 32;
assert(buffer == 0);
// (2 ** 0) * (w[0] + 977 * w[8]) +
buffer += (uint64_t)y[0] + 977 * overflow256;
y[0] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 32)* (w[1] + w[8]) +
buffer += (uint64_t)y[1] + overflow256;
y[1] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 64) * w[2] +
buffer += (uint64_t)y[2];
y[2] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 96) * w[3] +
buffer += (uint64_t)y[3];
y[3] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 128) * w[4] +
buffer += (uint64_t)y[4];
y[4] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 160) * w[5] +
buffer += (uint64_t)y[5];
y[5] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 192) * w[6] +
buffer += (uint64_t)y[6];
y[6] = buffer & 0xFFFFFFFF; buffer >>= 32;
// (2 ** 224) * w[7]
buffer += (uint64_t)y[7];
y[7] = buffer & 0xFFFFFFFF; buffer >>= 32;
}
// 512 bits are now reduced to 256 bits, however the result may exceed
// (2^256 - 2^32 - 977) and an additional subtraction may be necessary
}Тестовая программа, в которую зашита пара чисел — факториал 97 (число размером около 500 бит) и правильный результат вычисления (97! mod (2 ** 256 — 2 ** 32 — 977)), необходимый для проверки результатов работы функции:
#include
#include
void reduce_512(const uint32_t x[16], uint32_t y[8]) {
// put function code here
}
int main() {
// math.factorial(97), this is about 500 bits long
// least significant limbs go first
const uint32_t x[16] = {
0x00000000, 0x00000000, 0xc0000000, 0xc63bc975,
0xcb0e1818, 0xfe74c03b, 0x13559f1a, 0xca00bb56,
0xef9d44bc, 0xf57bf161, 0xf3e3d5c3, 0xab918234,
0xb69daa20, 0x4532ed8b, 0xafb0a77f, 0x01d62e2f
};
// math.factorial(97) % (2 ** 256 - 2 ** 32 - 977)
// least significant limbs go first
const uint32_t y_expected[8] = {
0x7999b163, 0xcf77a9bd, 0x7dfec23d, 0x80718b50,
0x6655e0fc, 0xcc6efc90, 0xd9b7c95d, 0x7c17a6d2
};
uint32_t y[8];
reduce_512(x, y); // stage 2 is ran once for this input
// print/verify
bool ok = true;
for (int i = 0; i < 8; ok &= (y[i] == y_expected[i]), i++) {
std::cout << std::hex << y[i] << std::endl;
}
std::cout << (ok ? "Ok" : "Failed") << std::endl;
return 0;
} p.s. Все примеры проверены в Microsoft Visual Studio 2022 (v17.5.4)
