{   This is a Pascal program which allows to search numerically            }
{   for local maxima of the packing                                        }
{   density of almost periodic sphere packings. This implementation works  }
{   for packings in dimensions 2,3,4,5                                     }

{   After having made experiments with this program, I conjectured         }
{   that local maxima are                                                  }
{   rational alpha parameters. The program makes a guess for the best      }
{   rational parameter. That critical points are indeed rational is proven }
{   in the paper:                                                          }

{   Oliver Knill, 'Maximizing the packing density on a class of almost     }
{   periodic sphere packings'                                              }
{   Expositiones Mathematicae, 14 (1996), p. 227-246, 1996                 }

{   This Pascal program was written by Oliver Knill in 1995 at Caltech     }
{   We had used more sophisticated programs when writing the paper         }
{   The other programs are however in kind of a mess and I don't post them }                               
{   If you don't change any of the parameters and compile the program      }
{   and run it, then the program will report critical points               }
{   of almost periodic packings with spheres on the integer lattice        }
{   in dimension 3. It will start with spheres of                          }
{   radius=Sqrt(25) and end the search with spheres of radius=sqrt(50).    }
{   By adjusting the waiting parameter (now waitingparameter=3000), the    }
{   scaling factor (now factor=5) and the number of resettings before in   }
{   creasing the radius of the spheres                                     }
{   (now numberoftries=3), one can search differently and more             }
{   persistently, and looking for better local minima.                     }

{   We still have the hope that a more sophisticated version of such a     }
{   program (when running on a more sophisticated computer) might find     }
{   in some dimension (probably bigger than 5) a packing which will have   }
{   record density. If you should decide to continue the search, I wish    }
{   good luck.                                                             }


program kepler(input,output); 
const pi=3.141592653589793;             { a well known constant            }
      np=20000;                         { arraylength                      }
      eee=1/10000000000;                { search treshold                  }
      dimension=3;                      { dimension must be 2,3,4 or 5     } 
      firstr=25;                        { radius^2 for initial search      } 
      lastr=50;                         { radius^2, where the search ends  }
      numberoftries=3;                  { number of times to start new     }
                                        { with same parameters and on      } 
                                        { initial search level             } 

      waitingparameter=3000;            { time until the search is refined }
      search1    = 1.0;                 { initial search factor            }
      factor     = 5;                   { scaling for adaptive search      }
      rand  = 0.54547523948761234;      { a 'random number', used because  }
                                        { because internal random number   }
                                        { generator is deterministic       }
 
type  realarray    = array[1..np] of real;
      vector2      = array[1..2] of integer;
      vector3      = array[1..3] of integer;
      vector4      = array[1..4] of integer;
      vector5      = array[1..5] of integer;
      realvector3  = array[1..3] of real; 
      realvector4  = array[1..4] of real; 
      realvector5  = array[1..5] of real; 
var   p2          :  array[1..np] of vector2;
      p           :  array[1..np] of vector3;
      p4          :  array[1..np] of vector4;
      p5          :  array[1..np] of vector5;
      N           :  integer;    {Number of balls }
      V           :  real;       {Scaling factor for computing    }
                                 {packing density from interval   }
                                 {length. Depends on the dimension}
      max         :  real;       {best density                    }
      r           :  real;       { radius of balls                }
      rn          :  integer;    { sqrt(rn)=r                     }
      a1,a2,a3,a4,a5 :  real;    { middle search point            }
      c1,c2,c3,c4,c5 :  real;    { search vectors                 }
      rec1,rec2,rec3 :  real;    { records                        }
      rec4,rec5      :  real;    { records                        }
      search      :  real;       { search size                    }
      display     :  real;       { density displayed if overcome  }
      u           :  realarray;  { the set R(r) alpha             }
      d           :  real;       { density                        }
      nosucc      :  integer;    { number of nosuccess in loop    }
      showintermed:  boolean;    { show intermediate steps        }
      den         :  integer;    { denominator in rational search }
      bestdens    :  real;       { best density for batch run     }

{  A horrible implementation of modulo which works in our case    }
{  Never bothered to make this nice                               }

function modulo(s: real):real; begin modulo:=s+30000-trunc(s+30000) end;

procedure initialisation2d;
var  ll           : integer;
     vec          : vector2;
     rr           : integer;
begin
  rr:=trunc(r)+1;
  ll:=1; vec[1]:=-rr; vec[2]:=-rr;
  while vec[1]<=r do begin
    while vec[2]<=r do begin
      if ((vec[1]*vec[1]+vec[2]*vec[2]) < rn) and
       ((vec[1]*vec[1]+vec[2]*vec[2]) > 0) then begin
        p2[ll]:=vec; ll:=ll+1;
      end; vec[2]:=vec[2]+1; 
    end; vec[1]:=vec[1]+1; vec[2]:=-rr;
  end;
  N:=ll-1; 
  V:=r*r*pi/4;
end;

procedure initialisation3d;
var  ll           : integer;
     vec          : vector3;
     rr           : integer;
begin
  rr:=trunc(r)+1;
  ll:=1; vec[1]:=-rr; vec[2]:=-rr; vec[3]:=-rr; 
  while vec[1]<=r do begin
    while vec[2]<=r do begin
      while vec[3]<=r do begin
        if ((vec[1]*vec[1]+vec[2]*vec[2]+vec[3]*vec[3]) < rn) and 
         ((vec[1]*vec[1]+vec[2]*vec[2]+vec[3]*vec[3]) > 0) then begin
          p[ll]:=vec; ll:=ll+1; 
        end; vec[3]:=vec[3]+1; 
      end; vec[2]:=vec[2]+1; vec[3]:=-rr;
    end; vec[1]:=vec[1]+1; vec[3]:=-rr; vec[2]:=-rr;
  end; 
  N:=ll-1; 
  V:=r*r*r*pi/6;
end;

procedure initialisation4d;
var  l            : integer;
     vec          : vector4;
     rr           : integer;
begin
  rr:=trunc(r);
  l:=1; vec[1]:=-rr; vec[2]:=-rr; vec[3]:=-rr; vec[4]:=-rr;
  while vec[1]<=r do begin
    while vec[2]<=r do begin
      while vec[3]<=r do begin
        while vec[4]<=r do begin
          if ((vec[1]*vec[1]+vec[2]*vec[2]
              +vec[3]*vec[3]+vec[4]*vec[4]) < rn) and
           ((vec[1]*vec[1]+vec[2]*vec[2]
            +vec[3]*vec[3]+vec[4]*vec[4]) > 0) then begin
            p4[l]:=vec; l:=l+1;
          end; vec[4]:=vec[4]+1;
        end; vec[3]:=vec[3]+1; vec[4]:=-rr;
      end; vec[2]:=vec[2]+1; vec[4]:=-rr; vec[3]:=-rr;
    end; vec[1]:=vec[1]+1; vec[4]:=-rr; vec[3]:=-rr; vec[2]:=-rr;
  end; N:=l-1;
  V:=r*r*r*r*pi*pi/32;
end;

procedure initialisation5d;
var  l            : integer;
     vec          : vector5;
     rr           : integer;
begin
  rr:=trunc(r);
  l:=1; vec[1]:=-rr; vec[2]:=-rr; vec[3]:=-rr; vec[4]:=-rr; vec[5]:=-rr; 
  while vec[1]<=r do begin
    while vec[2]<=r do begin
      while vec[3]<=r do begin
        while vec[4]<=r do begin
          while vec[5]<=r do begin
            if ((vec[1]*vec[1]+vec[2]*vec[2]
                +vec[3]*vec[3]+vec[4]*vec[4]+vec[5]*vec[5]) < rn) and
             ((vec[1]*vec[1]+vec[2]*vec[2]
              +vec[3]*vec[3]+vec[4]*vec[4]+vec[5]*vec[5]) > 0) then begin
              p5[l]:=vec; l:=l+1;
            end; vec[5]:=vec[5]+1;  
          end; vec[4]:=vec[4]+1; vec[5]:=-rr; 
        end; vec[3]:=vec[3]+1; vec[4]:=-rr; vec[5]:=-rr; 
      end; vec[2]:=vec[2]+1; vec[3]:=-rr; vec[4]:=-rr; vec[5]:=-rr; 
    end; vec[1]:=vec[1]+1; vec[2]:=-rr; vec[3]:=-rr; vec[4]:=-rr; vec[5]:=-rr; 
  end; N:=l-1; 
  V:=r*r*r*r*r*pi*pi/60;
end;

{  The following sorting procedure is taken from the Numerical Receipes  }
{  It is only used for coverings and not in this program.                }
procedure sort(n: integer; var ra: realarray);
label 99;
var  lsort,j,i,ir      : integer;
     rra               : real;
begin
  lsort:= (n div 2)+1; ir := n; 
  while true do begin
    if lsort>1 then begin lsort:=lsort-1; rra:=ra[lsort] end
           else begin rra:=ra[ir]; ra[ir]:=ra[1]; ir:=ir-1;
                  if ir=1 then begin ra[1]:=rra; goto 99 end
           end; i:=lsort; j:=lsort+lsort; 
    while j <= ir do begin 
      if j<ir then if ra[j]<ra[j+1] then j:=j+1;
                   if rra<ra[j] then 
                   begin ra[i]:=ra[j]; i:=j; j:=j+j
                   end
              else j:=ir+1 end; ra[i]:=rra
    end;  99:
end;

function min(n: integer; var ra: realarray):real;
var  g  : integer;
     m  : real;
begin
  m:=ra[1]; 
  for g:=2 to n do
  begin
    if ra[g]<m then begin m:=ra[g]; end;
  end;
  min:=m;
end;

function minimum(a,b:real):real; begin if a<b then minimum:=a else minimum:=b end;

procedure showparameter;
begin
  if showintermed=true then 
  begin
    if dimension=2 then 
      writeln('{',rec1:12:10,',', rec2:12:10,'}',d:12:10);
    if dimension=3 then  
      writeln('{',rec1:12:10,',', rec2:12:10,',',rec3:12:10,'}',d:12:10);
    if dimension=4 then 
      writeln('{',rec1:12:10,',', rec2:12:10,','
                 ,rec3:12:10,',', rec4:12:10,'}',d:12:10);
    if dimension=5 then 
      writeln('{',rec1:12:10,',', rec2:12:10,','
                 ,rec3:12:10,',', rec4:12:10,',',rec5:12:10,'}',d:12:10);
  end; 
end; 

procedure findrational;
var k: integer; 
    best,a : real; 
begin
  best:=1.0; 
  for k:=1 to 5000 do
  begin
    a:= minimum(modulo(k*a1),modulo(-k*a1))
       +minimum(modulo(k*a2),modulo(-k*a2)); 
    if dimension>2 then a:=a+minimum(modulo(k*a3),modulo(-k*a3));
    if dimension>3 then a:=a+minimum(modulo(k*a4),modulo(-k*a4));
    if dimension>4 then a:=a+minimum(modulo(k*a5),modulo(-k*a5));
    if a<best then begin best:=a;  den:=k; end; 
  end; 
end;  

procedure showresult;
begin
  if d>bestdens then
  begin
  bestdens:=d; 
  writeln('dimension= ', dimension, '    r*r= ',rn); 
  if dimension=2 then
  begin
    writeln('alpha={',a1:12:10,',', a2:12:10,'}  ');
    writeln('The density is ',d:12:10);
  end;
  if dimension=3 then
  begin
    writeln('alpha={',a1:12:10,',', a2:12:10,',',a3:12:10,'}  ');
    writeln('The density is ',d:12:10);
  end;
  if dimension=4 then
  begin
    writeln('alpha={',a1:12:10,',', a2:12:10,','
               ,a3:12:10,',', a4:12:10,'}  ');
    writeln('The density is ',d:12:10);
  end;
  if dimension=5 then
  begin
    writeln('alpha={',a1:12:10,',', a2:12:10,','
               ,a3:12:10,',', a4:12:10,',',a5:12:10,'}  ');
    writeln('The density is ',d:12:10);
  end;
  findrational;
  writeln('Best denominator guess:', den);
  writeln('The nominators are then integers near:');
  write((a1*den):12:10,' ', (a2*den):12:10,' '); 
  if dimension>2 then write((a3*den):12:10,' '); 
  if dimension>3 then write((a4*den):12:10,' '); 
  if dimension>4 then write((a5*den):12:10,' '); writeln;
  writeln('***************************************************************');
  end; 
end;

procedure run2d;
label  96,95;
var l        : integer;
begin
  c1:=modulo(a1+random(1)); c2:=modulo(a2+random(1)); 
  nosucc:=0;
  while true do begin
    c1:=modulo(c1+random(1)); c2:=modulo(c2+random(1)); 
    for l:=1 to N do 
    begin
      u[l]:=modulo((a1+(c1-1/2)*search)*p2[l,1]
                  +(a2+(c2-1/2)*search)*p2[l,2]);
      if (u[l]<max) then begin nosucc:=nosucc+1; goto 96; end;
    end; 
    d:=min(N,u)*V;
    if d>display then
    begin
      nosucc:=0;
      rec1:=modulo(a1+(c1-1/2)*search);
      rec2:=modulo(a2+(c2-1/2)*search);
      showparameter; 
      max:=d/V; display:=d;
    end;
    96:
    if nosucc>waitingparameter then
    begin
      nosucc:=0;
      a1:=rec1; a2:=rec2; search:=search/factor;
      if search<eee then goto 95; 
      if showintermed=true then writeln('.')
    end;
  end;
  95:
end;

procedure run3d;
label  96,95;
var m        : integer;
begin
  c1:=modulo(a1+random(1)); c2:=modulo(a2+random(1));
  c3:=modulo(a3+random(1));
  nosucc:=0; 
  while true do begin
    c1:=modulo(c1+random(1)); c2:=modulo(c2+random(1));
    c3:=modulo(c3+random(1));
    for m:=1 to N do begin
      u[m]:=modulo((a1+(c1-1/2)*search)*p[m,1]
                  +(a2+(c2-1/2)*search)*p[m,2]
                  +(a3+(c3-1/2)*search)*p[m,3]);
      if (u[m]<max) then begin nosucc:=nosucc+1; goto 96; end;
    end;
    d:=min(N,u)*V;
    if d>display then
    begin
      nosucc:=0;
      rec1:=modulo(a1+(c1-1/2)*search);
      rec2:=modulo(a2+(c2-1/2)*search);
      rec3:=modulo(a3+(c3-1/2)*search);
      showparameter; 
      max:=d/V; display:=d;
    end;
    96:
    if nosucc>waitingparameter then
    begin
      nosucc:=0;
      a1:=rec1; a2:=rec2; a3:=rec3;
      search:=search/factor;
      if search<eee then goto 95; 
      if showintermed=true then writeln('.')
    end;
  end;
  95:
end;

procedure run4d;
label  96,95;
var l        : integer;
begin
  c1:=modulo(a1+random(1)); c2:=modulo(a2+random(1));
  c3:=modulo(a3+random(1)); c4:=modulo(a4+random(1)); 
  nosucc:=0; 
  while true do begin
    c1:=modulo(c1+random(1)); c2:=modulo(c2+random(1));
    c3:=modulo(c3+random(1)); c4:=modulo(c4+random(1));
    for l:=1 to N do begin
      u[l]:=modulo((a1+(c1-1/2)*search)*p4[l,1]
                  +(a2+(c2-1/2)*search)*p4[l,2]
                  +(a3+(c3-1/2)*search)*p4[l,3]
                  +(a4+(c4-1/2)*search)*p4[l,4]);
      if (u[l]<max) then begin nosucc:=nosucc+1; goto 96; end;
    end;
    d:=min(N,u)*V;
    if d>display  then
    begin
      nosucc:=0;
      rec1:=modulo(a1+(c1-1/2)*search);
      rec2:=modulo(a2+(c2-1/2)*search);
      rec3:=modulo(a3+(c3-1/2)*search);
      rec4:=modulo(a4+(c4-1/2)*search);
      showparameter; 
      max:=d/V; display:=d;
    end;
    96:
    if nosucc>waitingparameter then
    begin
      nosucc:=0;
      a1:=rec1; a2:=rec2; a3:=rec3; a4:=rec4;
      search:=search/factor;
      if search<eee then goto 95;
      if showintermed=true then writeln('.')
    end;
  end;
  95:
end;

procedure run5d;
label  96,95;
var l        : integer;
begin 
  c1:=modulo(a1+random(1)); c2:=modulo(a2+random(1));
  c3:=modulo(a3+random(1)); c4:=modulo(a4+random(1));
  c5:=modulo(a5+random(1));
  nosucc:=0;
  while true do begin
    c1:=modulo(c1+random(1)); c2:=modulo(c2+random(1)); 
    c3:=modulo(c3+random(1)); c4:=modulo(c4+random(1));
    c5:=modulo(c5+random(1)); 
    for l:=1 to N do begin
      u[l]:=modulo((a1+(c1-1/2)*search)*p5[l,1]
                  +(a2+(c2-1/2)*search)*p5[l,2]
                  +(a3+(c3-1/2)*search)*p5[l,3]
                  +(a4+(c4-1/2)*search)*p5[l,4]
                  +(a5+(c5-1/2)*search)*p5[l,5]);
      if (u[l]<max) then begin nosucc:=nosucc+1; goto 96; end;
    end;
    d:=min(N,u)*V;
    if d>display then
    begin
      nosucc:=0;
      rec1:=modulo(a1+(c1-1/2)*search);
      rec2:=modulo(a2+(c2-1/2)*search);
      rec3:=modulo(a3+(c3-1/2)*search);
      rec4:=modulo(a4+(c4-1/2)*search);
      rec5:=modulo(a5+(c5-1/2)*search);
      showparameter; 
      max:=d/V; display:=d;
    end;
    96:
    if nosucc>waitingparameter then
    begin
      nosucc:=0;
      a1:=rec1; a2:=rec2; a3:=rec3; a4:=rec4; a5:=rec5; 
      search:=search/factor;
      if search<eee then goto 95;
      if showintermed=true then writeln('.')
    end;
  end;
  95:
end;

procedure initialisation;
begin 
  if dimension=2 then initialisation2d; 
  if dimension=3 then initialisation3d; 
  if dimension=4 then initialisation4d; 
  if dimension=5 then initialisation5d;
  writeln; 
  writeln('################################################################');
  writeln('  The initialisation for the new radius is done! The program    '); 
  writeln('  searches now for critical points with the parameters          ');
  writeln('  Square of the radius: r*r= ',rn,'                             ');
  writeln('  Dimension: d= ', dimension,'                                  ');  
  writeln('  Info: Number of lattice points in Ball(r) is: N= ',N,'        ');      
{ writeln('  Scaling factor: V= ',V:7:6);                                  ')}
  writeln('################################################################');
  max:=display/V;
end; 

procedure run;
begin 
  if dimension=2 then run2d;
  if dimension=3 then run3d;
  if dimension=4 then run4d;
  if dimension=5 then run5d;
end;

procedure batchrun;
var i   : integer; 
    irn : integer; 
  procedure shortinit; 
  begin
    a1:=modulo(random(1)+rand); a2:=modulo(random(1)+2*rand);
    a3:=modulo(random(1));      a4:=modulo(random(1)); 
    a5:=modulo(random(1));
    search:=search1; display:=0.0; max:=0.0; 
  end;
begin
  showintermed:=false;             { Do not display intermediate results }
  for irn:=firstr to lastr do 
  begin
    rn:=irn; r:=sqrt(rn); shortinit;
    initialisation; 
    bestdens:=0.0; 
    for i:=1 to numberoftries do begin 
                       shortinit; 
                       run; 
                       showresult; 
                     end; 
  end; 
end; 

begin { main program }
  batchrun; 
end.