---
title: "素数定理シミュレーション"
author: "宮部賢志"
date: "2015年9月29日"
output: html_document
---

#素数定理
**素数定理**とは，自然数の中に素数がどのくらいの「割合」で含まれるかを
述べる定理である．
**素数計数関数**(Prime-counting function)\(\pi(x)\)を，
正の実数\(x\)以下の素数の関数を表す関数とする．
素数定理は，
$$\pi(x)\sim\dfrac{x}{\ln x}$$
であることを主張する．
より精密な近似として，対数積分
$\displaystyle Li(x)=\int_2^x\dfrac{dt}{\ln t}$
を使って，
$$\pi(x)\sim Li(x)$$
が成立することが知られている．

#素数定理のシミュレーション
以下では，この素数定理のシミュレーションをしよう．
ここではライブラリを使うが，提出するレポートではできる限りライブラリを使わないで書くこと．
そうでなければプログラミングの学習にはならない．

まず，整数関連の関数のライブラリを読み込む．
```{r}
library(numbers)
```

次に，いくつまでの整数を数えるのかの変数を設定する．
```{r}
nmax <- 10^3
```

次に素数計数関数(prime-counting function)を次のように計算する．
```{r}
pcf <- rep(NA, nmax)
for(i in 1:nmax){
  pcf[i] <- length(Primes(i))
}
```

最初の10項は以下のようになる．
```{r}
pcf[1:10]
```

比較する関数を設定する．
```{r}
# 関数x/log(x)
f <- function(x) { return (x/log(x)) }

# Liはオイラーの対数積分
library(gsl)
Li <- function(x) {expint_Ei(log(x)) - expint_Ei(log(2))}
```

これらの関数を図にすると以下のようになる．
```{r}
ymax <- Li(nmax)
plot(pcf, xlim = c(1,nmax), ylim=c(1,ymax), xlab="n", ylab="prime counting function")
par(new=T)
plot(f, xlim = c(1,nmax), ylim=c(1,ymax), col="blue", xlab="", ylab="")
par(new=T)
plot(Li, xlim = c(1,nmax), ylim=c(1,ymax), col="red", xlab="", ylab="")
```

$Li$の方がより良い近似を与えていることが分かる．

上記のプログラムはとても遅い．
同じ計算を繰り返し行っているからである．
高速化の方法を各自考えてみよう．
