/// Read in inequivalent anomaly charges and test them for desired properties
#include <iostream>
#include <sstream>
#include <string>
#include <cstring>
#include <algorithm>
#include <vector>
#include "searchU1.h"

using namespace std;

int main() {

  header();
  int numberOfSolutions   = 0;
  int numberOfInteresting = 0;
  int paritySymmetric     = 0;
  int fnuPconditions      = 0;
  
  string line;
  while (getline(cin,line)) {
    istringstream input(line);
    string word1;
    input >> word1;

    /// The line doesn't start with a comment
    if (word1.find("#") == string::npos) {
      numberOfSolutions++;
      /// 3 families of SM fermions, defined in increasing order of charge
      solN a;
      inputsolN(line, a);

      /// Perform tests on a
      /// Firstly, does the solution have at least one non-zero charge for a Q?
      bool nonZeroQcharge = false;
      if (a.Q[0] != 0 || a.Q[2] != 0) nonZeroQcharge = true;
      /// Does the solution have only zero charges for d's
      bool zeroDcharges = false;
      if (a.d[0] == 0 && a.d[1] == 0 && a.d[2] == 0) zeroDcharges = true;
      /// Does the solution have at least one non-zero charge for an L?
      bool nonZeroLcharge = false;
      if (a.L[0] != 0 || a.L[2] != 0) nonZeroLcharge = true;
      /// We want at least one zero L charge and one e charge
      bool niceEcharge = false;
      if (a.L[1] == 0 && a.e[1] == 0) niceEcharge = true;

      /// These are to check Ellis/Fairbairn etc's assertions:
      bool q1Equalsq2 = false;
      if (a.Q[0] == a.Q[1] || a.Q[1] == a.Q[2]) q1Equalsq2 = true;

      solN neg = -a;
      if (neg == a) paritySymmetric++;
      
      /// Our requirements for an `interesting' solution
      if (nonZeroQcharge && zeroDcharges && q1Equalsq2 && niceEcharge) 	{
	outputsolN(a);
	numberOfInteresting++;
      }
      
      /// Check our Eq 3.26 applies
      int fnuPlus = a.n[0] + a.n[1] + a.n[2];
      if (fnuPlus == 0) {
	int fnuBar = a.n[2] + a.n[1] - 2 * a.n[0];
	if (fnuBar % 3 != 0) {
	  cout << "Found exception to 3.26: \n";
	  outputsolN(a);
	}
	int fuBar = a.u[2] + a.u[1] - 2 * a.u[0];
	int fQBar = a.Q[2] + a.Q[1] - 2 * a.Q[0];
	int fdBar = a.d[2] + a.d[1] - 2 * a.d[0];
	int feBar = a.e[2] + a.e[1] - 2 * a.e[0];
	int fLBar = a.L[2] + a.L[1] - 2 * a.L[0];

	int ru = 0;
	if (fuBar % 3 == 1) ru = 1;
	if (fuBar % 3 == 2) ru = -1;

	if ( (fuBar  % 3 == 0 && fQBar   % 3 == 0 && fdBar   % 3 == 0
	      && feBar % 3 == 0 && fLBar % 3 == 0) ||
	     ((fuBar-ru) % 3 == 0 && (fQBar-ru) % 3 == 0 && (fdBar-ru) % 3 == 0
	      && (feBar % 3) == 0 && fLBar % 3 == 0) ||
	     ((fuBar-ru) % 3 == 0 && (fQBar-ru) % 3 == 0 && (fdBar+ru) % 3 == 0
	      && (feBar+ru) % 3 == 0 && (fLBar-ru) % 3 == 0) ||
	     ((fuBar-ru) % 3 == 0 && (fQBar+ru) % 3 == 0 && (fdBar-ru) % 3 == 0
	      && (feBar+ru) % 3 == 0 && (fLBar-ru) % 3 == 0) ||
	     ((fuBar-ru) % 3 == 0 && (fQBar+ru) % 3 == 0 && (fdBar+ru) % 3 == 0
	      && (feBar-ru) % 3 == 0 && (fLBar+ru) % 3 == 0) ) fnuPconditions++;
	else {
	  cout << "Found Exception to 3.26: \n";
	  outputsolN(a);
	}
      }
    }
  }
  cout << "# Total number of solutions=" << numberOfSolutions << endl;
  cout << "# Parity symmetric=" << paritySymmetric << endl;
  cout << "# number of interesting solutions=" << numberOfInteresting << endl;
  cout << "# number satisfying (3.26) in correct conditions=" << fnuPconditions << endl;
}