{* Calculates the first digit distribution of the primes for N = a power of 10 *}
{* also calculates the chi-square value for the deviation from the expected occurence *}
{* of 1/9*100% of each first digit wore the prime numbers to be random. *}
{* Copyright Johan Gerard van der Galien 27 August 2002 johan.van.der.galien@satoconor.com *}

program prime9;
var N,a,b,c,d,e1,e2,e3,e4,e5,e6,e7,e8,e9,d10,d100,d1000,d10000,
    d100000,d1000000,d10000000,d100000000,f,prime:longint;
    CHI,h1,h2,h3,h4,h5,h6,h7,h8,h9,perc1,perc2,perc3,perc4,perc5,perc6,
    perc7,perc8,perc9,percTOT:real;
    R:text;
label 1,2,3;

begin
  assign(R,'D:\PASCAL\LOGDATA\PRIME9.TXT');
  rewrite(R);
1:N:=N+1;
  if N=1000000000 then goto 3;
  if N=1 then goto 1;
  if N=2 then
  begin
    prime:=2;
    goto 2;
  end;
  if N=3 then
  begin
    prime:=3;
    goto 2;
  end;
  c:=trunc(sqrt(N))+1;
  for a:=2 to c do
  begin
    b:=N mod a;
    if b=0 then goto 1;
  end;
2:d:=d+1;
  prime:=N;
  if N<10 then
  begin
    if prime=1 then e1:=e1+1;
    if prime=2 then e2:=e2+1;
    if prime=3 then e3:=e3+1;
    if prime=4 then e4:=e4+1;
    if prime=5 then e5:=e5+1;
    if prime=6 then e6:=e6+1;
    if prime=7 then e7:=e7+1;
    if prime=8 then e8:=e8+1;
    if prime=9 then e9:=e9+1;
  end;
  if (N>10) and (N<100) then
  begin
    d10:=N div 10;
    if d10=1 then e1:=e1+1;
    if d10=2 then e2:=e2+1;
    if d10=3 then e3:=e3+1;
    if d10=4 then e4:=e4+1;
    if d10=5 then e5:=e5+1;
    if d10=6 then e6:=e6+1;
    if d10=7 then e7:=e7+1;
    if d10=8 then e8:=e8+1;
    if d10=9 then e9:=e9+1;
  end;
  if (N>100) and (N<1000) then
  begin
    d100:=N div 100;
    if d100=1 then e1:=e1+1;
    if d100=2 then e2:=e2+1;
    if d100=3 then e3:=e3+1;
    if d100=4 then e4:=e4+1;
    if d100=5 then e5:=e5+1;
    if d100=6 then e6:=e6+1;
    if d100=7 then e7:=e7+1;
    if d100=8 then e8:=e8+1;
    if d100=9 then e9:=e9+1;
  end;
  if (N>1000) and (N<10000) then
  begin
    d1000:=N div 1000;
    if d1000=1 then e1:=e1+1;
    if d1000=2 then e2:=e2+1;
    if d1000=3 then e3:=e3+1;
    if d1000=4 then e4:=e4+1;
    if d1000=5 then e5:=e5+1;
    if d1000=6 then e6:=e6+1;
    if d1000=7 then e7:=e7+1;
    if d1000=8 then e8:=e8+1;
    if d1000=9 then e9:=e9+1;
  end;
  if (N>10000) and (N<100000) then
  begin
    d10000:=N div 10000;
    if d10000=1 then e1:=e1+1;
    if d10000=2 then e2:=e2+1;
    if d10000=3 then e3:=e3+1;
    if d10000=4 then e4:=e4+1;
    if d10000=5 then e5:=e5+1;
    if d10000=6 then e6:=e6+1;
    if d10000=7 then e7:=e7+1;
    if d10000=8 then e8:=e8+1;
    if d10000=9 then e9:=e9+1;
  end;
  if (N>100000) and (N<1000000) then
  begin
    d100000:=N div 100000;
    if d100000=1 then e1:=e1+1;
    if d100000=2 then e2:=e2+1;
    if d100000=3 then e3:=e3+1;
    if d100000=4 then e4:=e4+1;
    if d100000=5 then e5:=e5+1;
    if d100000=6 then e6:=e6+1;
    if d100000=7 then e7:=e7+1;
    if d100000=8 then e8:=e8+1;
    if d100000=9 then e9:=e9+1;
  end;
  if (N>1000000) and (N<10000000) then
  begin
    d1000000:=N div 1000000;
    if d1000000=1 then e1:=e1+1;
    if d1000000=2 then e2:=e2+1;
    if d1000000=3 then e3:=e3+1;
    if d1000000=4 then e4:=e4+1;
    if d1000000=5 then e5:=e5+1;
    if d1000000=6 then e6:=e6+1;
    if d1000000=7 then e7:=e7+1;
    if d1000000=8 then e8:=e8+1;
    if d1000000=9 then e9:=e9+1;
  end;
  if (N>10000000) and (N<100000000) then
  begin
    d10000000:=N div 10000000;
    if d10000000=1 then e1:=e1+1;
    if d10000000=2 then e2:=e2+1;
    if d10000000=3 then e3:=e3+1;
    if d10000000=4 then e4:=e4+1;
    if d10000000=5 then e5:=e5+1;
    if d10000000=6 then e6:=e6+1;
    if d10000000=7 then e7:=e7+1;
    if d10000000=8 then e8:=e8+1;
    if d10000000=9 then e9:=e9+1;
  end;
  if (N>100000000) and (N<1000000000) then
  begin
    d100000000:=N div 100000000;
    if d100000000=1 then e1:=e1+1;
    if d100000000=2 then e2:=e2+1;
    if d100000000=3 then e3:=e3+1;
    if d100000000=4 then e4:=e4+1;
    if d100000000=5 then e5:=e5+1;
    if d100000000=6 then e6:=e6+1;
    if d100000000=7 then e7:=e7+1;
    if d100000000=8 then e8:=e8+1;
    if d100000000=9 then e9:=e9+1;
  end;
  goto 1;
