#ifndef POLYTOPE_H
#define POLYTOPE_H
 
#include <iostream>
#include <giac/config.h>
#include <giac/giac.h>

#include <string>
#include <map>
#include <algorithm>
#include <ctime>

using namespace std;
using namespace giac;

#include "CMfield.h"
#include "Prism.h"
#include "Util.h"
#include "edge.h"
#include "ridge.h"
#include "sgraph.h"
#include "FiniteGroup.h"
#include "FiniteLinearGroup.h"

extern context ct;
class Polytope
{
  public:

  Polytope(string s);

  ~Polytope(){
    exit(1);
  };

  std::clock_t c_start, c_end;
  
  bool verbose;

  CMfield K;

  void createGroup();
  void startAlgo();
  void startFromFaceData();

  bool is_symmetric;

  bool iscompact, isarithmeticitychecked, isarithmetic, naindex;

  bool emb_dim1, emb_dim2, check_cycles, check_euler, check_sphere, perform_algo, draw_pics;

  gen p0, p0c, p0_rootof, p0c_rootof;
  bool is_p0_rootof_computed = false;
  
  gen H; // Matrix of the ambient Hermitian form
  gen H_rootof;
  gen R1, R2, R3, R1i, R2i, R3i, J, Ji, P, Pi, Ppower;

  vecteur Ppowers;

  bool needppower;

  gen Q; // change of basis to diagonal form, using eigenvectors of P

  // Matrix of the symmetry of the polytope (this may be (123) or (1J))
  bool issymmetric;

  int tnum, tden; // only used for Mostow groups
  int pval; // pval is the order of generators
  gen taugen, tau, tauc, mult, multconj; 

  WVector e1w, e2w, e3w;

  gen extramirror;

  int orderofp, orderofppower;

  gen euler_char, euler_char_top;

  gen ipmin;

  vector<int> word_p1, word_m1, word_123, word_321, w1, w2, w3;
  vector<Prism> faces;
  vector<Prism> faceReps;

  vector<WVector> vertices;
  vector<WVector> vertReps;

  vector<vector<int> > vertexPOrbits;
  vector<vector<int> > vertexOrbits;
  vector<int> currentVertexOrbit;
  vector<int> orbitsToVisit;

  vector<vector<int> > facePOrbits;
  vector<vector<int> > edgePOrbits;
  vector<vector<int> > ridgePOrbits;

  vector<vector<int> > vertexPairings;
  vector<vector<int> > vertexPOrbitPairings;
  vector<vector<int> > vertexPOrbitPairingMaps;
  vector<vector<int> > edgePairings;
  vector<vector<int> > ridgePairings;
  vector<int> facePairings;

  vector<edge> edges;
  vector<ridge> ridges;

  vector<vector<int> > relations;

  gen Proj(gen v, gen w);  
  gen Inn(gen v, gen w);
  gen SqNorm(gen v);
  gen BoxProd(gen v, gen w);

  gen refl(gen x, gen v, gen z);
  gen reflMat(gen v, gen z);

  bool testSame(WVector V, WVector W);
  bool isInList(WVector V, vector<WVector> li);
  int findInList(WVector V, vector<WVector> li);
  int findInList(gen V, vector<WVector> li);

  bool isInList(gen a, vector<gen> list); 

  void initSporadicGroup();
  void initThompsonGroup(); 
  void initMostowGroup();   

  Prism createPrism(WVector v0, WVector v1, WVector v2);
  void updateFreeRidges();
  void expand();

  bool isPrismNew(Prism A);
  bool isVertexNew(WVector V);
  bool comparePrisms(Prism A, Prism B);
  bool samePolygon(vector<WVector> V, vector<WVector> W);  
  bool samePolygon(vector<gen> V, vector<gen> W);  

  bool computeVertices(Prism &A);

  bool isOnPrism(gen X, Prism A);
  bool isOnMirror(gen X, Prism A);
  bool isOnTopFace(gen X, Prism A);
  gen reflectAcrossSpine(gen X, Prism A);

  Prism createPimage(Prism A);
  Prism createPinvimage(Prism A);

  Prism applyReflection(Prism A);
  Prism applyInverseReflection(Prism A);

  Prism applyMatrix(gen M, Prism A);

  void addFaceOrbit(Prism A);

  bool isVertexSafe(WVector V, Prism A);

  vector<int> findFaces(WVector v);
  vector<int> findFaces(gen V);
  vector<int> findBisectors(WVector v);
  vector<int> findBisectors(gen V);
  vector<WVector> findMirrors(WVector v);
  vector<WVector> findMirrors(gen V);

