平方剰余を理解したい

はじめに

某CTFでCryptoの問題に挑戦したのですが、$\rm{mod}$ $N$における平方根の求め方をよく理解できておらず、やや苦戦してしまいました。せっかくなので、周辺知識を再確認しておきます。
よくお世話になるこちらの書籍などを参考にして、Pythonで実装までしようと思います。
※計算することにフォーカスするため、証明は省略します。

以下、$Z^{*}_{N}=\{0, 1,\dots N-1\}$と表記します。

※注意
・🐺が勉強としてやっているだけなので、参考程度に見てください。
・この記事は予告なく加筆・修正される場合があります

平方剰余

まずはじめに、$\rm{mod}$ $N$で平方根が存在するかどうかという話です。

平方剰余とは

$x\in{Z^{*}_{N}}$に対し、$a=x^{2}\mod{N}$となる$a\in{Z^{*}_{N}}$を、法$N$における平方剰余といいます。逆に、$a=x^{2}\mod{N}$とならない$a\in{Z^{*}_{N}}$を、法$N$における平方非剰余といいます。

たとえば、
$1^2\equiv{1}\pmod{7}$
$2^2\equiv{4}\pmod{7}$
$3^2\equiv{2}\pmod{7}$
$4^2\equiv{2}\pmod{7}$
$5^2\equiv{4}\pmod{7}$
$6^2\equiv{1}\pmod{7}$
となるので、$1, 2, 4$は$\rm{mod}$ $7$における平方剰余、$\rm{mod}$ $7$における平方非剰余です。

ルジャンドル記号

$p$が奇素数で$a\in{Z^{*}_{p}}$のとき、$a$が法$p$における平方剰余か平方非剰余かを示すルジャンドル記号$(a/p)$というものがあり、次のように定義されています。

• $a$が法$p$における平方剰余のとき、$\left(\frac{a}{p}\right)=1$
• $a$が法$p$における平方非剰余のとき、$\left(\frac{a}{p}\right)=-1$
• ($a\equiv{0}{\pmod{p}}$のとき、$\left(\frac{a}{p}\right)=0$)

ルジャンドル記号は

$$\left(\frac{a}{p}\right)\equiv{a^{\frac{p-1}{2}}}\pmod{p}$$

で求めることができ、これをオイラーの規準といいます。

オイラーの規準は、Pythonで次のように書けます。

