おおきな数字に対して、素数リストを作るためのプログラム

※サイト運営にサーバーは必須です※
～このサイトもエックスサーバーを使用しています～

はじめに

素数のリストを作りたい場合、エラトスネスの篩に基づいたプログラムはよく使われる。

しかし、ある程度（というか、かなり）大きな数を取り扱いたいと考えた場合、ある問題に直面する。

それは、「メモリ不足」だ。

例えば、1000までの素数リストを作る際、1～1000の数字が素数か素数でないか覚えておくための配列を宣言する必要がある。

C#であれば、以下のようなコードで、素数か素数でないか真偽を入れておく配列が宣言される。

bool[] n_list = new bool[1000] ;

※ネット上では、boolでなくintで宣言するプログラムを書いている人もいるが、メモリ的には、真偽の2パターンしか入らないboolの方がメモリの観点で得なはず。

数字が1000までなら、問題は起きないが、これが1000でなく１億とかになると、問題になる。

例えばintを1つ宣言するには、4バイト必要。boolだと、(おそらく)1バイト。

ここでは、bool型で宣言すると考える。
1000までの数字に対応するメモリを用意するのに1KB(1000バイト)
100万までの数字に対応するメモリを用意するのに1MB(1000KB)
10億までの数字に対応するメモリを用意するのに1GB(1000MB)

普通のパソコンでメモリは、2~32GB程度。

実際には、全てのメモリをエラトスネスの篩のために全力投球することはできない。OSを動かすためなどにメモリが必要なはず。

それを考えると、1GB必要となる10億の数字をいっぺんに扱うのは、パソコンによってかなり厳しかったりする。

※タスクマネージャーを立ち上げておくと、どのぐらいのメモリを食っているかリアルタイムで見ることができる。私のパソコン(メモリ4GB)で、20億のエラトスネスの篩のプログラムは一応動いた(最後まで、実行が完了するかどうかは確認していない)。容量は2000MB近く確保されていたので、上の説明はそんなにずれてはいないだろう。

※このぐらいの数字になると、int型で扱える数字の最大値（およそ21億）に対しても注意を払う必要が出てくる。気を付けないとオーバーフローする危険がある。

この記事では、大きな数字でも安定して、エラトスネスの篩が動くようなプログラムを考えたい。

