MWC - Генератор PRNG и хэш функция

Я какое-то время назад занимался идеей реализации генератора псевдослучайных чисел для задач молекулярной динамики, с хорошим распределением ошибки, чтобы не давал водяных знаков. Искал алгоритмы, которые можно использовать на GPU и допускающие параллельный расчет с нескольких позиций.

Нашел методы MWC64 и MWC128 (Multiply With Carry) основанные на умножении по модулю простых чисел Ризеля. Сейчас в контексте распределенных вычислений мы рассматриваем методы хеширования для быстрого поиска тензоров в распределенных вычислениях и в частности хэши пригодные для параллельного вычисления и проверки на GPU. Для этого мы выбираем хеш функции на модульной арифметике и простых числах.

Методы пост-квантовой криптографии требуют, чтобы хеш функция строилась с использованием модульной арифметики и NTT-Friendly модулей - это простые числа Прота (Proth primes): $q=A2^s+1$. Функция сжатия данных может быть основана на использовании MDS матриц перестановок, которые можно считать с использованием аналога быстрыго фурье преобразования на кольце полиномов (NTT). Использование нескольких простых чисел позволяет реализовать Redundancy RNS непозиционную систему исчисления, основанную на теореме об остатках (CRT). Т.е. нас ко всему прочему интересует расчет нескольких хешей в виде вектора по взаимно простым числам.

Algorithm 1. MWC reduction

---

*Require*: $p=K\cdot 2^{s}-1$, где $K\lt 2^{w}$, $\beta=2^{s}$
*Require*: input $x=x_1\beta + x_0 \lt p2^{s}$
*Ensure*: output $r \equiv x\beta^{-1} \bmod p$

1. $r ← x_1 + x_0\cdot K$
2. $r ← r - p \text{ if } r\geqslant p$

---

Метод основан на операции сдвига на $\beta$ на пол- слова и тождестве: $$ r\beta^{-1}\pmod{p} \equiv (r \bmod \beta)\cdot\lfloor (p+1)/\beta \rfloor + \lfloor r/\beta \rfloor, \text{ где } \beta=2^{w/2} \leqslant 2^s~. $$ $$ (r\bmod \beta) + \lfloor r/\beta \rfloor\cdot \beta \equiv r $$ $$ \beta\cdot a \equiv 1 \pmod{a\beta -1} $$

Algorithm 2. NTT-Friendly reduction

---

*Require*: $p=K\cdot 2^{s}+1$, где $K\lt 2^{w}$, $\beta=2^{s}$
*Require*: input $x=x_1\beta + x_0 \lt p2^{s}$
*Ensure*: output $r \equiv x\beta^{-1} \bmod p$

1. $r ← x_1 - x_0\cdot K$
2. $r ← r + p \text{ if } r\lt 0$

---

Алгоритм основан на тождестве $$ r\beta^{-1}\pmod{p} \equiv \lfloor r/\beta \rfloor - (r \bmod \beta)\cdot\lfloor (p-1)/\beta \rfloor, \text{ где } \beta=2^w \leqslant 2^s~. $$ $$ \beta\cdot a \equiv -1 \pmod{a\beta +1} $$

// Функция MWC хэша без знака 
uint32_t mwc32_hash(uint32_t h, uint16_t d, uint32_t q, uint32_t a){
    h = h + d;
    h = (uint16_t)h*a + (h>>16);
    return h;
}
// Функция MWC хэша со знаком NTT-friendly
int32_t mwc32s_hash(int32_t h, int16_t d, int32_t q, int32_t a){
    h = h + d;
    h = (h>>16) - (uint16_t)h*a;
    return h;
}
// Множество простых чисел Прота 128 шт:
0x7ffe0001, 0x7ff80001, 0x7fea0001, 0x7fe90001, 0x7fd70001, 0x7fd20001, 0x7fcb0001, 0x7fbd0001, 
0x7fb90001, 0x7fb40001, 0x7fad0001, 0x7f7b0001, 0x7f4e0001, 0x7f440001, 0x7f420001, 0x7f3c0001, 
0x7f330001, 0x7f210001, 0x7f180001, 0x7f000001, 0x7efc0001, 0x7ee20001, 0x7ecf0001, 0x7eba0001, 
0x7eb50001, 0x7eb40001, 0x7ea90001, 0x7e780001, 0x7e640001, 0x7e550001, 0x7e520001, 0x7e4b0001, 

