k-smooth numberの分布を考える

そろそろ100桁くらいの数の素因数分解でもチャレンジしたいなと思ってるけど、general number field sieveは難しすぎてまだ無理。ちょうど手ごろな感じのquadratic sieveというのがあったので、これなら実装できそうである。110桁程度までだと、GNFSよりQSの方が効率がいいらしい。

　さて、このquadratic sieveだが、大抵の数は小さな素因数に分解できるというところが大きなポイントとなる。一般に、最大素因数の数がk以下の正の整数をk-smooth numberと呼びます。では実際に、smooth numberってどれくらいあるのでしょうか？

　以下の問題設定をします。

[1, N]までの正の整数に含まれるk-smooth numberの数がN/2を超えるのはkがいくらのときでしょうか？

このkが小さければ、小さいほどQSを実装する側としては嬉しいわけです。他にも面白い解法がありそうですが、取りあえずパソコンのパワーを信じてがりがりエラトステネスで篩にかけます。普通のエラトステネスの場合は、素数か合成数かをマークしますが、今回は最大の素因数をキャッシュしていきます。

view plainprint?

1. const int MAX_N = 1000;  
2. int fac[MAX_N+1];  
3. map<int, int> smooth_pop;  
4.   
5. int main() {  
6.     memset(fac, 0, sizeof(fac));  
7.     fac[1] = 1;  
8.   
9.     for (int i = 2; i <= MAX_N; i++) {  
10.         if (fac[i] > 0)  
11.             continue;  
12.         for (int j = i; j <= MAX_N; j += i)  
13.             fac[j] = i;  
14.     }  
15.   
16.     /* for debug 
17.     for (int i = 1; i <= MAX_N; i++) 
18.         cout << i << ": " << fac[i] << endl; 
19.     */  
20.   
21.     smooth_pop.clear();  
22.     for (int i = 1; i <= MAX_N; i++)  
23.         ++smooth_pop[fac[i]];  
24.   
25.     int cnt = 0;  
26.     for (map<int, int>::iterator itr = smooth_pop.begin(); itr != smooth_pop.end(); itr++) {  
27.         cnt += itr->second;  
28.         printf("%d %.3lf\n", itr->first, 1.*cnt/MAX_N);  
29.     }  
30.   
31.     return 0;  
32. }  

const int MAX_N = 1000;
int fac[MAX_N+1];
map<int, int> smooth_pop;

int main() {
    memset(fac, 0, sizeof(fac));
    fac[1] = 1;

    for (int i = 2; i <= MAX_N; i++) {
        if (fac[i] > 0)
            continue;
        for (int j = i; j <= MAX_N; j += i)
            fac[j] = i;
    }

    /* for debug
    for (int i = 1; i <= MAX_N; i++)
        cout << i << ": " << fac[i] << endl;
    */

    smooth_pop.clear();
    for (int i = 1; i <= MAX_N; i++)
        ++smooth_pop[fac[i]];

    int cnt = 0;
    for (map<int, int>::iterator itr = smooth_pop.begin(); itr != smooth_pop.end(); itr++) {
        cnt += itr->second;
        printf("%d %.3lf\n", itr->first, 1.*cnt/MAX_N);
    }

    return 0;
}

以下、Nの値と、k-smooth numberの数がN/2以上になるkの値です。

N	k
10^3	43
10^6	3257
10^8	53773

　10^8程度の数をquadratic sieveで分解したい場合は、53773-smooth numberまで調べれば、2個に1個は使える数が入っているということです。篩にかける際に、実際に必要となるメモリのサイズは、p(53773)の二乗程度です。p(x)はxが何番目の素数かを表します。p(x)はここでは素数定理を用いてx / log(x)で概算します。p(53773) = 4937です。4937 * 4937程度の行列なら特別なこと（非ゼロの要素だけもつようなデータ構造）をしなくても、普通にメモリに格納できます。この感覚だと100桁程度の数をQSで分解しようとすると、スパースな行列を扱うためのアルゴリズムが必要になりそうです・・・。

以下、k-smooth numberが{i | i = 1,2 ... N}の集合全体の何%を占めるかを示します。横軸がk、縦軸が占める割合です。予想どおりのグラフの形になりました。上から順にN=1,000、 1,000,000、 100,000,000のときです。