素数列挙について
Competitive Programming Advent Calendar 3日目は、数学っぽい話をしたいと思います。
こんな感じの問題を見たことがあると思います。また問題としてでなくても、解く過程にこのようなサブ問題を解かなければいけない場合もよくあると思います。素数については説明しなくてもいいですよね?
このような問題を素数列挙と呼ぶことにします。素数列挙ができれば、大きい数の素数判定や素因数分解をめっちゃ高速化したり、トーティエント関数、メビウス関数等、数学系のいろんな関数を求めたりできます。最近のもので素数列挙がほぼ必須のものだと
Codeforces Beta Round #86 (Div. 1 Only) C. Double Happiness
ICPC 国内予選 2011 A Chebyshev's Theorem
ProjectEuler Problem 351
ProjectEuler Problem 354
KUPC E Fox Number
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2286:RUPC F Farey Sequence
CodeChef April Cook-Off Challenge A Prime Conjecture
1000万個目の素数を超高速に出力せよ
interviewstreet EQUATIONS
あたりでしょうか。ProjectEulerではしょっちゅう使っている気がします。いろんなことの前提になるのであれば、当然高速に動くコードがあったほうが良いわけです。今回は素数列挙を高速に解く方法について語ります。10^8,10^9くらいまでの素数を高速に列挙orカウントすることを目標に。この野郎知ってるよという節はばんばん飛ばしてください。
今回有為さんはいっぱい○っぱいなので申し訳ないですがJavaで書きます。素数列挙コードは競技コーダーの方はよく書いている(蟻本にもちょろっとある)ので、ぐぐればそれなりに出てきます。
ぐちょく
まずナイーブな方法から書いて行きましょう。ある自然数Nが素数であることを愚直に判定するには2以上√N以下の整数で試し割って、どれも割り切れなければいいわけです。N以下の素数列挙はこれを1〜Nに対し行えばいいだけです。実行時間O(N√N) (定数係数は2/3程度), メモリO(N)です。もしもライブラリがなく、Nが10^5以下で早解き勝負になっているなら、僕は迷わずこれを書きます。素数の個数の上限はもっと抑えられますが、それは後述します。
static int[] sieveNaively(int n) { int[] primes = new int[n]; int p = 0; outer: for(int i = 2;i <= n;i++){ for(int j = 2;j*j <= i;j++){ if(i % j == 0)continue outer; } primes[p++] = i; } return Arrays.copyOf(primes, p); } static int countNaively(int n) { int ct = 0; outer: for(int i = 2;i <= n;i++){ for(int j = 2;j*j <= i;j++){ if(i % j == 0)continue outer; } ct++; } return ct; }
あんまりnが大きいとj*jのところでintの範囲を超えてしまう場合がありますがご愛嬌。
実行時間は以下のとおり。sieveもcountもほぼ同じ実行時間でした。
n | 実行時間(ideone) |
---|---|
10^4 | 0.05s |
10^5 | 0.11s |
10^6 | 1.46s |
10^7 | >15s |
Eratosthenesの篩(ふるい)
http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
ただの試し割りでは空振り率が高いので、既に求まった素数の倍数を選んで撃ち落そうというのがEratosthenesの篩です。2〜Nの整数を整列してフラグを立てます。これに、
列の先頭から見て行って、フラグが立っている数qに対し、q^2以上のqの倍数のフラグを全部消す。
という操作を行うと、残ったフラグの立っている数はすべて素数になります。先頭から見て行って出会うqは必ず素数になっていて、この素数で他の数を試し割りします。r*q(r
static int[] sieveEratos(int n) { BitSet primeset = new BitSet(); primeset.flip(2, n+1); for(int q = primeset.nextSetBit(2);q != -1 && q*q <= n;q = primeset.nextSetBit(q+1)){ for(int i = q*q;i <= n;i += q){ primeset.clear(i); } } int[] primes = new int[primeset.cardinality()]; int p = 0; for(int q = primeset.nextSetBit(2);q != -1;q = primeset.nextSetBit(q+1)){ primes[p++] = q; } return primes; }
n | 実行時間(ideone) |
---|---|
10^6 | 0.11s |
10^7 | 0.66s |
10^8 | 7.95s |
10^9 | MLE |
BitSetが活躍します。nextSetBitで立っているビットだけを走査したり、cardinalityで立っているビットの個数を数えたりできます。countはcardinalityをただ返せばいいですね。
Wheel factorization
http://en.wikipedia.org/wiki/Wheel_factorization
見に行く整数の個数の式で、2,3,5といった小さい素数のところの値が大きくなるのに気づいたでしょうか。実際10^9までの素数の逆数の和は3程度ですが、その半分がこの3つの素数で占められています。これらをあらかじめ除いてから篩にかけちゃおう、というのがWheel factorizationです。
たとえば2,3,5を除くとすると、残るのは30m+1,7,11,13,17,19,23,29の8個になります。30個から8個に減らせたのでそのぶん定数倍高速化が見込めますが、30m+dがqで割り切れるようなmを求めなくてはいけません。これは拡張ユークリッド互除法を使ってO(log q)で求めることができます。
ここでは昔の僕の2だけ消したしょぼい篩を書くにとどめておきます。
public static int[] sieveEratosthenesOld(int n) { int[] ret = new int[(int)(n / Math.log(n)) + (int)(n / (Math.log(n)*Math.log(n))*1.5) ]; ret[0] = 2; int pos = 1; // 0:3 1:5 // 2x+3=n BitSet bs = new BitSet(n/2+1); int sup = (n - 3) / 2; for(int i = bs.nextClearBit(0);i <= sup; i = bs.nextClearBit(i+1)){ int p = 2 * i + 3; ret[pos++] = p; for(int j = i + p;j <= sup;j += p){ bs.set(j); } } return Arrays.copyOf(ret, pos); }
n | 実行時間(ideone) |
---|---|
10^7 | 0.35s |
10^8 | 4.76s |
10^9 | MLE |
retの要素数のところにおぞましい式が出てきていますがこれは次節で。
素数の個数の概算
http://en.wikipedia.org/wiki/Prime_number_theorem
素数定理という有名なものがあります。素数の個数のだいたいの値はこれでわかりますが、あまり精度がよくないので、詳しい値は求められません。
対数積分Li(x)を1回部分積分すると、の項が出てきます。オーダー等を計算するにはだいたいこれで十分でしょう。前節のコードの式は2回部分積分したものの2項目に1ではなく1.5をかけて上限の個数を近似しています。1.5は実験で決めています。
1〜Nから一度に素数列挙しようとすると、やはりO(N/log N)以上のメモリは必要になってしまいます。N=10^9だと若干心配ですね。
自力高速化
BitSetではなくintでビットを持つことにすると、32bitを1単位として持つので、そのなかで複数出てきうるような小さい素数に対してはパターンが数種類程度で収まり、これをまとめて適用することで定数倍高速化できます。現在使っている篩はこのコードです。
public static int[] sieveEratosthenes(int n) { if(n <= 32){ int[] primes = {2,3,5,7,11,13,17,19,23,29,31}; for(int i = 0;i < primes.length;i++){ if(n < primes[i]){ return Arrays.copyOf(primes, i); } } return primes; } int u = n+32; double lu = Math.log(u); int[] ret = new int[(int)(u/lu+u/lu/lu*1.5)]; ret[0] = 2; int pos = 1; int[] isp = new int[(n+1)/32/2+1]; int sup = (n+1)/32/2+1; int[] tprimes = {3,5,7,11,13,17,19,23,29,31}; for(int tp : tprimes){ ret[pos++] = tp; int[] ptn = new int[tp]; for(int i = (tp-3)/2;i < tp<<5;i+=tp)ptn[i>>5] |= 1<<(i&31); for(int i = 0;i < tp;i++){ for(int j = i;j < sup;j += tp)isp[j] |= ptn[i]; } } // 3,5,7 // 2x+3=n int[] magic = {0, 1, 23, 2, 29, 24, 19, 3, 30, 27, 25, 11, 20, 8, 4, 13, 31, 22, 28, 18, 26, 10, 7, 12, 21, 17, 9, 6, 16, 5, 15, 14}; int h = n/2; for(int i = 0;i < sup;i++){ for(int j = ~isp[i];j != 0;j &= j-1){ int pp = i<<5|magic[(j&-j)*0x076be629>>>27]; int p = 2*pp+3; if(p > n)break; ret[pos++] = p; for(int q = pp;q <= h;q += p)isp[q>>5] |= 1<<(q&31); } } return Arrays.copyOf(ret, pos); }
n | 実行時間(ideone) |
---|---|
10^7 | 0.13s |
10^8 | 1.72s |
10^9 | MLE |
途中に出てくるmagicとかいうのは2^nからnを高速に得るための黒魔術です。
http://chessprogramming.wikispaces.com/De+Bruijn+sequence
区間篩
[L,L+B]の区間上での篩をしたいとき、√(L+B)までの素数が求まっていれば篩うことができます。Lが大きくて、Bが小さい場合に使えるテクニックです。区間篩自体はO(Bloglog(L+B))でしょうか。
public static BitSet sieveBySegment(long l, int b){ BitSet com = new BitSet(); com.set(0, l+1); for(int p : primes){ if((long)p*p > l+b)break; for(int i = (int)(p-1-(l-1)%p);i <= b;i+=p){ com.clear(i); } } return com; }
さあ、これまでは前座、前座あかりです。
Atkinの篩
http://en.wikipedia.org/wiki/Sieve_of_Atkin
Eratosthenesより計算量の良いものとしてAtkinの篩があります。実行時間はO(N/loglogN)です。
これは次の事実を利用しています。nはsquarefree(同じ素因数を2個以上含まない)として、
が素数⇔の正の解(x,y)が奇数個存在する。
が素数⇔の正の解(x,y)が奇数個存在する。
が素数⇔の解(x,y) (x>y>0)が奇数個存在する。
これの証明は単項イデアル整域等出てきてお見せ出来なかったので端折ります。
つまりは、2,3,5以外、この(x,y)をN以下について全部列挙して、squarefreeなやつだけ選べばいいのです。(x,y)の個数は高々O(N).
Spaghetti Sourceまんまですが、Atkinの篩にはあまり思い入れはないのでいいでしょう。
public static BitSet sieveAtkin(int n) { int sqn = (int)Math.sqrt(n); BitSet isp = new BitSet(n); int b; for (int z = 1; z <= 5; z += 4) { for (int y = z; y <= sqn; y += 6) { for (int x = 1; x <= sqn && (b = 4 * x * x + y * y) <= n; ++x)isp.flip(b); for (int x = y + 1; x <= sqn && (b = 3 * x * x - y * y) <= n; x += 2)isp.flip(b); } } for (int z = 2; z <= 4; z += 2) { for (int y = z; y <= sqn; y += 6) { for (int x = 1; x <= sqn && (b = 3 * x * x + y * y) <= n; x += 2)isp.flip(b); for (int x = y + 1; x <= sqn && (b = 3 * x * x - y * y) <= n; x += 2)isp.flip(b); } } for (int y = 3; y <= sqn; y += 6) { for (int z = 1; z <= 2; ++z) { for (int x = z; x <= sqn && (b = 4 * x * x + y * y) <= n; x += 3)isp.flip(b); } } for(int i = isp.nextSetBit(5);i != -1 && i <= sqn;i = isp.nextSetBit(i+2)){ int zz = i*i; for (int k = zz; k <= n; k += zz)isp.clear(k); } isp.set(2); isp.set(3); return isp; }
n | 実行時間(ideone) |
---|---|
10^7 | 0.22s |
10^8 | 3.39s |
10^9 | MLE |
Atkinの篩進化形
http://scholar.google.co.jp/scholar?cluster=17579060183460050316&hl=ja&as_sdt=0,5
前節のAtkinの篩は10^9でMLEしてしまいました。前にも書いたように1〜Nを一度に篩おうとするとそれなりのメモリが必要なのです。これから書くものは実行時間のオーダーそのままで、メモリがなんとO(√N)でできる篩です。メモリだけ減っても速くならないというわけではなく、既存のものもL1,L2キャッシュに乗りやすくなるので高速化が期待できます。
区間[60L,60(L+B))の素数をAtkinの篩で篩うことを考えます。当然、の解を求めることになりますが、xy平面上でみると、これらは2次曲線の軌跡からなる帯のような形をしています。この中を、特殊解から、なら(-15,0)または(0,30)のステップ幅で移動して解を列挙していきます。これで使用するメモリはO(B). B=60√Nとして選び、[60,60(1+B)),[60(1+B),60(1+2B)),..について区間篩を繰り返すと、全体での実行時間オーダーがちょうどO(N/loglogN)に、メモリオーダーがになるようです。
特殊解はあらかじめ計算できておくのでその辺は埋め込みです。でもあまり実行速度には影響ありませんでした。
また、squarefreeの判定のために√Nまでの素数が必要なので、別にsieveしておきます。これは自分自身を呼んでも構いません。
// 4x^2+y^2=d(mod 60)の特殊解(x,y)=(f,g) static final int[][] DFG1 = {{1,0,1},{1,0,11},{1,0,19},{1,0,29},{1,2,15},{1,3,5},{1,3,25},{1,5,9},{1,5,21},{1,7,15},{1,8,15},{1,10,9},{1,10,21},{1,12,5},{1,12,25},{1,13,15},{13,1,3},{13,1,27},{13,4,3},{13,4,27},{13,6,7},{13,6,13},{13,6,17},{13,6,23},{13,9,7},{13,9,13},{13,9,17},{13,9,23},{13,11,3},{13,11,27},{13,14,3},{13,14,27},{17,2,1},{17,2,11},{17,2,19},{17,2,29},{17,7,1},{17,7,11},{17,7,19},{17,7,29},{17,8,1},{17,8,11},{17,8,19},{17,8,29},{17,13,1},{17,13,11},{17,13,19},{17,13,29},{29,1,5},{29,1,25},{29,4,5},{29,4,25},{29,5,7},{29,5,13},{29,5,17},{29,5,23},{29,10,7},{29,10,13},{29,10,17},{29,10,23},{29,11,5},{29,11,25},{29,14,5},{29,14,25},{37,2,9},{37,2,21},{37,3,1},{37,3,11},{37,3,19},{37,3,29},{37,7,9},{37,7,21},{37,8,9},{37,8,21},{37,12,1},{37,12,11},{37,12,19},{37,12,29},{37,13,9},{37,13,21},{41,2,5},{41,2,25},{41,5,1},{41,5,11},{41,5,19},{41,5,29},{41,7,5},{41,7,25},{41,8,5},{41,8,25},{41,10,1},{41,10,11},{41,10,19},{41,10,29},{41,13,5},{41,13,25},{49,0,7},{49,0,13},{49,0,17},{49,0,23},{49,1,15},{49,4,15},{49,5,3},{49,5,27},{49,6,5},{49,6,25},{49,9,5},{49,9,25},{49,10,3},{49,10,27},{49,11,15},{49,14,15},{53,1,7},{53,1,13},{53,1,17},{53,1,23},{53,4,7},{53,4,13},{53,4,17},{53,4,23},{53,11,7},{53,11,13},{53,11,17},{53,11,23},{53,14,7},{53,14,13},{53,14,17},{53,14,23}}; // 3x^2+y^2=d(mod 60)の特殊解(x,y)=(f,g) static final int[][] DFG2 = {{7,1,2},{7,1,8},{7,1,22},{7,1,28},{7,3,10},{7,3,20},{7,7,10},{7,7,20},{7,9,2},{7,9,8},{7,9,22},{7,9,28},{19,1,4},{19,1,14},{19,1,16},{19,1,26},{19,5,2},{19,5,8},{19,5,22},{19,5,28},{19,9,4},{19,9,14},{19,9,16},{19,9,26},{31,3,2},{31,3,8},{31,3,22},{31,3,28},{31,5,4},{31,5,14},{31,5,16},{31,5,26},{31,7,2},{31,7,8},{31,7,22},{31,7,28},{43,1,10},{43,1,20},{43,3,4},{43,3,14},{43,3,16},{43,3,26},{43,7,4},{43,7,14},{43,7,16},{43,7,26},{43,9,10},{43,9,20}}; // 3x^2-y^2=d(mod 60)の特殊解(x,y)=(f,g) static final int[][] DFG3 = {{11,0,7},{11,0,13},{11,0,17},{11,0,23},{11,2,1},{11,2,11},{11,2,19},{11,2,29},{11,3,4},{11,3,14},{11,3,16},{11,3,26},{11,5,2},{11,5,8},{11,5,22},{11,5,28},{11,7,4},{11,7,14},{11,7,16},{11,7,26},{11,8,1},{11,8,11},{11,8,19},{11,8,29},{23,1,10},{23,1,20},{23,2,7},{23,2,13},{23,2,17},{23,2,23},{23,3,2},{23,3,8},{23,3,22},{23,3,28},{23,4,5},{23,4,25},{23,6,5},{23,6,25},{23,7,2},{23,7,8},{23,7,22},{23,7,28},{23,8,7},{23,8,13},{23,8,17},{23,8,23},{23,9,10},{23,9,20},{47,1,4},{47,1,14},{47,1,16},{47,1,26},{47,2,5},{47,2,25},{47,3,10},{47,3,20},{47,4,1},{47,4,11},{47,4,19},{47,4,29},{47,6,1},{47,6,11},{47,6,19},{47,6,29},{47,7,10},{47,7,20},{47,8,5},{47,8,25},{47,9,4},{47,9,14},{47,9,16},{47,9,26},{59,0,1},{59,0,11},{59,0,19},{59,0,29},{59,1,2},{59,1,8},{59,1,22},{59,1,28},{59,4,7},{59,4,13},{59,4,17},{59,4,23},{59,5,4},{59,5,14},{59,5,16},{59,5,26},{59,6,7},{59,6,13},{59,6,17},{59,6,23},{59,9,2},{59,9,8},{59,9,22},{59,9,28}}; static final int[] DALL = {1,7,11,13,17,19,23,29,31,37,41,43,47,49,53,59}; static final int[] D1 = {1,13,17,29,37,41,49,53}; static final int[] D2 = {7,19,31,43}; static final int[] D3 = {11,23,47,59}; static final int[] U60 = new int[]{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59}; static int[][] segs = new int[60][]; public static int[] sieveOfAtkinBernstain(int N) { if(N <= 60){ int ind = Arrays.binarySearch(U60, N); if(ind < 0)ind = -ind-2; return Arrays.copyOf(U60, ind+1); } int B = 60*(int)Math.sqrt(N); int[] primes = sieveOfAtkinBernstain((int)Math.sqrt(N)); int u = N+32; double lu = Math.log(u); int[] ret = new int[(int)(u/lu+u/lu/lu*1.5)]; int r = U60.length; System.arraycopy(U60, 0, ret, 0, U60.length); for(int d : DALL)segs[d] = new int[(B>>>5)+1]; for(int L = 1;L <= N/60;L+=B){ for(int d : DALL)Arrays.fill(segs[d], 0); // sieve! for(int[] D : DFG1)enum1(D[1], D[2], D[0], L, B); for(int[] D : DFG2)enum2(D[1], D[2], D[0], L, B); for(int[] D : DFG3)enum3(D[1], D[2], D[0], L, B); // sieve non-squarefree for(int p : primes){ if((long)p*p > 60*(L+B))break; if(p >= 7){ // kp^2=60(L+x)+d // A*p^2-X*60=60L+d long b = -exGCDSimple(p*p, 60); int p2 = p*p; if(b < 0)b += p2; for(int d : DALL){ int x = (int)(b*(60*L+d)%p2); for(;x < B;x += p2){ segs[d][x>>>5] &= ~(1<<(x&31)); } } } } // arrange primes inner: for(int j = 0;j < (B>>>5)+1;j++){ for(int x = 0;x < 32;x++){ long base = 60*(L+x+(j<<5)); for(int d : DALL){ if(base+d > N)break inner; if(segs[d][j]<<31-x<0){ ret[r++] = 60*(L+x+(j<<5))+d; } } } } } ret = Arrays.copyOf(ret, r); return ret; } public static int countOfAtkinBernstain(int N) { if(N <= 60){ int ind = Arrays.binarySearch(U60, N); return ind < 0 ? -ind-1 : ind; } int B = 60*(int)Math.sqrt(N); int[] primes = sieveOfAtkinBernstain((int)Math.sqrt(N)); int r = U60.length; for(int d : DALL)segs[d] = new int[(B>>>5)+1]; for(long L = 1;60*L <= N;L+=B){ for(int d : DALL)Arrays.fill(segs[d], 0); // sieve! for(int[] D : DFG1)enum1(D[1], D[2], D[0], L, B); for(int[] D : DFG2)enum2(D[1], D[2], D[0], L, B); for(int[] D : DFG3)enum3(D[1], D[2], D[0], L, B); // sieve non-squarefree for(int p : primes){ if((long)p*p > 60*(L+B))break; if(p >= 7){ // kp^2=60(L+x)+d // A*p^2-X*60=60L+d long b = -exGCD(p*p, 60); int p2 = p*p; if(b < 0)b += p2; for(int d : DALL){ int x = (int)(b*(60*L+d)%p2); for(;x < B;x += p2){ segs[d][x>>>5] &= ~(1<<(x&31)); } } } } // count primes for(int j = 0;j < (B>>>5)+1;j++){ long base = 60*(L+(j<<5)+31); for(int d : DALL){ if(base+d > N){ // 60(L+32j+x)+d <= N // x <= (N-d)/60-L-32j int sup = (int)((N-d)/60-L-(j<<5)+1); if(sup <= 0)break; r += Integer.bitCount(segs[d][j]&(1<<sup)-1); }else{ r += Integer.bitCount(segs[d][j]); } } } } return r; } public static long exGCDSimple(long a, long b) { long r = 0, s = 1; while(b > 0){ long c = a / b; long d; d = a; a = b; b = d % b; d = r; r = s; s = d - c * s; } return r; } static void enum1(int f, int g, int d, long L, int B) { int x = f, y0 = g; int k0 = (4*f*f+g*g-d)/60; while(k0<L+B){ k0 += 2*x + 15; x += 15; } while(true){ x -= 15; k0 -= 2*x+15; if(x <= 0)break; while(k0 < L){ k0 += y0+15; y0 += 30; } int k = k0, y = y0; while(k < L+B){ segs[d][(int)(k-L)>>>5] ^= 1<<((k-L)&31); k += y+15; y += 30; } } } static void enum2(int f, int g, int d, long L, int B) { int x = f, y0 = g; int k0 = (3*f*f+g*g-d)/60; while(k0<L+B){ k0 += x + 5; x += 10; } while(true){ x -= 10; k0 -= x+5; if(x <= 0)break; while(k0 < L){ k0 += y0+15; y0 += 30; } int k = k0, y = y0; while(k < L+B){ segs[d][(int)(k-L)>>>5] ^= 1<<((k-L)&31); k += y+15; y += 30; } } } static void enum3(int f, int g, int d, long L, int B) { int x = f, y0 = g; int k0 = (3*f*f-g*g-d)/60; outer: while(true){ while(k0 >= L+B){ if(x <= y0)break outer; k0 -= y0 + 15; y0 += 30; } int k = k0, y = y0; while(k>=L && y<x){ segs[d][(int)(k-L)>>>5] ^= 1<<((k-L)&31); k -= y+15; y += 30; } k0 += x+5; x += 10; } }
sieve
n | 実行時間(ideone) |
---|---|
10^7 | 0.15s |
10^8 | 0.76s |
10^9 | MLE |
digest(結果を格納せず、何も書き出さないWriterに逐次書きだす)
n | 実行時間(ideone) |
---|---|
10^7 | 0.15s |
10^8 | 0.75s |
10^9 | 6.76s |
count
n | 実行時間(ideone) |
---|---|
10^7 | 0.09s |
10^8 | 0.33s |
10^9 | 2.83s |
こんな感じです。論文を見て自前で実装したので、もっと高速化できると思います。
2次式の篩
呼び方がよくわからないので・・ProjectEulerで稀に、2次式の形の篩を要求されることがあります。これはO(NloglogN)でできます。考えてみてくださいね。
http://projecteuler.net/problem=216