/* Output from p2c, the Pascal-to-C translator */
/* From input file "kepler.pas" */


/*   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.                                                             */



#include <limits.h>                           /*  for INT_MAX         */
#include <stdlib.h>                           


#define pi              3.141592653589793
    /* a well known constant            */

#define np              20000   /* arraylength                      */

#define eee             (1 / 10000000000.0)
    /* search treshold                  */

#define dimension       5   /* dimension must be 2,3,4 or 5     */
#define firstr          25   /* radius^2 for initial search      */
#define lastr           50   /* radius^2, where the search ends  */
#define numberoftries   3    /* number of times to start new     */
/* with same parameters and on      */
/* initial search level             */

#define waitingparameter  300  /* time until the search is refined */

#define search1         1.0   /* initial search factor            */

#define factor          5   /* scaling for adaptive search      */

#define randseed        0.95532523948761234
/* a 'random number', used because  */
/* because internal random number   */
/* generator is deterministic       */

typedef double realarray[np];
typedef long vector2[2];
typedef long vector3[3];
typedef long vector4[4];
typedef long vector5[5];
typedef double realvector3[3];
typedef double realvector4[4];
typedef double realvector5[5];


 vector2 p2[np];
 vector3 p[np];
 vector4 p4[np];
 vector5 p5[np];
 long N;   /*Number of balls */
 double V;   /*Scaling factor for computing    */
/*packing density from interval   */
/*length. Depends on the dimension*/
 double max;   /*best density                    */
 double r;   /* radius of balls                */
 long rn;   /* sqrt(rn)=r                     */
 double a1, a2, a3, a4, a5;   /* middle search point            */
 double c1, c2, c3, c4, c5;   /* search vectors                 */
 double rec1, rec2, rec3;   /* records                        */
 double rec4, rec5;   /* records                        */
 double search;   /* search size                    */
 double display;   /* density displayed if overcome  */
 realarray u;   /* the set R(r) alpha             */
 double d;   /* density                        */
 long nosucc;   /* number of nosuccess in loop    */
 int showintermed;   /* show intermediate steps        */
 long den;   /* denominator in rational search */
 double bestdens;   /* best density for batch run     */


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

double modulo(s) double s; { return (s + 3000000 - (long)(s + 3000000)); }

void initialisation2d()
{
  long ll;
  vector2 vec;
  long rr;

  rr = (long)r + 1;
  ll = 1;
  vec[0] = -rr;
  vec[1] = -rr;
  while (vec[0] <= r) {
    while (vec[1] <= r) {
      if (vec[0] * vec[0] + vec[1] * vec[1] < rn &&
	  vec[0] * vec[0] + vec[1] * vec[1] > 0) {
	memcpy(p2[ll - 1], vec, sizeof(vector2));
	ll++;
      }
      vec[1]++;
    }
    vec[0]++;
    vec[1] = -rr;
  }
  N = ll - 1;
  V = r * r * pi / 4;
}


void initialisation3d()
{
  long ll;
  vector3 vec;
  long rr;

  rr = (long)r + 1;
  ll = 1;
  vec[0] = -rr;
  vec[1] = -rr;
  vec[2] = -rr;
  while (vec[0] <= r) {
    while (vec[1] <= r) {
      while (vec[2] <= r) {
	if (vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] < rn &&
	    vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] > 0) {
	  memcpy(p[ll - 1], vec, sizeof(vector3));
	  ll++;
	}
	vec[2]++;
      }
      vec[1]++;
      vec[2] = -rr;
    }
    vec[0]++;
    vec[2] = -rr;
    vec[1] = -rr;
  }
  N = ll - 1;
  V = r * r * r * pi / 6;
}


void initialisation4d()
{
  long l;
  vector4 vec;
  long rr;

  rr = (long)r;
  l = 1;
  vec[0] = -rr;
  vec[1] = -rr;
  vec[2] = -rr;
  vec[3] = -rr;
  while (vec[0] <= r) {
    while (vec[1] <= r) {
      while (vec[2] <= r) {
	while (vec[3] <= r) {
	  if (vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] +
	      vec[3] * vec[3] < rn && vec[0] * vec[0] + vec[1] * vec[1] +
				      vec[2] * vec[2] + vec[3] * vec[3] > 0) {
	    memcpy(p4[l - 1], vec, sizeof(vector4));
	    l++;
	  }
	  vec[3]++;
	}
	vec[2]++;
	vec[3] = -rr;
      }
      vec[1]++;
      vec[3] = -rr;
      vec[2] = -rr;
    }
    vec[0]++;
    vec[3] = -rr;
    vec[2] = -rr;
    vec[1] = -rr;
  }
  N = l - 1;
  V = r * r * r * r * pi * pi / 32;
}


