// Jonthan D. Hauenstein
// Copyright August 2011
// Developed for
// "Certifying solutions to square systems of polynomial-exponential equations"
// by J.D. Hauenstein and V. Levandovskyy

#include <stdio.h>
#include <time.h>

void increment(int *v)
{ // increment the vector
  int i, max[6];
  max[0] = max[2] = max[4] = 5;
  max[1] = max[3] = max[5] = 4;

  for (i = 0; i < 6; i++)
    if (v[i] < max[i])
    {
      v[i]++;
      break;
    }
    else
      v[i] = 1;

  return;
}

int checkVector(int *v)
{ // check the vector for consistency
  return !((v[0] > 1 && v[1] > 1) || (v[2] > 1 && v[3] > 1) || (v[4] > 1 && v[5] > 1));
}

int main()
{ // perform coefficient-parameter continuation
  int i, v[6], count = 0, run = 0, total = 7999, solns = 0;
  char ch;
  FILE *F = NULL;
  FILE *SOLNS = fopen("final_solutions", "w+");

  // initialize v
  v[0] = 2; 
  v[1] = v[2] = v[3] = v[4] = v[5] = 1;

  // setup the beginning of SOLNS
  fprintf(SOLNS, "              \n\n");
  F = fopen("start", "r");
  fscanf(F, "%d\n\n", &solns);
  while ((ch = fgetc(F)) != EOF)
    fprintf(SOLNS, "%c", ch);

  while (count < total)
  { // check the vector
    if (checkVector(v))
    { // setup the homotopy
      run++;
      F = fopen("input_middle", "w");

      fprintf(F, "L1 = L1%d;\n", v[0]);
      fprintf(F, "L2 = L2%d;\n", v[1]);
      fprintf(F, "L3 = L3%d;\n", v[2]);
      fprintf(F, "L4 = L4%d;\n", v[3]);
      fprintf(F, "L5 = L5%d;\n", v[4]);
      fprintf(F, "L6 = L6%d;\n", v[5]);

      fclose(F);

      // create the input file
      system("cat input_top input_middle input_bottom > input");

      // run Bertini
      system("./bertini input > /dev/null");

      // verify nonsingular_solutions
      F = fopen("nonsingular_solutions", "r");
      fscanf(F, "%d\n\n", &i);
      solns += i;
      // copy to SOLNS
      while ((ch = fgetc(F)) != EOF)
        fprintf(SOLNS, "%c", ch);

      // close F
      fclose(F);
    }

    // increment
    increment(v);
    count++;
  }

  // write the number of solutions
  rewind(SOLNS);
  fprintf(SOLNS, "%d", solns);

  // close SOLNS
  fclose(SOLNS);

  return 0;
}