r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

スコットランド独立賛成票が素数~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でも素数だった
前後2000くらいの素数率は7%だった。
まあ、そこそこ珍しいことが起こった模様。奇跡ってほどじゃない。

ちなみに、反対派は素数ではなく 2 * 571 * 1753 = 2001926で3つの
素数に分解された。
独立反対派にとっては、グレートブリテン島イングランドウェールズスコットランド
3つが合体して1つなのだと、こじつけて終了する。(アイルランドはこの際別の島と
いうことで無視w


#中澤さんから別解のご指摘あり追記

> 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の方は素因数分解まで行っているのに、十分速い。
マニュアルにも1000万以下の数なら、factorize()で十分速い
と書いてあった。餅は餅屋だなー。

しかも、
http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
http://leto.net/x/2009/01/primes-primality-and-psuedopri.html
isprimeが1(多分素数)であっても、ほとんど正しいらしい。
軽くwhile文で試したら142348791までは少なくとも正しかった。(これ以上回してもエラーは探せないと思われ。)
それもそのはず、

library(Rmpfr)
four <- mpfr(1,200000)
four^(40)
> print.mpfr(x4, n=30, m=10, deci=0, fill = "")
12089 2581961462 9174706176. 分の1で間違うってなんだそれ。
ほぼ、普通に使う限り正しいのね。

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が表示されない)