このページの記事一覧
- 2016/07/10 【Java】フィボナッチ数(数列)を求める
- 2016/06/07 【Java】素数判定③（アトキンの篩(ふるい)）
- 2016/06/06 【Java】素数判定③（サンダラムの篩(ふるい)）
- 2016/06/05 【Java】素数判定②（エラトステネスの篩(ふるい)）
« prev next »
【Java】フィボナッチ数(数列)を求める
2016/07/10 Sun [edit]
「エラトステネスの篩(ふるい)」同様、数列的なテーマとしてもよく取り上げられる「フィボナッチ数」。「黄金比(黄金数)」とも深い関係があり、実は日常的にもよく使われているとも言われている（自然現象からデザイン的なものまで）。
まぁ、それはさておき(笑)、プログラミングで計算するのは意外と簡単だったりする。せっかくなので、いくつかの方法で求めてみよう。
■動的計画法でフィボナッチ数列を求める
まずは漸化式を確認してみよう。Wikipedia にも載っているが、もう少しプログラミング用に変形した方がわかりやすいかも知れない。
(漸化式)
F0 = 0,
F1 = 1,
Fn = Fn-1 + Fn-2 (n ≧ 0)
F0 = 0,
F1 = 1,
Fn = Fn-1 + Fn-2 (n ≧ 0)
●動的計画法でフィボナッチ数列を作る
import java.util.Arrays;
int n = 20; //n 上限
//動的計画法
long[] dp = new long[n + 1];
dp[0] = 0;
dp[1] = 1;
for (int i = 2; i <= n; i++) {
dp[i] = dp[i - 1] + dp[i - 2];
}
//F(n)を表示
System.out.println("F(" + n + ") = " + dp[n]);
//数列を表示
System.out.println(Arrays.toString(dp));
F(20) = 6765
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
一般的には自然数（0 は含まない整数）で書かれているものも多いが、コンピュータで数列を作る場合は 0 も入れておいた方が何かと便利だろう。インデクスに対応できる利点もある。漸化式をそのままコーディングした感じになるので、説明はいらないだろう。
気をつけないといけないのは、フィボナッチ数は非常に大きな値となるので、long 型で実装した場合、n = 93 でオーバーフローすることだ。それ以上の値の場合は、後述する BigDecimal（任意精度実数）を使うのも良いだろう。ちなみに、BigInteger（任意精度整数）でも構わないが、BigDecimal の方がわずかに速いようだ。
■再帰でフィボナッチ数(数列)を求める
次に再帰で フィボナッチ数 F(n) を求めてみよう。再帰でもとてもシンプルに書けることがわかる。
●再帰でフィボナッチ数を求める
int n = 20;
//F(n)を表示
System.out.println("F(" + n + ") = " + fibonacci(n));
//再帰関数
long fibonacci(int n) {
if (n <= 1) {
return n;
}
return fibonacci(n - 1) + fibonacci(n - 2);
}
F(20) = 6765
ただ、上記の方法だと n = 40 ほどになると非常に重くなる。なぜなら、再帰で足しあわせている部分「fibonacci(n - 1) + fibonacci(n - 2)」は入れ子になっていて、同じ計算を何度も繰り返しているため、無駄が多いからだ。この関数は引数が同じなら戻値も同じになるので、メモ化することで高速化を図れる。またメモは数列にもなるので表示してみよう。
●メモ化再帰でフィボナッチ数列を作る
import java.util.Arrays;
int n = 20;
long[] fib = new long[n + 1]; //メモ(グローバルに置く)
//F(n)を表示
System.out.println("F(" + n + ") = " + fibonacci(n));
//数列を表示
System.out.println(Arrays.toString(fib));
//メモ化再帰
long fibonacci(int n) {
if (n <= 1) {
return fib[n] = n;
}
if (fib[n] > 0) {
return fib[n];
}
return fib[n] = fibonacci(n - 1) + fibonacci(n - 2);
}
F(20) = 6765
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
この方法なら n = 40 以上でも一瞬で求められる。
ただし、上記の動的計画法でも再帰でも、long で実装した場合、 n = 93 でオーバーフローになるので気をつけよう。
■ビネの公式でフィボナッチ数を求める
もう１つ、一般項を求める式（ビネの公式）を用いてフィボナッチ数を求めてみよう。
●ビネの公式でフィボナッチ数を求める
final double GOLDEN_RATIO = 1.618033988749895; //黄金比
final double SQRT5 = Math.sqrt(5); //√5
int n = 20;
//ビネの公式
//(上の式) n = 72 から誤差が出る
double ans1 = (Math.pow(1.0 + SQRT5, n) - Math.pow(1.0 - SQRT5, n)) / (Math.pow(2, n) * SQRT5);
System.out.println("F(" + n + ") = " + (long)ans1 + " (" + ans1 + ")");
//(下の式) n = 71 から誤差が出る
double ans2 = Math.floor(Math.pow(GOLDEN_RATIO, n) / SQRT5 + 0.5);
System.out.println("F(" + n + ") = " + (long)ans2 + " (" + ans2 + ")");
F(20) = 6765 (6765.000000000005)
F(20) = 6765 (6765.0)
F(20) = 6765 (6765.0)
ただし、これらは近似式のため誤差が出るので注意しよう。 double で実装した場合、上の式では n = 72 から、下の式では n = 71 から誤差が出る。式の解説は Wikipedia に載っているので参照して欲しいが、プログラミングするなら下の方の式が使いやすい。カッコの上の部分が無いような記号は床関数の記号で、小数点切り捨てと考えれば、Math.floor() でも整数型のキャストでも良いだろう（正の整数の場合）。よくわからない記号は「数学記号の表」などで調べると良い。
ビネの公式での誤差の出る範囲を抑える１つの方法は、BigDecimal のような任意精度実数を使って、桁の精度を上げることである。試しに黄金比や √5 の値を 50 桁まで上げて、先ほどの式（下の方の式、上の方の式は上手くいかない）を BigDecimal に書き換えてみると、n = 231 までは誤差が出ない。
●ビネの公式を BigDecimal で書いてみる(小数点以下 50桁)
import java.math.BigDecimal;
import java.math.RoundingMode;
int n = 231; //232 から誤差が出る
/** 任意精度実数の黄金比 (50桁) */
final BigDecimal BIG_DECIMAL_GOLDEN = new BigDecimal("1.61803398874989484820458683436563811772030917980576");
/** 任意精度実数の √5 (50桁) */
final BigDecimal BIG_DECIMAL_SQRT5 = new BigDecimal("2.23606797749978969640917366873127623544061835961152");
BigDecimal ans = BIG_DECIMAL_GOLDEN.pow(n).divide(BIG_DECIMAL_SQRT5, 50, RoundingMode.FLOOR).add(new BigDecimal(0.5));
System.out.println("F(" + n + ") = " + ans.setScale(0, RoundingMode.FLOOR) + " (" + ans + ")");
F(231) = 844617150046923109759866426342507997914076076194 (844617150046923109759866426342507997914076076194.15703878686668096410036906990688047595887460281556)
ビネの公式も黄金比や √5 の精度を上げれば、もっと大きな n でも誤差が出ないようにできると思うが、実際には計算が非常に重くなるのであまり実用的ではないかも知れない。その場合、はじめに書いた動的計画法を BigDecimal に書き換えた方が実行速度も速く、誤差もない。
●はじめの動的計画法を BigDecimal に書き換えたもの
import java.math.BigDecimal;
import java.util.Arrays;
int n = 100;
BigDecimal[] dp = new BigDecimal[n + 1];
dp[0] = BigDecimal.ZERO;
dp[1] = BigDecimal.ONE;
for (int i = 2; i <= n; i++) {
dp[i] = dp[i - 1].add(dp[i - 2]);
}
//F(n)を表示
System.out.println("F(" + n + ") = " + dp[n]);
//数列を表示
System.out.println(Arrays.toString(dp));
F(100) = 354224848179261915075
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073, 4807526976, 7778742049, 12586269025, 20365011074, 32951280099, 53316291173, 86267571272, 139583862445, 225851433717, 365435296162, 591286729879, 956722026041, 1548008755920, 2504730781961, 4052739537881, 6557470319842, 10610209857723, 17167680177565, 27777890035288, 44945570212853, 72723460248141, 117669030460994, 190392490709135, 308061521170129, 498454011879264, 806515533049393, 1304969544928657, 2111485077978050, 3416454622906707, 5527939700884757, 8944394323791464, 14472334024676221, 23416728348467685, 37889062373143906, 61305790721611591, 99194853094755497, 160500643816367088, 259695496911122585, 420196140727489673, 679891637638612258, 1100087778366101931, 1779979416004714189, 2880067194370816120, 4660046610375530309, 7540113804746346429, 12200160415121876738, 19740274219868223167, 31940434634990099905, 51680708854858323072, 83621143489848422977, 135301852344706746049, 218922995834555169026, 354224848179261915075]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437, 701408733, 1134903170, 1836311903, 2971215073, 4807526976, 7778742049, 12586269025, 20365011074, 32951280099, 53316291173, 86267571272, 139583862445, 225851433717, 365435296162, 591286729879, 956722026041, 1548008755920, 2504730781961, 4052739537881, 6557470319842, 10610209857723, 17167680177565, 27777890035288, 44945570212853, 72723460248141, 117669030460994, 190392490709135, 308061521170129, 498454011879264, 806515533049393, 1304969544928657, 2111485077978050, 3416454622906707, 5527939700884757, 8944394323791464, 14472334024676221, 23416728348467685, 37889062373143906, 61305790721611591, 99194853094755497, 160500643816367088, 259695496911122585, 420196140727489673, 679891637638612258, 1100087778366101931, 1779979416004714189, 2880067194370816120, 4660046610375530309, 7540113804746346429, 12200160415121876738, 19740274219868223167, 31940434634990099905, 51680708854858323072, 83621143489848422977, 135301852344706746049, 218922995834555169026, 354224848179261915075]
最後に long, double (ビネの公式：下の式), BigDecimal で計算した結果をダンプしてみよう。オーバーフローや誤差を確認しておけば、どれを使うかの目安になるかも知れない。値を確認するには以下のようなサイトを参考にすれば良い。
（参考）[数列]
・１００番目までのフィボナッチ数列
・１番目から１０００番目までのフィボナッチ数列
F(0)
long = 0
double = 0
BigDecimal = 0
F(1)
long = 1
double = 1
BigDecimal = 1
F(2)
long = 1
double = 1
BigDecimal = 1
F(3)
long = 2
double = 2
BigDecimal = 2
F(4)
long = 3
double = 3
BigDecimal = 3
F(5)
long = 5
double = 5
BigDecimal = 5
F(6)
long = 8
double = 8
BigDecimal = 8
F(7)
long = 13
double = 13
BigDecimal = 13
F(8)
long = 21
double = 21
BigDecimal = 21
F(9)
long = 34
double = 34
BigDecimal = 34
F(10)
long = 55
double = 55
BigDecimal = 55
F(11)
long = 89
double = 89
BigDecimal = 89
F(12)
long = 144
double = 144
BigDecimal = 144
F(13)
long = 233
double = 233
BigDecimal = 233
F(14)
long = 377
double = 377
BigDecimal = 377
F(15)
long = 610
double = 610
BigDecimal = 610
F(16)
long = 987
double = 987
BigDecimal = 987
F(17)
long = 1597
double = 1597
BigDecimal = 1597
F(18)
long = 2584
double = 2584
BigDecimal = 2584
F(19)
long = 4181
double = 4181
BigDecimal = 4181
F(20)
long = 6765
double = 6765
BigDecimal = 6765
F(21)
long = 10946
double = 10946
BigDecimal = 10946
F(22)
long = 17711
double = 17711
BigDecimal = 17711
F(23)
long = 28657
double = 28657
BigDecimal = 28657
F(24)
long = 46368
double = 46368
BigDecimal = 46368
F(25)
long = 75025
double = 75025
BigDecimal = 75025
F(26)
long = 121393
double = 121393
BigDecimal = 121393
F(27)
long = 196418
double = 196418
BigDecimal = 196418
F(28)
long = 317811
double = 317811
BigDecimal = 317811
F(29)
long = 514229
double = 514229
BigDecimal = 514229
F(30)
long = 832040
double = 832040
BigDecimal = 832040
F(31)
long = 1346269
double = 1346269
BigDecimal = 1346269
F(32)
long = 2178309
double = 2178309
BigDecimal = 2178309
F(33)
long = 3524578
double = 3524578
BigDecimal = 3524578
F(34)
long = 5702887
double = 5702887
BigDecimal = 5702887
F(35)
long = 9227465
double = 9227465
BigDecimal = 9227465
F(36)
long = 14930352
double = 14930352
BigDecimal = 14930352
F(37)
long = 24157817
double = 24157817
BigDecimal = 24157817
F(38)
long = 39088169
double = 39088169
BigDecimal = 39088169
F(39)
long = 63245986
double = 63245986
BigDecimal = 63245986
F(40)
long = 102334155
double = 102334155
BigDecimal = 102334155
F(41)
long = 165580141
double = 165580141
BigDecimal = 165580141
F(42)
long = 267914296
double = 267914296
BigDecimal = 267914296
F(43)
long = 433494437
double = 433494437
BigDecimal = 433494437
F(44)
long = 701408733
double = 701408733
BigDecimal = 701408733
F(45)
long = 1134903170
double = 1134903170
BigDecimal = 1134903170
F(46)
long = 1836311903
double = 1836311903
BigDecimal = 1836311903
F(47)
long = 2971215073
double = 2971215073
BigDecimal = 2971215073
F(48)
long = 4807526976
double = 4807526976
BigDecimal = 4807526976
F(49)
long = 7778742049
double = 7778742049
BigDecimal = 7778742049
F(50)
long = 12586269025
double = 12586269025
BigDecimal = 12586269025
F(51)
long = 20365011074
double = 20365011074
BigDecimal = 20365011074
F(52)
long = 32951280099
double = 32951280099
BigDecimal = 32951280099
F(53)
long = 53316291173
double = 53316291173
BigDecimal = 53316291173
F(54)
long = 86267571272
double = 86267571272
BigDecimal = 86267571272
F(55)
long = 139583862445
double = 139583862445
BigDecimal = 139583862445
F(56)
long = 225851433717
double = 225851433717
BigDecimal = 225851433717
F(57)
long = 365435296162
double = 365435296162
BigDecimal = 365435296162
F(58)
long = 591286729879
double = 591286729879
BigDecimal = 591286729879
F(59)
long = 956722026041
double = 956722026041
BigDecimal = 956722026041
F(60)
long = 1548008755920
double = 1548008755920
BigDecimal = 1548008755920
F(61)
long = 2504730781961
double = 2504730781961
BigDecimal = 2504730781961
F(62)
long = 4052739537881
double = 4052739537881
BigDecimal = 4052739537881
F(63)
long = 6557470319842
double = 6557470319842
BigDecimal = 6557470319842
F(64)
long = 10610209857723
double = 10610209857723
BigDecimal = 10610209857723
F(65)
long = 17167680177565
double = 17167680177565
BigDecimal = 17167680177565
F(66)
long = 27777890035288
double = 27777890035288
BigDecimal = 27777890035288
F(67)
long = 44945570212853
double = 44945570212853
BigDecimal = 44945570212853
F(68)
long = 72723460248141
double = 72723460248141
BigDecimal = 72723460248141
F(69)
long = 117669030460994
double = 117669030460994
BigDecimal = 117669030460994
F(70)
long = 190392490709135
double = 190392490709135
BigDecimal = 190392490709135
F(71)
long = 308061521170129
double = 308061521170130 //← ここから誤差が出る
BigDecimal = 308061521170129
F(72)
long = 498454011879264
double = 498454011879265
BigDecimal = 498454011879264
F(73)
long = 806515533049393
double = 806515533049395
BigDecimal = 806515533049393
F(74)
long = 1304969544928657
double = 1304969544928660
BigDecimal = 1304969544928657
F(75)
long = 2111485077978050
double = 2111485077978055
BigDecimal = 2111485077978050
F(76)
long = 3416454622906707
double = 3416454622906716
BigDecimal = 3416454622906707
F(77)
long = 5527939700884757
double = 5527939700884772
BigDecimal = 5527939700884757
F(78)
long = 8944394323791464
double = 8944394323791488
BigDecimal = 8944394323791464
F(79)
long = 14472334024676221
double = 14472334024676260
BigDecimal = 14472334024676221
F(80)
long = 23416728348467685
double = 23416728348467744
BigDecimal = 23416728348467685
F(81)
long = 37889062373143906
double = 37889062373144008
BigDecimal = 37889062373143906
F(82)
long = 61305790721611591
double = 61305790721611752
BigDecimal = 61305790721611591
F(83)
long = 99194853094755497
double = 99194853094755776
BigDecimal = 99194853094755497
F(84)
long = 160500643816367088
double = 160500643816367552
BigDecimal = 160500643816367088
F(85)
long = 259695496911122585
double = 259695496911123328
BigDecimal = 259695496911122585
F(86)
long = 420196140727489673
double = 420196140727490880
BigDecimal = 420196140727489673
F(87)
long = 679891637638612258
double = 679891637638614270
BigDecimal = 679891637638612258
F(88)
long = 1100087778366101931
double = 1100087778366105090
BigDecimal = 1100087778366101931
F(89)
long = 1779979416004714189
double = 1779979416004719360
BigDecimal = 1779979416004714189
F(90)
long = 2880067194370816120
double = 2880067194370824700
BigDecimal = 2880067194370816120
F(91)
long = 4660046610375530309
double = 4660046610375544800
BigDecimal = 4660046610375530309
F(92)
long = 7540113804746346429
double = 7540113804746369000
BigDecimal = 7540113804746346429
F(93)
long = -6246583658587674878 //← ここでオーバーフロー
double = 12200160415121914000
BigDecimal = 12200160415121876738
F(94)
long = 1293530146158671551
double = 19740274219868283000
BigDecimal = 19740274219868223167
F(95)
long = -4953053512429003327
double = 31940434634990200000
BigDecimal = 31940434634990099905
F(96)
long = -3659523366270331776
double = 51680708854858490000
BigDecimal = 51680708854858323072
F(97)
long = -8612576878699335103
double = 83621143489848690000
BigDecimal = 83621143489848422977
F(98)
long = 6174643828739884737
double = 135301852344707190000
BigDecimal = 135301852344706746049
F(99)
long = -2437933049959450366
double = 218922995834555900000
BigDecimal = 218922995834555169026
F(100)
long = 3736710778780434371
double = 354224848179263100000
BigDecimal = 354224848179261915075
（関連記事）
【Java】任意精度整数(多倍長整数)で演算する
■参考になる書籍
tb: 0
cm: --
【Java】素数判定③（アトキンの篩(ふるい)）
2016/06/07 Tue [edit]
「より現代的なアルゴリズム」(Wikipedia)と書いてあるので、たぶんそうなのだろう(笑)。残念ながら英語版しかないのだが、「エラトステネスの篩」のような篩系のアルゴリズムで、「アトキンの篩」というものもある。内容的には更なる奇数の分類わけという感じだが、説明は大変そうなので、こちらのサイトへ丸投げしよう(笑)。私も「ふむふむナルホド」と頷きつつ、たぶん本当の意味は理解していない(笑)。「奇数を分類し、いくつかの方程式の解が奇数個なら素数」という話が書いてあるね。とりあえず色々試してみたら確かに正しい答えが出たので、私としては問題ない(笑)。
（参考）
・アトキンのふるい
内容的にはいくつかのソースを試して、一番高速だったものを、具体的な値を返すように改造しただけのものだ。コードもほぼそのまま。
●素数のリストを作る②（アトキンの篩(ふるい)）
import java.util.ArrayList;
import java.util.List;
/**<h1>素数のリストを作る（アトキンの篩）</h1>
* <p>奇数を分類わけし、いくつかの方程式の解での判別と平方数を取り除く篩をかける。</p>
* @param n ： 取得する素数の最大の値 (>=0) [n を含む]
* @return<b>List<Integer></b> ： 素数のリスト
*/
public static final List<Integer> primeListAtkin(final int n) {
final List<Integer> list = new ArrayList<Integer>();
if (n < 2) { //n が2以上とわかっていれば、ここは不要
return list;
}
list.add(2);
if (n == 2) { //n が3以上とわかっていれば、ここは不要
return list;
}
final boolean[] sieve = new boolean[n + 1];
final int sqr = (int)Math.sqrt(n);
int k;
for (int i = 1; i <= 5; i += 4) {
for (int y = i; y <= sqr; y += 6) {
for (int x = 1; x <= sqr && (k = 4 * x * x + y * y) <= n; x++) {
sieve[k] = !sieve[k];
}
for (int x = y + 1; x <= sqr && (k = 3 * x * x - y * y) <= n; x += 2) {
sieve[k] = !sieve[k];
}
}
}
for (int i = 2; i <= 4; i += 2) {
for (int y = i; y <= sqr; y += 6) {
for (int x = 1; x <= sqr && (k = 3 * x * x + y * y) <= n; x += 2) {
sieve[k] = !sieve[k];
}
for (int x = y + 1; x <= sqr && (k = 3 * x * x - y * y) <= n; x += 2) {
sieve[k] = !sieve[k];
}
}
}
for (int i = 1; i <= 2; ++i) {
for (int y = 3; y <= sqr; y += 6) {
for (int x = i; x <= sqr && (k = 4 * x * x + y * y) <= n; x += 3) {
sieve[k] = !sieve[k];
}
}
}
for (int i = 5; i <= sqr; i += 2) {
if (sieve[i]) {
for (int j = i * i; j <= n; j += i * i) {
sieve[j] = false;
}
}
}
sieve[2] = sieve[3] = true;
for (int i = 3; i <= n; i += 2) {
if (sieve[i]) {
list.add(i);
}
}
return list;
}
//メインでは...
int n = 100;
List<Integer> primes = primeListAtkin(n);
for (int v : primes) {
System.out.print(v + " ");
}
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
今回も他のアルゴリズムと比較のため、n が０と１のときにも対応させてある。素数とは「１とそれ自身でしか割り切れない、２以上の整数」と定義されているので、不要なら０と１ははじめから除外しても良い（２以上または３以上を使うこと前提ならば、関数内冒頭の「if (n < 2)～」「if (n == 2)～」は削除しても構わない）。なるべく同じ条件で比較したかったので、「エラトステネスの篩」や「サンダラムの篩」のコードに合わせたら、関数冒頭が少し長くなってしまった。
実際にこれまでに作った「エラトステネスの篩」と「サンダラムの篩」とこの「アトキンの篩」を比較してみると、「アトキンの篩」「サンダラムの篩」「エラトステネスの篩」の順に速い。ただし、篩(ふるい)の真偽値（sieve[]）でだけなら、「アトキンの篩」「エラトステネスの篩」「サンダラムの篩」となった。実装の仕方もあるのだろうが、やはりコンピュータ時代に生み出された「アトキンの篩」のアルゴリズムの方が比較的速いようだ。とは言え「エラトステネスの篩」は改造しやすいため、出題に合わせて最適化しやすいメリットもある。逆に数式が複雑に絡むコードは改造がしづらく、手を加えればバグを埋め込むことにもなり兼ねないので、スピード特化で使うのも良い気がする。結局は複数のアルゴリズムを知っている方がメリットを活かせるのかも知れない。
ググっていて感じたのは「サンダラムの篩」はあまり例がないのに対し、「アトキンの篩」は思いのほか実装例がたくさんあったことだ。なので色々試したが、やはりなるべく平方根や除算を使わず、シンプルな実装の方が速かった。isPrime() のとき書いたポイントは非常に小さなことだけど、実は意外と効果は高い。以下にその実装例を挙げておくが、アルゴリズムはすべて同じなわけだから、その「ちょっとした違い」で実行速度が変わることを感じてみるのも勉強になるかも知れない。
（参考）
・アトキンのふるい
・エラトステネスの篩 vs アトキンの篩
・アトキンの篩
・ダメ出し：アトキンの篩
・Rubyで、「エラトステネスのふるい」を改良したと言われる「アトキンのふるい」をいまさら実装する
・Atkinの篩
（関連記事）
【Java】素数判定②（エラトステネスの篩(ふるい)）
【Java】素数判定③（サンダラムの篩(ふるい)）
【Java】素数判定①（試し割り法）
tb: 0
cm: --
【Java】素数判定③（サンダラムの篩(ふるい)）
2016/06/06 Mon [edit]
「エラトステネスの篩」とは別の方法で篩(ふるい)をかけていく方法で、「サンダラムの篩」というものもある。内容的には数学の証明のような感じで、素数と合成数の種類から分類するような感じになる。Wikipedia の説明ではわかりづらいので、こちらのサイトなどを見ると良いだろう。具体的でわかりやすく、C言語のソースも載っている。内容的にはこのソースをほぼそのままで、具体的な値を返すように改造しただけのものである。
●素数のリストを作る②（サンダラムの篩(ふるい)）
import java.util.ArrayList;
import java.util.List;
/**<h1>素数のリストを作る（サンダラムの篩）</h1>
* <p>整数 n > 1 は、素数または合成数で、合成数だけを篩い落せば素数が残される。</p>
* @param n ： 取得する素数の最大の値 (>=0) [n を含む]
* @return<b>List<Integer></b> ： 素数のリスト
*/
public static final List<Integer> primeListSundaram(final int n) {
final List<Integer> list = new ArrayList<Integer>();
if (n < 2) { //n が2以上とわかっていれば、ここは不要
return list;
}
final int h = ((n + 1) >>> 1); //÷2 の代わり
final boolean[] sieve = new boolean[h];
int k;
for (int i = 1; i * i <= n; i++) {
for (int j = i; (k = i + j + 2 * i * j) < h; j++) {
sieve[k] = true; //k = 4,7,10,12,13,16,17,…
}
}
list.add(2);
for (int i = 1; i < h; i++) {
if (!sieve[i]) {
list.add(2 * i + 1); //素数の数列
}
}
return list;
}
//メインでは...
int n = 100;
List<Integer> primes = primeListSundaram(n);
for (int v : primes) {
System.out.print(v + " ");
}
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
今回も他のアルゴリズムと比較も兼ねて、n が０と１のときにも対応させてある。素数とは「１とそれ自身でしか割り切れない、２以上の整数」と定義されているので、不要なら０と１ははじめから除外しても良い（２以上を使うこと前提ならば、関数内冒頭の「if (n < 2)～」は削除しても構わない）。
サンダラムの篩ではエラトステネスの篩のときのように、篩だけ（boolean[] 配列だけ）の例は掲載しなかったが、この篩は「偶数の合成数と奇数の合成数を篩った数列」であって、素数のフラグそのものではないためだ。エラトステネスの篩と同じようにフラグの配列が欲しければ、関数内でもう１つ配列を作り、最後のリストを作っている部分（list.add(2 * i + 1)）を配列に変えれば良い。ただし、２つ配列を作りオーダーがO(n/2)のためか、篩だけならエラトステネスの篩の方が速い（エラトステネスは配列１つでO(√n)）。またリストを作る場合はO(n/2)で良いためか、サンダラムの篩の方が速い（エラトステネスはO(n)）。実装にもよると思うが、使い分けしても良いかも知れない。
コードは比較のため、前回に合わせた書き方という感じもあるが、他に参考になったサイトなども載せておこう。配列やリスト、言語の特質なんかもあるかも知れないが、色々試してみると「なるほどな～」と思ったりする(笑)。
（参考）
・サンダラムの篩 (MEMORANDUM)
・サンダラムの篩 (Trivial Contents)
・サンダラムの篩 (私的数学塾)
もう１つ、現代風のアルゴリズム(?)の「アトキンの篩」というものもある。Wikipedia では英語版しかないが、他のサイトで用例があったので次回やってみよう。
（関連記事）
【Java】素数判定②（エラトステネスの篩(ふるい)）
【Java】素数判定④（アトキンの篩(ふるい)）
【Java】素数判定①（試し割り法）
tb: 0
cm: --
【Java】素数判定②（エラトステネスの篩(ふるい)）
2016/06/05 Sun [edit]
前回の isPrime() は１つの数値が素数かどうかを調べる方法だったが、今回の「エラトステネスの篩(ふるい)」は素数かどうかの一覧のような感じのもの。Wikipedia にも書いてあるように、コンピュータがある以前からあるアルゴリズムである。「ユークリッドの互除法」に並んで最古のアルゴリズムとも言われているらしい。故にわかりやすく、色々応用しやすい長所がある。その内容はググればいくらでも出てくるので他のサイトに譲るが、Wikipedia の右端にあるアニメーションを見ているだけでもその方法は十分にわかると思う。
今回は素数かどうかのフラグを列挙した配列と実際の素数の値を列挙したリストを作ってみよう。
●素数かどうかを判定する②(エラトステネスの篩(ふるい))
/**<h1>エラトステネスの篩（素数フラグ）を作る</h1>
* <p>偶数以外で、3 から √n までに対して順番に倍数をふるい落としていく。</p>
* @param n ： 調べる素数の最大の値 (>=0) [n を含む]
* @return<b>boolean[]</b> ： 素数である数のインデクスが true となる配列
*/
public static final boolean[] eratosthenesSieve(final int n) {
final boolean[] sieve = new boolean[n + 1]; //素数フラグ
if (n < 2) { //n が2以上とわかっていれば、ここは不要
return sieve;
}
sieve[2] = true;
for (int i = 3; i <= n; i += 2) { //2の倍数以外(奇数)を候補にする
sieve[i] = true;
}
for (int i = 3; i * i <= n; i += 2) { //奇数のみ [O(√n[/2])で可]
if (sieve[i]) {
for (int j = i * 2; j <= n; j += i) { //倍数を候補から除外する
sieve[j] = false;
}
}
}
return sieve;
}
//メインでは...
int n = 10;
boolean[] sieve = eratosthenesSieve(n);
for (int i = 0; i < sieve.length; i++) {
System.out.print(sieve[i] + " "); //素数である i が true になる
}
false false true true false true false true false false false
真偽値は素数となるインデクスの数が true となる（上の例の場合、[2], [3], [5], [7] が true となっているので、2～10(n) までの数の中で 2, 3, 5, 7 が素数ということを表している）。
素数とは「１とそれ自身でしか割り切れない、２以上の整数」と定義されているので、０と１ははじめから除外しても良いが、前回の isPrime() が０と１にも対応していたので、比較のため合わせておいた。２以上を使うこと前提ならば、関数内冒頭の「if (n < 2)～」は削除しても構わない。
コード自体は isPrime() とあまり変わらないので、必要なら前回の解説を見て欲しい。
もう１つ、今度は具体的な素数の値のリストを返す関数を作ろう。作り方は簡単で、上記のエラトステネスの篩をしながら、true となっている値をリストに追加してくだけだ。
●素数のリストを作る①（エラトステネスの篩(ふるい)）
import java.util.ArrayList;
import java.util.List;
/**<h1>素数のリストを作る（エラトステネスの篩）</h1>
* <p>偶数以外で、3 から √n までに対して順番に倍数をふるい落としていく。</p>
* @param n ： 取得する素数の最大の値 (>=0) [n を含む]
* @return<b>List<Integer></b> ： 素数のリスト
*/
public static final List<Integer> primeListEratosthenes(final int n) {
final List<Integer> list = new ArrayList<Integer>();
if (n < 2) { //n が2以上とわかっていれば、ここは不要
return list;
}
final boolean[] sieve = new boolean[n + 1]; //素数フラグ
sieve[2] = true;
list.add(2);
for (int i = 3; i <= n; i += 2) { //2の倍数以外(奇数)を候補にする
sieve[i] = true;
}
//エラトステネスのふるい
for (int i = 3; i <= n; i += 2) { //奇数のみ[リスト作成のため O(n[/2])]
if (sieve[i]) {
list.add(i);
for (int j = i * 2; j <= n; j += i) { //倍数を候補から除外する
sieve[j] = false;
}
}
}
return list;
}
//メインでは...
int n = 100;
List<Integer> primes = primeListEratosthenes(n);
for (int v : primes) {
System.out.print(v + " ");
}
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
今回は簡単にするため、List を用いているが、配列を戻値にしても良いだろう。ただ、配列の場合はあらかじめ要素数を決めないといけないので、篩（sieve[i]）が true になっているものをカウントしておいて（上記の例の場合、リストに追加（list.add()）している場所）、カウント後に配列のインスタンスを生成する必要がある。または取り得る最大の個数（上限値）が決まっているなら、引数に配列の参照を受け取り、値を順に追加した上で、最後に追加したインデクスを戻値として返す方法もある（メインではその戻値で [0]～[idx] の間だけ使うようにする）。具体的な値のリストが必要ないなら、篩の真偽値だけの方がもちろん速い。
ちなみに「n までの素数の個数」は「素数定理」や「ベルトランの仮説」にある通り、「n が十分大きいときには n と 2n の間の素数の個数は n / log n に近いことが言え、特にベルトランの仮説によって保証されている1つの素数の存在よりもより強く、より多くの素数が n と 2n の間に存在していることが分かる」(Wikipedia)と言われているので、「n / log n」で近似してみるのも良いかも知れない。または「ルジャンドルの素数定理」で「n 以下の素数の個数は n/(log n - 1.08366) で近似できる」とか「リーマンゼータ関数」など色々あるようなので、覗いてみると面白いかも知れない。ただ、私もいくつか試してみたが、元々「与えられた数より小さい素数の個数」なので、配列の要素数定義には使えなかった（n + 1 したり、関数の答えを +1 してみたが、素数の個数が必要な配列の要素数を上回る個数になることもある[ArrayIndexOutOfBoundsException]）。いい加減な個数で良いのなら「n * 1.25505871293248 / log n」あたりにしておくと、n までの素数の個数が上回ることはないようだ（この係数は全探索で算出したものなので、Java 以外で適用できるかわからないが…（→実数計算は言語によって微妙に誤差が出るため））。まぁ、無駄は多くなるが n 個で要素数を作っても良いと思う。参考までに素数の個数に関する資料を載せておこう。
（参考）
・素数定理の歴史
・素数個数関数の近似、素数定理、グラフ
・平方数間の素数の個数(1)、平方数間の素数の個数(2)
ネットに落ちているサンプルコードも同じアルゴリズムでありながら、実装によって実行速度は随分違うので、配列やリスト、真偽値のみなど、用途に合わせて使い分けた方が良いだろう。特にエラトステネスの篩は他のアルゴリズムと違って改造しやすい利点があるので、例えば「開始～終了インデクスの間だけ値を取り出す」とか、n を上限数ではなく「n 個の素数を取り出す」など簡単に作れる。オンラインジャッジでの出題でも利用頻度は高いので、いくつかのパターンをライブラリとして作っておくとかなり役に立つだろう。
他にも同じ篩系のアルゴリズムとして「サンダラムの篩」というものもある。せっかくなので、比較のためにも少しだけやってみよう。
（関連記事）
【Java】素数判定③（サンダラムの篩(ふるい)）
【Java】素数判定④（アトキンの篩(ふるい)）
【Java】素数判定①（試し割り法）
■参考になる書籍
tb: 0
cm: --
| h o m e |
