몽고메리 곱셈¶

소수 판별이나 소인수분해 같은 정수론의 많은 알고리즘과 RSA 같은 암호학 알고리즘에서는 큰 수를 법으로 하는 연산이 많이 필요합니다. $x y \bmod{n}$과 같은 곱셈은 곱에서 $n$을 몇 번 빼야 하는지 알기 위해 나눗셈이 필요하므로, 일반적인 알고리즘으로 계산하면 상당히 느립니다. 그리고 나눗셈은 특히 큰 수에서 매우 비용이 큰 연산입니다.

몽고메리 (모듈러) 곱셈은 이러한 곱셈을 더 빠르게 계산할 수 있게 해 주는 방법입니다. 곱을 나누고 $n$을 여러 번 빼는 대신, $n$의 배수를 더해 하위 비트를 상쇄한 다음 하위 비트를 그냥 버립니다.

몽고메리 표현¶

하지만 몽고메리 곱셈을 공짜로 쓸 수 있는 것은 아닙니다. 이 알고리즘은 몽고메리 공간에서만 동작합니다. 그리고 곱셈을 시작하기 전에 수를 그 공간으로 변환해야 합니다.

이 공간을 위해서는 $n$과 서로소인, 즉 $\gcd(n, r) = 1$인 양의 정수 $r \ge n$이 필요합니다. 실제로는 항상 양의 정수 $m$에 대해 $r$을 $2^m$으로 선택하는데, 그러면 곱셈, 나눗셈, $r$을 법으로 하는 연산을 시프트와 기타 비트 연산으로 효율적으로 구현할 수 있기 때문입니다. 짝수는 소인수분해하기 어렵지 않으므로, 거의 모든 응용에서 $n$은 홀수입니다. 따라서 $2$의 모든 거듭제곱은 $n$과 서로소입니다.

몽고메리 공간에서 수 $x$의 대표값 $\bar{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$이므로, $0 < r^{-1}, n^{\prime} < n$이면서 다음을 만족하는 두 수 $r^{-1}$과 $n^{\prime}$이 존재함을 알 수 있습니다

$$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$에 대해 성립합니다. 이는 $x \cdot n^{\prime}$에 $r$의 임의의 배수를 더하거나 뺄 수 있다는 뜻이며, 다시 말해 $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}$$

따라서 $2^1$을 법으로 하는 $a$의 역원으로 $x = 1$에서 시작하여 이 기법을 몇 번 적용하면 되며, 반복할 때마다 $x$의 올바른 비트 수가 두 배로 늘어납니다.

구현¶

GCC 컴파일러는 __int128과 __uint128 타입으로 128비트 정수를 지원하므로, 세 수가 모두 64비트 정수일 때에도 $x \cdot y \bmod n$을 여전히 효율적으로 계산할 수 있습니다.

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

그러나 256비트 정수를 위한 타입은 없습니다. 따라서 여기서는 128비트 곱셈에 대한 구현을 보이겠습니다.

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;
};

빠른 변환¶

현재의 방법으로 수를 몽고메리 공간으로 변환하는 것은 꽤 느립니다. 더 빠른 방법이 있습니다.

다음 관계에 주목할 수 있습니다:

$$\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;
};