/******************************************************************
 ******************************************************************
 ***                                                            ***
 *** Computation of the Seifert Matrix and Alexander Polynomial ***
 *** of a knot or link with a given pretzel representation.     ***
 ***                                                            ***
 *** Algorithm by Julia Collins, Edinburgh, 2007                ***
 *** Ported to C++ by Thomas Koeppe                             ***
 *** Bugfix by Lukas Lewark                                     ***
 *** Version: 2016-11-10                                        ***
 ***                                                            ***
 ***                                                            ***
 *** Usage: The program reads a braid or pretzel representation ***
 ***        from the standard input, which is either a space-   ***
 ***        separated list of integers or a sequence of letters ***
 ***        and integers.                                       ***
 ***                                                            ***
 ***        It computes the strand permutations in the braid,   ***
 ***        the adjacency of homology cycles, the Seifert       ***
 ***        matrix and the Alexander polynomial.                ***
 ***                                                            ***
 ***  Requirements: GMP library required. The attached          ***
 ***        "matrix.h" provides the necessary matrix routines.  ***
 ***        Compilation with GCC for example via                ***
 ***       g++ -o pretzeljulia -O2 pretzel-v2.cpp -lgmpxx -lgmp ***
 ***                                                            ***
 ***                                                            ***
 ******************************************************************
 ******************************************************************/



#include <vector>
#include <string>
#include <sstream>
#include <gmpxx.h>
#include "matrix-v2.h"


/** The type "Pretzel" is a vector holding the pretzel information,
 ** which consists of pairs <unsigned int, int> of strand number and twisting.
 ** A consistent pretzel has only pairs <positive int, odd int>.
 ** The type "IntVector" is just a list of signed integers.
 **/
typedef std::pair<unsigned int, int> PretzelItem;
typedef std::vector<PretzelItem>     Pretzel;
typedef std::vector<int>             IntVector;

/** We also implement a convenient stream output for vectors and pairs.
 **/
template <typename T, typename U>
std::ostream & operator<<(std::ostream & os, const std::pair<T, U> & p)
{
  return (os << "(" << p.first << "," << p.second << ")");
}
template <typename T>
std::ostream & operator<<(std::ostream & os, const std::vector<T> & v)
{
  if (!v.size()) { os << "[]"; return os; }
  os << "[";
  for (size_t i = 0; i < v.size() - 1; i++)
    os << v[i] << ", ";
  os << v[v.size() - 1] << "]";
  return os;
}



/** The signum function for integers.
 **/
int sign(const int n)
{
  if (n < 0) return -1;
  else if (n == 0) return 0;
  else return 1;
}


/** This auxilliary routine reads the braid data from the
 ** standard input (the console).
 **/
void readInput(Pretzel & pretzel)
{
  int tmp;
  std::string instr, tmpstr;
  std::cerr << "Please enter a braid: ";
  std::getline(std::cin, instr);
  std::cerr << std::endl;
  
  if (instr.length() && ((int(instr[0]) >= 48 && int(instr[0]) <= 57) || instr[0] == '-'))
  {
    std::istringstream is(instr);
    while (is >> tmp) 
      pretzel.push_back(PretzelItem(abs(tmp), sign(tmp)));
    std::cerr << "Interpreting your input \"" << instr << "\" as numerical braid notation as " << pretzel << std::endl;
    return;
  }


  bool reading_number = false;
  bool found_letter = false;
  bool item_ready = false;
  size_t index = 0;
  PretzelItem tmpitem;
  int thesign = 1;
  for (std::string::iterator i = instr.begin(); i != instr.end(); )
  {
    // Next character is a letter
    if ( (int(*i) >=65 && int(*i) <= 90) || (int(*i) >= 97 && int(*i) <= 122) )
    {
      if (item_ready)
      {
        if (reading_number)
        {
          if ( atoi(tmpstr.c_str()) % 2 != 1 && atoi(tmpstr.c_str()) % 2 != -1 )
            throw "Error while parsing number!";

          tmpitem.second = atoi(tmpstr.c_str()) * thesign;
          reading_number = false;
        }
        else
          tmpitem.second = thesign;

        pretzel.push_back(tmpitem);
        thesign = 1;
        item_ready = false;
      }

      if (int(*i) < 97) tmpitem.first = int(*i) - 64;
      else            { tmpitem.first = int(*i) - 96; thesign = -1; }
      tmpstr.clear();
      item_ready = true;
    }

    // Next character is a number or "-"
    else if ( (int(*i) >= 48 && int(*i) <= 57) || *i == '-' )
    {
      tmpstr += *i;
      reading_number = true;
    }

    else if (*i != ' ') std::cerr << "Skipping unrecognized character '" << *i << "'." << std::endl;

    i++;
  }

  // Finally, after no more input is available, add the last item
  if (item_ready)
  {
    if (reading_number)
    {
      if ( atoi(tmpstr.c_str()) % 2 != 1 && atoi(tmpstr.c_str()) % 2 != -1 )
        throw "Error while parsing number!";
      tmpitem.second = atoi(tmpstr.c_str()) * thesign;
    }
    else
      tmpitem.second = thesign;
    pretzel.push_back(tmpitem);
  }
  std::cerr << "Interpreting your input \"" << instr << "\" as alpha-numeric pretzel notation as " << pretzel << std::endl;
}


