lea 23 апр 2023 в 19:07

Быстрое нахождениe остатка от деления больших чисел для делителей специального вида

Средний

11 мин

Туториал

В этой статье я расскажу об одном способе вычисления 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_8f564d667ec7eb3c

p.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 <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)

Теги:

Хабы:

Быстрое нахождениe остатка от деления больших чисел для делителей специального вида

Зачем это нужно?

Теория

Генератор коэффициентов на Python

Игрушечные примеры

Серьёзный пример: 97! mod (2 256 - 2 32 - 977)

Публикации

Истории

Работа

Ближайшие события

Быстрое нахождениe остатка от деления больших чисел для делителей специального вида

Зачем это нужно?

Теория

Генератор коэффициентов на Python

Игрушечные примеры

Серьёзный пример: 97! mod (2 ** 256 - 2 ** 32 - 977)

Серьёзный пример: 97! mod (2 256 - 2 32 - 977)