初等整数論の基本の復習

2016/8/20誤記修正

基本的な数学がいまいち使いこなせてないので復習しておく。

参考にしたのは蟻本と以下のサイト。
http://www2.cc.niigata-u.ac.jp/~takeuchi/tbasic/BackGround/

ユークリッドの互除法

整数aとbについて、a > b > 0 かつ a % b を r としたとき、 ${ GCD(a,\,b) = GCD(b,\,r) }$ が成り立つ。

    static int Gcd(int a, int b)
    {
        if (b == 0) return a;
        var p = a > b ? a : b;
        return Gcd(b, p % b);
    }

毎回、半分以下になっていくので計算量は ${ O\left ( log \, max (a, b) \right ) }$ 程度になる。

拡張ユークリッドの互除法

${ ax + by = GCD(a,b) }$ の整数解 ${(x, y)}$ は常に存在し、計算することができる。
この方程式の解を求める関数を int ExtGcd(int a, int b, ref int x, ref int y)とすると、次のように求められる。返り値はGCD(a, b)。

    static int ExtGcd(int a, int b, ref int x, ref int y)
    {
        int d = a;
        if (b != 0)
        {
            d = ExtGcd(b, a % b, ref y, ref x);
            y -= (a / b) * x;
        }
        else
        {
            x = 1;
            y = 0;
        }
        return d;
    }

手でやってみると、なぜこうなるかよく分かる。計算量はGcd()と同じくらい。

エラトステネスの篩

整数iが素数かどうか判定する。また、j番目の素数を求める。これは簡単。

    static void Sieve(int[] prime, bool[] isPrime)
    {
        for (int i = 0; i < prime.Length; i++) prime[i] = -1;
        for (int i = 0; i < isPrime.Length; i++) isPrime[i] = true;
        isPrime[0] = isPrime[1] = false;

        var idx = 0;
        for (int i = 2; i < isPrime.Length; i++)
        {
            if (isPrime[i])
            {
                prime[++idx] = i;
                for (int j = 2 * i; j < isPrime.Length; j += i) isPrime[j] = false;
            }
        }
    }

蟻本によると、計算量は ${ O(n \, log \, log \, n) }$ 程度で、通常の範囲ならほぼ線形と考えてよいらしい。

整数の合同

整数aとbのmで割った余りが等しい（a-bがmで割り切れる）とき、
${ a \equiv b \, (mod\: m) }$
と表し、aとbはmを法にして合同、という。

線形の合同には次のような性質がある。
${ a \equiv a \, (mod\: m) }$
${ a \equiv b \, (mod\: m) }$ 　ならば　 ${ b \equiv a \, (mod\: m) }$
${ a \equiv b \, (mod\: m) }$ 　かつ　 ${ b \equiv c \, (mod\: m) }$ 　ならば　 ${ a \equiv c \, (mod\: m) }$

${ a \equiv c \, (mod\: m) }$ 、 ${ b \equiv d \, (mod\: m) }$ 　とすると、
${ a + b \equiv c + d (mod\: m) }$
${ a - b \equiv c - d (mod\: m) }$
${ a \times b \equiv c \times d (mod\: m) }$
が成り立つ（除算は成り立たない）。

${ a \times c \equiv b \times c(mod\: m) }$ の場合、
${ a \equiv b(mod\:\, m\, /\, gcd(m,\, c)) }$
が成り立つ。

べき乗法

以下の単純なやり方で普通は十分。

    static ulong ModPow(ulong x, ulong n, ulong mod)
    {
        ulong ret = 1;
        while (n > 0)
        {
            if ((n & 1) == 1) ret = ret * x % mod;
            x = x * x % mod;
            n >>= 1;
        }
        return ret;
    }

計算量は ${ O(log\, n) }$ 。
modが大きいとオーバーフローする（途中でmod^2程度の数が表れる）ので、改善したアルゴリズムがここに紹介されていた。
http://www2.cc.niigata-u.ac.jp/~takeuchi/tbasic/BackGround/power.html

中国の剰余定理

${ m_{1}, m_{2} }$ が互いに素な整数なら、任意に与えられる整数 ${a_1, a_2}$ に対し
${ x \equiv a_{1} (mod\: m_{1}) }$
${ x \equiv a_{2} (mod\: m_{2}) }$
を満たす整数xが、 ${ m_{1}m_{2} }$ を法にしてただ一つ存在する。

これを拡張して、整数 ${ m_{1}, m_{2}, ...m_k }$ がどの二つも互いに素なら、任意に与えられる整数 ${a_1, a_2, ... a_k}$ に対し
${ x \equiv a_{1} (mod\: m_{1}) }$
${ x \equiv a_{2} (mod\: m_{2}) }$
...
${ x \equiv a_{k} (mod\: m_{k}) }$
を満たす整数xが、 ${ m_{1}m_{2}...m_k }$ を法にしてただ一つ存在する。

これを利用して、TopcoderのSRM657(Div2Hard)を解いてみる。

整数a、b、cが与えられたとき、 ${ ax^{2} + bx + c \equiv 0\: mod\: 10^{9} }$ を満たすxを求める。
（a、b、cは0～999,999,999）

xの範囲が0,...999999999と大きいので、単純にループさせては時間がかかりすぎる。

${10^{9}}$ は ${ 2^{9} }\times{ 5^{9} }$ かつ ${ 2^{9} }$ と ${ 5^{9} }$ は互いに素なので、中国の剰余定理を利用して

${ ax^{2} + bx + c \equiv 0\: mod\: 2^{9} }$
${ ax^{2} + bx + c \equiv 0\: mod\: 5^{9} }$

を満たすxを求めればよい。

それぞれの式から、

${ x \equiv x_{2}\: mod \: 2^{9} }$
${ x \equiv x_{5}\: mod \: 5^{9} }$

が導き出されるので、片方を満たす最小非負代表を求めてから、もう片方の式を満足できるまで法の数だけ足していけばよい。
下の例では、 ${ x \equiv x_{5}\: mod \: 5^{9} }$ を満たす最小値を算出し、 ${ x \equiv x_{2}\: mod \: 2^{9} }$ が満たされるまで ${ 5^{9} }$ を足していっている。

public class PolynomialRemainder {
    int POW2 = 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2 * 2;
    int POW5 = 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5;

    ulong GetAx2BxC(ulong mod, ulong x, ulong a, ulong b, ulong c)
    {
        return ((a * (x * x % mod) % mod) + (b * x) % mod + c) % mod;
    }

    int GetRoot(int mod, int a, int b, int c)
    {
        for (int x = 0; x < mod; x++)
        {
            var ret = GetAx2BxC((ulong)mod, (ulong)x, (ulong)a, (ulong)b, (ulong)c);
            if (ret % (ulong)mod == 0)
            {
                return x;
            }
        }

        return -1;
    }

    public int findRoot(int a, int b, int c)
    {
        var x2 = GetRoot(POW2, a, b, c);
        var x5 = GetRoot(POW5, a, b, c);

        if (x2 == -1 || x5 == -1) return -1;

        int x = x5;
        while (x % POW2 != x2)
            x += POW5;

        return x;
    }
}

ちゃんとした教科書なしで復習したせいか、かえってフラストレーションが溜まった気がする。
以下をポチったので、今月は数学月間にしよう。
http://www.amazon.co.jp/dp/4062573865/ref=pe_492632_159100282_TE_itemwww.amazon.co.jp
http://www.amazon.co.jp/dp/4062575957/ref=pe_492632_159100282_TE_itemwww.amazon.co.jp