2013-10-09

素数の数の近似

1からxまでに存在する素数の数に関する話
- この件に関してウェブサーチをするけれど、リーマン予想の一般的な話とかばかりで、今日のテレビ録画の「ゼータ関数の非自明な解を使った素数の数の推定精度の向上」のための情報がなかなかわかりやすく書いてあるところに行きつかない。こちらが『経緯』などからしてわかりやすい。
numbersパッケージを使って、地道に数えてみよう

library(numbers)
xs <- 1:1000
y1 <- rep(0,length(xs))
for(i in 1:length(xs)){
	y1[i] <- length(Primes(xs[i]))
}
y1

いくつかの近似公式が知られているので、それをx=10^(1:12)についてみてみる
- xが大きいと、地道な計算は大変なので、既に計算されている数を使おう

xs <- 10^(1:12)
piN <- c(4, 25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534, 455052511, 4118054813, 37607912018)

まずは、ガウスの近似、 $\frac{x}{\ln{x}}$

y.g <- xs/log(xs)
cbind(xs, piN, y.g=round(y.g))

> cbind(xs, piN, y.g=round(y.g))
         xs         piN         y.g
 [1,] 1e+01           4           4
 [2,] 1e+02          25          22
 [3,] 1e+03         168         145
 [4,] 1e+04        1229        1086
 [5,] 1e+05        9592        8686
 [6,] 1e+06       78498       72382
 [7,] 1e+07      664579      620421
 [8,] 1e+08     5761455     5428681
 [9,] 1e+09    50847534    48254942
[10,] 1e+10   455052511   434294482
[11,] 1e+11  4118054813  3948131654
[12,] 1e+12 37607912018 36191206825

次に、いわゆる関数。
- gslパッケージにある対数積分(logarithmic integral)関数expint_Eiを用いて

Li <- function(x) expint_Ei(log(x)) - expint_Ei(log(2))
y.L <- Li(xs)
cbind(xs, piN, y.g=round(y.g),y.L=round(y.L))

> cbind(xs, piN, y.g=round(y.g),y.L=round(y.L))
         xs         piN         y.g         y.L
 [1,] 1e+01           4           4           5
 [2,] 1e+02          25          22          29
 [3,] 1e+03         168         145         177
 [4,] 1e+04        1229        1086        1245
 [5,] 1e+05        9592        8686        9629
 [6,] 1e+06       78498       72382       78627
 [7,] 1e+07      664579      620421      664917
 [8,] 1e+08     5761455     5428681     5762208
 [9,] 1e+09    50847534    48254942    50849234
[10,] 1e+10   455052511   434294482   455055614
[11,] 1e+11  4118054813  3948131654  4118066400
[12,] 1e+12 37607912018 36191206825 37607950280

次に、リーマン関数が、なる関係を持つが、この第１項のもよい近似式なので、これを使ってみる。第二項については後述(だけれど、これが本題)
- $R(x)=1+\sum_{n=1}^{\infty}(\frac{1}{n\zeta(n+1)}\frac{(\ln{x})^n}{n!})$ (グラム級数)
- ここの $\zeta(x)$ はxが１以上なのでVGAMパッケージのzeta関数を使えるから(もちろんこちらの方が前記事の自前関数より速いから)、そちらを使って

my.R.fx <- function(x,n=10000){
	#term1 <- my.zeta.2((1:n)+1)
	term1 <- zeta((1:n)+1)
	term2 <- (1:n)*log(log(x)) - cumsum(log(1:n))
	sum(1/((1:n)*term1)*exp(term2))+1
}
y.R <- rep(0,length(xs))
for(i in 1:length(xs)){
	y.R[i] <- my.R.fx(xs[i])
}
cbind(xs, piN, y.g=round(y.g),y.L=round(y.L),y.R=round(y.R))

> cbind(xs, piN, y.g=round(y.g),y.L=round(y.L),y.R=round(y.R))
         xs         piN         y.g         y.L         y.R
 [1,] 1e+01           4           4           5           5
 [2,] 1e+02          25          22          29          26
 [3,] 1e+03         168         145         177         168
 [4,] 1e+04        1229        1086        1245        1227
 [5,] 1e+05        9592        8686        9629        9587
 [6,] 1e+06       78498       72382       78627       78527
 [7,] 1e+07      664579      620421      664917      664667
 [8,] 1e+08     5761455     5428681     5762208     5761552
 [9,] 1e+09    50847534    48254942    50849234    50847455