0x7e370001, 0x7e360001, 0x7e250001, 0x7e100001, 0x7e0f0001, 0x7e040001, 0x7e000001, 0x7dfd0001, 
0x7df20001, 0x7de90001, 0x7de60001, 0x7de30001, 0x7dbf0001, 0x7dbe0001, 0x7db90001, 0x7db20001, 
0x7d9e0001, 0x7d980001, 0x7d970001, 0x7d8e0001, 0x7d6b0001, 0x7d560001, 0x7d4c0001, 0x7d3a0001, 
0x7d290001, 0x7d200001, 0x7d1f0001, 0x7d190001, 0x7d160001, 0x7d110001, 0x7d0d0001, 0x7d010001, 

0x7cfc0001, 0x7cf30001, 0x7cea0001, 0x7ce00001, 0x7cda0001, 0x7cd80001, 0x7cd50001, 0x7cd40001, 
0x7ccb0001, 0x7cb10001, 0x7cae0001, 0x7c9c0001, 0x7c980001, 0x7c800001, 0x7c5c0001, 0x7c480001, 
0x7c450001, 0x7c3b0001, 0x7c390001, 0x7c1d0001, 0x7c1a0001, 0x7c0e0001, 0x7c030001, 0x7c020001, 
0x7bff0001, 0x7bf30001, 0x7bd90001, 0x7bd50001, 0x7bd00001, 0x7bcd0001, 0x7bc10001, 0x7bb40001, 

0x7baf0001, 0x7b9f0001, 0x7b8b0001, 0x7b850001, 0x7b700001, 0x7b6c0001, 0x7b640001, 0x7b5a0001, 
0x7b3c0001, 0x7b3a0001, 0x7b340001, 0x7b310001, 0x7b270001, 0x7b250001, 0x7b070001, 0x7af40001, 
0x7ada0001, 0x7a8c0001, 0x7a880001, 0x7a860001, 0x7a7d0001, 0x7a770001, 0x7a650001, 0x7a550001, 
0x7a4f0001, 0x7a470001, 0x7a460001, 0x7a380001, 0x7a310001, 0x7a2b0001, 0x7a220001, 0x7a100001, 

Модульная арифметика со знаком для PQC

Наша задача реализация методов модульной арифметики для ZKP и PQC-пост-квантовой криптографии. Методов оказалось много. Раньше я ссылался на учебник [GECC] (сокр. Guide to Elliptic Curve Criptography), где разобраны разнообразные алгоритмы для криптографии, включая модульную арифметику и арифметику Галуа. В настоящее время появился ряд методов, к ним относятся, оптимизированные методы умножения Шоупа (Shoup) и Плантарда(Plantard, Thomas), а также новая тенденция в использовании модульной арифметики со знаком.

PQC - пост квантовая криптография основана на использование модульной арифметики в кольце полиномов с добавлением квантового шума Ring-LWE. Коэффициенты полиномов считаются в форме векторов длиной 256 элементов и выше.

ZKP - Zero-Knowledge Proof - концепция построения доказательства с нулевым разглашением, применяется например, для доказательства использования сложной функции без раскрытия информации о самой функции; для доказательства использования приватных данных без раскрытия данных. Применяется в распределенных вычислительных системах MPC - Multi-Party Computation.

Какой метод модульной арифметики может быть наиболее эффективным?
При выборе метода мы ориентируемся на систему команд целевой платформы: векторные инструкции CPU или GPU. И в настоящее время это должны быть 16/32/64 битные числа и вектора. Алгоритмы должны выражаться в векторных инструкциях сохраняющих разрядность элементов, таких как `mul_hi` - high product и `mul_lo` - low product, а также операции типа dot product с отложенным редуцированием по модулю могут дать прирост производительности. Операции могут быть с расширением знака или без знака. Заметим, что знаковая и беззнаковая операции умноженения идентичны signed low product и в целевой системе команд одна из операций может быть не представлена и выполняться через приведение типа. В системе команд могут быть представлены операции типа `madd` - умножения с накоплением и `msub` - умножение с вычитанием, в двух вариантах - старшая и младшая часть. Отдельно хочется упомянуть инструкции со странной разрядностью 52, которая рождается на месте мантиссы от double, в системе команд Intel AVX-IFMA.

