オイラーの定理についてちょっとだけ

Google Code Jam Japan決勝の問題B「バクテリア」が非常におもしろかったので、オイラーの定理について以下の本で復習したことをまとめておきます。（数論の入門書としてかなりお勧めです。オイラーの定理、フェルマーの小定理、拡張ユークリッド互除法、中国の剰余定理などプロコンに頻出の定理はほぼ網羅しています。）

※一部自分の考えも入っているため、間違えている箇所がありましたらご指摘ねがいます。

　a^N を nで割った余りについて考えます。n = p1^k1 * p2^k2* .... pm^kmと素因数分解できるとします。このとき、a = A * p1^a1 * p2^a2 * .... pm^amとおきます。(Aはnと互いに素な数)。考えられるパターンは、以下の3つです。

1. a1 = a2 = .... am = 0のとき（aとnが互いに素なとき）
2. a1 >0, a2 > 0, ... am >0のとき(aがnの素因数をすべて含むとき)
3. その他のとき(aがいくつかのnの素因数を含むとき)

　1. はオイラーの定理より、a^φ(n) = 1なので、φ(n)はa^Nをnで割った余りの列の周期になります。

　2. のとき。a1 * N >= k1, a2 * N >= k2, ..., am * N >= kmとなるようなNに対してa^Nをnで割ったあまりは0となります。十分に大きなNに対しては、φ(n)はa^Nをnで割った余りの列の周期になります。"十分に大きな"とはどれくらいか？式を見た感じだと、かなり悪く見積もっても最悪でN = nとすればよさそうです。

　3. のとき。これは長くなるので書けませんが、

φ(n) = φ(p1^k1) * φ(p2^k2) * .... φ (pm^km)

となるのを利用して、mod(pi^ki), i = 1,2, ..., mにおいてa^Nを評価して結果を中国の剰余定理で統合するようなことをしています。この場合もφ(n)はa^Nをnで割った余りの列の周期になります。

断言はできませんが、こちらも式を追う限りだと、最悪でもN=nとすれば周期的になりそうです。

　プログラムでテストしてみます。

view plainprint?

1. const int MAX_N = 1000;  
2.   
3. int gcd(int a, int b) {  
4.     if (b == 0)  
5.         return a;  
6.     return gcd(b, a%b);  
7. }  
8.   
9. int euler(int x) {  
10.     int ret = x;  
11.   
12.     for (int i = 2; i * i <= x; i++) {  
13.         if (x % i == 0)  
14.             ret = ret * (i-1) / i;  
15.         while (x % i == 0)  
16.             x /= i;  
17.     }  
18.     if (x > 1)  
19.         ret = ret * (x-1) / x;  
20.   
21.     return ret;  
22. }  
23.   
24. int modPow(int x, int p, int mod) {  
25.     int ret = 1;  
26.     while (p > 0) {  
27.         if (p & 1)  
28.             ret = ret * x % mod;  
29.         x = x * x % mod;  
30.         p >>= 1;  
31.     }  
32.     return ret;  
33. }  
34.   
35. bool cover(int x, int n) {  
36.     for (int i = 2; i * i <= n; i++) {  
37.         if (n % i == 0 && x % i > 0)  
38.             return false;  
39.         while (x % i == 0)  
40.             x /= i;  
41.         while (n % i == 0)  
42.             n /= i;  
43.     }  
44.     if (n > 1 && n != x)  
45.         return false;  
46.     return true;  
47. }  
48.   
49. bool assert1(int x, int mod) {  
50.     int v = euler(mod);  
51.   
52.     if (modPow(x, v, mod) != 1) {  
53.         fprintf(stderr, "Error1 %d, %d", x, mod);  
54.         return false;  
55.     }  
56.     return true;  
57. }  
58.   
59. int assert2(int x, int mod) {  
60.     int y = 1;  
61.     int i = 0;  
62.     for (; y > 0; i++)  
63.         y = y * x % mod;  
64.   
65.     if (i > mod) {  
66.         fprintf(stderr, "Error2 %d, %d", x, mod);  
67.             return -1;  
68.     }  
69.     return i;  
70. }  
71.   
72. int assert3(int x, int mod) {  
73.     set<int> used;  
74.     int y = x % mod;  
75.   
76.     int i = 0;  
77.     for (; !used.count(y);i++) {  
78.         used.insert(y);  
79.         y = y * x % mod;  
80.     }  
81.     if (i > mod) {  
82.         fprintf(stderr, "Error3 %d, %d", x, mod);  
83.             return -1;  
84.     }  
85.   
86.     return i;  
87. }  
88.   
89. int main() {  
90.     // consider sequences i^k mod j  
91.   
92.     int passed1 = 0;  
93.     int passed2 = 0;  
94.     int passed3 = 0;  
95.   
96.     int loop2 = 0;  
97.     int loop3 = 0;  
98.     for (int i = 1; i < MAX_N; i++) {  
99.         for (int j = 2; j < MAX_N; j++) {  
100.             int g = gcd(i, j);  
101.   
102.             if (g == 1)  
103.                 passed1 += assert1(i, j);    // eurler's theorem  
104.             else if (cover(g, j)) {          // j has all prime factor of i  
105.                 int loop = assert2(i, j);  
106.                 passed2 +=  loop > 0;  
107.                 loop2 += loop;  
108.             }  
109.             else {  
110.                 int loop = assert3(i, j);  
111.                 passed3 += loop > 0;  
112.                 loop3 += loop;  
113.             }  
114.         }  
115.     }  
116.   
117.     int total = passed1 + passed2 + passed3;  
118.     printf("case1 event prob: %.2lf\n", 1.*passed1/total);  
119.     printf("case2 event prob: %.2lf average loop: %.2lf\n", 1.*passed2/total, 1.*loop2/passed2);  
120.     printf("case3 event prob: %.2lf average loop: %.2lf\n", 1.*passed3/total, 1.*loop3/passed3);  
121.   
122.     return 0;  
123. }  

