Множення Монтгомері

Багато алгоритмів теорії чисел, як-от перевірка на простоту чи факторизація цілих чисел, а також алгоритми криптографії, на кшталт RSA, потребують безлічі операцій за модулем великого числа. Множення на кшталт $x y \bmod{n}$ обчислюється доволі повільно за допомогою звичних алгоритмів, оскільки воно вимагає ділення, щоб дізнатися, скільки разів треба відняти $n$ від добутку. А ділення — це справді дорога операція, особливо коли йдеться про великі числа.

Множення (за модулем) Монтгомері — це метод, який дозволяє обчислювати такі множення швидше. Замість того, щоб ділити добуток і багаторазово віднімати $n$, він додає кратні $n$, щоб занулити молодші біти, а потім просто відкидає ці молодші біти.

Коли підходить цей алгоритм?

• Чи виконується багато множень за тим самим модулем $n$ (наприклад, у тестах простоти чи факторизації)? Лише так одноразова вартість переходу в простір Монтгомері окупається.
• Чи модуль $n$ непарний, тобто $\gcd(n, 2^m) = 1$? Для парного $n$ вибір $r = 2^m$ не працює.
• Чи вузьке місце саме в операції за модулем (ділення), яку хочеться замінити зсувами? Якщо множень мало, простіше обійтися звичайним %.

Представлення Монтгомері

Однак множення Монтгомері не дається задарма. Алгоритм працює лише в просторі Монтгомері. І нам потрібно перетворити наші числа в цей простір, перш ніж починати множити.

Для цього простору нам потрібне додатне ціле число $r \ge n$, взаємно просте з $n$, тобто $\gcd(n, r) = 1$. На практиці ми завжди вибираємо $r$ рівним $2^m$ для додатного цілого $m$, оскільки тоді множення, ділення та операції за модулем $r$ можна ефективно реалізувати за допомогою зсувів та інших бітових операцій. $n$ майже в усіх застосуваннях буде непарним числом, бо парне число факторизувати неважко. Тож кожен степінь $2$ буде взаємно простим з $n$.

Представник $\bar{x}$ числа $x$ у просторі Монтгомері означається як:

$$\bar{x} := x \cdot r \bmod n$$

Зауважте, що це перетворення насправді і є тим самим множенням, яке ми хочемо оптимізувати. Тож це все ще дорога операція. Однак перетворити число в простір потрібно лише один раз. Щойно ви опинилися в просторі Монтгомері, ви можете ефективно виконувати скільки завгодно операцій. А наприкінці ви перетворюєте остаточний результат назад. Тож поки ви виконуєте багато операцій за модулем $n$, це не буде проблемою.

Усередині простору Монтгомері ви все ще можете виконувати більшість операцій як звичайно. Ви можете додавати два елементи ($x \cdot r + y \cdot r \equiv (x + y) \cdot r \bmod n$), віднімати, перевіряти рівність і навіть обчислювати найбільший спільний дільник числа з $n$ (оскільки $\gcd(n, r) = 1$). Усе це за допомогою звичних алгоритмів.

Однак для множення це не так.

Ми очікуємо, що результатом буде:

$$\bar{x} * \bar{y} = \overline{x \cdot y} = (x \cdot y) \cdot r \bmod n.$$

Але звичайне множення дасть нам:

$$\bar{x} \cdot \bar{y} = (x \cdot y) \cdot r \cdot r \bmod n.$$

Тому множення в просторі Монтгомері означається так:

$$\bar{x} * \bar{y} := \bar{x} \cdot \bar{y} \cdot r^{-1} \bmod n.$$

Редукція Монтгомері

Множення двох чисел у просторі Монтгомері вимагає ефективного обчислення $x \cdot r^{-1} \bmod n$. Ця операція називається редукцією Монтгомері і також відома як алгоритм REDC.

Оскільки $\gcd(n, r) = 1$, ми знаємо, що існують два числа $r^{-1}$ та $n^{\prime}$ з $0 < r^{-1}, n^{\prime} < n$ такі, що

$$r \cdot r^{-1} + n \cdot n^{\prime} = 1.$$

