#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/*****************************************************************************/
/*                                                                           */
/*  Equidistribution of roots of an n-th degree congruence over prime moduli */
/*      by Leo Goldmakher, Princeton University Mathematics Dept., 2003      */
/*                                                                           */
/*  This program is to test numerically the conjecture that any n-th degree  */
/*  polynomial f(x) satisfies the following relation:                        */
/*                                                                           */
/*  #{(p,n) s.t. p <= x, f(n) = 0 (mod p), a <= (n/p) < b} ~ (b - a) * pi(x) */
/*                                                                           */
/*  for any a, b between 0 and 1. In other words, we wish to show that the   */
/*  roots of f(x) modulo p, normalized via division by p, are                */
/*  equidistributed on the interval [0,1).                                   */
/*                                                                           */
/*  A brief overview of the variables used:                                  */
/*                                                                           */
/*  "lim" shall denote the upper limit on primes, ie we're counting roots    */
/*  modulo all p <= lim.                                                     */
/*                                                                           */
/*  "size" is an estimate on the number of such roots over all p <= lim.     */
/*  It would be an interesting problem to develop a good theoretic estimate. */
/*                                                                           */
/*  "counter" is what the size should have been, ie the actual number of     */
/*  roots. The program outputs counter at the end of the program; so far     */
/*  it appears that counter is approximately lim / 5.                        */
/*                                                                           */
/*  "dist[]" is the array storing all the normalized roots. Since we don't   */
/*  at present know how to predict the exact number of such roots, we just   */
/*  let the size of dist[] be "size". In all the extra entries of dist[] we  */
/*  store 2's.                                                               */
/*                                                                           */
/*  "r" is the number of elements of dist[] between a and b                  */
/*                                                                           */
/*  "pi" is the the number of primes <= lim (ie pi(lim))                     */
/*                                                                           */
/*  "distribution" is the left side of our principal relation, divided by    */
/*  the right side. Our conjecture is that distribution tends to 1 as lim    */
/*  tends to infinity.                                                       */
/*                                                                           */
/*  The function f(x) is inputted by the user above the main()               */
/*                                                                           */
/*  "index[]" is essentially an array of pointers; we use it to keep track   */
/*  of where entries are once we sort dist[] (using qsort)                   */
/*                                                                           */
/*  "eprdat.fil" is a data file containing, in order, the first 79000 primes */
/*  We input and store these primes into an array "prime[]".                 */
/*                                                                           */
/*  For theoretical work on this problem, we refer the reader to the paper   */
/*  "Equidistribution of roots of a quadratic congruence to prime moduli",   */
/*  by W. Duke, J. B. Friedlander, H. Iwaniec, printed in Ann. of Math. 141  */
/*  (1995), 423-441.                                                         */
/*                                                                           */
/*****************************************************************************/

#define lim 10000
#define size 5000

FILE *fp;
double dist[size];
double poisson[size];

int compfunc(const void *x, const void *y)
{
  if (dist[*(int *)x] < dist[*(int *)y]) return -1;
  else
    if (dist[*(int *)x] == dist[*(int *)y]) return 0;
    else
      return 1;
}

int compfunc1(const void *x, const void *y)
{
  if (poisson[*(int *)x] < poisson[*(int *)y]) return -1;
  else
    if (poisson[*(int *)x] == poisson[*(int *)y]) return 0;
    else
      return 1;
}

long f(long c,long d,long e,long g,long x)
{
  long f = c*x*x*x + d*x*x + e*x + g;
  return f;
}

int main()
{
  long c, d, e, g, p, x, y, r, i, j=0, counter=0, pi=0;
  long A[size], B[size];
  int index[size], index1[size];
  double distribution;
  double a=0.1, b=0.2;
  long eprsize=79000;        /* number of primes in the file eprdat.fil */
  long prime[eprsize];       /* array in which we'll store primes       */
  FILE *fpr, *fnr;

  fpr = fopen("eprdat.fil", "r");
  for (i = 0; i < eprsize; ++i)
    {
      fscanf(fpr, "%ld ", &prime[i]);    /* stores list of first 79000  */
    }                                    /* primes in array prime[]     */
  fclose(fpr);

  fp = fopen("testdata.fil", "w");

  for (c = 1; c <= 5; c++)
    {
      for (d = 0; d < 5; d++)
	{
	  for (e = 0; e < 5; e++)
	    {
	      for (g = 1; g <= 5; g++)
		{
		  j=0;
		  counter=0;
		  pi=0;
		  p = prime[0];
		  for (i=0;i<size;i++)
		    {
		      A[i]=0;
		      B[i]=0;
		      index[i]=0;
		    }

		  while (p <= lim)
		    {
		      for (x = 0; x < p; x++)
			{
			  if ((f(c,d,e,g,x) % p) == 0.0)
			    {
			      B[counter] = p;
			      A[counter] = x;
			      counter = counter + 1;
			    }
			}
		      
		      j = j + 1;
		      pi = pi + 1;
		      p = prime[j];
		    }
		  
		  for (i = 0; i < counter; i++)
		    {
		      dist[i] = 1.0 * A[i] / B[i];
		      index[i] =  i;
		    }
		  
		  for (i = counter; i < size; i++)
		    {
		      dist[i] = 2;
		      index[i] = i;
		    }
		  
		  qsort(index, size, sizeof(int), compfunc);
		  
		  i = 0;
		  while (dist[index[i]] < a) i++;
		  
		  r = 0;
		  while (dist[index[i]] < b)
		    {
		      i++;
		      r++;
		    }
		  
		  distribution = 1.0 * r / (pi * (b - a));

		  if (distribution < 0.5)
		    fprintf(fp, "%ld x^3 + %ld x^2 + %ld x + %ld\t %lf **\n",c,d,e,g,distribution);
		  if (distribution >= 0.5)
		    fprintf(fp, "%ld x^3 + %ld x^2 + %ld x + %ld\t %lf\n",c,d,e,g,distribution);
		}
	    }
	}
    }


  //  /* Now we study whether diffences between adjacent elements */
  //  /* of the dist[] array follow a poissonian distribution.    */
  //  /* For this to work effectively, set a=0, b=1               */

  //  for (i = 0; i < counter - 1; i++)
  //    {
  //      poisson[i] = 1.0 * (dist[index[i + 1]] - dist[index[i]]);
  //      index1[i] = i;
  //    }

  //  for (i = counter - 1; i < size; i++)
  //    {
  //      poisson[i] = 2;
  //      index1[i] = i;
  //    }

  //  qsort(index1, size, sizeof(int), compfunc1);

  /* Now print all those elements of the (hopefully!) equidistributed     */
  /* set, which are located between a and b (ie, if you want ALL elements */
  /* in the set, make sure a=0 and b=1!)                                  */

  //  for (i = 0; i < counter; i++)
  //    {
  //      printf("dist[%ld] = %lf\n", i, dist[index[i]]);
  //    }

  /* Or, print elements of the poisson[] to see whether it IS poisson */

  
  //  fp = fopen("testdata.fil", "w");    /* opens for writing */

  //  for (i = 0; i < counter - 1; i++)
  //    {
  //      printf("poisson[%ld] = %lf\n", i, poisson[index1[i]]);
  //      fprintf(fp, "%ld\t %lf\n", index1[i], poisson[index1[i]]);

  //    }

  //  printf("lim = %ld\n", lim);
  //  printf("counter = %ld\n", counter);
  //  printf("distribution = %lf\n", distribution);

  
  fclose(fp);     

  return 0;
}