콘텐츠로 이동

모듈러 곱셈 역원

정의

정수 $a$의 모듈러 곱셈 역원은 어떤 법 $m$에 대해 $a \cdot x$가 $1$과 합동이 되게 하는 정수 $x$입니다. 이를 형식적으로 쓰면, 다음을 만족하는 정수 $x$를 찾고자 하는 것입니다.

$$a \cdot x \equiv 1 \mod m.$$

또한 $x$를 간단히 $a^{-1}$로 표기하겠습니다.

모듈러 역원이 항상 존재하는 것은 아니라는 점에 유의해야 합니다. 예를 들어, $m = 4$, $a = 2$라고 합시다. $m$을 법으로 하는 가능한 모든 값을 확인해 보면, 위 방정식을 만족하는 $a^{-1}$을 찾을 수 없다는 것이 분명해질 것입니다. 모듈러 역원은 $a$와 $m$이 서로소일 때(즉, $\gcd(a, m) = 1$일 때), 그리고 오직 그때에만 존재한다는 것을 증명할 수 있습니다.

이 글에서는 모듈러 역원이 존재할 때 이를 구하는 두 가지 방법과, 모든 수에 대한 모듈러 역원을 선형 시간에 구하는 한 가지 방법을 소개합니다.

확장 유클리드 알고리즘을 이용해 모듈러 역원 구하기

미지수 $x$, $y$에 대한 다음 방정식을 생각해 봅시다:

$$a \cdot x + m \cdot y = 1$$

이것은 이변수 선형 디오판토스 방정식입니다. 링크된 글에서 보인 것처럼, $\gcd(a, m) = 1$일 때 이 방정식은 해를 가지며, 그 해는 확장 유클리드 알고리즘으로 구할 수 있습니다. $\gcd(a, m) = 1$은 모듈러 역원이 존재하기 위한 조건이기도 하다는 점에 유의하십시오.

이제 양변에 $m$을 법으로 취하면 $m \cdot y$를 없앨 수 있으며, 방정식은 다음과 같이 됩니다:

$$a \cdot x \equiv 1 \mod m$$

따라서 $a$의 모듈러 역원은 $x$입니다.

구현은 다음과 같습니다.

int x, y;
int g = extended_euclidean(a, m, x, y);
if (g != 1) {
    cout << "No solution!";
}
else {
    x = (x % m + m) % m;
    cout << x << endl;
}

x를 수정하는 방식에 주목하십시오. 확장 유클리드 알고리즘으로 얻은 x는 음수일 수 있으므로 x % m 역시 음수일 수 있으며, 따라서 먼저 m을 더해 양수로 만들어야 합니다.

이진 거듭제곱을 이용해 모듈러 역원 구하기

모듈러 역원을 구하는 또 다른 방법은 오일러 정리를 이용하는 것입니다. 오일러 정리에 따르면 $a$와 $m$이 서로소일 때 다음 합동식이 성립합니다.

$$a^{\phi (m)} \equiv 1 \mod m$$

$\phi$는 오일러 피 함수입니다. 여기서도 $a$와 $m$이 서로소라는 것이 모듈러 역원이 존재하기 위한 조건이었다는 점에 유의하십시오.

$m$이 소수이면 이는 페르마의 소정리로 단순화됩니다.

$$a^{m - 1} \equiv 1 \mod m$$

위 방정식의 양변에 $a^{-1}$을 곱하면 다음을 얻습니다:

  • 임의의 (단, 서로소인) 법 $m$에 대해: $a ^ {\phi (m) - 1} \equiv a ^{-1} \mod m$
  • 소수인 법 $m$에 대해: $a ^ {m - 2} \equiv a ^ {-1} \mod m$

이 결과로부터 $O(\log m)$ 시간에 동작하는 이진 거듭제곱 알고리즘을 이용하여 모듈러 역원을 쉽게 구할 수 있습니다.

이 방법은 앞 문단에서 설명한 방법보다 이해하기 쉽지만, $m$이 소수가 아닌 경우에는 오일러 피 함수를 계산해야 하며, 이를 위해서는 $m$을 소인수분해해야 하는데 이는 매우 어려울 수 있습니다. $m$의 소인수분해를 알고 있다면 이 방법의 복잡도는 $O(\log m)$입니다.

유클리드 나눗셈을 이용해 소수 법에 대한 모듈러 역원 구하기

소수 모듈러스 $m > a$가 주어졌다고 합시다(그렇지 않더라도 모듈러 연산을 한 번 적용하여 더 작게 만들 수 있습니다). 유클리드 나눗셈에 따르면

