メモ用紙2.0

日記 - 2013-10

まどか☆マギカ

劇場版まどか☆マギカ観てきました。感想。とても良い話だった。

色紙は杏さやでした。杏子かわいい！

ここにはかつてコメントが表示されていました

2013-10-25

研究

前回の記事にツッコミが入ったので更新しました。

なんかやりたいことは溜まってるんですが、研究してPOJ解いてアニメ観ると1日が終わっているのでなかなか進捗のない1週間でした。まあ、研究している分だけ良いとも言いますが。でも今やってる作業は、自分の研究に関わってくるかというと半々くらいのとこなので、手放しで喜べない感じもします。楽しいけど。

ここにはかつてコメントが表示されていました

2013-10-17

高速素因数分解

知ってしまえば当たり前の話ですが、意外と蟻本にも載っていない話なので。

しばらく前に、Project Eulerの問題で素因数分解をする必要に迫られました。素因数分解をするには、当然ながら素数を求めなくてはいけません。これはエラトステネスの篩などの簡単なアルゴリズムで、比較的高速に列挙することができます。

問題はこの後、ある整数Nが与えられたとき、これを素因数分解するにはどうするか。ぱっと思いつくのは、さっき求めた素数リストで順にNを割っていく方法です。 N以下の素数の個数は、素数定理より N/ln(N) 個くらいで、素因数の個数は高々ln(N)個であることを合わせて考えると O(N/ln(N) + ln(N)) で素因数分解ができます。しかしよく考えると、素数を列挙しておかなくても2から順番に試し割りしていけば O(N + ln(N)) で素因数分解ができるので、期待するほど速くもありません……。（訂正:2013/10/07 04:42）また、実際には小さい素数の積で表せてしまう数も多いので、ループ回数自体は N/ln(N) 回よりもずっと少ないことが期待できます。

それでも、大量の数を素因数分解しようとすると結構な時間がかかってしまいます。そこで、エラトステネスの篩の途中で、それぞれの数について最小の素因数を覚えておくようにします。こうすると、任意のNについて最小の素因数が O(1) で求められるので、Nをこの素因数で割って、また最小の素因数を探して……と繰り返していくと、さっきのアルゴリズムで「ハズレ」の素数を引かないのと同じことになるので O(ln(N)) くらいで素因数分解ができるようになります。爆速！

というわけで、実際に計測してみました。 2からMAXまでの整数を素因数分解して、全ての素因数を合計するコードです。

 1#include <iostream>
 2#include <ctime>
 3#include <vector>
 4#include <algorithm>
 5
 6using namespace std;
 7
 8const int MAX = 10000000;
 9
10bool is_prime[MAX+1];
11int min_factor[MAX+1];
12vector<int> primes;
13
14int main() {
15    // Sieve of Eratosthenes
16    fill_n(is_prime, MAX+1, true);
17    is_prime[0] = is_prime[1] = false;
18    min_factor[0] = 0;
19    min_factor[1] = 1;
20    for(int i = 2; i <= MAX; ++i) {
21        if(is_prime[i]) {
22            min_factor[i] = i;
23            primes.push_back(i);
24            for(long long j = static_cast<long long>(i)*i; j <= MAX; j += i) {
25                is_prime[j] = false;
26                if(min_factor[j] == 0) min_factor[j] = i;
27            }
28        }
29    }
30
31    ////////// Factoring by division loop //////////
32    {
33        long long sum = 0;
34        clock_t start = clock();
35        for(int i = 2; i <= MAX; ++i) {
36            int tmp = i;
37            for(int p : primes) {
38                while(tmp % p == 0) {
39                    sum += p;
40                    tmp /= p;
41                }
42                if(tmp == 1) break;
43                if(is_prime[tmp]) {
44                    sum += tmp;
45                    break;
46                }
47            }
48        }
49        clock_t end = clock();
50        cout << "Factoring by division loop:" << endl;
51        cout << "\tValue: " << sum << endl;
52        cout << "\tTime: " << static_cast<double>(end-start) / CLOCKS_PER_SEC << endl;
53    }
54    ////////// Factoring by min_factor //////////
55    {
56        long long sum = 0;
57        clock_t start = clock();
58        for(int i = 2; i <= MAX; ++i) {
59            int tmp = i;
60            while(tmp > 1) {
61                const int fact = min_factor[tmp];
62                sum += fact;
63                tmp /= fact;
64            }
65        }
66        clock_t end = clock();
67        cout << "Factoring by min_factor:" << endl;
68        cout << "\tValue: " << sum << endl;
69        cout << "\tTime: " << static_cast<double>(end-start) / CLOCKS_PER_SEC << endl;
70    }
71    return 0;
72}