void initialisation5d()
{
  long l;
  vector5 vec;
  long rr;

  rr = (long)r;
  l = 1;
  vec[0] = -rr;
  vec[1] = -rr;
  vec[2] = -rr;
  vec[3] = -rr;
  vec[4] = -rr;
  while (vec[0] <= r) {
    while (vec[1] <= r) {
      while (vec[2] <= r) {
	while (vec[3] <= r) {
	  while (vec[4] <= r) {
	    if (vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] +
		vec[3] * vec[3] + vec[4] * vec[4] < rn &&
		vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] +
		vec[3] * vec[3] + vec[4] * vec[4] > 0) {
	      memcpy(p5[l - 1], vec, sizeof(vector5));
	      l++;
	    }
	    vec[4]++;
	  }
	  vec[3]++;
	  vec[4] = -rr;
	}
	vec[2]++;
	vec[3] = -rr;
	vec[4] = -rr;
      }
      vec[1]++;
      vec[2] = -rr;
      vec[3] = -rr;
      vec[4] = -rr;
    }
    vec[0]++;
    vec[1] = -rr;
    vec[2] = -rr;
    vec[3] = -rr;
    vec[4] = -rr;
  }
  N = l - 1;
  V = r * r * r * r * r * pi * pi / 60;
}


/*  The following sorting procedure is taken from the Numerical Receipes  */
/*  It is only used for coverings and not in this program.                */
void sort(n, ra)
long n;
double *ra;
{
  long lsort, j, i, ir;
  double rra;

  lsort = n / 2 + 1;
  ir = n;
  while (1==1) {
    if (lsort > 1) {
      lsort--;
      rra = ra[lsort - 1];
    } else {
      rra = ra[ir - 1];
      ra[ir - 1] = ra[0];
      ir--;
      if (ir == 1) {
	ra[0] = rra;
	goto _L99;
      }
    }
    i = lsort;
    j = lsort + lsort;
    while (j <= ir) {
      if (j < ir) {
	if (ra[j - 1] < ra[j])
	  j++;
      }
      if (rra < ra[j - 1]) {
	ra[i - 1] = ra[j - 1];
	i = j;
	j += j;
      } else
	j = ir + 1;
    }
    ra[i - 1] = rra;
  }
_L99: ;
}


 double min(n, ra)
long n;
double *ra;
{
  long g;
  double m;

  m = ra[0];
  for (g = 1; g < n; g++) {
    if (ra[g] < m)
      m = ra[g];
  }
  return m;
}


 double minimum(a, b)
double a, b;
{
  if (a < b)
    return a;
  else
    return b;
}


void showparameter()
{
  if (showintermed != 1)
    return;
  if (dimension == 2)
    printf("{%12.10f,%12.10f}%12.10f\n", rec1, rec2, d);
  if (dimension == 3)
    printf("{%12.10f,%12.10f,%12.10f}%12.10f\n", rec1, rec2, rec3, d);
  if (dimension == 4)
    printf("{%12.10f,%12.10f,%12.10f,%12.10f}%12.10f\n",
	   rec1, rec2, rec3, rec4, d);
  if (dimension == 5)
    printf("{%12.10f,%12.10f,%12.10f,%12.10f,%12.10f}%12.10f\n",
	   rec1, rec2, rec3, rec4, rec5, d);
}


void findrational()
{
  long k;
  double best, a;

  best = 1.0;
  for (k = 1; k <= 5000; k++) {
    a = minimum(modulo(k * a1), modulo(-k * a1)) + minimum(modulo(k * a2),
							   modulo(-k * a2));
    if (dimension > 2)
      a += minimum(modulo(k * a3), modulo(-k * a3));
    if (dimension > 3)
      a += minimum(modulo(k * a4), modulo(-k * a4));
    if (dimension > 4)
      a += minimum(modulo(k * a5), modulo(-k * a5));
    if (a < best) {
      best = a;
      den = k;
    }
  }
}


void showresult()
{
  if (d <= bestdens)
    return;
  bestdens = d;
  printf("dimension= %12ld    r*r= %12ld"); printf("%d\n",rn);
  if (dimension == 2) {
    printf("alpha={%12.10f,%12.10f}  \n", a1, a2);
    printf("The density is %12.10f\n", d);
  }
  if (dimension == 3) {
    printf("alpha={%12.10f,%12.10f,%12.10f}  \n", a1, a2, a3);
    printf("The density is %12.10f\n", d);
  }
  if (dimension == 4) {
    printf("alpha={%12.10f,%12.10f,%12.10f,%12.10f}  \n", a1, a2, a3, a4);
    printf("The density is %12.10f\n", d);
  }
  if (dimension == 5) {
    printf("alpha={%12.10f,%12.10f,%12.10f,%12.10f,%12.10f}  \n",
	   a1, a2, a3, a4, a5);
    printf("The density is %12.10f\n", d);
  }
  findrational();
  printf("Best denominator guess:%12ld\n", den);
  printf("The nominators are then integers near:\n");
  printf("%12.10f %12.10f ", a1 * den, a2 * den);
  if (dimension > 2)
    printf("%12.10f ", a3 * den);
  if (dimension > 3)
    printf("%12.10f ", a4 * den);
  if (dimension > 4)
    printf("%12.10f ", a5 * den);
  printf("\n***************************************************************\n");
}