def legendre_symbol(a: int, p: int) -> int:
	legendre_symbol = pow(a, (p - 1) // 2, p)
	if legendre_symbol == 1:
		return 1
	elif legendre_symbol == p - 1:
		return -1
	else:
		return 0

平方根

$p$が奇素数のとき、$k\in{Z}$を使って$p=4k+3$または$p=4k+1$と表すことができます。それぞれの場合で、法$p$における平方剰余$a$の平方根を求める方法を示します。

$p=4k+3$のとき

$a$が法$p$における平方剰余であれば、オイラーの規準より

$$a^{\frac{p-1}{2}}\mod{p}=1$$

となります。よって、

$$(\pm a^{(p+1)/4})^2\equiv a^{(p+1)/2}\equiv a^{(p-1)/2}\cdot a=a\pmod{p}$$

となるので、$\pm a^{(p+1)/4}\mod{p}$が法$p$における$a$の平方根です。

Pythonで次のように書けます。

def mod_sqrt(a:int, p:int) -> list:
	a %= p
	x = pow(a, (p + 1) // 4, p)
	return [x, p - x]

具体的な値で試してみましょう。奇素数$p$と$(a/p)=1$な$a$を適当に用意しました。

>>> p = 1122278419650073878178964596327
>>> p % 4
3
>>> a = 299321295248724911916559376201
>>> legendre_symbol(a, p)
1
>>> print(mod_sqrt(a, p))
[239651537938722371683492056230, 882626881711351506495472540097]
>>> pow(239651537938722371683492056230, 2, p) == a
True
>>> pow(882626881711351506495472540097, 2, p) == a
True

$p=4k+1$のとき

$p=4k+1$のときは、少し複雑な方法が必要となります。以下に示すTonelli-Shanksアルゴリズムで求められます。

1. $(z/p)=-1$($z$は法$p$における平方非剰余)となる$z\in{Z_{p}^{*}}$を選ぶ
2. $p-1=2^{s}Q$ ($Q$は奇数)となる$s, Q$を求める
3. 以下のように初期値を設定する
   • $M=s$
   • $c=z^{Q}\mod{p}$
   • $t=a^{Q}\mod{p}$
   • $R=a^{(Q+1)/2}\mod{p}$
4. 以下を繰り返す
   • $t=0$のとき、$r=0$
   • $t=1$のとき、$r=R$
   • それ以外のとき、$t^{2i}=1 (0<i<M)$となる最小の$i$を求めて$b=c^{2^{M-i-1}}\mod{p}$を計算し、以下のように値を更新する
     • $M\leftarrow i$
     • $c\leftarrow{b^{2}}\mod{p}$
     • $t\leftarrow{tb^{2}}\mod{p}$
     • $R\leftarrow{Rb}\mod{p}$

この結果として得られる$\pm{r}$が法$p$における$a$の平方根です。

Pythonで次のように書けます。

def mod_sqrt(a:int, p:int) -> list:
    a %= p
	
    # 1. 平方非剰余を選ぶ
    z = 2
    while legendre_symbol(z, p) != -1:
        z += 1
	
    # 2. s, Qを求める
    s = 0
    q = p - 1
    while q % 2 == 0:
        q //= 2
        s += 1
		
    # 3. M, c, t, Rの初期値を設定する
    m = s
    c = pow(z, q, p)
    t = pow(a, q, p)
    r = pow(a, (q + 1) // 2, p)
	
    # 4. M, c, t, Rの値を更新する
    while t != 1 :
        if t == 0: # t = 0のときは0を返す
                return [0, 0]
        i = 1
        while pow(t, pow(2, i), p) != 1:
                i += 1
        b = pow(c, pow(2, m - i - 1, p))
        m = i
        c = pow(b, 2, p)
        t = (t * b**2) % p
        r = (r * b) % p
    return [r, p - r]

具体的な値で試してみましょう。こちらも、奇素数$p$と$(a/p)=1$な$a$を適当に用意しました。

>>> p = 822985351289460529207697791361
>>> p % 4
1
>>> a = 635837859252342443813240556561
>>> legendre_symbol(a, p)
1
>>> print(mod_sqrt(a, p))
[249110913840436119894350236958, 573874437449024409313347554403]
>>> pow(249110913840436119894350236958, 2, p) == a
True
>>> pow(573874437449024409313347554403, 2, p) == a
True

まとめ

完成したスクリプトはこんな感じです。

def legendre_symbol(a: int, p: int) -> int:
    legendre_symbol = pow(a, (p - 1) // 2, p)
    if legendre_symbol == 1:
        return 1
    elif legendre_symbol == p - 1:
        return -1
    else:
        return 0

def mod_sqrt(a:int, p:int) -> list:
    a %= p
    
    if p % 4 == 3:
        a %= p
        x = pow(a, (p + 1) // 4, p)
        return [x, p - x]
	
    elif p % 4 == 1:
    # 1. 平方非剰余を選ぶ
        z = 2
        while legendre_symbol(z, p) != -1:
            z += 1
        
        # 2. s, Qを求める
        s = 0
        q = p - 1
        while q % 2 == 0:
            q //= 2
            s += 1
            
        # 3. M, c, t, Rの初期値を設定する
        m = s
        c = pow(z, q, p)
        t = pow(a, q, p)
        r = pow(a, (q + 1) // 2, p)
        
        # 4. M, c, t, Rの値を更新する
        while t != 1 :
            if t == 0: # t = 0のときは0を返す
                    return [0, 0]
            i = 1
            while pow(t, pow(2, i), p) != 1:
                    i += 1
            b = pow(c, pow(2, m - i - 1, p))
            m = i
            c = pow(b, 2, p)
            t = (t * b**2) % p
            r = (r * b) % p
        return [r, p - r]
    
def main():
    print(mod_sqrt(3, 11))
    # [5, 6]

if __name__ == "__main__":
     main()

おまけ

$\rm{mod}$ $N$における平方根を求めるpythonスクリプトをごちゃごちゃと書きましたが、実はSageMathでsqrtメソッドとして実装されています。

>>> from sage.all import *
>>> p = 822985351289460529207697791361
>>> a = 635837859252342443813240556561
>>> mod(a, p).sqrt(all=True)
[249110913840436119894350236958, 573874437449024409313347554403]