結果：

 1N = 10000000 (1千万)
 2Factoring by division loop:
 3    Value: 5495532284181
 4    Time: 1.21331
 5Factoring by min_factor:
 6    Value: 5495532284181
 7    Time: 0.384788
 8
 9N = 100000000 (1億)
10Factoring by division loop:
11    Value: 476028748414934
12    Time: 19.7341
13Factoring by min_factor:
14    Value: 476028748414934
15    Time: 4.85803

素因数テーブルを使うと、~~ケタ違いに（訂正:2013/10/17 04:42）~~かなり速くなっていることが分かります。ということで、知らなかった人はぜひこの手法を覚えましょう。

2013-10-25追記

素数だけで試し割りをするときにはMAXまでの素数表を持つ必要がないので、上記のプログラムでは不公平だとコメントで指摘されました。ループ回数はsqrt(N)回くらいしか違わないのであまり効いてこなさそうですが、一応お互いにベストな条件で計測してみます。

 1#include <iostream>
 2#include <ctime>
 3#include <vector>
 4#include <algorithm>
 5
 6using namespace std;
 7
 8const int MAX = 10000000;
 9const int LIM = sqrt(MAX);
10
11bool is_prime[MAX+1];
12int min_factor[MAX+1];
13vector<int> primes;
14
15void division_loop() {
16    fill_n(is_prime, MAX+1, true);
17    is_prime[0] = is_prime[1] = false;
18    for(int i = 2; i <= LIM; ++i) {
19        if(is_prime[i]) {
20            primes.push_back(i);
21            for(long long j = static_cast<long long>(i)*i; j <= MAX; j += i) {
22                is_prime[j] = false;
23            }
24        }
25    }
26
27    ////////// Factoring by division loop //////////
28    {
29        long long sum = 0;
30        for(int i = 2; i <= MAX; ++i) {
31            int tmp = i;
32            for(int p : primes) {
33                while(tmp % p == 0) {
34                    sum += p;
35                    tmp /= p;
36                }
37                if(tmp == 1) break;
38                if(is_prime[tmp]) {
39                    sum += tmp;
40                    break;
41                }
42            }
43        }
44        cout << "division_loop: " << sum << endl;
45    }
46}
47
48void min_factor_loop() {
49    fill_n(is_prime, MAX+1, true);
50    is_prime[0] = is_prime[1] = false;
51    min_factor[0] = 0;
52    min_factor[1] = 1;
53    for(int i = 2; i <= MAX; ++i) {
54        if(is_prime[i]) {
55            min_factor[i] = i;
56            for(long long j = static_cast<long long>(i)*i; j <= MAX; j += i) {
57                is_prime[j] = false;
58                if(min_factor[j] == 0) min_factor[j] = i;
59            }
60        }
61    }
62
63    ////////// Factoring by min_factor //////////
64    {
65        long long sum = 0;
66        for(int i = 2; i <= MAX; ++i) {
67            int tmp = i;
68            while(tmp > 1) {
69                const int fact = min_factor[tmp];
70                sum += fact;
71                tmp /= fact;
72            }
73        }
74        cout << "min_factor_loop: " << sum << endl;
75    }
76}
77
78void (*fp[])() = {division_loop, min_factor_loop};
79
80int main(int argc, char **argv) {
81    if(argc < 2) {
82        cerr << "Usage: " << argv[0] << " number" << endl;
83        return 0;
84    }
85    int i = atoi(argv[1]);
86    clock_t start = clock();
87    fp[i]();
88    clock_t end = clock();
89    cout << static_cast<double>(end - start) / CLOCKS_PER_SEC << " sec." << endl;
90    return 0;
91}

結果：

 1N = 1000万
 2$ ./factor 0
 3division_loop: 5495532284181
 41.26641 sec.
 5$ ./factor 1
 6min_factor_loop: 5495532284181
 70.527824 sec.
 8
 9N = 1億
10$ ./factor 0
11division_loop: 476028748414934
1219.4185 sec.
13$ ./factor 1
14min_factor_loop: 476028748414934
156.51549 sec.

結果としてはやはりたいして変わらないようです。

ところでmin_factorが前回と比べて遅くなっていますが、どうも前回はキャッシュが効いていたみたいですね……。 RubyのbenchmarkみたいなのがC++にも欲しいところです。

ここにはかつてコメントが表示されていました

日記 - 2013-10

2013-10-27

まどか☆マギカ

2013-10-25

研究

2013-10-17

高速素因数分解

2013-10-25追記