Day 132-Q4 | Cipollaのアルゴリズム（mod p 平方根）

問題

素数$p$（奇素数）と整数$n$（$0\le n最小の非負整数解（$x$と$p-x$のうち小さい方）を出力せよ。存在しない場合は-1を出力せよ。

入力形式

p n

制約

$3 \le p \le 10^9+7$（素数）

$0 \le n < p$

入出力例

入力例1

1000000007 241383831

出力例1

入力例2

13 2

出力例2

-1

$2$は$13$を法とする平方剰余ではない：$1^2,\dots,12^2 \bmod 13$の集合は$\{1,4,9,3,12,10\}$で$2$を含まない。

概念図

ヒント（段階的開示）

ヒント1（方向性）

$p\equiv3\pmod4$なら$x=n^{(p+1)/4}\bmod p$の1発公式が使えるが、一般の$p$には使えない。Tonelli-Shanks法は乗法群の位数を$2^s\cdot q$に分解して反復するが実装がやや煩雑。ここでは全く異なる発想のCipollaのアルゴリズムを使う。

ヒント2（アプローチ）

オイラーの規準$n^{(p-1)/2}\bmod p$で平方剰余判定する。ランダムに$a$を選び$w=a^2-n\bmod p$が非剰余になるものを探す（確率1/2）。$w$が非剰余なら拡大体$GF(p^2)=\{x+y\sqrt w\}$を考えられ、その中で$(a+\sqrt w)^{(p+1)/2}$を計算すると虚部が0になり実部が$n$の平方根になる。

ヒント3（誘導）

def mul(A, B, w, p):
    x1, y1 = A
    x2, y2 = B
    return ((x1*x2 + y1*y2 % p * w) % p, (x1*y2 + x2*y1) % p)

def power(base, e, w, p):
    result = (1, 0)
    while e > 0:
        if e & 1:
            result = mul(result, base, w, p)
        base = mul(base, base, w, p)
        e >>= 1
    return result

模範解答 (Python)

import sys
import random

def legendre(a, p):
    return pow(a, (p - 1) // 2, p)

def cipolla(n, p):
    n %= p
    if n == 0:
        return 0
    if legendre(n, p) != 1:
        return -1

    a = random.randint(1, p - 1)
    while legendre((a * a - n) % p, p) != p - 1:
        a = random.randint(1, p - 1)
    w = (a * a - n) % p

    def mul(A, B):
        x1, y1 = A
        x2, y2 = B
        return ((x1 * x2 + y1 * y2 % p * w) % p, (x1 * y2 + x2 * y1) % p)

    result = (1, 0)
    base = (a, 1)
    e = (p + 1) // 2
    while e > 0:
        if e & 1:
            result = mul(result, base)
        base = mul(base, base)
        e >>= 1

    x = result[0]
    return min(x, p - x)

def solve():
    p, n = map(int, sys.stdin.read().split())
    print(cipolla(n, p))

solve()

Step-by-Step 解説

1平方剰余判定
$n^{(p-1)/2}\bmod p$が1なら平方剰余、$p-1$なら非剰余。これを先にやらないと次のwhileが無限ループする。

2非剰余wを見つける
$a$をランダムに選び$w=a^2-n\bmod p$が非剰余になるまで繰り返す。期待2回程度で見つかる。

3拡大体GF(p)[√w]での計算
$w$が非剰余なので$x+y\sqrt w$の全体は乗算について閉じた体（GF(p²)と同型）になる。

4(a+√w)^((p+1)/2)が答えになる理由
フロベニウス写像より$(a+\sqrt w)^p=a-\sqrt w$。よって$(a+\sqrt w)^{p+1}=a^2-w=n$。半分乗すれば√nの候補。虚部は理論上必ず0になる。

5計算量
反復2乗法でO(log p)回の拡大体乗算（各O(1)）。乱択の期待反復回数は定数。全体O(log p)。

よくあるミス

ミス	原因	正しい書き方
非剰余nでも無限ループする	事前のLegendre記号判定を省略	`legendre(n,p)!=1`なら即座に-1を返す分岐を最初に入れる
w=0のケースでの挙動を誤解する	legendre(0,p)=0で非剰余判定条件を満たさず再試行され続けると誤解	whileループはそのまま`!=p-1`条件で正しく再試行される（0は非剰余として扱われない）
出力がxとp-xで安定しない	2つの根のどちらが出るか固定していない	`min(x, p-x)`のように正規化してから出力する
p=2特有の分岐を気にしてコードを複雑化	制約に反するケースを想定してしまう	制約でpは奇素数と保証されているので不要

次のステップ

発展: Tonelli-Shanks法で同じ問題を解き、Cipollaとの実装量・速度を比較する
発展2: mod p^k（素数冪）やmod n（合成数、Garnerのアルゴリズムで統合）への一般化
次回予告: 次回のテーマは今後のセッションで選定

自己評価

完全理解おおむね理解要復習

Day 132-Q4 — Cipollaのアルゴリズム（mod p 平方根の別解法）

問題