  void useRelations(vector<int> &w);
  void useRelations(Prism &A);

  void computeStabilizer(int k);
  void studyVertexStabilizers();
  void findNextVertexInCycle(int k);

  void addEdge(edge e);
  bool isSame(edge e1, edge e2);
  bool isFlipped(edge e1, edge e2);
  bool isEdgeNew(edge e);

  void addRidge(ridge e);
  bool isSame(ridge e1, ridge e2);
  bool isRotated(ridge e1, ridge e2);
  bool isRotatedFlip(ridge e1, ridge e2);
  bool isRidgeNew(ridge e);

  int applyPtoVertex(int k);
  vector<int> getPOrbitOfVertex(int k);
  vector<int> getPOrbitOfFace(int k);
  int findNextInFacePOrbit(int k);
  void addEdgeOld(edge e);

  void createBraidRelation(vector<int> w1, vector<int> w2);

  gen evalWord(vector<int> w);
  gen evalWordAsRootOf(vector<int> w);

  int applyPairing(ridge e);
  void checkCycles();

  int applyPairingToVertex(int jface, int j);
  int applyPairingToEdge(int jface, int j);
  int applyPairingToRidge(int jface, int j);

  int applyPToVertex(int pow, int j);
  int applyPToEdge(int pow, int j);
  int applyPToRidge(int pow, int j);
  int applyPinvToVertex(int pow, int j);
  int applyPinvToEdge(int pow, int j);
  int applyPinvToRidge(int pow, int j);

  void drawFace(int j); 
  void drawFaceSector(int j, int k1, int k2); 
  void drawProjection(int axis); 

  bool checkSide();
  void checkEdge(int j);
  void checkRidge(int j);

  void checkPairings();
  void computeAllVertices();
  void computeAllEdges();
  void computeAllRidges();

  bool checkt1val(gen rur, gen par, gen t1v, int &res);
  bool checkt1valAttempt(gen rur, gen par, gen t1v, int &res);
  bool checkEqualt1val(gen rur, gen par, gen t1v);

  bool checkt2val(gen rur, gen par, gen t1v, int &res);
  bool checkt2valAttempt(gen rur, gen par, gen t1v, int &res);
  bool checkEqualt2val(gen rur, gen par, gen t1v);

  bool checkt3val(gen rur, gen par, gen t3v, gen vl, gen wl, gen vt, int &res);
  bool checkt3valAttempt(gen rur, gen par, gen t3v, gen vl, gen wl, gen vt, int &res);
  bool checkEqualt3val(gen rur, gen par, gen t3v, gen vl, gen wl, gen vt);

  bool checkEqSign(gen rur, gen par, gen eq, int &res);
  bool checkEqSignAttempt(gen rur, gen par, gen eq, int &res);
  bool checkEqZero(gen rur, gen par, gen eq);

  void drawRidgePic(int i);

  void computeEulerCharacteristic();
  vecteur eulerdim0();
  vecteur eulerdim1();
  vecteur eulerdim2();
  vecteur eulerdim3();

  bool isStabilizedByPpower(int j);
  bool isRidgeComplex(int j, gen &mi);

  void exportFaceData();
  void importFaceData();

  void checkSphere();
  void checkLinks();

  gen computeAngle(int j, gen X);
  gen computeAngleNum(int j, gen X);

  bool expressAsWord(gen M, vecteur list, gen fixed_pt, gen mirror, vector<int> &inds);
  
  bool isArithmetic();
  bool isArithmeticAlt();
  string checkSignature(gen A);

  bool checkTraceField(gen &res);

  bool findEigenVectors_P();

  bool checkPairOfPrisms(Prism A1, Prism A2, bool &res);

  gen lowerDegree(gen g);

  bool isPPower(gen M);

  int testPoint(gen V);

  bool findMirror(gen M, gen &res);

  bool my_jordan(gen M, vecteur &bas, vecteur &evs);
  bool my_jordan_knowing_linear_order(gen M, int o, vecteur &bas, vecteur &evs);
  bool my_jordan_knowing_linear_order_alt(gen M, int o, vecteur &bas, vecteur &evs);
 
  bool identify_in_field(gen x, gen pol, gen &res);

  void reduce_degree(gen &x);
  void reduce_degree_vect(gen &x);
  void reduce_degree_mat(gen &x);

  bool my_giac_kernel(gen M, gen &res);

  bool my_giac_are_dep(gen v, gen w);

  void normalize_vector(gen &v);

  string convert_time(long double x);
  
};

#endif