読者です 読者をやめる 読者になる 読者になる

初等整数論の基本の復習

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