Быстрое нахождени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).

Принцип работы:

1. Определяем количество битовых слов размера s, необходимых для представления числа из m бит;

2. Находим степень двойки (2 ** n), по границе которой будем разбивать каждый коэффициент на «младшие» (1..n) и «старшие» (n+1..m) биты;

3. Находим степени двойки, являющиеся коэффициентами битовых слов в разложении делимого;

4. Пока хотя бы один из коэффициентов превышает (2 ** n) — дробим такие коэффициенты на «младшие» и «старшие» биты и пересобираем из них новый коэффициент: («младшие биты» + «старшие биты» * omega).

Примеры работы генератора коэффициентов:

1. Коэффициенты для сокращения 32-битного числа до 8-битного по модулю (2 ** 8 — 17) = 239; размер слова — 8 бит:

>>> r = build_reducer(32, 8, 8, 17) >>> print_reducer(r, 8) 01 11 32 85

2. Коэффициенты для сокращения 32-битного числа до 16-битного по модулю (2 ** 16 — 666) = 64870; размер слова — 8 бит:

>>> r = build_reducer(32, 16, 8, 666) >>> print_reducer(r, 16) 0001 0100 029a 9f34

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

4. Коэффициенты для сокращения 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

5. Коэффициенты для сокращения 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

6. Коэффициенты для сокращения 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_8f564d667ec7eb3c

p.s. (2 ** 256 — 2 ** 32 — 977) — параметр p эллиптической кривой SECP256K1; (2 ** 256 — 432420386565659656852420866394968145599) — параметр N эллиптической кривой SECP256K1. (См. https://en.bitcoin.it/wiki/Secp256k1)

Что можно сказать, глядя на эти коэффициенты:

1. Чем меньше единичных бит в представлении числа omega — тем проще выглядят коэффициенты;

2. Чем меньше размер битового слова — тем быстрее происходит перенос в младшие разряды и тем самым сокращается длина битового представления x mod p;

3. Коэффициенты при младших битовых словах содержат по одному единичному биту на позициях, соответствующим началам битовых слов. Т.е. первые 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 <iostream>  // 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 <iostream> #include <cassert>  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)