$$m = k \cdot a + r$$

여기서 $k = \left\lfloor \frac{m}{a} \right\rfloor$이고 $r = m \bmod a$이며, 그러면 다음이 성립합니다.

$$ \begin{align*} & \implies & 0 & \equiv k \cdot a + r & \mod m \\ & \iff & r & \equiv -k \cdot a & \mod m \\ & \iff & r \cdot a^{-1} & \equiv -k & \mod m \\ & \iff & a^{-1} & \equiv -k \cdot r^{-1} & \mod m \end{align*} $$

$m$이 소수가 아니면 이 추론이 성립하지 않는다는 점에 유의하십시오. 일반적인 경우에는 $a^{-1}$이 존재한다고 해서 $r^{-1}$도 존재하는 것은 아니기 때문입니다. 이를 확인하기 위해 위 공식으로 $12$를 법으로 하는 $5^{-1}$을 계산해 봅시다. $5 \cdot 5 \equiv 1 \bmod 12$이므로 결과로 $5$를 얻어야 합니다. 그러나 $12 = 2 \cdot 5 + 2$이므로 $k=2$, $r=2$가 되는데, $2$는 $12$를 법으로 역원을 갖지 않습니다.

하지만 모듈러스가 소수라면 $0 < a < m$인 모든 $a$는 $m$을 법으로 역원을 가지며, $m$에 대한 수 $a$의 모듈러 역원을 계산하는 다음과 같은 재귀 함수(C++)를 작성할 수 있습니다.

int inv(int a) {
  return a <= 1 ? a : m - (long long)(m/a) * inv(m % a) % m;
}

이 재귀의 정확한 시간 복잡도는 알려져 있지 않습니다. 시간 복잡도는 $O(\frac{\log m}{\log\log m})$와 $O(m^{\frac{1}{3} - \frac{2}{177} + \epsilon})$ 사이 어딘가에 있습니다. On the length of Pierce expansions를 참고하십시오. 실제로 이 구현은 빠릅니다. 예를 들어 모듈러스가 $10^9 + 7$이면 항상 50번 미만의 반복 안에 끝납니다.

이 공식을 적용하면 $[1, m-1]$ 범위의 모든 수에 대한 모듈러 역원을 $O(m)$에 미리 계산할 수도 있습니다.

inv[1] = 1;
for(int a = 2; a < m; ++a)
    inv[a] = m - (long long)(m/a) * inv[m%a] % m;

$m$을 법으로 하는 수 배열의 모듈러 역원 구하기

배열이 주어지고, 그 안의 모든 수에 대해 모듈러 역원을 구하려 한다고 가정합시다(모든 수는 역원을 가집니다). 각 수마다 역원을 계산하는 대신, 분수에 접두사 곱(자기 자신 제외)과 접미사 곱(자기 자신 제외)을 곱해 확장하면 결국 역원을 단 한 번만 계산하면 됩니다.

$$ \begin{align} x_i^{-1} &= \frac{1}{x_i} = \frac{\overbrace{x_1 \cdot x_2 \cdots x_{i-1}}^{\text{prefix}_{i-1}} \cdot ~1~ \cdot \overbrace{x_{i+1} \cdot x_{i+2} \cdots x_n}^{\text{suffix}_{i+1}}}{x_1 \cdot x_2 \cdots x_{i-1} \cdot x_i \cdot x_{i+1} \cdot x_{i+2} \cdots x_n} \\ &= \text{prefix}_{i-1} \cdot \text{suffix}_{i+1} \cdot \left(x_1 \cdot x_2 \cdots x_n\right)^{-1} \end{align} $$

코드에서는 누적 곱 배열(자기 자신은 제외하고, 항등 원소에서 시작)을 만들고, 모든 수의 곱에 대한 모듈러 역원을 계산한 다음, 여기에 누적 곱과 접미사 곱(자기 자신은 제외)을 곱하기만 하면 됩니다. 접미사 곱은 뒤에서 앞으로 순회하며 계산합니다.

std::vector<int> invs(const std::vector<int> &a, int m) {
    int n = a.size();
    if (n == 0) return {};
    std::vector<int> b(n);
    int v = 1;
    for (int i = 0; i != n; ++i) {
        b[i] = v;
        v = static_cast<long long>(v) * a[i] % m;
    }
    int x, y;
    extended_euclidean(v, m, x, y);
    x = (x % m + m) % m;
    for (int i = n - 1; i >= 0; --i) {
        b[i] = static_cast<long long>(x) * b[i] % m;
        x = static_cast<long long>(x) * a[i] % m;
    }
    return b;
}

연습 문제