#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <unistd.h>

//------------------------------------------------------------------------
// The program counts points on a (singular) double octic for the
// first 25 primes (between 2 and 97). It then compares the results
// modulo p with the coefficients of weight four newforms. It displays
// the newform which is most likely to occur in the middle cohomology
// of the double octic.
// It is possible to play around with double octics by trying different
// equations for branchlocus = ...
// You could also add a row with 26 zeroes to w4modforms.dat. If this
// "modular form" is evaluated to be the most likely one then probably
// your double octic is rational.
//------------------------------------------------------------------------

int
main()
{
  //--------------------------------------------------------------
  // determine the prime numbers < 100 by simple Eratosthenes
  //--------------------------------------------------------------

  cout << "Determining prime numbers..." << endl;

  bool isprime[100];
  for( long i = 1; i < 100; i++ )
    isprime[i] = true;
  isprime[1] = false;
  for( long i = 2; i < 100; i++ )
    for( long j = i*i; j < 100; j+=i )
      isprime[j] = false;

  //--------------------------------------------------------------
  // store the computed primes into the primes[] array
  //--------------------------------------------------------------

  long primes[25];   // the first 25 primes
  long j = 1;
  for( long i = 0; i < 25; i++ )
    {
      while( !isprime[j] )
	j++;
      primes[i] = j;
      j++;
    }

  //--------------------------------------------------------------
  // load table with weight 4 modular forms
  //--------------------------------------------------------------

  cout << "Loading table of modular forms..." << endl;

  long thislevel;
  long level[4000];
  long coeff[4000][100];
  long nummodforms = 0;

  ifstream modforms( "w4modforms.dat", ios::in );
  if( !modforms ) {
    cerr << "'w4modforms.dat' not found." << endl;
    exit( 1 );
  }

  while( modforms >> thislevel )
    {
      level[nummodforms] = thislevel;
      for( long i = 0; i < 25; i++ )
	modforms >> coeff[nummodforms][i];
      nummodforms++;
    }
  modforms.close();

  //--------------------------------------------------------------
  // count points on the double octic
  //--------------------------------------------------------------

  cout << "----------------------------------------------" << endl;
  cout << "Counting points on the double octic..." << endl;

  bool issquare_p[100];   // issquare_p[i] is true if i is a square mod p
  long branchlocus;       // evaluated equation of the octic (i.e., the branch locus)
  long numpoints[25];     // number of points on the double octic for the first 25 primes
  long x,y,z,t;           // coordinates (x:y:z:t) on P3
  long p;                 // current prime

  // count points for prime after prime
  for( long i = 0; i < 25; i++ )
    {
      numpoints[i] = 0;
      p = primes[i];

      // Create table of squares mod p
      for( long j = 0; j < p; j++ )
	issquare_p[j] = false;

      for( long j = 0; j < p; j++ )
	issquare_p[j*j%p] = true;

      // main loop over the points in P3(x:y:z:t)
      for( x = 0; x <= 1; x++ )
	for( y = 0; y <= ( x == 0 ? 1 : p-1 ); y++ )
	  for( z = 0; z <= ( x == 0 && y == 0 ? 1 : p-1 ); z++ )
	    for( t = ( x == 0 && y == 0 && z == 0 ? 1 : 0 );
		 t <= ( x == 0 && y == 0 && z == 0 ? 1 : p-1 ); t++ )
	      {
		// evaluate the equation for the octic surface
		// the sample equation produces a rigid double octic of level 8
		// CAUTION: Be sure to reduce mod p in time to avoid overflow
		branchlocus = ( x *
				y % p *
				z *
				t % p *
				( x + y ) *
				( y + z ) % p *
				( z + t ) *
				( t + x ) ) % p;

		// add to the number of points depending on the value of the branch locus
		if( branchlocus == 0 )
		  numpoints[i]++;
		else
		  {
		    // CAUTION: in C++, "%" reduces negative numbers to negative numbers
		    if( branchlocus < 0 )
		      branchlocus += p;
		    if( issquare_p[branchlocus] )
		      numpoints[i] += 2;
		  }
	      }
      cout << "p = " << p;
      if( p < 10 )
	cout << " ";
      cout << " : " << numpoints[i] << " points on the double octic." << endl;
    }

  //--------------------------------------------------------------
  // compare numbers of points with coefficients of modular forms
  //--------------------------------------------------------------

  cout << "----------------------------------------------" << endl;
  cout << "Comparing with coefficients of modular forms..." << endl;

  long oldlevel = 0;
  long numofthislevel = 1;

  long test[100];
  long sum_modp;
  long bestlevel = -1;
  long besttest = -1;
  long bestnum = -1;
  long maxtest = -1;
  long bestindex = -1;

  for( long thismodform = 0; thismodform < nummodforms; thismodform++ )
    {
      if( level[thismodform] == oldlevel )
	numofthislevel++;
      else
	{
	  numofthislevel = 1;
	  oldlevel = level[thismodform];
	}

      for( long i = 0; i < 100; i++ )
	test[i] = 0;

      maxtest = -1;

      // compare for all 25 primes
      for( long i = 0; i < 25; i++ )
	{
	  p = primes[i];
	  sum_modp = ( numpoints[i] + coeff[thismodform][i] ) % p;
	  while( sum_modp < 0 )
	    sum_modp += p;
	  test[sum_modp]++;

	  // With this extra check, the program will not determine the actual
	  // newform associated with the double octic, but the twist of minimal level.
	  // If you do not like this just comment out these 7 lines.
	  if( coeff[thismodform][i] % p != 0 && p > 2 )
	    {
	      sum_modp = ( numpoints[i] - coeff[thismodform][i] ) % p;
	      while( sum_modp < 0 )
		sum_modp += p;
	      test[sum_modp]++;
	    }
	}

      for( long i = 0; i < 100; i++ )
	if( maxtest == -1 || test[i] > maxtest )
	  maxtest = test[i];

      if( maxtest > besttest )
	{
	  bestlevel = level[thismodform];
	  besttest = maxtest;
	  bestnum = numofthislevel;
	  bestindex = thismodform;
	}
    }

  //--------------------------------------------------------------
  // print the most likely modular form
  // CAUTION: the coefficients should fit for all good primes
  //--------------------------------------------------------------

  cout << "The most likely newform seems to be " << bestlevel << "/" << bestnum << "." << endl;
  cout << "It fits for " << besttest << " of 25 primes." << endl;
}