{* First digit distribution of all factors of natural numbers including *}
{* 1 and n, squares double. With Linear Regression and Chi-Square analysis*}
{* Copyright Johan Gerard van der Galien 8-11-2003 johan.van.der.galien@satoconor.com *}

program FACBEN21LR;

var N,a:longword;
    b,c,d,e1,e2,e3,e4,e5,e6,e7,e8,e9,d10,d100,d1000,d10000,
    d100000,d1000000,d10000000,d100000000,f1,f2,f3,f4,f5,f6,
    f7,f8,f9,g,h,i,j:int64;
    x:byte;
    CHI,h1,h2,h3,h4,h5,h6,h7,h8,h9,part1,part2,part3,part4,part5,part6,
    part7,part8,part9,partTOT,sumlogx,sumlogy,sumlogxlogy,sumsquarelogx,
    sumsquarelogy,slope,intersection,sigmax,sigmay,correlation:real;
    R:textfile;

begin
  assign(R,'H:\PASCAL\FACTORIAL\LOGDATA\LOGDATA1.TXT');
  rewrite(R);
  e1:=0;
  e2:=0;
  e3:=0;
  e4:=0;
  e5:=0;
  e6:=0;
  e7:=0;
  e8:=0;
  e9:=0;
  d:=0;
  sumlogx:=0;
  sumsquarelogx:=0;
  for x:= 1 to 9 do
  begin
    sumlogx:=sumlogx+ln(x)/ln(10);
    sumsquarelogx:=sumsquarelogx+sqr(ln(x)/ln(10));
  end;
  writeln(R,'First digit distribution of all factors of natural numbers');
  writeln(R,'including 1 and n, squares double. With Lineair Regression and Chi-Square analysis.');
  for N:=1 to 1000 do
  begin
    h:=N;
    g:=N mod 1000;
    if g=0 then writeln(N);
    if N<10 then
    begin
      if N=1 then e1:=e1+1;
      if N=2 then e2:=e2+1;
      if N=3 then e3:=e3+1;
      if N=4 then e4:=e4+1;
      if N=5 then e5:=e5+1;
      if N=6 then e6:=e6+1;
      if N=7 then e7:=e7+1;
      if N=8 then e8:=e8+1;
      if N=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;
    e1:=e1+1;
    d:=d+2;
    {writeln(R,'N = ',N);
    writeln(R,'factor number ',(d-1),' = ',1);
    writeln(R,'N = ',N);
    writeln(R,'factor number ',d,' = ',N);}
    for a:=2 to trunc(sqrt(N+1)) do
    begin
      b:=N mod a;
      if b=0 then
      begin
        d:=d+1;
        {writeln(R,'N = ',N);
        writeln(R,'factor number ',d,' = ',a);}
        if a<10 then
        begin
          if a=1 then e1:=e1+1;
          if a=2 then e2:=e2+1;
          if a=3 then e3:=e3+1;
          if a=4 then e4:=e4+1;
          if a=5 then e5:=e5+1;
          if a=6 then e6:=e6+1;
          if a=7 then e7:=e7+1;
          if a=8 then e8:=e8+1;
          if a=9 then e9:=e9+1;
        end;
        if (a>=10) and (a<100) then
        begin
          d10:=a 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 (a>=100) and (a<1000) then
        begin
          d100:=a 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 (a>=1000) and (a<10000) then
        begin
          d1000:=a 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 (a>=10000) and (a<100000) then
        begin
          d10000:=a 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 (a>=100000) and (a<1000000) then
        begin
          d100000:=a 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 (a>=1000000) and (a<10000000) then
        begin
          d1000000:=a 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 (a>=10000000) and (a<100000000) then
        begin
          d10000000:=a 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 (a>=100000000) and (a<1000000000) then
        begin
          d100000000:=a 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;
        j:= N div a;
        d:=d+1;
        {writeln(R,'N = ',N);
        writeln(R,'factor number ',d,' = ',j);}
        if j<10 then
        begin
          if j=1 then e1:=e1+1;
          if j=2 then e2:=e2+1;
          if j=3 then e3:=e3+1;
          if j=4 then e4:=e4+1;
          if j=5 then e5:=e5+1;
          if j=6 then e6:=e6+1;
          if j=7 then e7:=e7+1;
          if j=8 then e8:=e8+1;
          if j=9 then e9:=e9+1;
        end;
        if (j>=10) and (j<100) then
        begin
          d10:=j 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 (j>=100) and (j<1000) then
        begin
          d100:=j 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 (j>=1000) and (j<10000) then
        begin
          d1000:=j 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 (j>=10000) and (j<100000) then
        begin
          d10000:=j 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 (j>=100000) and (j<1000000) then
        begin
          d100000:=j 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 (j>=1000000) and (j<10000000) then
        begin
          d1000000:=j 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 (j>=10000000) and (j<100000000) then
        begin
          d10000000:=j 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 (j>=100000000) and (j<1000000000) then
        begin
          d100000000:=j 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;
      end;
    end;
    if (h=100) or (h=1000) or (h=10000) or (h=100000) or (h=1000000)
    or (h=10000000) or (h=100000000) or (h=1000000000) then
    begin
      writeln(R);
      writeln(R,'Cummulative amount of factors under ',h,'=',d);
      part1:=e1/d;
      part2:=e2/d;
      part3:=e3/d;
      part4:=e4/d;
      part5:=e5/d;
      part6:=e6/d;
      part7:=e7/d;
      part8:=e8/d;
      part9:=e9/d;
      partTOT:=part1+part2+part3+part4+part5+part6+part7+part8+part9;
      sumlogy:=ln(part1)/ln(10)+ln(part2)/ln(10)+ln(part3)/ln(10)+
      ln(part4)/ln(10)+ln(part5)/ln(10)+ln(part6)/ln(10)+ln(part7)/ln(10)+
      ln(part8)/ln(10)+ ln(part9)/ln(10);
      sumsquarelogy:=sqr(ln(part1)/ln(10))+sqr(ln(part2)/ln(10))+sqr(ln(part3)/ln(10))+
      sqr(ln(part4)/ln(10))+sqr(ln(part5)/ln(10))+sqr(ln(part6)/ln(10))+sqr(ln(part7)/ln(10))+
      sqr(ln(part8)/ln(10))+sqr(ln(part9)/ln(10));
      sumlogxlogy:=(ln(1)/ln(10))*(ln(part1)/ln(10))+
      (ln(2)/ln(10))*(ln(part2)/ln(10))+(ln(3)/ln(10))*(ln(part3)/ln(10))+
      (ln(4)/ln(10))*(ln(part4)/ln(10))+(ln(5)/ln(10))*(ln(part5)/ln(10))+
      (ln(6)/ln(10))*(ln(part6)/ln(10))+(ln(7)/ln(10))*(ln(part7)/ln(10))+
      (ln(8)/ln(10))*(ln(part8)/ln(10))+(ln(9)/ln(10))*(ln(part9)/ln(10));
      slope:=(sumlogxlogy-(sumlogx*sumlogy/9))/(sumsquarelogx-(sqr(sumlogx)/9));
      intersection:=(sumlogy-slope*sumlogx)/9;
      sigmax:=sqrt((sumsquarelogx-(sqr(sumlogx)/9))/8);
      sigmay:=sqrt((sumsquarelogy-(sqr(sumlogy)/9))/8);
      correlation:=slope*sigmax/sigmay;
      writeln(R,'e(n) is the observed cumulative amount of factors per digit');
      writeln(R,'part1=',part1,' e1=',e1);
      writeln(R,'part2=',part2,' e2=',e2);
      writeln(R,'part3=',part3,' e3=',e3);
      writeln(R,'part4=',part4,' e4=',e4);
      writeln(R,'part5=',part5,' e5=',e5);
      writeln(R,'part6=',part6,' e6=',e6);
      writeln(R,'part7=',part7,' e7=',e7);
      writeln(R,'part8=',part8,' e8=',e8);
      writeln(R,'part9=',part9,' e9=',e9);
      writeln(R,'Total=',partTOT);
      f1:=round((ln(1+1/1)/ln(10))*d);
      f2:=round((ln(1+1/2)/ln(10))*d);
      f3:=round((ln(1+1/3)/ln(10))*d);
      f4:=round((ln(1+1/4)/ln(10))*d);
      f5:=round((ln(1+1/5)/ln(10))*d);
      f6:=round((ln(1+1/6)/ln(10))*d);
      f7:=round((ln(1+1/7)/ln(10))*d);
      f8:=round((ln(1+1/8)/ln(10))*d);
      f9:=round((ln(1+1/9)/ln(10))*d);
      h1:=(e1-f1)*(e1-f1)/f1;
      h2:=(e2-f2)*(e2-f2)/f2;
      h3:=(e3-f3)*(e3-f3)/f3;
      h4:=(e4-f4)*(e4-f4)/f4;
      h5:=(e5-f5)*(e5-f5)/f5;
      h6:=(e6-f6)*(e6-f6)/f6;
      h7:=(e7-f7)*(e7-f7)/f7;
      h8:=(e8-f8)*(e8-f8)/f8;
      h9:=(e9-f9)*(e9-f9)/f9;
      writeln(R,'Expected cumulative amount of factors digit 1=',f1);
      writeln(R,'Expected cumulative amount of factors digit 2=',f2);
      writeln(R,'Expected cumulative amount of factors digit 3=',f3);
      writeln(R,'Expected cumulative amount of factors digit 4=',f4);
      writeln(R,'Expected cumulative amount of factors digit 5=',f5);
      writeln(R,'Expected cumulative amount of factors digit 6=',f6);
      writeln(R,'Expected cumulative amount of factors digit 7=',f7);
      writeln(R,'Expected cumulative amount of factors digit 8=',f8);
      writeln(R,'Expected cumulative amount of factors digit 9=',f9);
      writeln(R,'CHI h1=',h1);
      writeln(R,'CHI h2=',h2);
      writeln(R,'CHI h3=',h3);
      writeln(R,'CHI h4=',h4);
      writeln(R,'CHI h5=',h5);
      writeln(R,'CHI h6=',h6);
      writeln(R,'CHI h7=',h7);
      writeln(R,'CHI h8=',h8);
      writeln(R,'CHI h9=',h9);
      CHI:=h1+h2+h3+h4+h5+h6+h7+h8+h9;
      writeln(R,'CHI-Square=',CHI);
      writeln(R,'Slope =',slope);
      writeln(R,'Intersection =',intersection);
      writeln(R,'Correlation coefficient =',correlation);
    end;
  end;
  close(R);
end.