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


/* 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                                  */
/* cc -lm density.c                                        */
/* 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  */


#define pi              3.141592653589793

#define np              300000L


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 packing and covering */
 double r;   /* radius of balls */
 double a1, a2, a3, a4, a5;   /* alpha  */
 long rn;   /* sqrt(rn)=r      */
 realarray u;   /* the set R(r) alpha  */
 double d;   /* density             */
 long dimension;   /* dimension           */


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


void writearray(n, a)
long n;
double *a;
{
  long i;

  for (i = 0; i < n; i++)
    printf("%7.3f ", a[i]);
  putchar('\n');
}


 double squareroot(s)
double s;
{
  return exp(log(s) / 2);
}


 void rationalreadln(a)
double *a;
{
  long u, v;

  scanf("%ld%*[^\n]", &u);
  getchar();
  putchar('/');
  scanf("%ld%*[^\n]", &v);
  getchar();
  *a = (double)u / v;
}


 void datainput()
{
  printf("\nGive the dimension: \n");
  scanf("%ld%*[^\n]", &dimension);
  getchar();
  printf("Give the radius:  \n");
  scanf("%lg%*[^\n]", &r);
  getchar();
  if (r == 0) {
    printf("Give the radius squared: \n");
    scanf("%ld%*[^\n]", &rn);
    getchar();
    r = squareroot((double)rn);
  } else
    rn = (long)floor(r * r + 0.5);
  printf("Give alpha1:  \n");
  scanf("%lg%*[^\n]", &a1);
  getchar();
  if (a1 == 0)
    rationalreadln(&a1);
  printf("Give alpha2:  \n");
  scanf("%lg%*[^\n]", &a2);
  getchar();
  if (a2 == 0)
    rationalreadln(&a2);
  if (dimension > 2) {
    printf("Give alpha3:  \n");
    scanf("%lg%*[^\n]", &a3);
    getchar();
    if (a3 == 0)
      rationalreadln(&a3);
  }
  if (dimension > 3) {
    printf("Give alpha4:  \n");
    scanf("%lg%*[^\n]", &a4);
    getchar();
    if (a4 == 0)
      rationalreadln(&a4);
  }
  if (dimension <= 4)
    return;


  printf("Give alpha5:  \n");
  scanf("%lg%*[^\n]", &a5);
  getchar();
  if (a5 == 0)
    rationalreadln(&a5);
}


 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;
}


 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;
}


 void showresult()
{
  printf("--------------------------------------------------------------\n");
  printf("dimension= %12ld    r*r=  %12ld   radius= %5.3f\n",
	 dimension, rn, r);
  printf("---------------------------------------------------------------\n");
  if (dimension == 2)
    printf("alpha={%12.10f,%12.10f}  %12.10f\n", a1, a2, d);
  if (dimension == 3)
    printf("alpha={%12.10f,%12.10f,%12.10f}  %12.10f\n", a1, a2, a3, d);
  if (dimension == 4)
    printf("alpha={%12.10f,%12.10f,%12.10f,%12.10f}  %12.10f\n",
	   a1, a2, a3, a4, d);
  if (dimension == 5)
    printf("alpha={%12.10f,%12.10f,%12.10f,%12.10f,%12.10f}  %12.10f\n",
	   a1, a2, a3, a4, a5, d);
  printf("--------------------------------------------------------------\n");

}


 void run()
{
  long l, FORLIM;

  if (dimension == 2) {
    FORLIM = N;
    for (l = 0; l < FORLIM; l++)
      u[l] = modulo(a1 * p2[l][0] + a2 * p2[l][1]);
  }
  if (dimension == 3) {
    FORLIM = N;
    for (l = 0; l < FORLIM; l++)
      u[l] = modulo(a1 * p[l][0] + a2 * p[l][1] + a3 * p[l][2]);
  }
  if (dimension == 4) {
    FORLIM = N;
    for (l = 0; l < FORLIM; l++)
      u[l] = modulo(a1 * p4[l][0] + a2 * p4[l][1] + a3 * p4[l][2] + a4 * p4[l]
								    [3]);
  }
  if (dimension == 5) {
    FORLIM = N;
    for (l = 0; l < FORLIM; l++)
      u[l] = modulo(a1 * p5[l][0] + a2 * p5[l][1] + a3 * p5[l][2] + a4 * p5[l]
		      [3] + a5 * p5[l][4]);
  }
  d = min(N, u) * V;
}


 void initialisation()
{
  if (dimension == 2)
    initialisation2d();
  if (dimension == 3)
    initialisation3d();
  if (dimension == 4)
    initialisation4d();
  if (dimension == 5)
    initialisation5d();
  printf(" Number of lattice points:  \n");
  printf("%12ld\n", N);
}


main(argc, argv)
int argc;
char *argv[];
{  /* main program */
  datainput();
  initialisation();
  run();
  showresult();
}




/* End. */