素数のこさを予想する
このワークシートはMath by Codeの一部です。
＜素数の減り方＞
前回、双子素数の分布を調べてみて、
最初は急に減りますが、そのあとは相当なだらかな減り方になっていました。
素数自体はどうでしょう。
幅で切って調べるとどうしても変動がでるので、
平均的な値、こさを調べる。
だから、区間で切るのでなく、累積で素数個数をたして、全体の範囲で割ろう。
まずは、数値化する。
最後に1/ln(x)と比べてグラフ化してみよう。
ご存知のように、底をeにした指数関数e^xの逆関数をln(x)とか、log(x)とかきます。
単調に増加するイメージはわきますね。
最初だけ急にふえてあとならだかに増えるあの曲線です。
その逆数ですから、
1/ln(x)は、最初だけ急に減少してあとなだらかに減少する曲線になりますね。
[IN]julia
#=====================================================
function isPrime(n)
lim = Int(round(n^0.5))
if n<2
return false
end
for num in 2:lim
if n % num ==0
return false
end
end
return true
end
limit=5000
size=limit/50
println(limit,"までの素数のこさ")
clxs=[]
rates=[]
ps=0
for cluster in 0:50
starting=0+cluster*size
ending=size-1+cluster*size
clx=size*(cluster+1)
N= [x for x in starting:ending]
P= filter(n->(isPrime(n)==true),N)
ps +=length(P)
rate = ps/clx
push!(rates,rate)
push!(clxs,clx)
end
println(clxs)
println(rates)
#=====================================================
limit=5000のとき、
[OUT]
5000までの素数のこさ
Any[100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0, 1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0, 2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0, 3600.0, 3700.0, 3800.0, 3900.0, 4000.0, 4100.0, 4200.0, 4300.0, 4400.0, 4500.0, 4600.0, 4700.0, 4800.0, 4900.0, 5000.0, 5100.0]
Any[0.25, 0.23, 0.20666666666666667, 0.195, 0.19, 0.18166666666666667, 0.17857142857142858, 0.17375, 0.1711111111111111, 0.168, 0.16727272727272727, 0.16333333333333333, 0.16230769230769232, 0.15857142857142856, 0.15933333333333333, 0.156875, 0.1564705882352941, 0.15444444444444444, 0.15263157894736842, 0.1515, 0.15095238095238095, 0.14863636363636365, 0.14869565217391303, 0.14875, 0.1468, 0.1453846153846154, 0.14555555555555555, 0.14535714285714285, 0.14448275862068966, 0.14333333333333334, 0.14258064516129032, 0.14125, 0.1403030303030303, 0.14058823529411765, 0.1397142857142857, 0.13972222222222222, 0.13945945945945945, 0.13894736842105262, 0.1382051282051282, 0.1375, 0.1378048780487805, 0.13666666666666666, 0.1372093023255814, 0.13613636363636364, 0.13555555555555557, 0.13521739130434782, 0.13489361702127659, 0.13458333333333333, 0.13346938775510203, 0.1338, 0.13352941176470587]
100までに25個、200までに46個、…..... という結果が、
x=100で、rates=0.25, x=200で、rate=0.23.,,,,,,,,というデータになっていることがわかるね。
もっと、ｘの桁数を上げて調べるとどうなるだろうか？
limit=1000000のとき、
[OUT]
1000000までの素数のこさ
Any[20000.0, 40000.0, 60000.0, 80000.0, 100000.0, 120000.0, 140000.0, 160000.0, 180000.0, 200000.0, 220000.0, 240000.0, 260000.0, 280000.0, 300000.0, 320000.0, 340000.0, 360000.0, 380000.0, 400000.0, 420000.0, 440000.0, 460000.0, 480000.0, 500000.0, 520000.0, 540000.0, 560000.0, 580000.0, 600000.0, 620000.0, 640000.0, 660000.0, 680000.0, 700000.0, 720000.0, 740000.0, 760000.0, 780000.0, 800000.0, 820000.0, 840000.0, 860000.0, 880000.0, 900000.0, 920000.0, 940000.0, 960000.0, 980000.0, 1.0e6,
1.02e6]
Any[0.1131, 0.105075, 0.10095, 0.0979625, 0.09592, 0.094175, 0.09292857142857143, 0.09176875, 0.09078888888888889, 0.08992, 0.08917272727272728, 0.08842083333333334, 0.08783461538461539, 0.08725714285714285, 0.08665666666666667, 0.086275, 0.08582941176470588, 0.08543611111111112, 0.085, 0.08465, 0.08426190476190476, 0.08395681818181819, 0.08360434782608696, 0.08334375, 0.083076, 0.08280961538461538, 0.08254074074074073, 0.08227142857142856, 0.08204827586206896, 0.08183, 0.08163225806451613, 0.081365625, 0.08115757575757576, 0.080975, 0.08077571428571428, 0.08059583333333334, 0.0804472972972973, 0.08023421052631578, 0.08008717948717949, 0.07993875, 0.07977560975609756, 0.07963095238095239, 0.07946744186046512, 0.07934431818181818, 0.07919333333333334, 0.07905869565217391, 0.07892234042553191, 0.07876875, 0.07863979591836735, 0.078498, 0.07840196078431373]
limit=10000000にすると
[OUT]
10000000までの素数のこさ
Any[200000.0, 400000.0, 600000.0, 800000.0, 1.0e6, 1.2e6, 1.4e6, 1.6e6, 1.8e6, 2.0e6, 2.2e6, 2.4e6, 2.6e6, 2.8e6, 3.0e6, 3.2e6, 3.4e6, 3.6e6, 3.8e6, 4.0e6, 4.2e6, 4.4e6, 4.6e6, 4.8e6, 5.0e6, 5.2e6, 5.4e6, 5.6e6, 5.8e6, 6.0e6, 6.2e6, 6.4e6, 6.6e6, 6.8e6, 7.0e6, 7.2e6, 7.4e6, 7.6e6, 7.8e6, 8.0e6, 8.2e6, 8.4e6, 8.6e6, 8.8e6, 9.0e6, 9.2e6, 9.4e6, 9.6e6, 9.8e6, 1.0e7, 1.02e7]
Any[0.08992, 0.08465, 0.08183, 0.07993875, 0.078498, 0.07744833333333333, 0.07651857142857142, 0.075704375, 0.07504, 0.0744665, 0.07393727272727273, 0.07345916666666667, 0.07303076923076923, 0.07262928571428572, 0.072272, 0.0719403125, 0.07162911764705883, 0.07131277777777778, 0.07104921052631578, 0.0707865, 0.07055095238095238, 0.0703034090909091, 0.07009586956521739, 0.069883125, 0.0697026, 0.06950134615384615, 0.0693262962962963, 0.06914321428571428, 0.06896431034482758, 0.06880816666666667, 0.06865290322580646, 0.0685015625, 0.06835742424242425, 0.0682164705882353, 0.06809257142857143, 0.06796097222222222, 0.0678327027027027, 0.06770592105263158, 0.06758384615384615, 0.067472125, 0.06735597560975609, 0.0672472619047619, 0.06714406976744186, 0.06704613636363636, 0.06694322222222222, 0.06683880434782609, 0.06674468085106383, 0.06665114583333333, 0.06655765306122449, 0.0664579, 0.06636960784313725]
質問：素数のこさをグラフにしてみましょう。どんな特徴が見えますか？
視覚化のためにパッケージをインストールします。
jupyter notebook のコードとして次の一行を入力して実行しましょう。
import Pkg; Pkg.add("PyPlot")
しばらくすると、パッケージがインストールされて、PyPlotが使えるようになります。
#julia
#======================
using PyPlot
xs= clxs
# 描画
ys1= rates
ys2= 1 ./(log.(xs))
plt.ylim(0, 0.12)
plt.plot(xs,ys1)
plt.plot(xs,ys2)
xs[end],ys2[end],ys1[end]
#======================
limit=1000000のとき
[OUT]
(1.02e6, 0.07227881195095438, 0.07840196078431373)
limit=10000000のとき
[OUT]
(1.02e7, 0.06196593774218297, 0.06636960784313725)
整数ｘまでの範囲での素数のこさは約1/ln(x)。
これはグラフから予想がつくね。
整数ｘまでの素数の個数をπ(x)とすると、素数のこさはπ(x)/xですね。
だから
π(x)/x≒1/ln(x)
という式ができそうです。これが素数定理です。
これをかきかえると、
ln(x)・π(x)/x≒1
π(x)/(x/ln(x))≒1
つまり、
素数の個数はx/ln(x)
に近づくだろうということだね。
質問：rustを使って、もっと高速に素数の濃さが１/ln(x)に近づくことを視覚化するにはどうしますか。
juliaでは、ｎの平方根までの整数の中に約数がなければ素数だという判定を毎回使ってましたが、事前に素数表を作ってしまいましょう。
エラトステネスのふるいを使います。
素数表の範囲をlimitまでとします。すべてtrueにしておき、0，1をfalseとします。2はtrueだから素数という意味です。２が素数とわかったら、そのあとの2の倍数をすべてfalseにすべて変えます。
kが素数trueだとわかったらkの2乗未満のすでにkより小さい素数の倍数としてfalseになっているはずなので、kの2乗以上のkの倍数だけをfalseにします。こうして、kをlimitまで動かし、trueのまま残った数が素数です。素数表の名前をprimesとします。
素数表のlimitまでの範囲を50等分ではなく100等分して、
一区切りの長さstep_size=limit/100とします。
区切りに0番から99番まで番号をつけ、その間のcnt番目までの平均濃度の求め方は、
基本的にjuliaとほぼ同じです。
区切りの開始start=cnt * step_sizeから
終了end=(cnt+1) * step_sizeの手前までの素数個数がsegment_primes_countです。
let segment_primes_count = (start..end)
.filter(|&x| primes[x])
.count();
count_sumにsegment_primes_countを累積してたしていき、
endまでの範囲での濃度を求めます。
それを視覚化するためのdataに追加していきます。
count_sum += segment_primes_count;
let x_pos = end as f64;
let density = count_sum as f64 / end as f64;
data.push((end as f64, density));
dataの横軸がendで、たて軸がdensityですね。
実際の素数濃度のグラフを青、１/ln(x)のグラフを緑にして
共存させることで理論値との差が減ることを視覚化します。
[IN]rust in Cargo.toml
[dependencies]
plotters = "0.3"
[IN]rust in main.rs
use plotters::prelude::*;
fn main() {
let limit = 100_000_000;
let step_size = limit / 100;
let primes = create_prime_table(limit);
let mut data: Vec<(f64, f64)> = Vec::new();
let mut count_sum = 0;
println!("{}までの素数濃度の計算中", limit);
for cnt in 0..100 {
let start = cnt * step_size;
let end = (cnt + 1) * step_size;
let segment_primes_count = (start..end)
.filter(|&x| primes[x])
.count();
count_sum += segment_primes_count;
let density = count_sum as f64 / end as f64;
data.push((end as f64, density));
}
println!("{}までのグラフ作成中", limit);
draw_chart(data, limit as i64);
}
//あの有名なエラトステネスのふるいをつかます。
fn create_prime_table(limit: usize) -> Vec {
let mut table = vec![true; limit + 1];
table[0] = false;table[1] = false;// 0,1は素数ではない。
let sqrt_limit = (limit as f64).sqrt() as usize;
for p in 2..=sqrt_limit {
if table[p] {
// p*pから開始してpずつ飛ばして効率化をはかる。
for i in (p * p..=limit).step_by(p) {
table[i] = false;
}
}
}
table
}
// 座標データdataから折れ線グラフ(.png)を直下に作ります。
fn draw_chart(data: Vec<(f64, f64)>, limit_val: i64) { // 引数名をlimit_valに変更
if data.is_empty() { return; }
let xs_data: Vec = data.iter().map(|(x, _)| *x).collect();
let filename = format!("primes_density_{}.png", limit_val); // ファイル名変更
let title = format!("Primes Density to {} with 1/ln(x)", limit_val); // タイトル変更
let root = BitMapBackend::new(&filename, (1000, 700)).into_drawing_area();
root.fill(&WHITE).unwrap();
let mut chart = ChartBuilder::on(&root)
.caption(title, ("sans-serif", 50).into_font())
.margin(15)
.x_label_area_size(40)
.y_label_area_size(40)
.build_cartesian_2d(
*xs_data.first().unwrap()..*xs_data.last().unwrap(),
0.0..0.3 // 縦軸の最大値を0.3で固定してグラフ比較しやすくする。
)
.unwrap();
chart.configure_mesh()
.x_desc("Range")
.y_desc("Density")
.draw()
.unwrap();
// 実測値 (青色)
chart.draw_series(LineSeries::new(
data.iter().map(|(x, y)| (*x, *y)),
&BLUE,
)).unwrap()
.label("Actual Density") // 凡例の追加
.legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &BLUE));
// 理論値 (緑色)
chart.draw_series(LineSeries::new(
(2..=limit_val as usize)
.filter_map(|x| {
let x_f64 = x as f64;
if x_f64 > 1.0 { Some((x_f64, 1.0 / x_f64.ln())) } else { None }
}),
&GREEN,
)).unwrap()
.label(" 1 / ln(x)") // 凡例を追加
.legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &GREEN));
// 凡例を描画
chart.configure_series_labels()
.background_style(&WHITE.mix(0.8))
.border_style(&BLACK)
.draw()
.unwrap();
root.present().expect("ファイルの書き込みに失敗しました。");
println!("グラフ: {}ができました。", filename);
}
２つの素数の型
素数には4n+1,4n+3の２つのタイプがあります。
おもしろいことに、はじめのうちは4n+3型の方がかなり多く現れます。
しかし、２つのタイプを割合で素数に占める割合でくらべると、徐々に0.5に近づいていくようです。
質問：この3型優位の偏り（チェビシェフの偏り）を視覚化するにはどうしたらよいでしょうか。
やはり、先に素数表を作り、
素数が見つかるたびに、4n+3型が素数に占める割合data_rate_3_1にデータを追加していきます。
あとで比較するために、半分の割合のデータhalfも追加します。
この２つのデータdata_rate_3_1、halfを赤、青のグラフにして描画するなどで視覚化できるでしょう。
最初は、3型がとても多いのがわかりますね。
[IN]
//rust in cargo.toml
[dependencies]
plotters = "0.3"
//rust in main.rs
use plotters::prelude::*;
fn main() {
let limit = 10_000_000;
let primes = create_prime_table(limit);
let mut half=Vec::new();
let mut data_rate_3_1 = Vec::new();
let mut count_4n_3 = 0;
let mut sum_primes =0;
for x in 2..=limit {
if primes[x] {
if x % 4 == 3 {
count_4n_3 += 1;
}
sum_primes +=1;
let diff_rate:f64= count_4n_3 as f64/sum_primes as f64;
half.push((x as f64, 0.5 as f64));
data_rate_3_1.push((x as f64, diff_rate as f64));
}
}
draw_comparison_chart(half, data_rate_3_1, limit as i64);
}
fn create_prime_table(limit: usize) -> Vec {
let mut table = vec![true; limit + 1];
table[0] = false; table[1] = false;
let sqrt_limit = (limit as f64).sqrt() as usize;
for p in 2..=sqrt_limit {
if table[p] {
for i in (p * p..=limit).step_by(p) {
table[i] = false;
}
}
}
table
}
fn draw_comparison_chart(d1: Vec<(f64, f64)>, d3: Vec<(f64, f64)>, limit: i64) {
let filename = format!("prime_types_mod4_{}.png", limit);
let title = format!("prime_types_mod4_{}", limit);
let root = BitMapBackend::new(&filename, (1000, 700)).into_drawing_area();
root.fill(&WHITE).unwrap();
// 縦軸の最大値
let max_y = 1.0;
let mut chart = ChartBuilder::on(&root)
.caption(title, ("sans-serif", 40).into_font())
.margin(20)
.x_label_area_size(40)
.y_label_area_size(50)
.build_cartesian_2d(0.0..limit as f64, 0.0..max_y)
.unwrap();
chart.configure_mesh().draw().unwrap();
// 半分 (青)
chart.draw_series(LineSeries::new(d1, &BLUE)).unwrap()
.label("Half")
.legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &BLUE));
// 4n+3 型の割合 (赤)
chart.draw_series(LineSeries::new(d3, &RED)).unwrap()
.label("3(mod4)type_rate")
.legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &RED));
chart.configure_series_labels()
.background_style(&WHITE.mix(0.8))
.border_style(&BLACK)
.draw()
.unwrap();
root.present().expect("Unable to save file");
println!("グラフ作成完了: {}", filename);
}