/** Given a braid, this function computes its strands
 ** and their permutations, i.e. the size of the result
 ** equals the number of strands of the given braid.
 **/
void computeStrandPermutations(const Pretzel & pretzel, IntVector & strand_permutation, unsigned int & MaximalElement)
  /* The algorithm follows each strand in turn through the braid to see
   * where its final position is.  For example, if we are following strand
   * 2 and we see a crossing labelled '1', then we know that strand 2 will
   * switch places with strand 1. If the crossing is labelled '2' then strand
   * 2 will switch places with strand 3.  Any other number will leave strand 2
   * as strand 2.
   */
{
  MaximalElement = pretzel[0].first;
  for (size_t i = 1; i < pretzel.size(); i++)
    if (pretzel[i].first > MaximalElement) MaximalElement = pretzel[i].first;

  for (int n = 1; n <= MaximalElement + 1; n++)
  {
    int m = n;
    for (int i = 0; i < pretzel.size(); i++)
    {
      if (pretzel[i].first == m)
        m++;
      else if (pretzel[i].first == m - 1)
        m--;
    }
    strand_permutation.push_back(m);
  }
}


/** This function assumes that "strand_permutation" is a vector whose
 ** entries are 1-based indexes of itself; it returns the number of
 ** cycles in the permutation.
 **/
int countCycles(const IntVector & strand_permutation)
  /* The algorithm is best illustrated with an example.  If our permutation
   * is (3 1 2 4) this means that strand 1 ends up as strand 3; strand 2 ends 
   * up as strand 1; strand 3 ends up as strand 2 and strand 4 stays as strand 4.  
   * This gives us a sequence 1 -> 3 -> 2 -> 1 which is a cycle.  The algorithm 
   * follows these sequences until it reaches the number where it started, and each 
   * time it does this it increases the count of the cycles.  The numbers through 
   * which the sequence goes (in this case 2 and 3) are recorded so that these don't 
   * have to be checked as potential starting points of new cycles.
   */
{
  if (!strand_permutation.size()) { std::cerr << "Error: No strand permutations!" << std::endl; exit(4); }
  std::vector<bool> strand_item_checked(strand_permutation.size(), false);

  int count(0);

  for (size_t i = 0; i < strand_permutation.size(); i++)
  {
    if (strand_item_checked[i]) continue;
    if (strand_permutation[i] == 0) { std::cerr << "Error: Invalid strand permutation!" << std::endl; exit(4); }

    strand_item_checked[i] = true;
    int k = strand_permutation[i];
    while (k - 1 != i)
    {
      strand_item_checked[k-1] = true;
      k = strand_permutation[k - 1];
    }
    count++;
  }
  return count;
}


/** If a braid does not mention a particular strand, it decomposes into disjoint groups of strands.
 ** These are calculated here and stored in separate vectors.
 **/