Обидва числа $r^{-1}$ та $n^{\prime}$ можна обчислити за допомогою розширеного алгоритму Евкліда.

Користуючись цією тотожністю, ми можемо записати $x \cdot r^{-1}$ як:

$$\begin{aligned} x \cdot r^{-1} &= x \cdot r \cdot r^{-1} / r = x \cdot (-n \cdot n^{\prime} + 1) / r \\ &= (-x \cdot n \cdot n^{\prime} + x) / r \equiv (-x \cdot n \cdot n^{\prime} + l \cdot r \cdot n + x) / r \bmod n\\ &\equiv ((-x \cdot n^{\prime} + l \cdot r) \cdot n + x) / r \bmod n \end{aligned}$$

Ці еквівалентності виконуються для будь-якого цілого $l$. Це означає, що ми можемо додавати чи віднімати довільне кратне $r$ до $x \cdot n^{\prime}$, або, іншими словами, ми можемо обчислити $q := x \cdot n^{\prime}$ за модулем $r$.

Це дає нам такий алгоритм для обчислення $x \cdot r^{-1} \bmod n$:

function reduce(x):
    q = (x mod r) * n' mod r
    a = (x - q * n) / r
    if a < 0:
        a += n
    return a

Оскільки $x < n \cdot n < r \cdot n$ (навіть якщо $x$ є добутком множення) і $q \cdot n < r \cdot n$, ми знаємо, що $-n < (x - q \cdot n) / r < n$. Тому остаточна операція за модулем реалізується за допомогою однієї перевірки та одного додавання.

Як ми бачимо, ми можемо виконати редукцію Монтгомері без жодних важких операцій за модулем. Якщо ми оберемо $r$ як степінь $2$, то операції за модулем та ділення в алгоритмі можна обчислити за допомогою бітових масок і зсувів.

Другим застосуванням редукції Монтгомері є перенесення числа назад із простору Монтгомері у звичайний простір.

Трюк зі швидким оберненням

Щоб ефективно обчислити обернений елемент $n^{\prime} := n^{-1} \bmod r$, ми можемо скористатися таким трюком (натхненним методом Ньютона):

$$a \cdot x \equiv 1 \bmod 2^k \Longrightarrow a \cdot x \cdot (2 - a \cdot x) \equiv 1 \bmod 2^{2k}$$

Це легко довести. Якщо ми маємо $a \cdot x = 1 + m \cdot 2^k$, то маємо:

$$\begin{aligned} a \cdot x \cdot (2 - a \cdot x) &= 2 \cdot a \cdot x - (a \cdot x)^2 \\ &= 2 \cdot (1 + m \cdot 2^k) - (1 + m \cdot 2^k)^2 \\ &= 2 + 2 \cdot m \cdot 2^k - 1 - 2 \cdot m \cdot 2^k - m^2 \cdot 2^{2k} \\ &= 1 - m^2 \cdot 2^{2k} \\ &\equiv 1 \bmod 2^{2k}. \end{aligned}$$

Це означає, що ми можемо почати з $x = 1$ як оберненого до $a$ за модулем $2^1$, застосувати трюк кілька разів, і на кожній ітерації ми подвоюватимемо кількість правильних бітів числа $x$.

Реалізація

За допомогою компілятора GCC ми все ще можемо ефективно обчислити $x \cdot y \bmod n$, коли всі три числа є 64-бітними цілими, оскільки компілятор підтримує 128-бітні цілі через типи __int128 та __uint128.

C++

long long result = (__int128)x * y % n;

Python

# Python має нативні цілі довільної точності, тож проміжний добуток
# x * y не переповнюється, і додатковий 128-бітний тип не потрібен.
result = x * y % n

TypeScript

// BigInt має довільну точність, тож x * y не переповнюється;
// усі операнди мають бути bigint (літерали з суфіксом n).
const result: bigint = (x * y) % n;

Go

// math/bits.Mul64 повертає повний 128-бітний добуток (hi, lo) двох
// 64-бітних чисел, а bits.Div64 ділить 128-бітне (hi, lo) на n.
// Так ми відтворюємо (__int128)x * y % n без типу __int128.
hi, lo := bits.Mul64(x, y)
_, result := bits.Div64(hi%n, lo, n)