※単純にエラトスネスの篩を実装した場合のプログラムがどんなもんか確認したい場合は、以下の記事を参照
→エラトステネスの篩の素数判定プログラム(c#)

※追記：boolはプログラムによって容量が違っていたりする。0と1だけ扱うのだから1ビットで済みそうな気もするが、実際にはプログラムの言語やバージョンでbool型に必要な容量に違いが存在する。(1byte=8bit。8ビットが集まって1バイト) 。

理由は、ちょっとよくわからないが、プログラムで扱う最低容量が、1バイトだから、bool に対しても、1バイトは最低でも確保されていることがあるらしい。プログラムの言語によって、2バイトだったり、場合によっては、intと同じ４バイトあったりする。4バイト確保されているパターンの場合は、intへの型変換の関係上で確保されていたりするためと思われる。

確実に1バイトで済ませたい場合は、0～255の数字を扱うbyte型を採用するのが一番無難な気がする。

プログラムのアルゴリズム

単純に1万（10000）までの素数を求めたいと考えたとする(n_ini=10000)。

この1万を分割し、100の単位に分ける(dn=100)。

最初の１回目は1～100までの数字を考える(配列はn_list[0]~ n_list[99])。

100の平方根（√）は10なので、2～10までの数字で割り切れなったら、その数字は素数と考える。

※プログラム上、2を特別扱いした方が早くなる。2で割る前処理が終わった後は、3、5、7、9の奇数だけを取り出して、割りきれるかどうか判定している。

2回目は、101～200を考える。

2回目以降から一部処理が変わる。1回目で得られた、素数リスト(p_list)をもとに、素数の探索を行うようになる。

まず、2だけ先に処理する。その後、200の平方根（√）の整数部に対応する14以下の数字で、かつ、素数の数をピックアップする。ここでは、3、5、7、11、13の数字を使って探索をする。

※プログラム上は、101～200がn_list[0]~ n_list[99]に対応。offsetとして100だけゲタをはかして、１回目と同様の処理をする。

これが終わったら、101～200で素数判定を受けた数字だけを、素数のリストに加えていく。

これを１万まで繰り返す。

以上のような感じで計算している。

確保するメモリは、1～10000でなく、1～100までの分があれば問題ない。
※正確には、素数を記録するためなどにもメモリは使われている。

プログラム

プログラムは以下のようになる。
言語はC#。VisualStudio2015のコンソールアプリケーションとして作成。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using static System.Console;
using System.Diagnostics;
using static System.Math;
using System.IO;

//using System.Collections.Generic;

namespace prime_list
{
    class Program
    {
        static void Main(string[] args)
        {

            Stopwatch st = new Stopwatch();
            st.Start();



            int n_ini = 1000000000;//この数字まで、素数を調べる
            int dn = 10000;//刻み幅（分割単位）
            int a = n_ini/dn; //a*dn+b(配列の大きさに相当)
            int b = n_ini % dn;//b余り（bが0でないときの処理はそのうち実装。このプログラムではn_ini/dnが割り切れる場合を想定。つまりb=0）
            int n_root;//ルート
            int index = 1;//「何」番目の素数か記録する

            List<int> p_list = new List<int>();//素数を覚えておくリスト
            p_list.Add(2);//2を素数リストにあらかじめ追加しておく

            bool[] n_list = new bool[dn];//bool の初期はfalse

            //出力準備
            string file_name = "prime_list" + n_ini + ".txt";
            FileStream fs = new FileStream(file_name, FileMode.Create);//上書き保存する
            StreamWriter sw = new StreamWriter(fs);





            //最初の1回目
            //WriteLine(2 + "\r\n");
            sw.Write(index + "\t" +2 + "\r\n");
            index += 1;

            for (int j = 2; 2 * j - 1 < dn; ++j)
            {
                n_list[2* j - 1] = true;
            }


            n_root = (int)Sqrt(dn);

            for (int i=2; (i+1)<= n_root;i+=2)
            {
                if(n_list[i] == false)
                {
                    for (int j = (i + 1); (i + 1) * j - 1 < dn; j+=2)
                    {
                        n_list[(i + 1) * j - 1] = true;
                    }
                }

            }



            for (int k = 2; k < dn; k+=2)
            {
                if (n_list[k] == false)
                {
                    //WriteLine(k+1 + "\r\n");
                    sw.Write(index+"\t"+ Convert.ToString(k + 1) + "\r\n");
                    index += 1;
                    p_list.Add(k+1);
                }
            }

            int offset;
            int max_num;

            //2回目以降
            for (int t = 1; t < a; ++t)
                {

                offset = dn * t;
                max_num= dn * (t+1);
                n_root = (int)Sqrt(max_num);

                //n_list = new bool[dn];//配列の初期化(これで初期化するとメモリ的に損しないか？)
                //n_listの初期化
                for (int j=0; j < dn; ++j)
                {
                    n_list[j] = false;
                }


                //2だけ前処理
                for (int j = (offset/2) + 1; 2 * j - 1  < max_num; ++j)
                {
                    n_list[2 * j - 1 - offset] = true;
                }

                //3以上の素数の処理
                for (int i = 1; p_list[i] <= n_root; ++i)
                {
                    int m;
                    m= (int)Ceiling(((double)(offset+1)/(double)(p_list[i])));
                    if (m % 2 == 0)
                    {
                        m += 1;
                    }
                    for (int j = m; p_list[i] * m - 1 < max_num; m += 2)
                    {
                        n_list[p_list[i] * m - 1 - offset] = true;
                    }
                }

                //結果の処理
                for (int k = 0; k < dn; k += 2)
                    {
                        if (n_list[k] == false)
                        {
                        //WriteLine(k + 1 + offset + "\r\n");
                        sw.Write(index + "\t" + Convert.ToString(k + 1+ offset) + "\r\n");
                        index += 1;
                        p_list.Add(k + 1 + offset);
                        }
                    }

            }            //2回目以降の処理終了


            //最後にbが0でなかったときの処理(未実装)




            st.Stop();
            //WriteLine(st.Elapsed);
            sw.Write("//測定時間は" + st.Elapsed + "\r\n");
            sw.Close();//ファイルを閉じる



        }
    }
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Threading.Tasks;

using static System.Console;

using System.Diagnostics;

using static System.Math;

using System.IO;

//using System.Collections.Generic;

namespace prime_list

{

class Program

{

static void Main(string[] args)

{

Stopwatch st = new Stopwatch();

st.Start();

int n_ini = 1000000000;//この数字まで、素数を調べる

int dn = 10000;//刻み幅（分割単位）

int a = n_ini/dn; //a*dn+b(配列の大きさに相当)

int b = n_ini % dn;//b余り（bが0でないときの処理はそのうち実装。このプログラムではn_ini/dnが割り切れる場合を想定。つまりb=0）

int n_root;//ルート

int index = 1;//「何」番目の素数か記録する

List<int> p_list = new List<int>();//素数を覚えておくリスト

p_list.Add(2);//2を素数リストにあらかじめ追加しておく

bool[] n_list = new bool[dn];//bool の初期はfalse

//出力準備

string file_name = "prime_list" + n_ini + ".txt";

FileStream fs = new FileStream(file_name, FileMode.Create);//上書き保存する

StreamWriter sw = new StreamWriter(fs);

//最初の1回目

//WriteLine(2 + "\r\n");

sw.Write(index + "\t" +2 + "\r\n");

index += 1;

for (int j = 2; 2 * j - 1 < dn; ++j)

{

n_list[2* j - 1] = true;

}

n_root = (int)Sqrt(dn);

for (int i=2; (i+1)<= n_root;i+=2)

{

if(n_list[i] == false)

{

for (int j = (i + 1); (i + 1) * j - 1 < dn; j+=2)

{

n_list[(i + 1) * j - 1] = true;

}

for (int k = 2; k < dn; k+=2)

{

if (n_list[k] == false)

{

//WriteLine(k+1 + "\r\n");

sw.Write(index+"\t"+ Convert.ToString(k + 1) + "\r\n");

index += 1;

p_list.Add(k+1);

}

int offset;

int max_num;

//2回目以降

for (int t = 1; t < a; ++t)

{

offset = dn * t;

max_num= dn * (t+1);

n_root = (int)Sqrt(max_num);

//n_list = new bool[dn];//配列の初期化(これで初期化するとメモリ的に損しないか？)

//n_listの初期化

for (int j=0; j < dn; ++j)

{

n_list[j] = false;

}

//2だけ前処理

for (int j = (offset/2) + 1; 2 * j - 1 < max_num; ++j)

{

n_list[2 * j - 1 - offset] = true;

}

//3以上の素数の処理

for (int i = 1; p_list[i] <= n_root; ++i)

{

int m;

m= (int)Ceiling(((double)(offset+1)/(double)(p_list[i])));

if (m % 2 == 0)

{

m += 1;

}

for (int j = m; p_list[i] * m - 1 < max_num; m += 2)

{

n_list[p_list[i] * m - 1 - offset] = true;

}

//結果の処理

for (int k = 0; k < dn; k += 2)

{

if (n_list[k] == false)

{

//WriteLine(k + 1 + offset + "\r\n");

sw.Write(index + "\t" + Convert.ToString(k + 1+ offset) + "\r\n");

index += 1;

p_list.Add(k + 1 + offset);

}

} //2回目以降の処理終了

//最後にbが0でなかったときの処理(未実装)

st.Stop();

//WriteLine(st.Elapsed);

sw.Write("//測定時間は" + st.Elapsed + "\r\n");

sw.Close();//ファイルを閉じる

}

注意点

スピードを稼ぐために、2を別処理している。
なるべく処理するスピードを上げるために、色々工夫をしているが、その分、プログラムはわかりづらくなっている。プログラムを書いた本人が言うのもなんだが、あまり見たくない。
Falseの場合が素数で、Trueなら素数でないと判定している。篩をかける前は、全てが素数である(False)と考えている
探索したい最大の数(n_ini)に対して、分割単位(dn)がきちんと割り切れる状況を考えている。つまり、n_ini/dnが整数で、bがゼロである。bが0でなかったときの処理は実装していない。
このプログラムに、バグが残っていないと保証できない(一応、他のサイトと結果を比較してずれてはいないことを確認)。