3:writeln(R,'Amount of prime numbers under ',N,'=',d);
  perc1:=100*e1/d;
  perc2:=100*e2/d;
  perc3:=100*e3/d;
  perc4:=100*e4/d;
  perc5:=100*e5/d;
  perc6:=100*e6/d;
  perc7:=100*e7/d;
  perc8:=100*e8/d;
  perc9:=100*e9/d;
  percTOT:=perc1+perc2+perc3+perc4+perc5+perc6+perc7+perc8+perc9;
  writeln(R,'e(n) is the observed amount of prime numbers per digit');
  writeln(R,'perc1=',perc1,' e1=',e1);
  writeln(R,'perc2=',perc2,' e2=',e2);
  writeln(R,'perc3=',perc3,' e3=',e3);
  writeln(R,'perc4=',perc4,' e4=',e4);
  writeln(R,'perc5=',perc5,' e5=',e5);
  writeln(R,'perc6=',perc6,' e6=',e6);
  writeln(R,'perc7=',perc7,' e7=',e7);
  writeln(R,'perc8=',perc8,' e8=',e8);
  writeln(R,'perc9=',perc9,' e9=',e9);
  writeln(R,'Total=',percTOT);
  if N=10 then f:=1 else f:=round(d/9);
  h1:=SQR(e1-f)/f;
  h2:=SQR(e2-f)/f;
  h3:=SQR(e3-f)/f;
  h4:=SQR(e4-f)/f;
  h5:=SQR(e5-f)/f;
  h6:=SQR(e6-f)/f;
  h7:=SQR(e7-f)/f;
  h8:=SQR(e8-f)/f;
  h9:=SQR(e9-f)/f;
  writeln(R,'Expected amount of prime numbers per digit=',f);
  CHI:=h1+h2+h3+h4+h5+h6+h7+h8+h9;
  writeln(R,'CHI-square=',CHI);
close(R);
end.