void computeGroups(const Pretzel & pretzel, std::vector<Pretzel> & groups)
{
  IntVector missing; // store the missing numbers here.
  std::vector<unsigned int> abses(pretzel.size());
  for (size_t i = 0; i < pretzel.size(); i++)
    abses[i] = pretzel[i].first;

  // "abses" will contain a unique, sorted list of absolute values in the braid:
  std::sort(abses.begin(), abses.end());
  abses.erase(std::unique(abses.begin(), abses.end()), abses.end());

  // "missing" will contain a list of all values not in "abses":
  int diff = 1;
  for (size_t i = 0; i < abses.size(); i++)
    while(abses[i] != i + diff)
    {
      missing.push_back(i + diff);
      diff++;
    }


  // Now remove from "missing" all consecutive numbers, keeping track of the number of removals in "diff":
  diff = 0;
  for (IntVector::iterator i = missing.begin(); i != missing.end(); )
  {
     if (*i + 1 == *(i+1)) { i = missing.erase(i); diff++; }
     else i++;
  }
  
  // Now extract the independent strands from "pretzel" into "groups":
  for (size_t i = 0; i < missing.size(); i++)
  {
    groups.push_back(Pretzel());
    int below = ( (i == 0) ? -1 : missing[i - 1] );
    for (size_t j = 0; j < pretzel.size(); j++)
      if (pretzel[j].first < missing[i] && int(pretzel[j].first) > below)
        groups[i].push_back(pretzel[j]);
  }
  groups.push_back(Pretzel());
  if (missing.size())
    for (size_t j = 0; j < pretzel.size(); j++)
      if (pretzel[j].first > missing[missing.size() - 1])
        groups[groups.size() - 1].push_back(pretzel[j]);

  // Lastly, each reduction in "diff" corresponds to an unlinked unknot, so we add empty subbraids:
  while (diff--) groups.push_back(Pretzel());
}


/** Given a braid, this function finds the homology generators:
 ** the crossings braid[i] and braid[homology[i]-1] are adjacent.
 ** homology[i] = 0 means there is no adjacency.
 **/
std::vector<size_t> computeHomology(const Pretzel & pretzel)
  /* The algorithm takes each crossing in turn and looks through
   * the braid to find the next crossing with the same modulus.
   * This is because the modulus of the crossing tells us between
   * which strands it lies.
   */
{
  if (!pretzel.size()) return std::vector<size_t>();

  std::vector<size_t> homology(pretzel.size() - 1, 0);
  for (size_t i = 0; i < homology.size(); i++)
  {
    size_t j = i;
    while (homology[i] == 0 && j < homology.size())
    {
      j++;
      if (pretzel[i].first == pretzel[j].first)
        homology[i] = j + 1;
    }
  }
  return homology;
}

/* Type of the entries of the matrices used below to compute the Alexander polynomial.
 * In a previous version, double was used, which led to overflows for some braids.
 * mpq_class is a rational number of arbitrary precision provided by the GMP library.
 */
typedef mpq_class AlexMatrixEntriesType;

mpq_class pow(const mpq_class& r, int e) {
    mpq_class result(1);
    for (int i = 0; i < e; ++i)
        result *= r;
    return result;
}