const int MAX_N = 1000;

int gcd(int a, int b) {
    if (b == 0)
        return a;
    return gcd(b, a%b);
}

int euler(int x) {
    int ret = x;

    for (int i = 2; i * i <= x; i++) {
        if (x % i == 0)
            ret = ret * (i-1) / i;
        while (x % i == 0)
            x /= i;
    }
    if (x > 1)
        ret = ret * (x-1) / x;

    return ret;
}

int modPow(int x, int p, int mod) {
    int ret = 1;
    while (p > 0) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

bool cover(int x, int n) {
    for (int i = 2; i * i <= n; i++) {
        if (n % i == 0 && x % i > 0)
            return false;
        while (x % i == 0)
            x /= i;
        while (n % i == 0)
            n /= i;
    }
    if (n > 1 && n != x)
        return false;
    return true;
}

bool assert1(int x, int mod) {
    int v = euler(mod);

    if (modPow(x, v, mod) != 1) {
        fprintf(stderr, "Error1 %d, %d", x, mod);
        return false;
    }
    return true;
}

int assert2(int x, int mod) {
    int y = 1;
    int i = 0;
    for (; y > 0; i++)
        y = y * x % mod;

    if (i > mod) {
        fprintf(stderr, "Error2 %d, %d", x, mod);
            return -1;
    }
    return i;
}

int assert3(int x, int mod) {
    set<int> used;
    int y = x % mod;

    int i = 0;
    for (; !used.count(y);i++) {
        used.insert(y);
        y = y * x % mod;
    }
    if (i > mod) {
        fprintf(stderr, "Error3 %d, %d", x, mod);
            return -1;
    }

    return i;
}

int main() {
    // consider sequences i^k mod j

    int passed1 = 0;
    int passed2 = 0;
    int passed3 = 0;

    int loop2 = 0;
    int loop3 = 0;
    for (int i = 1; i < MAX_N; i++) {
        for (int j = 2; j < MAX_N; j++) {
            int g = gcd(i, j);

            if (g == 1)
                passed1 += assert1(i, j);    // eurler's theorem
            else if (cover(g, j)) {          // j has all prime factor of i
                int loop = assert2(i, j);
                passed2 +=  loop > 0;
                loop2 += loop;
            }
            else {
                int loop = assert3(i, j);
                passed3 += loop > 0;
                loop3 += loop;
            }
        }
    }

    int total = passed1 + passed2 + passed3;
    printf("case1 event prob: %.2lf\n", 1.*passed1/total);
    printf("case2 event prob: %.2lf average loop: %.2lf\n", 1.*passed2/total, 1.*loop2/passed2);
    printf("case3 event prob: %.2lf average loop: %.2lf\n", 1.*passed3/total, 1.*loop3/passed3);

    return 0;
}

以下、実験結果。ランダムに2つの正の整数を取ったとき、半分以上の確率でお互いに素になるようです。iがjの素因数をすべてカバーするケースはまれでMAX_Nが大きくなるに従って減っていくように見えます。また、i^nがmod jにおいて周期的になるのはn < MAX_Nと見積もっておけば十分だといえます。（case 2については、式と結果から考えてlogオーダーっぽい。）

MAX_N = 100

case1 event prob: 0.61

case2 event prob: 0.09 average loop: 2.07

case3 event prob: 0.30 average loop: 7.57

MAX_N = 1000

case1 event prob: 0.61

case2 event prob: 0.02 average loop: 2.70

case3 event prob: 0.37 average loop: 46.79