Векторные инструкции не используют признаки переполнения и переноса разряда - это сильно осложняет реализацию методов для работы с большими числами произвольной разрядности. В этих случаях можно прибегнуть к непозиционной системе исчисления, основанной на CRT-китайской теореме об остатках, см. RNS-Residue Number System.

Оптимизированные алгоритмы опираются на сдвиг 16-, 32-, 52- или 64 бит и отложенное редуцирование Баррета при сложении и вычитании в ядре, при использовании модульных операций в цикле.

Algorithm 1.1. Signed Barrett modular reduction

---

*Require*: $0\lt q\lt\beta/2$, $-\beta q\leqslant z = z_1 \beta + z_0 \lt 2^{31}q$, где $-2^{31} \leqslant z_0 \lt 2^{31}$
*Require*: precomputed $u = \lfloor (2^{\lfloor \log(q)\rfloor -1} \beta)/q \rceil$
*Ensure*: Returns $r \equiv a \bmod^{\pm} q$ ; $-q\lt r \leqslant q$

1. $t ← \lfloor au/\beta \rfloor \gg 2^{\lfloor \log(q)\rfloor -1}$ -- signed low product and arithmetic shift
2. $t ← t \cdot q \bmod^{\pm} \beta$ -- signed low product
3. $r ← z - t$ -- signed remainder

---

Algorithm 2. Montgomery's modular multiplication

---

*Require*: $b,w \in [0,p)$; $p\lt\beta/2$
*Require*: precomputed $J = p^{-1} \mod \beta$, $w' = \beta w \mod p$
*Ensure*: Returns $r = b\cdot w \mod p$

1. $U = u_1\beta + u_0 ← bw'$ -- (wide product)
2. $q ← u_0J \mod \beta$ -- (low product)
3. $h ← \lfloor qp/\beta\rfloor$ -- (high product)
4. $r ← u_1 - h$
5. $r ← r + p$ if $r \lt 0$

---

последнюю строчку можно исключить, если принять $r \in (-p, p)$. Две последовательные операции _low product_ в строчке (1) и (2) можно объединить.

Algorithm 2.1. Signed Montgomery's modular multiplication

---

*Require*: $b,w \in [0,p)$; $p\lt\beta/2$
*Require*: precompute $J = p^{-1} \bmod^{\pm} \beta~$; $w' = w\beta \bmod p~$;
*Require*: precompute $\tilde{w} = w'\cdot J \bmod^{\pm} \beta$ -- (signed low product)
*Ensure*: Returns $r = b\cdot w \bmod^{\pm} p$

1. $u_1 ← \lfloor bw'/\beta \rfloor$ -- (signed high product)
2. $q ← b\tilde{w} \bmod^{\pm} \beta$ -- (signed low product)
3. $t ← \lfloor qp/\beta\rfloor$ -- (signed high product)
4. $r ← u_1 - t$

---

Для алгоритма важен метод вычисления обратного значения $p^{-1}$. Мы его вычисляем в 32 битных числах, потом приводим к знаковому типу $p^{-1} \bmod^{\pm} \beta$.

*Listing* Функция для вычисления $q^{-1} \bmod 2^{s}$, для $s\lt 32$:

uint32_t mod_inverse(uint32_t q, int s) {
    uint32_t q_inv = 1;
    for (int i = 1; i < s; i++) {
        if (((q * q_inv) & ((~0uL)>>(31-i))) != 1) {
            q_inv += (1u<<i);
        }
    }
    return q_inv;
}

Алгоритм -авторский, так и родился в виде программы. Базовая идея - инверсия по модулю должна давать единицу для любых s. Аналогичный метод может быть использован для нахождения константы при замене деления. Подобный метод возникал в конечных полях Галуа и все они, вероятно, могут быть выражены через расширенный алгоритм Евклида.

Система остаточных классов (RNS)

Китайская теорема об остатках (CRT)

