#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <gmpxx.h>

// This program will test the representability of
// a quaternary quadratic form and output the
// results to a .log file.
// It takes a single argument, which is
// the output file of a Magma script
// which does all necessary preprocessing.

// This script uses the output of prep6.txt
// to check up to the bound computed. It
// generated squarefree integers less than the bound
// and checks them, 500000 at a time. 

// This program uses the C++ datatype bool
// but is otherwise essentially written in C.
// It relies on the GMP library for large integer arithmetic
// To compile it, you should run something like:
// g++ represent2.c -lgmpxx

// Function prototypes

int KroneckerSymbol(long a, long p);
double F4Prime(long p, long charnum, long N, unsigned short anisos[10]);
void bitset(unsigned long index, unsigned long terntheta[]);
bool bitread(unsigned long index, unsigned long terntheta[]);

// The Legendre symbol (a/p). It uses Euler's criterion
// in a fairly naive way, but it is efficient enough.

int KroneckerSymbol(long a, long p) {
  if (p==2) {
    if (a % 2 == 0) {
      return 0;
    }
    else {
      if (abs(a) % 8 == 1)
          return 1;
      else
      return -1;
    }
  }
  else {
    long a1 = a % p;
    if (a1 < 0) {
      a1 = a1 + p;
    }
    long b = 1;
    long i;
    for(i = 1; i<=((p-1)/2); i++)
      {
    b = (b * a1) % p;
      }
    if(b==(p-1)) 
      return -1;
    else
      return b;
  }
}

// A bound on the F_4 function evaluated at a prime.

double F4Prime(long p, long charnum, long N, unsigned short anisos[10])
{
  double F4;
  double pp;
  bool done;
  pp = p;
  if(pp < 10000) {
    // Compute exact
    F4 = sqrt(pp)/2;
    if(KroneckerSymbol(charnum,p) == -1) {
      if(N % p != 0) {
        F4 = F4*(pp-1)/(pp+1);
      }
    }
    // Check to see if it is anisotropic
    for(int i=0; i<10; i++)
      {
    done = false;
        if(done==false) {
          if(anisos[i]==0) { done = true; }
          if(anisos[i]>p) { done = true; }
          if(anisos[i]==p) {
            F4 = F4*(1/pp);
      }
    }      
      }
  }
  if(pp > 10000) {
    F4 = (sqrt(pp)/2)*(1 - 1/pp);
  }
  return F4;
}

// This sets a bit in an array of long integer. This is used
// in the ternary theta function computation.

void bitset(unsigned long index, unsigned long terntheta[])
{
  unsigned short whichbit;
  unsigned long whichlong;
  unsigned long orlong;
  
  whichlong = index/32;
  whichbit = index % 32;
  orlong = 1;
  orlong = orlong << whichbit;
  terntheta[whichlong] = terntheta[whichlong] | orlong;
}

// This reads the value of a bit in an array of longs.

bool bitread(unsigned long index, unsigned long terntheta[])
{
  unsigned short whichbit;
  unsigned long whichlong, orlong;
  bool ret;
  
  whichlong = index/32;
  whichbit = index % 32;
  orlong = 1;
  orlong = orlong << whichbit;
  if((terntheta[whichlong] & orlong) != 0) {
    ret = true;
  }
  else
  { ret = false;
  }
  return ret;
}

