素数の不思議：チェビシェフの偏り

Wikipediaにもあるのだがあまり詳細な数値が書かれていないのでここでまとめておく。
素数Pが4k+1型（４で割って１余る数）ならば２つの平方数の和で表せる。これは簡単に説明すると、P = a^2 + b^2 のようにかける、そのとき４で割ると１余る、ということだ。Pは２以外は奇数なので奇数足す偶数の形になる（偶数＋偶数＝偶数、奇数＋奇数＝偶数なので奇数＋偶数になる）。よって、a = 2n+1（奇数）, b = 2m（偶数）と書いてもよい。a^2 + b^2 = 4n^2+4n+1 + 4m^2 = 4(n^2+m^2+n)+1 なので４で割ると１余ることになる。
さてチェビシェフの偏りだが、4k+3タイプの素数の数は4k+1タイプの素数の数以上になっている区間が非常に長い（素数Pが数万程度の小さい数の時）という現象のことだ。等号が成立するのは5, 17, 41, 461 で次に等号が成立するのは26833 (1471個ずつ）、4k+1タイプのほうが多くなるのは、26861 (1473個＞1472個）でその次ぐらいから等号が成立した後、ずっと不等号がさらに続く。N1を4k+1タイプ、N3を4k+3タイプの素数の数とすると以下のようになっている。

Prime    N1    N3
5        1     1
17       3     3
41       6     6
461     44     44
26833   1471   1471
26849   1472   1472
26861   1473 > 1472
26863   1473   1473
26881   1474   1474
26893   1475   1475
26921   1476   1476

この後は N1 < N3がかなり続く

詳しくは小山先生の文献を参照してください。

AZ-prologによるコードは以下のようなソースです（is_p(201)などと２００以上の数を入れてください。素数表をつかっているので（２００以下の素数で30,000以下の素数を発生させています）。

p_l([2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
        101,103,107,109,113,127, 131, 137, 139, 149, 151, 157, 163, 167, 173,179,181,191,193,197,199]).

is_prime(X) :- p_l(L),is_prime_aux(X,L),write(X),write(' is prime number'),nl,!.
is_prime(_) :- !, fail.
is_prime_aux(X,[]) :- !.
is_prime_aux(X,[Y|Rest]) :-
    Z is X mod Y, Z =:= 0, !, fail.
is_prime_aux(X,[Y|Rest]) :- is_prime_aux(X,Rest).

disp(P,X,Y) :- X >= Y, 
    write(P),write(', '),write(X:Y),nl,!.
disp(_,_,_) :- !.

is_p(N) :- N > 30000,cnt_n1(N1),cnt_n3(N3),write(N1:N3),!.
is_p(N) :- is_prime(N), N1 is N + 2, update(N),cnt_n1(NN),cnt_n3(MM),
           write(NN:MM),nl,disp(N,NN,MM),is_p(N1).
is_p(N) :- N1 is N + 2, is_p(N1).

cnt_n1(21).
cnt_n3(24).

up_n1 :- cnt_n1(N),N1 is N + 1, retract(cnt_n1(_)),asserta(cnt_n1(N1)).
up_n3 :- cnt_n3(N),N1 is N + 1, retract(cnt_n3(_)),asserta(cnt_n3(N1)).
update(X) :- P is X mod 4, P =:= 1, up_n1.
update(X) :- P is X mod 4, P =:= 3, up_n3.


do :- p_l(L),do_aux(L,0,0).

do_aux([],N,M) :- write('4n+1:'),write(N),write(', 4n+3:'),write(M),nl,!.
do_aux([X|L],N,M) :- P is X mod 4, P =:= 1, N1 is N + 1,
    disp(X,N1,M),
    do_aux(L,N1,M).
do_aux([X|L],N,M) :- P is X mod 4, P =:= 3, M1 is M + 1,
    do_aux(L,N,M1).
do_aux([X|L],N,M) :- do_aux(L,N,M).