Система остаточных классов (Residue Number System, RNS) является непозиционной системой целых чисел, основанной на китайской теореме об остатках (CRT). В такой системе целое число $x$ представляется его остатками $x_i = x \mod p_i$ по базису взаимно простых чисел $\mathcal{B} = \{p_0, \ldots, p_{k-1}\}$. Множество $\mathcal{B} = \{p_0, \ldots, p_{k-1}\}$ формирует базис RNS, состоящий из $k$ каналов. Модули $p_i$ обычно выбираются с учетом ширины слова $w$, которая соответствует целевой архитектуре. Важным преимуществом такой системы является то, что операции сложения, вычитания и умножения могут выполняться параллельно в каждом канале: $$ z_i = x_i \circ y_i \mod p_i, \text{ где } \circ \in \{+, -, \times\} $$ Традиционно рассматриваются системы $\{2^n-2^s\pm1, 2^n, 2^n-1\}$

Ряд работ по использованию RNS в доказательствах ZKP и FHE:

• [2016/510] A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes
• [2018/117] An Improved RNS Variant of the BFV Homomorphic Encryption Scheme
• [2018/931] A Full RNS Variant of Approximate Homomorphic Encryption

• Исходники на GitHub
• Полный текст включает описания алгоритмов расширения и сжатия базиса, преобразования из RNS в позиционную систему

Обозначения

Для целого числа $q \geq 2$ мы отождествляем кольцо $\mathbb{Z}_q$ с его отображением на симметричном интервале $\mathbb{Z} \cap [-q/2, q/2)$. Для произвольного действительного числа $x$ мы обозначаем через $[x]_q$ отображение $x$ на этот интервал, а именно, действительное число $x' \in [-q/2, q/2)$, такое что $x' - x$ делится на $q$. Мы также обозначаем через $\lfloor x \rfloor$, $\lceil x \rceil$ и $\lfloor x\rceil$ округление $x$ вниз, вверх и до ближайшего целого числа, соответственно. Векторы мы обозначаем жирным шрифтом, и расширяем нотации $\lfloor \mathbf{x} \rfloor$, $\lceil \mathbf{x} \rceil$, $\lfloor\mathbf{x}\rceil$ на векторы поэлементно.

Мы выбираем множество из $k$ взаимно простых модулей $\{p_0, \dots, p_{k-1}\}$, где все числа целые больше 1, и пусть их произведение равно $P = \prod_{i=0}^{k-1} p_i$.

Для всех $i \in \{0, \dots, k-1\}$ мы также обозначаем $$ \hat{p}_i = P / p_i \in \mathbb{Z} \quad \text{и} \quad \tilde{p}_i = \hat{p}_i^{-1} \pmod{q_i} \in \mathbb{Z}_{q_i}, $$ где $\tilde{p}_i \in [-p_i/2, p_i/2)$ и $\hat{p}_i \cdot \tilde{p}_i = 1 \pmod{p_i}$.

Теорема об остатках (CRT)

Обозначим представление целого числа $x \in \mathbb{Z}_p$ относительно базиса RNS $\{p_0, \dots, p_{k-1}\}$ через $x \sim (x_0, \dots, x_{k-1})$, где $x_i = [x]_{p_i} \in \mathbb{Z}_{p_i}$. Формула, выражающая $x$ через $x_i$, имеет вид $$ x = \sum_{i=0}^{k-1} x_i \cdot \tilde{p}_i \cdot \hat{p}_i \pmod{P}~. $$

Эта формула может быть использована более чем одним способом для «реконструкции» значения $x \in \mathbb{Z}_p$ из $[x]_\mathcal{B}$. В данной работе мы используем следующие два факта: $$ x = \sum_{i=0}^{k-1} [x_i \cdot \tilde{p}_i]_{p_i} \cdot \hat{p}_i - e \cdot P \quad \text{для некоторого } e \in \mathbb{Z}, $$ и $$ x = \sum_{i=0}^{k-1} x_i \cdot \tilde{p}_i \cdot \hat{p}_i - e' \cdot P \quad \text{для некоторого } e' \in \mathbb{Z}, $$ где сумма во второй формуле берётся по $x_i \cdot \tilde{q}_i \cdot \hat{q}_i \in \big[-\cfrac{q_i q}{4}, \cfrac{q_i q}{4}\big)$.