Однак для 256-бітних цілих типу немає. Тому тут ми покажемо реалізацію для 128-бітного множення.

C++

using u64 = uint64_t;
using u128 = __uint128_t;
using i128 = __int128_t;

struct u256 {
    u128 high, low;

    static u256 mult(u128 x, u128 y) {
        u64 a = x >> 64, b = x;
        u64 c = y >> 64, d = y;
        // (a*2^64 + b) * (c*2^64 + d) =
        // (a*c) * 2^128 + (a*d + b*c)*2^64 + (b*d)
        u128 ac = (u128)a * c;
        u128 ad = (u128)a * d;
        u128 bc = (u128)b * c;
        u128 bd = (u128)b * d;
        u128 carry = (u128)(u64)ad + (u128)(u64)bc + (bd >> 64u);
        u128 high = ac + (ad >> 64u) + (bc >> 64u) + (carry >> 64u);
        u128 low = (ad << 64u) + (bc << 64u) + bd;
        return {high, low};
    }
};

struct Montgomery {
    Montgomery(u128 n) : mod(n), inv(1) {
        for (int i = 0; i < 7; i++)
            inv *= 2 - n * inv;
    }

    u128 init(u128 x) {
        x %= mod;
        for (int i = 0; i < 128; i++) {
            x <<= 1;
            if (x >= mod)
                x -= mod;
        }
        return x;
    }

    u128 reduce(u256 x) {
        u128 q = x.low * inv;
        i128 a = x.high - u256::mult(q, mod).high;
        if (a < 0)
            a += mod;
        return a;
    }

    u128 mult(u128 a, u128 b) {
        return reduce(u256::mult(a, b));
    }

    u128 mod, inv;
};

Python

# Python має цілі довільної точності, тож окремий тип u256 не потрібен:
# просто маскуємо результат до 128 та 256 бітів, де C++ покладається на
# переповнення u128/u256.
MASK128 = (1 << 128) - 1


class Montgomery:
    def __init__(self, n: int) -> None:
        self.mod = n
        self.inv = 1
        # обернене n^{-1} mod 2^128 трюком Ньютона (7 ітерацій)
        for _ in range(7):
            self.inv = (self.inv * (2 - n * self.inv)) & MASK128

    def init(self, x: int) -> int:
        x %= self.mod
        for _ in range(128):
            x <<= 1
            if x >= self.mod:
                x -= self.mod
        return x

    def reduce(self, x: int) -> int:
        # x — повне 256-бітне число; беремо молодші/старші 128 бітів
        x_low = x & MASK128
        x_high = x >> 128
        q = (x_low * self.inv) & MASK128
        # старші 128 бітів добутку q * mod
        a = x_high - ((q * self.mod) >> 128)
        if a < 0:
            a += self.mod
        return a

    def mult(self, a: int, b: int) -> int:
        return self.reduce(a * b)

TypeScript

// BigInt дає довільну точність, тож окремий тип u256 непотрібний:
// маскуємо до 128 та 256 бітів там, де C++ спирається на переповнення.
const MASK128 = (1n << 128n) - 1n;

class Montgomery {
  mod: bigint;
  inv: bigint;

  constructor(n: bigint) {
    this.mod = n;
    this.inv = 1n;
    // обернене n^{-1} mod 2^128 трюком Ньютона (7 ітерацій)
    for (let i = 0; i < 7; i++) {
      this.inv = (this.inv * (2n - n * this.inv)) & MASK128;
    }
  }

  init(x: bigint): bigint {
    x %= this.mod;
    for (let i = 0; i < 128; i++) {
      x <<= 1n;
      if (x >= this.mod) x -= this.mod;
    }
    return x;
  }

  reduce(x: bigint): bigint {
    // x — повне 256-бітне число; розділяємо на молодші/старші 128 бітів
    const xLow = x & MASK128;
    const xHigh = x >> 128n;
    const q = (xLow * this.inv) & MASK128;
    // старші 128 бітів добутку q * mod
    let a = xHigh - ((q * this.mod) >> 128n);
    if (a < 0n) a += this.mod;
    return a;
  }