[10,] 1e+10   455052511   434294482   455055614   455050683
[11,] 1e+11  4118054813  3948131654  4118066400  4118052495
[12,] 1e+12 37607912018 36191206825 37607950280 37607910542

かなり良い近似になっていることが確かめられた
さて、ここからが本題
ゼータ関数で実部が0.5のところにの解があるんですよ、というのがリーマン予想で、このような解を負の偶整数を自明解と呼ぶのに対して、非自明解とか複素零点とかいうのだけれど、その複素零点全体に渡って、足し合わせるのが、先のリーマン関数に出てきた第二項
- この近似式の残差項を複素零点を使ってさらに補正してみる
  - 残差項を入れない方がよっぽど成績が良いようだが…どこかがおかしいのか、「もっと大きい値の世界」では違うのか…

zeta.zero <- 0.5+c(14.134725,21.022040,25.010858,30.424876,32.935062,37.586178)*1i
xs <- seq(from=100,to=1000,by=50)
y1 <- rep(0,length(xs))
for(i in 1:length(xs)){
	y1[i] <- length(Primes(xs[i]))
}
y.R <- rep(0,length(xs))
for(i in 1:length(xs)){
	y.R[i] <- my.R.fx(xs[i],n=1000)
}
y.R.0 <- matrix(0,length(xs),length(zeta.zero)+1)
y.R.0[,1] <- y.R
for(i in 1:length(zeta.zero)){
	tmp <- rep(0,length(xs))
	for(j in 1:length(xs)){
		print(my.R.fx(xs[j]^zeta.zero[i]))
		tmp[j] <- -my.R.fx(xs[j]^zeta.zero[i],n=1000)
	}
	y.R.0[,i+1] <- y.R.0[,i] + tmp
}
cbind(xs, y1, y.R.0=round(Re(y.R.0)))
cbind(xs, y1, y.R.0=round(Re(y.R.0))-y1)

零点数を増やしても…(30くらいでＯＫなはずなのだが)

# ゼータ関数解の虚部のリストを取ってくる
non.Trivial.Zero.10000 <- read.table("http://www.dtc.umn.edu/~odlyzko/zeta_tables/zeros1")
zeta.zero <- unlist(non.Trivial.Zero.10000)
zeta.zero <- 0.5+zeta.zero*1i
library(numbers)
library(pracma)

my.R.fx <- function(x,n=10000){
	#term1 <- my.zeta.2((1:n)+1)
	term1 <- zeta((1:n)+1)
	term2 <-  (1:n)*log(log(x))- cumsum(log(1:n))
	sum(1/((1:n)*term1)*exp(term2))+1
}
#xs <- seq(from=100,to=1000,by=50)
xs <- 1:100
y1 <- rep(0,length(xs))
for(i in 1:length(xs)){
	y1[i] <- length(Primes(xs[i]))
}
y.R <- rep(0,length(xs))
for(i in 1:length(xs)){
	y.R[i] <- my.R.fx(xs[i],n=10)
}
n.k <- 30

y.R.0 <- matrix(0,length(xs),n.k+1)
y.R.0[,1] <- y.R
for(i in 1:n.k){
	tmp <- rep(0,length(xs))
	for(j in 1:length(xs)){
		#print(my.R.fx(xs[j]^zeta.zero[i]))
		tmp[j] <- my.R.fx(xs[j]^zeta.zero[i],n=10)
	}
	y.R.0[,i+1] <- y.R.0[,i] - tmp
}
#cbind(xs, y1, y.R.0=round(Re(y.R.0)))
#cbind(xs, y1, y.R.0=round(Re(y.R.0))-y1)

matplot(y.R.0[,1:n.k],type="l")
matplot(cbind(y1,y.R,y.R.0[,n.k]),type="l")
matplot(cbind(y1,y.R),type="l")
 matplot(cbind(y1,y.R,Re(y.R.0[,n.k+1])),type="l")

Mathematicaサイト