Представление RNS в кольце

Пусть $\mathcal{B} = \{p_0, \ldots, p_{k-1}\}$ — это базис взаимно простых чисел, и пусть $P = \prod_{i=0}^{k-1} p_i$. Обозначим через $[\cdot]_\mathcal{B}$ отображение из $\mathbb{Z}_P \mapsto \mathbb{Z}_{p_0} \times \mathbb{Z}_{p_1} \times \cdots \times \mathbb{Z}_{p_{k-1}}$, определённое как $a \mapsto [a]_\mathcal{B} = ([a]_{p_i})_{0 \leq i \le k}$ -- отображение из множества целых чисел на множество остатков в базисе взаимно простых чисел. Это изоморфизм кольца по Теореме об остатках (CRT), и $[a]_\mathcal{B}$ называется представлением числа $a \in \mathbb{Z}_P$ в системе остаточных классов (RNS). Основное преимущество представления RNS заключается в возможности выполнения компонентных арифметических операций в малых кольцах $\mathbb{Z}_{p_i}$, что снижает асимптотическую и практическую вычислительную сложность. Этот изоморфизм кольца над целыми числами может быть естественным образом расширен до изоморфизма в кольце полиномов $[ \cdot ]_\mathcal{B} : \mathcal{R}_P \to \mathcal{R}_{p_0} \times \cdots \times \mathcal{R}_{p_{k-1}}$ путём пересчета коэффициентов над циклическими кольцами.

Расширение базиса CRT

Пусть $x \in \mathbb{Z}_P$ задано в представлении CRT $(x_0, \dots, x_{k-1})$, и предположим, что мы хотим расширить базис CRT, вычислив $[x]_q \in \mathbb{Z}_q$ для некоторого другого модуля $q > 1$. Используя уравнение (2), мы хотели бы вычислить $$ [x]_q = \left[ \left( \sum_{i=0}^{k-1} [x_i \cdot \tilde{p}_i]_{p_i} \cdot \hat{p}_i \right) - e \cdot P \right]_q~. $$ Основная сложность здесь заключается в вычислении $e$, которое является целым числом в $\mathbb{Z}_k$. Формула для $e$ выглядит так: $$ e = \left\lfloor \frac{\sum_{i=0}^{k-1} [x_i \cdot \tilde{p}_i]_{p_i} \cdot \hat{p}_i}{P} \right\rceil = \left\lfloor \sum_{i=0}^{k-1} [x_i \cdot \tilde{p}_i]_{p_i} \cdot \frac{\hat{p}_i}{P} \right\rceil = \left\lfloor \sum_{i=0}^{k-1} \frac{[x_i \cdot \tilde{p}_i]_{p_i}}{p_i} \right\rceil. $$

Чтобы получить $e$, мы вычисляем для каждого $i \in \{0, \dots, k-1\}$ элемент $\xi_i := [x_i \cdot \tilde{p}_i]_{p_i}$ используя арифметику целых чисел, а затем рациональное число $z_i := \xi_i / p_i$ в формате с плавающей точкой одинарной точности. Затем суммируем все $z_i$ и округляем результат, чтобы получить $e$. {округление к меньшему для чисел без знака, проверить}: $$ e+{x\over P} = \sum_{i=0}^{k-1} \frac{\xi_i}{p_i}, \quad e = \left\lfloor\sum_{i=0}^{k-1} \frac{\xi_i}{p_i}\right\rfloor, $$ -- в такой форме должно быть справедливо для модульной арифметики без знака [23].

После того как мы получили значение $e$, а также все $\xi_i$, мы можем напрямую вычислить уравнение (2) по модулю $q$, чтобы получить $$ [x]_q = \left[ \left( \sum_{i=0}^{k-1} \xi_i \cdot [\hat{p}_i]_q \right) - e \cdot [P]_q \right]_q~. $$ Заметим, если все значения $[\hat{p}_i]_q$ представить в качестве элементов вектора, то вычисление сводится к операции скалярного произведения векторов размерности $k$ по модулю $q$.