スコットランド独立賛成票が素数～Rで素数をチェックする関数色々
#2014.09.23中澤さんの別解など追加しました
http://business.nikkeibp.co.jp/article/topics/20140919/271500/
目的
＃本当に素数か確かめる。
＃その前後にどのくらい素数があるか。
＃ついでに、素数を判定するパッケージを示す。（デスマコロシアムでさんざん検索した遺産www
yes <- 1617989 #賛成票 no <- 2001926 #反対票 yes / (yes + no) #賛成の得票率 binom.test(yes, (yes + no)) #意味ないけど、50%から有意にずれているか diff(qbinom(c(0.05, 0.95), (yes + no), yes / (yes + no))) #yes率から計算した95％CI。 #幅3000くらいなので、余裕を持って±2000くらい計算することにする library(gmp) isprime(yes) #1なので多分素数。確定できず library(matlab) isprime(yes) #素数確定 is.sosu <- function(x){ #素数判定の自作関数作成。（2は判定できない） sum(x %% 2:sqrt(x) == 0) == 0 } #Rの様々な素数判定関数たち system.time(replicate(1000, is.sosu(yes))) system.time(replicate(10, matlab::isprime(yes))) #library(matlab)のis.prime桁違いに遅い system.time(replicate(1000, gmp::isprime(yes))) #高速だけど確実ではない system.time(replicate(1000, DescTools::IsPrime(yes))) system.time(replicate(1000, pracma::isprime(yes))) system.time(replicate(1000, schoolmath::is.prim(yes))) #自作が一番早いのでこれを使う。2は判定できないが今回の目的には関係ない is.sosu(yes) #素数！ n <- 1616000:1620000 res <- sapply(n, is.sosu) res2 <- n[res == TRUE] x5 <- seq(5, length(res2), by = 5) res2[x5] <- paste(res2[x5], "\n", sep ="") cat("", res2, sep = " ") #1617989プラマイ2000くらいの間の素数リスト length(res2) / length(n) #プラマイ2000くらいの間の素数率は7％程度 res2[which(res2 == yes) - 1] #1個前の素数 res2[which(res2 == yes) - 1] - yes #12人少なくても素数だった res2[which(res2 == yes) + 1] #1個あとの素数 res2[which(res2 == yes) + 1] - yes #14人多くても素数だった > yes <- 1617989 #賛成票 > no <- 2001926 #反対票 > yes / (yes + no) #賛成の得票率 [1] 0.4469688 > > binom.test(yes, (yes + no)) #意味ないけど、50%から有意にずれているか Exact binomial test data: yes and (yes + no) number of successes = 1617989, number of trials = 3619915, p-value < 2.2e-16 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.4464565 0.4474811 sample estimates: probability of success 0.4469688 > > diff(qbinom(c(0.05, 0.95), (yes + no), yes / (yes + no))) #yes率から計算した95％CI。 [1] 3112 > #幅3000くらいなので、余裕を持って±2000くらい計算することにする > > library(gmp) > isprime(yes) #1なので多分素数。確定できず [1] 1 > > library(matlab) > isprime(yes) #素数確定 [1] 1 > > is.sosu <- function(x){ #素数判定の自作関数作成。 + sum(x %% 2:sqrt(x) == 0) == 0 + } > #Rの様々な素数判定関数たち > system.time(replicate(1000, is.sosu(yes))) ユーザ システム 経過 0.17 0.00 0.19 > system.time(replicate(10, matlab::isprime(yes))) #library(matlab)のis.prime桁違いに遅い ユーザ システム 経過 1.54 0.00 1.61 > system.time(replicate(1000, gmp::isprime(yes))) #高速だけど確実ではない ユーザ システム 経過 0.04 0.00 0.05 > system.time(replicate(1000, DescTools::IsPrime(yes))) ユーザ システム 経過 0.34 0.00 0.34 > system.time(replicate(1000, pracma::isprime(yes))) ユーザ システム 経過 0.67 0.00 0.67 > system.time(replicate(1000, schoolmath::is.prim(yes))) ユーザ システム 経過 1.17 0.00 1.21 > > #自作が一番早いのでこれを使う > > is.sosu(yes) #素数！ [1] TRUE > > n <- 1616000:1620000 > > res <- sapply(n, is.sosu) > res2 <- n[res == TRUE] > x5 <- seq(5, length(res2), by = 5) > res2[x5] <- paste(res2[x5], "\n", sep ="") > cat("", res2, sep = " ") #1617989プラマイ2000くらいの間の素数リスト 1616009 1616029 1616033 1616039 1616047 1616057 1616063 1616077 1616099 1616113 1616119 1616161 1616171 1616183 1616201 1616221 1616227 1616231 1616269 1616281 1616291 1616297 1616347 1616359 1616401 1616429 1616437 1616443 1616453 1616473 1616491 1616497 1616519 1616533 1616543 1616551 1616569 1616597 1616603 1616609 1616611 1616617 1616621 1616623 1616627 1616633 1616639 1616651 1616669 1616677 1616687 1616689 1616711 1616723 1616749 1616801 1616803 1616807 1616809 1616821 1616827 1616833 1616851 1616861 1616891 1616897 1616899 1616939 1616947 1616963 1616983 1617019 1617029 1617037 1617043 1617047 1617079 1617103 1617137 1617139 1617149 1617211 1617247 1617251 1617269 1617277 1617283 1617289 1617311 1617347 1617349 1617373 1617391 1617433 1617437 1617439 1617443 1617463 1617493 1617503 1617509 1617523 1617541 1617547 1617557 1617563 1617569 1617589 1617619 1617647 1617661 1617689 1617691 1617697 1617727 1617739 1617743 1617757 1617767 1617769 1617773 1617779 1617809 1617817 1617827 1617871 1617883 1617893 1617923 1617929 1617943 1617949 1617971 1617977 1617989 1618003 1618007 1618033 1618039 1618049 1618051 1618079 1618081 1618087 1618091 1618093 1618129 1618139 1618153 1618181 1618187 1618189 1618207 1618217 1618223 1618241 1618261 1618271 1618277 1618291 1618307 1618319 1618327 1618333 1618367 1618369 1618373 1618387 1618399 1618411 1618433 1618453 1618457 1618459 1618471 1618481 1618489 1618501 1618517 1618531 1618537 1618549 1618559 1618601 1618613 1618619 1618627 1618637 1618663 1618679 1618681 1618703 1618739 1618741 1618769 1618777 1618807 1618817 1618823 1618829 1618831 1618849 1618853 1618891 1618909 1618931 1618937 1618943 1618957 1618963 1618973 1618979 1619021 1619053 1619069 1619071 1619087 1619113 1619119 1619153 1619159 1619171 1619179 1619207 1619209 1619227 1619239 1619243 1619249 1619257 1619281 1619287 1619311 1619327 1619329 1619339 1619341 1619353 1619381 1619383 1619417 1619419 1619473 1619507 1619531 1619549 1619551 1619557 1619561 1619593 1619599 1619603 1619633 1619647 1619663 1619669 1619671 1619677 1619687 1619689 1619699 1619713 1619741 1619747 1619753 1619759 1619773 1619791 1619831 1619837 1619857 1619861 1619887 1619899 1619903 1619909 1619929 1619941 1619957 1619983 1619987 > length(res2) / length(n) #プラマイ2000くらいの間の素数率は7％程度 [1] 0.07023244 > res2[which(res2 == yes) - 1] #1個前の素数 [1] 1617977 > res2[which(res2 == yes) - 1] - yes #12人少なくても素数だった [1] -12 > res2[which(res2 == yes) + 1] #1個あとの素数 [1] 1618003 > res2[which(res2 == yes) + 1] - yes #14人多くても素数だった [1] 14
結論。
たしかに1617989は素数だった。
12人少なくて1617977でも素数だった
14人多くて1618003でも素数だった
前後２０００くらいの素数率は７％だった。
まあ、そこそこ珍しいことが起こった模様。奇跡ってほどじゃない。
ちなみに、反対派は素数ではなく 2 * 571 * 1753 = 2001926で３つの
素数に分解された。
独立反対派にとっては、グレートブリテン島はイングランド、ウェールズ、スコットランドの
３つが合体して１つなのだと、こじつけて終了する。（アイルランドはこの際別の島と
いうことで無視ｗ
#中澤さんから別解のご指摘あり追記
> library(gmp); > n <- 1616000:1620000 > system.time(replicate(100,n[sapply(sapply(n, factorize), length)==1])) ユーザ システム 経過 14.59 0.00 14.68 > > n <- 1616000:1620000 > system.time(replicate(100,n[sapply(n, is.sosu)])) ユーザ システム 経過 8.91 0.02 8.94 > > n <- 1616000:1620000 > system.time(replicate(100,n[sapply(n, isprime) != 0])) ユーザ システム 経過 2.56 0.00 2.59
gmpのisprime爆速。しかも多分素数と言いつつ全部正しかった。
gmpのfactorizeの方は素因数分解まで行っているのに、十分速い。
マニュアルにも１０００万以下の数なら、factorize()で十分速い
と書いてあった。餅は餅屋だなー。
しかも、
http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
http://leto.net/x/2009/01/primes-primality-and-psuedopri.html
isprimeが１（多分素数）であっても、ほとんど正しいらしい。
軽くwhile文で試したら142348791までは少なくとも正しかった。（これ以上回してもエラーは探せないと思われ。）
それもそのはず、
library(Rmpfr)
four <- mpfr(1,200000)
four^(40)
> print.mpfr(x4, n=30, m=10, deci=0, fill = "")
12089 2581961462 9174706176. 分の１で間違うってなんだそれ。
ほぼ、普通に使う限り正しいのね。
http://r-statistics-fan.hatenablog.com/entry/2014/09/20/160224多倍長出力をキレイに桁を揃えて出力する（改良版 - r-statistics-fanの日記
なお、print.mpfrは、数が小さいときにe-25みないな
指数表示になる場合に対応させた。
たとえば、4^(-40)だと
__________ __________ __________ __________ _________8.
2718061255 3027674871 4086920699 6285356581 2110900878
9062500000 0000000000 0000000000 0000000000 0000000000
0000000000 0000000000 0000000000 0000000000 0000000000
0000000000 0000000000 0000000000 0000000000 0000000000
e -25
みたいに出力されるように改良した。
（これまではe-25が表示されない）