void run2d()
{
  long l, FORLIM;

  c1 = modulo(a1 +  rand()*1.0*pi/INT_MAX);
  c2 = modulo(a2 +  rand()*1.0*pi/INT_MAX);
  nosucc = 0;
  while (1==1) {
    c1 = modulo(c1 + rand()*1.0*pi/INT_MAX);
    c2 = modulo(c2 + rand()*1.0*pi/INT_MAX);
    FORLIM = N;
    for (l = 0; l < FORLIM; l++) {
      u[l] = modulo((a1 + (c1 - 1.0 / 2) * search) * p2[l][0] +
		    (a2 + (c2 - 1.0 / 2) * search) * p2[l][1]);
      if (u[l] < max) {
	nosucc++;
	goto _L96;
      }
    }
    d = min(N, u) * V;
    if (d > display) {
      nosucc = 0;
      rec1 = modulo(a1 + (c1 - 1.0 / 2) * search);
      rec2 = modulo(a2 + (c2 - 1.0 / 2) * search);
      showparameter();
      max = d / V;
      display = d;
    }
_L96:
    if (nosucc <= waitingparameter)
      continue;
    nosucc = 0;
    a1 = rec1;
    a2 = rec2;
    search /= factor;
    if (search < eee)
      goto _L95;
    if (showintermed == 1)
      printf(".\n");
  }
_L95: ;
}


void run3d()
{
  long m, FORLIM;

  c1 = modulo(a1 +  rand()*1.0*pi/INT_MAX);
  c2 = modulo(a2 +  rand()*1.0*pi/INT_MAX);
  c3 = modulo(a3 +  rand()*1.0*pi/INT_MAX);
  nosucc = 0;
  while (1==1) {
    c1 = modulo(c1 +  rand()*1.0*pi/INT_MAX);
    c2 = modulo(c2 +  rand()*1.0*pi/INT_MAX);
    c3 = modulo(c3 +  rand()*1.0*pi/INT_MAX);
    FORLIM = N;
    for (m = 0; m < FORLIM; m++) {
      u[m] = modulo((a1 + (c1 - 1.0 / 2) * search) * p[m][0] +
		    (a2 + (c2 - 1.0 / 2) * search) * p[m][1] +
		    (a3 + (c3 - 1.0 / 2) * search) * p[m][2]);
      if (u[m] < max) {
	nosucc++;
	goto _L96;
      }
    }
    d = min(N, u) * V;
    if (d > display) {
      nosucc = 0;
      rec1 = modulo(a1 + (c1 - 1.0 / 2) * search);
      rec2 = modulo(a2 + (c2 - 1.0 / 2) * search);
      rec3 = modulo(a3 + (c3 - 1.0 / 2) * search);
      showparameter();
      max = d / V;
      display = d;
    }
_L96:
    if (nosucc <= waitingparameter)
      continue;
    nosucc = 0;
    a1 = rec1;
    a2 = rec2;
    a3 = rec3;
    search /= factor;
    if (search < eee)
      goto _L95;
    if (showintermed == 1)
      printf(".\n");
  }
_L95: ;
}