/* The main function */
int main(int argc, char * argv[])
{
  Pretzel              pretzel;            // The pretzel data given by the user
  IntVector            strand_permutation; // The permutation of the strands
  unsigned int         MaximalElement;     // The maximal element of "pretzel[].first"
  int                  components;         // The number of components of the pretzel (= cycles in the above)
  std::vector<Pretzel> groups;             // Sets of independent pretzels (e.g. [1 1 3 3] -> { [1 1], [3 3] }
  int                  genus1, genus2;     // The genus of the Seifert surface (using two different formulae)


  /* Read input into the vector "braid" */
  try { readInput(pretzel); }
  catch(...) { std::cerr << "An error occurred while reading your input." << std::endl
                         << "Make sure twist numbers are odd integers." << std::endl
                         << "Goodbye." << std::endl; exit(1); }
  if (!pretzel.size()) { std::cout << "No valid data entered! Goodbye." << std::endl; exit(1); }


  /* Establish the permutation of strands; stored in "strand_permutation" */
  computeStrandPermutations(pretzel, strand_permutation, MaximalElement);
  

  /* Next we calculate how many components the link has; i.e. the number of cycles in "strand_permutation" */
  components = countCycles(strand_permutation);


  /* If a strand does not get menitoned in the braid, the braid decomposes into disjoint groups of strands. */
  computeGroups(pretzel, groups);


  /* We reorder the braid silently to be partitioned along the groups */
  if (pretzel.size() > 1)
  {
    Pretzel::iterator k = pretzel.begin();
    for (size_t i = 0; i < groups.size(); i++)
      k = std::copy(groups[i].begin(), groups[i].end(), k);
  }


  /* Find the homology generators: the crossings pretzel[i] and pretzel[homology[i]-1] are adjacent.
     homology[i] = 0 means there is no adjacency. */
  std::vector<size_t> homology = computeHomology(pretzel);


  /* Compute the Seifert matrix "sm" from the homology basis. */
  SquareMatrix<int> sm(homology.size(), 0);
  for (size_t i = 0; i < homology.size(); i++)
  {
    if (!homology[i]) continue;

    for (size_t j = i; j < homology.size(); j++)
    {
      /* Self-linking */
      if (i == j)  sm(i, j) = -(pretzel[i].second + pretzel[homology[i] - 1].second) / 2;

      /* See Section 3.3 case 1 */
      else if (homology[i] > homology[j]) ;

      /* See Section 3.3 case 2 */
      else if (homology[i] < j + 1)       ;

      /* See Section 3.3 case 3 */
      else if (homology[i] == j + 1)
      {
        sm(i, j) = (pretzel[j].second - 1) / 2;
        sm(j, i) = (pretzel[j].second + 1) / 2;
      }

      /* See Section 3.3 case 4 */
      else if (abs(int(pretzel[i].first) - int(pretzel[j].first)) > 1) ;
      
      /* See Section 3.3 case 5 */
      else if (int(pretzel[i].first) - int(pretzel[j].first) ==  1)    sm(j, i) = -1;
      else if (int(pretzel[i].first) - int(pretzel[j].first) == -1)    sm(i, j) =  1;

      /* Debug: There should never be a 555 */
      else
      {
        sm(i, j) = sm(j, i) = 555;
        std::cerr << "Warning: Error in Seifert Matrix Algorithm! (i = " << i
                  << ", j = " << j << ", hom[i] = " << homology[i]
                  << ", hom[j] = " << homology[j] << ")" << std::endl;
      }

    }
  }


  /* Output some auxilliary data to stderr. */
  std::cerr << "You entered the pretzel " << pretzel << "." << std::endl
  //        << "The permutations are: " << strand_permutation << "." << std::endl
            << "This has " << components << " component(s), i.e. it is a ";
  if (components == 1) std::cerr << "knot.";
  else if (components > 1) std::cerr << "link.";
  else if (components == 0) std::cerr << "not very interesting object.";
  else { std::cerr << std::endl << "Error: Number of cycles is negative!" << std::endl; exit(2); }
  std::cerr << std::endl
  //        << "The homology data is: " << homology << "." << std::endl
  //        << "The unpruned Seifert data is " << sm << std::endl
  ;


  /* Remove zero-rows/columns. Start from the back so it can be done in a single pass.
   * Here we use index type "int" to allow for negative values in the stop condition.
   */
  for (int i = homology.size() - 1; i >= 0 ; i--)
    if (!homology[i]) { sm.remove_row(i); sm.remove_column(i); }


  /* Calculate the genus of the Seifert surface(s). */
  if (components == 0) { std::cerr << "The link has no components. Goodbye." << std::endl; exit(3); }
  genus1 = (int(sm.dim()) + int(groups.size()) - components) / 2;
  genus2 = (2 * int(groups.size()) - components - int(MaximalElement + 1) + int(pretzel.size())) / 2;

  /* The final output to stdout. */
  std::cout << "The Seifert matrix is: " << sm << std::endl;
  if (groups.size() > 1)
  {
    std::cout << "The pretzel is a disjoint union of unrelated sub-pretzels: { ";
    for (size_t i = 0; i < groups.size() - 1; i++)
      std::cout << groups[i] << ", ";
    std::cout << groups[groups.size() - 1] << " }." << std::endl
              << "The Seifert surface that was constructed has " << groups.size()
              << " disjoint components." << std::endl;
  }
  std::cout << "The genus of the Seifert surface is " << genus2;
  if (groups.size() > 1)
  {
    std::cout << " = ";
    for (size_t i = 0; i < groups.size(); i++)
    {
      if (groups[i].size() == 0)
        std::cout << "0";
      else
      {
        unsigned int MaxEl, MinEl;
        IntVector sub_strand_perms;
        Pretzel sub_pretzel(groups[i]);
        IntVector sub_pretzel_abs(groups[i].size());
        for (size_t j = 0; j < sub_pretzel.size(); j++) sub_pretzel_abs[j] = sub_pretzel[j].first;
        MinEl = *min_element(sub_pretzel_abs.begin(), sub_pretzel_abs.end()) - 1;
        for (size_t j = 0; j < sub_pretzel.size(); j++)
          sub_pretzel[j].first -= MinEl;
        //std::cerr << "Subbraid: " << sub_braid << std::endl;
        computeStrandPermutations(sub_pretzel, sub_strand_perms, MaxEl);
        std::cout << (2 - (int(MaxEl) + 1 - int(sub_pretzel.size())) - countCycles(sub_strand_perms)) / 2;
      }
      if (i < groups.size() - 1) std::cout << " + ";
    }
  }
 std::cout << "." << std::endl;

/* If the link has more than one connected component then we define the Alexander
   * polynomial to be zero, and end the program here. */
  if (groups.size() > 1)
  { std::cout << "The Alexander polynomial is 0." << std::endl; return 0; }

  /* Get the Alexander polynomial by fitting sm.dim() + 1 points:
   * "vandat" contains the data points 0..dim().
   * "vansol" contains the Vandermonde matrix augmentent by the solution vector p(0)..p(dim()).
   * The polynomial coefficients are in the last column of the ReducedRowEchelon(vansol).
   */
  std::vector<AlexMatrixEntriesType> poly_points(sm.dim() + 1);
  for (int i = 0; i < poly_points.size(); i++) poly_points[i] = i;
  Matrix<AlexMatrixEntriesType> vandat(vandermonde(poly_points));
  Matrix<AlexMatrixEntriesType> vansol(vandat.rows(), vandat.cols() + 1);
  for (int i = 0; i < vandat.rows(); i++)
  {
    for (int j = 0; j < vandat.cols(); j++)
      vansol(i,j) = vandat(i,j);

    // Here we make the Alexander polynomial!
    SquareMatrix<AlexMatrixEntriesType> am(sm);
    SquareMatrix<AlexMatrixEntriesType> AlexanderMatrix = am + am.transpose() * (-poly_points[i]);
    vansol(i, vandat.cols()) = AlexanderMatrix.determinant();  // fill in the solutions
    /* std::cout <<  "Alexander matrix for t = " << i << " is " << AlexanderMatrix
                 << " with determinant " << AlexanderMatrix.determinant() << "." << std::endl; */
  }
  Matrix<AlexMatrixEntriesType> Solution(vansol.gauss_jordan());

  std::cout << "The Alexander polynomial is p(t) = ";

  bool found_nonzero_term(false);
  for (int i = Solution.rows() - 1; i >= 0; i--)
  {
    if (Solution(i, Solution.cols() - 1) != 0)
    {
      if (found_nonzero_term) std::cout << " + ";
      found_nonzero_term = true;
      std::cout << Solution(i, Solution.cols() - 1);
      if (i > 1) std::cout << " * t^" << i;
      else if (i == 1) std::cout << " * t";
    }
  }
  if (!found_nonzero_term) std::cout << "0";
  std::cout << "." << std::endl;

  return 0;
}