平方剰余の根

ある整数xに対して

x^2 = a (mod p)

となるようなaを「pの平方剰余」という。

以下ではaが与えられたときに

x^2 = a (mod p)

を満たすxを求めることを考える。ただし今回は簡単のためpが素数の場合のみを考える。

解の存在

ルジャンドル記号(a/p)を以下のように定義する。

a = 0 (mod p)ならば, (a/p) = 0
aが法pで平方剰余ならば、(a/p) = 1
aが法pで平方剰余でなければ、(a/p) = -1

pが奇素数の場合は、オイラーの基準により

(a/p) = a^{(p-1)/2}  (mod p)

が成り立つ。よって繰り返し二乗法を用いることで解の存在を簡単に判定出来る。

long long modpow(long long x, long long p, long long mod) {
    long long ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

// a is a quadratic residue modulo p?
int Legendre(long long a, long long p) {

    long long ret = modpow(a, (p-1)/2, p);

    if (ret == p-1)
        ret = -1;

    return ret;
}

解の公式

ルジャンドル記号が1の場合は、

x^2 = a (mod p)

を満たすxが存在する。

p=2の場合はxを全列挙すればよいため、以下ではpは奇素数とする。
pを8で割ったときの余りが3, 5, 7の場合は解の公式が存在する。

i) p = 3 (mod 8)、またはp = 7 (mod 8) の場合

x = ± a^{(p+1)/4} (mod p)

が解となる。

ii) p = 5 (mod 8)の場合

a^{(p-1)/4} = 1なら x = ± a^{(p+3)/8} (mod p)
a^{(p-1)/4} = -1ならx = ± 2^{(p-1)/4} a^{(p+3)/8} (mod p)

が解になる。

// calculate x such that x^2 = a (mod p)
// p is an odd prime that satisfies p = 3, 5, 7 (mod 8)
long long modsqrt(long long a, long long p) {
    
    // case i)
    if (p % 4 == 3)
        return modpow(a, (p+1)/4, p);

    // case ii)
    if (p % 8 == 5) {
        if (modpow(a, (p-1)/4, p) == 1)
            return modpow(a, (p+3)/8, p);
        else
            return modpow(2, (p-1)/4, p) * modpow(a, (p+3)/8, p) % p;
    }
        
    throw;
}

これで75%の場合は公式を使えば解けることが分かる。
しかし残りの25%の場合（p = 1 (mod 8)のとき）は公式が存在しないので、以下で示すアルゴリズムを使って計算するとよい。

非決定性アルゴリズム

Tonelli–Shanks algorithmを使うと、任意の奇素数pについて

x^2 = a (mod p)

を解くことが出来る。

以下に実装例を載せておく。

#include <iostream>
#include <random>

using namespace std;

long long modpow(long long x, long long p, long long mod) {
    long long ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

// a is a quadratic residue modulo p?
int Legendre(long long a, long long p) {

    long long ret = modpow(a, (p-1)/2, p);

    if (ret == p-1)
        ret = -1;

    return ret;
}

// Tonelli–Shanks algorithm
// calculate R such that R^2 = n (mod p)
// p is an odd prime
// n is a positive integer that satisfies (n/p) = 1
// See http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
long long modsqrt(long long n, long long p) {

    // step 1
    long long Q = p-1;
    long long S = 0;
    while (Q % 2 == 0) {
        ++S;
        Q /= 2;
    }

    if (S == 1) {
        return modpow(n, (p+1)/4, p);
    }

    // step 2
    default_random_engine generator;
    uniform_int_distribution<long long> distribution(0, p);
    
    long long z = 1;
    while (Legendre(z, p) != -1) {
        z = distribution(generator);
    }

    long long c = modpow(z, Q, p);
    
    // step 3
    long long R = modpow(n, (Q+1)/2, p);
    long long t = modpow(n, Q, p);
    long long M = S;

    // step 4
    while (t != 1) {
        long long i;
        long long t2 = t;
        for (i = 1; i < M; i++) {
            t2 = t2 * t2 % p;
            if (t2 == 1)
                break;
        }

        long long b = modpow(c, 1LL<<(M-i-1), p);
        R = R * b % p;
        t = (t * b % p) * b % p;
        c = b * b % p;
        M = i;
    }
    
    return R;
}

// test
bool isPrime(long long n) {
    if (n < 2)
        return false;
    
    for (long long i = 2; i * i <= n; i++)
        if (n % i == 0)
            return false;
    
    return true;           
}

#include <cassert>

int main(int argc, char **argv) {
    
    for (int i = 3; i < 10000; i++) {
        if (!isPrime(i))
            continue;

        for (int j = 0; j < i; j++) {
            if (Legendre(j, i) != 1)
                continue;

            long long x = modsqrt(j, i);
            assert (x * x % i == j);
        }

    }

    return 0;
}

参考URL

[1] Euler's criterion - Wikipedia, the free encyclopedia
[2] Computing square roots mod p
[3] Tonelli–Shanks algorithm - Wikipedia, the free encyclopedia