  mult(a: bigint, b: bigint): bigint {
    return this.reduce(a * b);
  }
}

Go

// Go не має ані 128-бітних, ані 256-бітних цілих, а ручний u256 на
// парах uint64 був би громіздким. Тому використовуємо math/big: метод
// big.Int дає довільну точність, а маскування до 128/256 бітів
// відтворює переповнення u128/u256 з оригіналу.
import (
	"math/big"
)

var mask128 = new(big.Int).Sub(new(big.Int).Lsh(big.NewInt(1), 128), big.NewInt(1))

type Montgomery struct {
	mod *big.Int
	inv *big.Int
}

func NewMontgomery(n *big.Int) *Montgomery {
	m := &Montgomery{mod: new(big.Int).Set(n), inv: big.NewInt(1)}
	// обернене n^{-1} mod 2^128 трюком Ньютона (7 ітерацій)
	for i := 0; i < 7; i++ {
		t := new(big.Int).Mul(n, m.inv)         // n * inv
		t.Sub(big.NewInt(2), t)                 // 2 - n * inv
		m.inv.Mul(m.inv, t).And(m.inv, mask128) // inv * (2 - n*inv) mod 2^128
	}
	return m
}

func (m *Montgomery) Init(x *big.Int) *big.Int {
	x = new(big.Int).Mod(x, m.mod)
	for i := 0; i < 128; i++ {
		x.Lsh(x, 1)
		if x.Cmp(m.mod) >= 0 {
			x.Sub(x, m.mod)
		}
	}
	return x
}

func (m *Montgomery) Reduce(x *big.Int) *big.Int {
	// x — повне 256-бітне число; розділяємо на молодші/старші 128 бітів
	xLow := new(big.Int).And(x, mask128)
	xHigh := new(big.Int).Rsh(x, 128)
	q := new(big.Int).Mul(xLow, m.inv)
	q.And(q, mask128)
	// старші 128 бітів добутку q * mod
	prodHigh := new(big.Int).Mul(q, m.mod)
	prodHigh.Rsh(prodHigh, 128)
	a := new(big.Int).Sub(xHigh, prodHigh)
	if a.Sign() < 0 {
		a.Add(a, m.mod)
	}
	return a
}

func (m *Montgomery) Mult(a, b *big.Int) *big.Int {
	return m.Reduce(new(big.Int).Mul(a, b))
}

Швидке перетворення

Поточний метод перетворення числа в простір Монтгомері доволі повільний. Є швидші способи.

Можна помітити таке співвідношення:

$$\bar{x} := x \cdot r \bmod n = x \cdot r^2 / r = x * r^2$$

Перетворення числа в простір — це просто множення всередині простору цього числа на $r^2$. Тому ми можемо заздалегідь обчислити $r^2 \bmod n$ і просто виконати одне множення замість того, щоб зсувати число 128 разів.

У наведеному нижче коді ми ініціалізуємо r2 значенням -n % n, що еквівалентно $r - n \equiv r \bmod n$, і зсуваємо його 4 рази, щоб отримати $r \cdot 2^4 \bmod n$. Це число можна інтерпретувати як $2^4$ у просторі Монтгомері. Якщо ми піднесемо його до квадрата $5$ разів, то отримаємо $(2^4)^{2^5} = (2^4)^{32} = 2^{128} = r$ у просторі Монтгомері, що є саме $r^2 \bmod n$.

struct Montgomery {
    Montgomery(u128 n) : mod(n), inv(1), r2(-n % n) {
        for (int i = 0; i < 7; i++)
            inv *= 2 - n * inv;

        for (int i = 0; i < 4; i++) {
            r2 <<= 1;
            if (r2 >= mod)
                r2 -= mod;
        }
        for (int i = 0; i < 5; i++)
            r2 = mul(r2, r2);
    }

    u128 init(u128 x) {
        return mult(x, r2);
    }

    u128 mod, inv, r2;
};