Pythonで巨大なべき乗のmodを高速に計算する：((a^b)^(c^d)) mod n を求めるプログラム

5つの整数 a, b, c, d, n が与えられたとき、((a^b)^(c^d)) mod n の値を求めることを考えます。出力は整数値です。指数が何重にも重なった巨大な数をそのまま計算すると現実的な時間では処理できないため、数論のテクニックを活用して効率よく求めるのがポイントになります。

問題の例

たとえば、入力が a = 2、b = 3、c = 2、d = 4、n = 10 の場合を考えてみましょう。

2^3 = 8
2^4 = 16
8^16 = 281474976710656
281474976710656 mod 10 = 6

この場合の出力は 6 となります。

解法のポイント

((a^b)^(c^d)) は a^(b・c^d) と変形できますが、この指数は天文学的な大きさになる可能性があります。そこでオイラーの定理を利用します。

オイラーの定理とは、「a と n が互いに素のとき、a^φ(n) ≡ 1 (mod n) が成り立つ」というものです。ここで φ(n) は オイラーのトーティエント関数（1以上n以下の整数のうち、nと互いに素である数の個数）を表します。この性質を使うと、巨大な指数を φ(n) を法として縮約でき、膨大なべき乗計算を回避できます。

アルゴリズムの手順

1. helper() 関数：オイラー関数 φ(n) を計算する

• p := n、i := 2 と初期化します。
• i × i ≤ n の間、次の処理を繰り返します。
  • n が i で割り切れるなら、p := p − ⌊p / i⌋ とします。
  • n が i で割り切れる限り、n := ⌊n / i⌋ を繰り返して因数を除去します。
  • i が 2 以外なら i := i + 2（奇数のみを試す）、そうでなければ i := i + 1 とします。
• ループ後に n > 1 が残っていれば、それは素因数なので p := p − ⌊p / n⌋ とします。
• p を返します。

2. solve() 関数：特別なケースを順に判定する

• b = 0、または「c = 0 かつ d ≠ 0」の場合：指数全体が 0 になるため、a^0 mod n を返します。
• c = 1、または d = 0 の場合：c^d = 1 となるため、a^b mod n を返します。
• a = 0、または a mod n = 0 の場合：結果は 0 です。
• d = 1 の場合：(a^b)^c = a^(b・c) となるため、a^(b・c) mod n を返します。
• 上記以外の場合：
  • p := helper(n) で φ(n) を求めます。
  • e := (c^d mod p) + p とします。「+ p」を加えることで、縮約後の指数でも元の値と合同であることが保証され、境界ケースでの誤りを防げます。
  • ((a^b mod n)^e) mod n を返します。

Pythonによる実装例

それでは、実際のコードを見てみましょう。Pythonの組み込み関数 pow(a, b, m) は3引数形式で使うと、剰余を取りながら高速にべき乗計算（冪乗の二進法）を行ってくれるため、この種の問題に最適です。

def helper(n):
    # オイラーのトーティエント関数 φ(n) を計算する
    p = n
    i = 2
    while i * i <= n:
        if n % i == 0:
            p -= p // i
        while n % i == 0:
            n = n // i
        if i != 2:
            i += 2
        else:
            i += 1
    if n > 1:
        p -= p // n
    return p

def solve(a, b, c, d, n):
    if b == 0 or (c == 0 and d != 0):
        return pow(a, 0, n)
    if c == 1 or d == 0:
        return pow(a, b, n)
    if a == 0 or a % n == 0:
        return 0
    if d == 1:
        return pow(a, b * c, n)
    p = helper(n)
    e = pow(c, d, p) + p
    return pow(pow(a, b, n), e, n)

print(solve(2, 3, 2, 4, 10))

入力

2, 3, 2, 4, 10

出力

6

コードの解説

helper() の動作

helper() は試し割り法による素因数分解をもとに、オイラー関数の公式 φ(n) = n・∏(1 − 1/pᵢ) を実装しています。n をその素因数で順に割りながら、p から p/i ずつ差し引いていくことで φ(n) を求めます。ループ変数 i を 2 の次は奇数だけ増やしている点も、計算量を抑える工夫です。

solve() の動作

solve() では、まず指数が自明になるケース（0 や 1 になるパターン）をすべて先に処理することで、無駄な計算を省いています。その後に残る一般ケースに対してのみ、オイラーの定理による指数の縮約を適用します。e = pow(c, d, p) + p という式により、c^d がどんなに巨大でも、その値を φ(n) 以下の範囲に落とし込みつつ正しい合同関係を維持できます。最終的に、内側の pow(a, b, n) で底を先に縮約し、外側の pow(..., e, n) で仕上げのべき乗剰余を計算しています。

この手法を使えば、指数が数十桁・数百桁になったとしても、瞬時に答えを求めることができます。競技プログラミングで頻出の「べき乗タワー型の剰余計算」における定番テクニックなので、ぜひマスターしておきましょう。