// Elementary quadrics calculations

#ifndef _qi_elem_h_
#define _qi_elem_h_

#include <sstream>

/** RPL */
#include <libqi/rpl/bigint_matrix.h>

/** QI */
#include <libqi/kernel/QIHompoly.h>
#include <libqi/kernel/QIComponentTypes.h>
#include <libqi/kernel/QINumber.h>

using namespace std;

// Enter namespace QI
namespace QI {

// Find inertia of quadric using Descartes' rule
math_vector <int> inertia(const bigint_matrix &q);

// Signed inertia
math_vector <int> signed_inertia(const bigint_matrix &q);

// Find inertia of quadric using Descartes' rule when we know its rank is r
math_vector <int> inertia_known_rank(const bigint_matrix &q, 
					      const rpl_size_t &r);

// Find the inertia of a non-singular non-rational conic Q whose 3x3 matric is 
// q1+sqrt(D).q2. Being non-singular, the inertia is [3,0] or [2,1]. 
// Output 0 if inertia [3,0] and 1 otherwise 
bool inertia_non_rational_conic(const bigint_matrix &q1, const bigint_matrix &q2, 
				const bigint &D);

// Descartes: count number of sign changes in sequence
rpl_size_t descartes(const polynomial <bigint> &p);

// Compute the determinantal equation of the pencil generated by (q1,q2), 4x4 matrices
hom_polynomial <bigint> det_pencil(const bigint_matrix &q1, const bigint_matrix &q2);

// Compute the determinantal equation of the pencil generated by (q1,q2), 3x3 matrices
hom_polynomial <bigint> det_pencil3(const bigint_matrix &q1, const bigint_matrix &q2);

// Compute the determinantal equation of the pencil generated by (q1,q2), 2x2 matrices
hom_polynomial <bigint> det_pencil2(const bigint_matrix &q1, const bigint_matrix &q2);

// Compute the gcd of diff(D,lambda) and diff(D,mu)
hom_polynomial <bigint> gcd_pencil(const hom_polynomial <bigint> &pol);

// Nice priting of plane equation
void print_plane(const bigint_matrix &m, ostream &s);

// Nice printing of equation of quadric
void print_quadric(const bigint_matrix &m, ostream &s);

// Convert list of coeffs (x^2, xy, xz, xw, y^2, yz, yw, z^2, zw, w^2) to matrix form
bigint_matrix vectoquad(const math_vector <bigint> &qvec);

// $$ JC added
/** A proper way to designate the following function */
#define vector2matrix vectoquad

string quad2string (const bigint_matrix quad, bool homogeneous);

/** Returns the Euclidean type of the matrix */
InputEuclideanType getEuclideanType (const bigint_matrix &mat);

// Convert matrix to list of coeffs
math_vector <bigint> quadtovec(const bigint_matrix &q);

// Kill row k and column l of matrix q
bigint_matrix mat_minor(const bigint_matrix &q, const rpl_size_t &k, const rpl_size_t &l);

// Cofactor matrix
bigint_matrix cofactors_mat(const bigint_matrix &q);

// Compute the singular locus of matrix q - Output is a matrix of vectors forming
// a basis of singular locus
bigint_matrix singular(const bigint_matrix &q);

// Intersect two linear spaces
bigint_matrix linear_intersect(const bigint_matrix &m1, const bigint_matrix &m2);

// Compute a 4x4 projective transformation sending infinity (0 0 0 1) to the point
bigint_matrix send_to_infinity(const math_vector <bigint> &point);

// Compute a 4x4 projective transformation sending the line point1-point2 to z = w = 0
bigint_matrix send_to_zw(const math_vector <bigint> &point1, 
			 const math_vector <bigint> &point2);

// Nice printout of non-rational roots
void print_root(const math_vector <bigint> &root, ostream &s, const int sig);

// Compare two (projective) points
bool are_equal(const math_vector <bigint> &v, const math_vector <bigint> &w);

// Tell if k1+sqrt(d)*k2 = (a+sqrt(d)*b)*(k1'+sqrt(d)*k2')
bool are_equal(const bigint_matrix &k1, const bigint_matrix &k2, 
	       const bigint_matrix &k1p, const bigint_matrix &k2p, const bigint &d);

// Output a quadric of the pencil through the given rational point
bigint_matrix pencil_quadric_through_ratpoint(const bigint_matrix &q1, 
					      const bigint_matrix &q2,
					      const math_vector <bigint> &pt,
					      math_vector <bigint> &l);

// Sign of the algebraic number a+b*sqrt(c)
int sign(const bigint &a, const bigint &b, const bigint &c);

// Gauss reduction of quadratic form (used only in qi_inter_with_22)
void gauss(const bigint_matrix &q, bigint_matrix &tm, bigint_matrix &can);

// Compute a random matrix of bigints
bigint_matrix rand_mat(const bigint &size);

// Compute a random bigint (including negative numbers)
bigint rand_bigint(const bigint &size);

// Solve the equation q v = l dir in integers
math_vector <bigint> solve_proj(const bigint_matrix &q, const math_vector <bigint> &dir, 
				const math_vector <bigint> &sing);

} // end of namespace QI

#endif