void run4d()
{
  long l, FORLIM;

  c1 = modulo(a1 +  rand()*1.0*pi/INT_MAX);
  c2 = modulo(a2 +  rand()*1.0*pi/INT_MAX);
  c3 = modulo(a3 +  rand()*1.0*pi/INT_MAX);
  c4 = modulo(a4 +  rand()*1.0*pi/INT_MAX);
  nosucc = 0;
  while (1==1) {
    c1 = modulo(c1 +  rand()*1.0*pi/INT_MAX);
    c2 = modulo(c2 +  rand()*1.0*pi/INT_MAX);
    c3 = modulo(c3 +  rand()*1.0*pi/INT_MAX);
    c4 = modulo(c4 +  rand()*1.0*pi/INT_MAX);
    FORLIM = N;
    for (l = 0; l < FORLIM; l++) {
      u[l] = modulo((a1 + (c1 - 1.0 / 2) * search) * p4[l][0] +
		    (a2 + (c2 - 1.0 / 2) * search) * p4[l][1] +
		    (a3 + (c3 - 1.0 / 2) * search) * p4[l][2] +
		    (a4 + (c4 - 1.0 / 2) * search) * p4[l][3]);
      if (u[l] < max) {
	nosucc++;
	goto _L96;
      }
    }
    d = min(N, u) * V;
    if (d > display) {
      nosucc = 0;
      rec1 = modulo(a1 + (c1 - 1.0 / 2) * search);
      rec2 = modulo(a2 + (c2 - 1.0 / 2) * search);
      rec3 = modulo(a3 + (c3 - 1.0 / 2) * search);
      rec4 = modulo(a4 + (c4 - 1.0 / 2) * search);
      showparameter();
      max = d / V;
      display = d;
    }
_L96:
    if (nosucc <= waitingparameter)
      continue;
    nosucc = 0;
    a1 = rec1;
    a2 = rec2;
    a3 = rec3;
    a4 = rec4;
    search /= factor;
    if (search < eee)
      goto _L95;
    if (showintermed == 1)
      printf(".\n");
  }
_L95: ;
}


void run5d()
{
  long l, FORLIM;

  c1 = modulo(a1 +  rand()*1.0*pi/INT_MAX);
  c2 = modulo(a2 +  rand()*1.0*pi/INT_MAX);
  c3 = modulo(a3 +  rand()*1.0*pi/INT_MAX);
  c4 = modulo(a4 +  rand()*1.0*pi/INT_MAX);
  c5 = modulo(a5 +  rand()*1.0*pi/INT_MAX);
  nosucc = 0;
  while (1==1) {
    c1 = modulo(c1 +  rand()*1.0*pi/INT_MAX);
    c2 = modulo(c2 +  rand()*1.0*pi/INT_MAX);
    c3 = modulo(c3 +  rand()*1.0*pi/INT_MAX);
    c4 = modulo(c4 +  rand()*1.0*pi/INT_MAX);
    c5 = modulo(c5 +  rand()*1.0*pi/INT_MAX);
    FORLIM = N;
    for (l = 0; l < FORLIM; l++) {
      u[l] = modulo((a1 + (c1 - 1.0 / 2) * search) * p5[l][0] +
		    (a2 + (c2 - 1.0 / 2) * search) * p5[l][1] +
		    (a3 + (c3 - 1.0 / 2) * search) * p5[l][2] +
		    (a4 + (c4 - 1.0 / 2) * search) * p5[l][3] +
		    (a5 + (c5 - 1.0 / 2) * search) * p5[l][4]);
      if (u[l] < max) {
	nosucc++;
	goto _L96;
      }
    }
    d = min(N, u) * V;
    if (d > display) {
      nosucc = 0;
      rec1 = modulo(a1 + (c1 - 1.0 / 2) * search);
      rec2 = modulo(a2 + (c2 - 1.0 / 2) * search);
      rec3 = modulo(a3 + (c3 - 1.0 / 2) * search);
      rec4 = modulo(a4 + (c4 - 1.0 / 2) * search);
      rec5 = modulo(a5 + (c5 - 1.0 / 2) * search);
      showparameter();
      max = d / V;
      display = d;
    }
_L96:
    if (nosucc <= waitingparameter)
      continue;
    nosucc = 0;
    a1 = rec1;
    a2 = rec2;
    a3 = rec3;
    a4 = rec4;
    a5 = rec5;
    search /= factor;
    if (search < eee)
      goto _L95;
    if (showintermed == 1)
      printf(".\n");
  }
_L95: ;
}


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


void run()
{
  if (dimension == 2)
    run2d();
  if (dimension == 3)
    run3d();
  if (dimension == 4)
    run4d();
  if (dimension == 5)
    run5d();
}


void shortinit()
{
  a1 = modulo( rand()*1.0*pi/INT_MAX + randseed);
  a2 = modulo( rand()*1.0*pi/INT_MAX + 2 * randseed);
  a3 = modulo((double) rand()*1.0*pi/INT_MAX);
  a4 = modulo((double) rand()*1.0*pi/INT_MAX);
  a5 = modulo((double) rand()*1.0*pi/INT_MAX);
  search = search1;
  display = 0.0;
  max = 0.0;
}


void batchrun()
{
  long i, irn;

  showintermed = 0;  
  for (irn = firstr; irn <= lastr; irn++) {
    rn = irn;
    r = sqrt((double)rn);
    shortinit();
    initialisation();
    bestdens = 0.0;
    for (i = 1; i <= numberoftries; i++) {
      shortinit();
      run();
      showresult();
    }
  }
}


main(argc, argv)
int argc;
char *argv[];
{  
  batchrun();
}



/* End. */