int main(int argc, char *argv[])
{
  // Parse input
  // 

  FILE *in, *out;
  char form_num[80], level[80], charac[80], aniso_list[80];
  char F4_bound[80], d_str[80];
  char split[80];
  char curnum[10], tempst[5];
  char temp;
  unsigned long formnum, levelnum, charnum, d;
  signed short matrix[10];
  unsigned short anisos[10];
  unsigned short iter, iter2, maxaniso;
  unsigned short tempnum;
  double F4bound, dtemp, maxprime;
  // Eligible prime bound
  double epb; 
  // Number of primes in the prime file.
  const long PRIMESIZE = 21500000; 

  // Initizialize squarefree up front
  mpz_t squarefree[500000];

  mpz_array_init(squarefree[0],500000,80);  

  mpz_t exceptions[10000];
  long curexcep;
  mpz_array_init(exceptions[0], 9999, 80);
  
  //printf("Number of bytes in a long is %d.\n",sizeof(curexcep));
  if(argc != 2) {
    //printf("Program for determining which numbers are\n");
    //printf("represented by a quaternary quadratic form.\n");
    //printf("Usage: represent prepfile\n");
    //printf("Here, prepfile is the output of the MAGMA script\n");
    //printf("which does the necessary pre-processing.\n");
    exit(1);
  }
  
  if((in=fopen(argv[1], "r"))==NULL) {
    //printf("Cannot open file!.\n");
    exit(1);
  }
  printf("Input file successfully opened.\n");

  // Read the number of the form
  fgets(form_num, 79, in);
  // Read the level
  fgets(level, 79, in);
  // Read the character
  fgets(charac, 79, in);
  // Read the list of anisotropic primes
  fgets(aniso_list, 79, in);
  // Read the F4 bound
  fgets(F4_bound, 79, in);
  // Read the value of d
  fgets(d_str, 79, in);
  // Read the ternary lattice
  fgets(split, 79, in);

  printf("Data read.\n");

  //Initizialize anisos

  for(iter = 0; iter < 10; iter++)
  {
    anisos[iter] = 0;
    matrix[iter] = 32767;
  }

  // Parse the string form_num

  formnum = atoi(form_num);

  // Parse the character

  charnum = atoi(charac);

  // Parse the level

  levelnum = atoi(level);

  // Parse the list of anisotropic primes

  curnum[0] = '\0';
  for(iter = 1; iter < strlen(aniso_list)-2; iter++)
  {
    temp = aniso_list[iter];
    tempst[0] = temp;
    tempst[1] = '\0';
    if(temp!=',') {
      strcat(curnum,tempst);
    }

    if(temp==',') {
      tempnum = atoi(curnum);
      iter2 = 0;
      while(anisos[iter2]!=0)
      {
        iter2++;
      }
      anisos[iter2] = tempnum;
      curnum[0] = '\0';
    }
  }
  tempnum = atoi(curnum);
  iter2 = 0;
  while(anisos[iter2]!=0)
  {
    iter2++;
  }
  anisos[iter2] = tempnum;

  // Parse the string F4_bound
  
  F4bound = atof(F4_bound);

  // Parse the string d_str

  d = atoi(d_str);

  // Parse the ternary matrix

  curnum[0] = '\0';
  for(iter = 1; iter < strlen(split)-2; iter++)
  {
    temp = split[iter];
    tempst[0] = split[iter];
    tempst[1] = '\0';
    if(temp!=',') {
      strcat(curnum,tempst);      
    }

    if(temp==',') {
      tempnum = atoi(curnum);
      iter2 = 0;
      while(matrix[iter2]!=32767)
      {
        iter2++;
      }
      matrix[iter2] = tempnum;
      curnum[0] = '\0';
    }
  }  
  tempnum = atoi(curnum);
  iter2 = 0;
  while(matrix[iter2]!=32767)
  {
    iter2++;
  }
  matrix[iter2] = tempnum;
  
  printf("The form number is %lu.\n",formnum);
  printf("The level is %lu.\n",levelnum);
  printf("The character is %lu.\n",charnum);
  printf("The list of anisotropic primes is [");
  iter2 = 0;
  while(anisos[iter2]!=0)
  {
    iter2++;
  }
  if(iter2==0) { maxaniso = 0; }
  if(iter2>0) { maxaniso = anisos[iter2-1]; }
  if(maxaniso > 9999) {
    printf("Error! Large anisotropic prime!\n");
    exit(1);
  }
  
  for(iter = 0; iter < iter2; iter++)
  {
    printf("%u ",anisos[iter]);
  }
  printf("].\n");
    
  printf("The F4 bound is %f.\n",F4bound);
  printf("The value of d is %lu.\n",d);
  printf("The matrix of the split local cover is [");
  for(iter = 0; iter < 9; iter++)
  {
    printf("%d ",matrix[iter]);
  }
  printf("].\n");
  

  // Parse the precomputed list of primes

  FILE *prime;
  printf("Opening file for primes.\n");
  if((prime=fopen("primelist.txt", "r"))==NULL) {
    printf("Cannot open prime file primelist.txt!.\n");
    exit(1);
  }

  long *primelist;
  long primecount, curprime;
  printf("Allocating memory for prime list.\n");
  if((primelist = (long *)malloc(sizeof(curprime)*PRIMESIZE))==NULL) {
    printf("Cannot allocate space for list of primes!.\n");
    exit(1);
  }
  char ch;

  // Read the first PRIMESIZE primes.

  printf("Parsing prime file.\n");
 
  primecount = 0;
  curnum[0] = '\0';
  bool done;
  done = false;  
  do {
    ch = getc(prime);
    if(ch!='[' && ch!=' ' && ch!=',') {
      tempst[0] = ch;
      tempst[1] = '\0';
      strcat(curnum,tempst);
    }
    if(ch==',') {
      curprime = atol(curnum);
      primelist[primecount] = curprime;
      primecount = primecount + 1;
      curnum[0] = '\0';      
    }
    if(ch==']') {
      curprime = atol(curnum);
      primelist[primecount] = curprime;
      primecount = primecount + 1;
      done = true;
    }       
  } while(done==false);
  fclose(prime);
  printf("%lu primes parsed.\n",primecount);

  // Compute the largest possible eligible prime

  epb = F4bound;

  dtemp = F4Prime(2,charnum,levelnum,anisos);
  if(dtemp < 1) {
    epb = epb/dtemp;
  }
  dtemp = F4Prime(5,charnum,levelnum,anisos);
  if(dtemp < 1) {
    epb = epb/dtemp;
  }
  dtemp = F4Prime(7,charnum,levelnum,anisos);
  if(dtemp < 1) {
    epb = epb/dtemp;
  }  
  iter2 = 0;
  while(anisos[iter2]!=0) { iter2++; }
  for(iter = 0; iter < iter2; iter++)
    {
      dtemp = F4Prime(anisos[iter],charnum,levelnum,anisos);
      if(dtemp < 1 && anisos[iter] > 7) {
        epb = epb/dtemp;
      }
    }  
  
  maxprime = 2*epb*epb + 2*epb*sqrt(epb*epb+1) + 1;
  printf("We need primes up to %f.\n",maxprime);
  if(maxprime < maxaniso) { maxprime = maxaniso; }
  if(maxprime > 400000000) {
    printf("Not enough precomputed primes!\n");
    exit(1);
  }
  printf("F4 prime bound = %f.\n",epb);
  printf("Max prime = %f.\n",maxprime);

  // Create list of eligible primes

  // Allocate space for list of eligible primes and F4 bounds

  long *eligiblep;
  float *F4list;
  float F4ret;
  long approxnump;

  if(maxprime > 17) {
    approxnump = (long)floor(maxprime/(log(maxprime)-2));
  }
  if(maxprime <= 17) {
    approxnump = 7;
  }
  
  if((eligiblep = (long *)malloc(sizeof(curprime)*approxnump))==NULL) {
    printf("Cannot allocate space for list of primes!.\n");
    exit(1);
  }
  if((F4list = (float *)malloc(sizeof(F4ret)*approxnump))==NULL) {
    printf("Cannot allocate space for list of F4 numbers!\n");
    exit(1);
  }

  // Create list of primes and F4 bounds

  long count, primeindex;

  done = false;
  primeindex = 0;
  count = 0;
  do {
    curprime = primelist[primeindex];
    if(curprime > maxprime) {
      done = true;
    }
    else
      {
        primeindex++;
        F4ret = F4Prime(curprime,charnum,levelnum,anisos);
        if (F4ret < epb) {
          //printf("Prime %lu is eligible with F4 bound %f.\n",curprime,F4ret);
      eligiblep[count] = curprime;
          F4list[count] = F4ret;
          count++;
        }
      }

  } while (done==false);

  // Sort first entries of list

  bool reallydone;

  reallydone = false;
  printf("Done computing eligible primes.\n");
  printf("There were %lu.\n",count);
  if(count==0) {
    // Nothing to do!
    reallydone = true;
  }

  long maxsort;
  maxsort = 0;
  if(reallydone==false) {
    printf("The last is %lu.\n",eligiblep[count-1]);
    if(eligiblep[count-1] < 10000) {
      maxsort = count-1;
    }
    else
    {
      maxsort = 1227;
    }
  }

  //printf("Maxsort is %lu.\n",maxsort);

  // Output number of eligible primes

  FILE *outfile;
  char outfilename[80];
  char suffix[10] = ".log";
  char writedata[1000];  

  sprintf(outfilename, "%lu",formnum);
  strcat(outfilename,suffix);

  if((outfile=fopen(outfilename, "a"))==NULL) {
    printf("Cannot open log file %s!\n",outfilename);
    exit(1);
  }  
  printf("Output %s file opened.\n",outfilename);

  sprintf(writedata,"Verification\n");
  fputs(writedata,outfile);
  sprintf(writedata,"Number of eligible primes = %lu.\n",primeindex);
  fputs(writedata,outfile);
  sprintf(writedata,"F4 prime bound = %f.\n",epb);
  fputs(writedata,outfile);


  if(reallydone==true) {
    sprintf(writedata,"No verification necessary.\n");
    fputs(writedata,outfile);
    sprintf(writedata,"No new exceptions found.\n");
    fputs(writedata,outfile);

    double timing;
    timing = clock();
    timing = timing/CLOCKS_PER_SEC;
    sprintf(writedata,"Time required for verification: %f sec.\n",timing);
    fputs(writedata,outfile);
    fclose(outfile);
    exit(0);  
  }
  
  // bubble sort first maxsort entries of eligiblep and F4list

  printf("Bubble sorting.\n");
  
  bool changed;
  done = false;
  do {
    changed = false;
  for(iter = 0; iter < maxsort; iter++)
    {
      if(F4list[iter+1] < F4list[iter]) 
    {
          curprime = eligiblep[iter];
          F4ret = F4list[iter];
          eligiblep[iter] = eligiblep[iter+1];
          F4list[iter] = F4list[iter+1];
          eligiblep[iter+1] = curprime;
          F4list[iter+1] = F4ret;
          changed = true;
    }
    }
  if(changed==false) { done = true; }
  } while(done==false);

  printf("Done bubble sorting.\n");

  // Determine the maximum number of prime factors

  short maxfac;

  maxfac = 0;
  done = false;
  iter = 0;
  F4ret = F4list[0];
  do {
    //printf("Current max number of prime factors = %lu.\n",iter);
    //printf("F4 bound = %f.\n",F4ret);
    iter = iter + 1;
    if(F4ret > F4bound) {
      done = true;
    }

    if(iter>count) {
      done = true;
      iter = count + 3; 
    }
    if(iter==count) {
      done = true;
      iter = count + 2;
    }

    if(done==false) { F4ret = F4ret * F4list[iter]; }
  } while(done==false);  
  maxfac = iter-2;
  if(F4bound==0) { maxfac = 0; }

  printf("Determined maximum number of prime factors: %u\n",maxfac);
  sprintf(writedata,"Eligible numbers have at most %u prime factors.\n",
  maxfac);
  fputs(writedata,outfile);

  // Determine an upper bound for the largest eligible number

  double adjustment, enumbound;
  short maxprimecheck;
  maxprimecheck = 1227;
  if(maxprimecheck > primeindex) {
    maxprimecheck = primeindex;
  }
  adjustment = 1;
  
  for(iter = 0; iter < maxprimecheck; iter++)
  {
    curprime = primelist[iter];
    if(levelnum % curprime != 0 && KroneckerSymbol(charnum,curprime) == -1)
      {
    adjustment = adjustment*(double(curprime-1)/double(curprime+1));
      }
  }
  iter2 = 0;
  while(anisos[iter2]!=0)
  {
    adjustment = adjustment*(1/double(anisos[iter2]));
    iter2++;
  }

  enumbound = pow(2.0, (double)maxfac) * F4bound / adjustment;
  enumbound = enumbound*enumbound;

  printf("Determined upper bound for largest eligible number:%f\n",enumbound);
  sprintf(writedata,"Largest eligible number is <= %f.\n",enumbound);
  fputs(writedata,outfile);

  // Determine ternary theta precision

  // Hanke's theta precision function takes the single argument 
  // d. The ternary theta precision is selected to be
  // 2*sqrt(Upper bound for eligible numbers)*d*C
  // where C is some constant. Hanke sometimes takes C = 5 (or C = 9).

  // The ranges used in computing are the following.
  // The z max is sqrt(D* discrim of (x,y) subform / ternary discrim)
  // D = 2*theta precision + 1
  // x and y range from 0 to E, E = 500 (sometimes 800)

  long thetaprec, thetaprec2, thetaprec3;
  double xysubformdisc, ternarydisc;
  thetaprec = (long)ceil(2*sqrt(enumbound)*d*16); // Sometimes last const = 5
  
  xysubformdisc = matrix[0]*matrix[4] - matrix[1]*matrix[3];
  ternarydisc = matrix[0]*matrix[4]*matrix[8] + matrix[1]*matrix[5]*matrix[6];
  ternarydisc = ternarydisc + matrix[2]*matrix[3]*matrix[7];
  ternarydisc = ternarydisc - matrix[6]*matrix[4]*matrix[2];
  ternarydisc = ternarydisc - matrix[7]*matrix[5]*matrix[0];
  ternarydisc = ternarydisc - matrix[8]*matrix[3]*matrix[1];

  thetaprec2 = 
  (long)ceil(sqrt(2*thetaprec+1)*sqrt(xysubformdisc/ternarydisc));
  thetaprec3 = 800; // Sometimes this is 800

  sprintf(writedata,
  "Necessary precision of boolean theta function = %lu.\n",thetaprec);
  fputs(writedata,outfile);
  sprintf(writedata,
  "Maximum z value to be used for boolean theta function = %lu.\n",
  thetaprec2);
  fputs(writedata,outfile);
  sprintf(writedata,
  "Maximum x and y values to be used for boolean theta function = %lu.\n"
  ,thetaprec3);
  fputs(writedata,outfile);

  printf("Determined theta precision.\n");

  // Compute approximate boolean ternary theta function

  unsigned long *terntheta;
  unsigned long thetaindex;
  unsigned long thetasize;
  
  printf("Initializing ternary theta function memory.\n");
  
  thetasize = (long)ceil(thetaprec/32); 
  printf("Need %lu bytes for theta function.\n",thetasize*4);
  if((terntheta = (unsigned long *)malloc(sizeof(curprime)*thetasize))==NULL) {
    printf("Cannot allocate space for ternary theta function!.\n");
    exit(1);
  }
  for(thetaindex = 0; thetaindex < thetasize; thetaindex++)
  {
    terntheta[thetaindex] = 0;
  }
  signed long x, y, z;
  unsigned long q11, q12, q13, q22, q23, q33;
  unsigned long zterm, yzterm, val;
  q11 = matrix[0]/2;
  q12 = matrix[1];
  q13 = matrix[2];
  q22 = matrix[4]/2;
  q23 = matrix[5];
  q33 = matrix[8]/2;

  printf("Computing ternary theta function.\n");
  printf("Maximum z precision is %lu.\n",thetaprec2);

  for(z = thetaprec2; z >= 0; z--)
  {
    if(z%100==0) printf("Working on z = %ld.\n",z);
    zterm = q33*z*z;
    for(y = thetaprec3; y >= -thetaprec3; y--)
    {
      yzterm = zterm + q23*y*z + q22*y*y;
      for(x = thetaprec3; x >= -thetaprec3; x--)
      {
         val = yzterm + q11*x*x + q12*x*y + q13*x*z;
         if(val < thetaprec) {
           bitset(val,terntheta);
         }
      }
    }
  }
  printf("Approximate boolean theta function computation complete.\n");

  // Verify

  mpz_t temp2;
  // This fixes a runtime error with realloc
  // and allows temp2 to be as large as 2^200.
  mpz_init2(temp2,200);
  
  iter = maxfac;
  long ptrarray[maxfac];
  int lastptr, carryptr;
  double tempF4, checknumcount;
  long sqrfreecount, exceptioncount;
  long sqrfreeindex, depth, difflong;
  signed long w, startw;
  bool carry, repn, badflow;
  sqrfreecount = 0;
  checknumcount = 0;
  exceptioncount = 0;
  // We know that our forms represent 1, so we only have
  // to check numbers with a positive number of prime factors
  for(iter = maxfac; iter >= 1; iter--)
  {
    printf("Working on numbers with %d prime factors.\n",iter);
    lastptr = iter - 1;
    carryptr = iter - 1;
    for(iter2 = 0; iter2 <= lastptr; iter2++)
    {
      ptrarray[iter2] = iter2;
    }
    done = false;
    do {
      carry = false;
      if(ptrarray[lastptr] < count - 1)
      {
        tempF4 = 1;
        for(iter2 = 0; iter2 <= lastptr; iter2++)
        {
          tempF4 = tempF4*F4list[ptrarray[iter2]];
        }
        if(tempF4 < F4bound) {
          mpz_set_ui(temp2,1);
          for(iter2 = 0; iter2 <= lastptr; iter2++)
          {
            mpz_mul_ui(temp2,temp2,eligiblep[ptrarray[iter2]]);
          }
          mpz_set(squarefree[sqrfreecount],temp2);
          sqrfreecount = sqrfreecount + 1;
          ptrarray[lastptr] = ptrarray[lastptr] + 1;
          carryptr = lastptr;
        }
        else { carry = true; }
      }
      else {carry = true; }
      
      if(carry == true)
      {
        if(carryptr > 0) {
          carryptr = carryptr - 1;
          ptrarray[carryptr] = ptrarray[carryptr] + 1;
          for(iter2 = carryptr + 1; iter2 <= lastptr; iter2++)
          {
            ptrarray[iter2] = ptrarray[carryptr] + (iter2 - carryptr);
          }
        }
        else
        { done = true; }
      }
      if(sqrfreecount == 500000 || done==true) {
        //printf("Checking %lu numbers.\n",sqrfreecount);
        // Check numbers in squarefree
        checknumcount = checknumcount + sqrfreecount;
        for(sqrfreeindex = 0; sqrfreeindex < sqrfreecount; sqrfreeindex++)
        {
          repn = false;
          badflow = false;
          depth = 0;
          w = (signed long)floor(sqrt(mpz_get_d(squarefree[sqrfreeindex])/d));
          startw = w;
          do {
            if(w < 0) {
              w = 0;
              badflow = true;
            }
            mpz_set_ui(temp2, d);
            mpz_mul_ui(temp2, temp2, w);
            mpz_mul_ui(temp2, temp2, w);
            mpz_sub(temp2, squarefree[sqrfreeindex], temp2);
            difflong = mpz_get_ui(temp2);
            if(mpz_cmp_ui(temp2, thetaprec) > 0)
            {
              badflow = true;
              difflong = 0;
            }
            if(bitread(difflong,terntheta) == true)
            {
              repn = true;
            }
            else
            { 
              w--; 
            }
          } while((repn == false) && (badflow == false));
          // if badflow = true - add to exceptions
          if(badflow == true)
          {
            mpz_set(exceptions[exceptioncount],squarefree[sqrfreeindex]);
            exceptioncount++;
          }
        }
        sqrfreecount = 0;
      }
      
    } while(done==false);
  }

  sprintf(writedata, "A total of %f numbers were checked.\n",
  checknumcount);
  fputs(writedata, outfile);
  sprintf(writedata, "There were a total of %lu exceptions.\n",
  exceptioncount);
  fputs(writedata, outfile);
  sprintf(writedata, "They are the following:");
  fputs(writedata, outfile);

  bool large_excep;
  large_excep = false;
  unsigned long long int reallyLarge = 18446744073709551615; 
  //printf("Starting output of exceptions.\n");
  for(iter = 0; iter < exceptioncount; iter++)
  {
    //printf("Checking size of exception %u.\n",iter);
    //printf("Copying value.\n");
    if(mpz_cmp_ui(exceptions[iter],reallyLarge) > 0)
      {
	printf("One of the exceptions is larger than 2^64-1!\n");
        sprintf(writedata, "Really large exception found!\n");
	fputs(writedata, outfile);
        printf("Outputting exception %lu\n",mpz_get_ui(exceptions[iter])); 
      }
    curexcep = mpz_get_ui(exceptions[iter]);
    //printf("Done!\n");
    if(curexcep > 5000)
      {
	large_excep = true;
      }
    //printf("Outputing exception %u.\n",iter);
    //gmp_printf("It is %Zd.\n",curexcep);
    if(iter == 0) 
    {
      sprintf(writedata, "%lu",curexcep);
    }
    if(iter > 0)
    {
      sprintf(writedata, ",%lu",curexcep);
    }
    //printf("Writing to file.\n");
    fputs(writedata, outfile);
  }
  sprintf(writedata,"\n");
  fputs(writedata,outfile);

  if(large_excep==false)
    {
      sprintf(writedata,"No new exceptions found.\n");
      fputs(writedata,outfile);
    }

  if(large_excep==true)
    {
      sprintf(writedata,"Some new large exceptions found!\n");
      fputs(writedata,outfile);
    }

  // Write time required  


  double timing;
  timing = clock();
  timing = timing/CLOCKS_PER_SEC;
  sprintf(writedata,"Time required for verification: %f sec.\n",timing);
  fputs(writedata,outfile);
  fclose(outfile);
  
}