program density(input,output); 

{ Pascal program to determine the density of a quasi-     }
{ periodic sphere packing.                                }

{ AUTHOR :  Oliver Knill, 1995 Written at Caltech         }
{ See my paper                                            }
{ Expositiones Mathematicae, 14 (1996), p. 227-246, 1996  } 
{ for the mathematics                                     }

{ DESCRIPTION                                             }
{ This program was used as a helper program to play       }
{ around, check other programs. It is not suited          }
{ to find new dense packings.                             }

{ USAGE                                                   }
{ After you compile the program  with a Pascal compiler   }
{ i.e. in UNIX by typing                                  }
{ pc density.pas                                          }
{ you can run it                                          }
{ i.e. in UNIX by typing                                  }
{ ./a.out                                                 }

{ The use of the program is explained by an example       }
{ which shows how the dialogue can go.                    }
{ The dimension is either 2,3,4, or 5                     }

{   Give the dimension:                                   }
{   3                                                     }
{   Give the radius:                                      } 
{   12                                                    }
{   Give alpha1:                                          }
{   0.0472440944882                                       }
{   Give alpha2:                                          }
{   0.496062992126                                        }
{   Give alpha3:                                          }
{   0.616535433071                                        }
{   Number of lattice points:                             }
{   7122                                                  }
{  ----------------------------------------------------   }
{  dimension=      3    r*r=     144   radius= 12.0       }
{  ----------------------------------------------------   }
{  alpha= 0.0472441,0.4960630,0.616535   0.712424         }
{  ----------------------------------------------------   }

{ This density  0.712424 can be found in table 1 in the   }
{ article Oliver Knill, Maximizing the packing density... }
{ Expositiones Mathematicae, 14 (1996), p. 227-246, 1996  } 


const pi=3.141592653589793;   
      np=300000;
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 packing and covering }
      r           :  real; { radius of balls }
      a1,a2,a3,a4,a5 : real;  { alpha  } 
      rn          :  integer; { sqrt(rn)=r      }
      u           :  realarray;  { the set R(r) alpha  }
      d           :  real;       { density             }
      dimension   :  integer;    { dimension           }

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

procedure writearray(n:integer; var a: realarray);
var i : integer; begin for i:=1 to n do write(a[i]:7:3,' '); writeln end;

function squareroot(s: real):real; begin squareroot:=exp(ln(s)/2); end; 

procedure rationalreadln(var a: real); 
var u,v :  integer; 
begin
  readln(u); write('/'); readln(v); a:=u/v;
end;  

procedure datainput; 
begin   
  writeln; 
  writeln('Give the dimension: '); readln(dimension); 
  writeln('Give the radius:  '); readln(r); 
  if r=0 then 
  begin 
    writeln('Give the radius squared: '); readln(rn); r:=squareroot(rn); 
  end else rn:=round(r*r); 
  writeln('Give alpha1:  '); readln(a1); if a1=0 then rationalreadln(a1); 
  writeln('Give alpha2:  '); readln(a2); if a2=0 then rationalreadln(a2); 
  if dimension>2 then begin 
    writeln('Give alpha3:  '); readln(a3); if a3=0 then rationalreadln(a3) end; 
  if dimension>3 then begin 
    writeln('Give alpha4:  '); readln(a4); if a4=0 then rationalreadln(a4) end;
   if dimension>4 then begin 
    writeln('Give alpha5:  '); readln(a5); if a5=0 then rationalreadln(a5) end;


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;

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;

procedure showresult;
begin
  writeln('--------------------------------------------------------------'); 
  writeln('dimension= ', dimension, '    r*r=  ',rn, '   radius= ', r:5:3); 
 writeln('---------------------------------------------------------------'); 
  if dimension=2 then
    writeln('alpha={',a1:12:10,',', a2:12:10,'}  ',d:12:10);
  if dimension=3 then
    writeln('alpha={',a1:12:10,',', a2:12:10,',',a3:12:10,'}  ',d:12:10);
  if dimension=4 then
    writeln('alpha={',a1:12:10,',', a2:12:10,','
               ,a3:12:10,',', a4:12:10,'}  ',d:12:10);
  if dimension=5 then
    writeln('alpha={',a1:12:10,',', a2:12:10,','
               ,a3:12:10,',', a4:12:10,',',a5:12:10,'}  ',d:12:10);
  writeln('--------------------------------------------------------------'); 

end;

procedure run;
var l     : integer;
begin
  if dimension=2 then 
    for l:=1 to N do u[l]:=modulo(a1*p2[l,1]+a2*p2[l,2]);
  if dimension=3 then 
    for l:=1 to N do u[l]:=modulo(a1*p[l,1]+a2*p[l,2]+a3*p[l,3]);
  if dimension=4 then 
    for l:=1 to N do u[l]:=modulo(a1*p4[l,1]+a2*p4[l,2]+a3*p4[l,3]+a4*p4[l,4]);
  if dimension=5 then 
    for l:=1 to N do u[l]:=modulo(a1*p5[l,1]+a2*p5[l,2]
                                 +a3*p5[l,3]+a4*p5[l,4]+a5*p5[l,5]);
  d:=min(N,u)*V;
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(' Number of lattice points:  '); writeln(N); 
end; 

begin { main program }
  datainput; 
  initialisation;
  run; 
  showresult; 
end.