素数定理

私は素数定理に興味があり、実際に成り立つのか確かめてみたいと思いこの題材を選んだ。

素数定理とは\(π(x)\)を\(x\)以下の素数の数として

\[π(x)～x/log(x)\] と近似できることである。 \[π(x)～ \int_2^x dt/log(t) \] とするとより正確に近似出来ることも知られている。

n以下の素数の個数

近似を確認するためにまず必要な関数を作成する。

n以下の素数を数える関数を作るため、数xを与えたときにその数が素数か判定する関数is.primeを次のように作成した。

#xが素数か判定
is.prime  <- function(x){
  
  # xが整数ではない場合
  if (x != as.integer(x)) return("xが整数でなければならない")
  
  # x<2の場合
  if (x < 2){
    return("x≧2でなければならない")
  } 
  
  # xが2または3
  else if (x == 2 || x == 3) {
    return(TRUE)
    #それ以外
  } else {
    for (i in 2:(x-1) ){
      if ((x %% i) == 0) {
        return(FALSE)
      }
    }
    return(TRUE)
  }
}

次に上の関数is.primeを使いn以下の素数を全て表示る関数get.primeを作成した。

#n以下の素数を表示
get.prime <- function(n){
  
  #n<2の場合  
  if (n< 2) return("n以下の素数は無い")
  
  else {
    vec <- NULL
    for (i in 2:n){
      if(is.prime(i) == TRUE) vec <- append(vec,i)
    }
    
    return(vec)
  }
}

例えば100以下の素数を表示すると

get.prime(100)
##  [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83
## [24] 89 97

となる。そして表示された素数を数える関数count.primeを作成した。

#n以下の素数の個数
count.prime <- function(n){
  
  #またはn<2の場合  
  if (n < 2) return(0)
  
  else {
    vec <- NULL
    for (i in 2:n){
      if(is.prime(i) == TRUE) vec <- append(vec,i)
    }
    
    return(length(vec))
  }  
}

使用してみると

count.prime(100)
## [1] 25
count.prime(1000)
## [1] 168
count.prime(10000)
## [1] 1229

と正しく素数の数が表示されていることが分かる。

近似

n以下の素数を数える関数を作ることが出来たので、ここでは実際に近似を確かめてみようと思う。 まず、ここでは1~1000までの素数を数えることにする。

n <- 1000
x <- 1:n
y <- numeric(n)
for(i in x){
  y[i] <- count.prime(i)
}

次に比較する2つの関数をf,gと定義する

f <- function(x) x/log(x)
library(gsl)
g <- function(x) {expint_Ei(log(x)) - expint_Ei(log(2))}

準備が出来たので、図にして近似を確かめる。 \(π(x)\),\(x/log(x)\),\(\int_2^x dt/log(t)\)をそれぞれ赤、緑、青で表すとする。

plot(x,y,xlim=c(0,1000),col=2 ,ylim=c(0,200),xlab="n", ylab="") 
par(new=T)
plot(f(x),col=3, xlim=c(0,1000) ,ylim=c(0,200),ann=F)
par(new=T)
plot(g(x),col=4, xlim=c(0,1000) ,ylim=c(0,200),ann=F)

図を見ると近似されていることが確認することが出来た。

\(π(x)/f(x)\) と　\(π(x)/g(x)\)についても確かめてみると

plot(x,y/f(x),col=2,xlim=c(0,1000),ylim=c(0,2),xlab="n", ylab="")
par(new=T)
plot(x,y/g(x),col=3,xlim=c(0,1000),ylim=c(0,2),ann=F)

やはり1に収束していくことが確認できた。

まとめ

自ら素数を数える関数を作成し、π(x)が\(x/log(x)\)と\(\int_2^xdt/log(t)\)で近似出来ることを確認することが出来た。

図から\(\int_2^xdt/log(t)\)の方がより正確に近似出来ている事もわかった。

手作業では相当時間がかかることをRにより確かめることが出来て、素数定理が実際に成り立つことが確認できてよかった。

参考文献

http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html 関数の作成、plot関係についてはこのサイトを参考にした。

https://ja.wikipedia.org/wiki/%E7%B4%A0%E6%95%B0%E5%AE%9A%E7%90%86 素数定理、n以下の素数の個数などはここを参考にした。

http://seekr.jp/ 対数積分についてはこのサイトで検索しgslパッケージを使用すればいい事がわかった。