#  Copyright (c) 1997-2015
#  Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
#  http://www.polymake.org
#
#  This program is free software; you can redistribute it and/or modify it
#  under the terms of the GNU General Public License as published by the
#  Free Software Foundation; either version 2, or (at your option) any
#  later version: http://www.gnu.org/licenses/gpl.txt.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#-------------------------------------------------------------------------------

# @topic category property_types/Algebraic Types
# This category contains all "algebraic" types, such as matrices, vectors, polynomials, rings, ...

# @category Algebraic Types
# A type for vectors with entries of type //Element//.
# 
# You can perform algebraic operations such as addition or scalar multiplication.
# 
# You can create a new Vector by entering its elements, e.g.:
#	$v = new Vector<Int>(1,2,3);
# or
#	$v = new Vector<Int>([1,2,3]);
# @tparam Element
declare property_type Vector<Element=Rational> : c++ (include => ["polymake/Vector.h", "polymake/Matrix.h"], operators => '@arith:wary |:anchor |= @compare:wary') {

   method construct(Int) : c++;

   type_method init {
      my ($proto)=@_;
      if ($proto->params->[0] == typeof Float) {
         $proto->equal=\&compare_float_sequences;
      }
   }

   # The length of the vector
   # @return Int
   user_method dim() : c++;

   # Returns a [[Vector]] containing only the entries after //i//.
   # @param Int i
   # @return Vector
   user_method slice(*:anchor) : lvalue_opt : wary : c++;

   # Returns a [[Vector]] containing only the entries from //i// to //j//-1.
   # @param Int i
   # @param Int j
   # @return Vector
   user_method slice($$) : lvalue_opt : wary : c++;

   # Divides every entry by //a// (assuming that every entry is divisible by //a//).
   # @param Int a
   # @return Vector
   user_method div_exact(*) : lvalue : c++;
}

function permuted(Vector, *) : c++;

function permuted_inv(Vector, *) : c++;


# @category Arithmetic
# Compute the __greatest common divisor__ of the elements of the given vector.
# @param Vector<Element> v
# @return Element
# @example > $v = new Vector<Int>(3,6,9);
# > print gcd($v);
# | 3

user_function gcd(Vector) : c++ (include => "polymake/linalg.h");

# @category Arithmetic
# Compute the least common multiple of the elements of the given vector.
# @param Vector<Element> v
# @return Element
# @example > $v = new Vector<Integer>(1,3,6,9);
# > print lcm($v);
# | 18

user_function lcm(Vector) : c++ (include => "polymake/linalg.h");

##################################################################################

# @category Algebraic Types
# @tparam Scalar default: [[Rational]]
declare property_type Plucker<Scalar=Rational> : c++ \
   (operators => '+ *', include => "polymake/Plucker.h") {

    method construct(Vector) : c++;

    method construct(Int, Int, Vector) : c++;

    method construct(Int, Int) : c++;

    method point() : c++;

    user_method permuted(*) : c++; 

    user_method coordinates() : c++; 
}

##################################################################################

# @category Algebraic Types
# @tparam Element default: [[Rational]]
# @tparam Sym default: [[NonSymmetric]]
declare property_type Matrix<Element=Rational, Sym=NonSymmetric> : \
   c++ (name => "Matrix<%1>", include => "polymake/Matrix.h", operators => '@arith_nodiv:wary /:wary:anchor |:wary:anchor /=:wary |=:wary @compare:wary') {

   method construct(Int,Int) : c++;

   method construct(Vector+) {
      my ($proto, $vectors)=@_;
      my $M=$proto->construct->(scalar(@$vectors), $vectors->[0]->dim);
      my $i=0;
      $M->[$i++]=$_ for @$vectors;
      $M;
   }

   type_method init {
      my ($proto)=@_;
      if ($proto->params->[0] == typeof Float) {
         $proto->equal=sub : method {
            (undef, my $m1, my $m2)=@_;
            my $vec_proto=typeof Vector<Float>;
            if ($m1->rows==$m2->rows) {
               my $i=0;
               foreach my $v1 (@$m1) {
                  $vec_proto->equal->($v1, $m2->row($i++)) or return 0;
               }
               1
            }
         };
      }
   }

   # Returns the number of rows.
   # @return Int
   user_method rows() : c++;

   # Returns the number of columns.
   # @return Int
   user_method cols() : c++;

   # Returns the //i//-th row.
   # @param Int i
   # @return Vector<Element>
   user_method row($) : lvalue_opt : wary : c++;

   # Returns the //i//-th column.
   # @param Int i
   # @return Vector<Element>
   user_method col($) : lvalue_opt : wary : c++;

   # Returns a __minor__ of the matrix containing the rows in //r// and the columns in //c//. You can pass [[all_rows_or_cols|All]] if you want all rows or columns and ~ for the complement of a set. E.g. 
   #	$A->minor(All, ~[0]);
   # will give you the minor of a matrix containing all rows and all but the 0-th column.
   # @param Set r the rows
   # @param Set c the columns
   # @return Matrix
   user_method minor(*:anchor,*:anchor) : lvalue_opt : wary : c++;

   # Returns an element of the matrix.
   # The return value is an `lvalue', that is, it can be modified if the matrix object is mutable.
   # @param Int r the row index
   # @param Int c the column index
   # @return Element
   user_method elem($,$) : lvalue_opt : wary : c++(name => '()');

   # Returns the __diagonal__ of the matrix.
   # @param Int i //i//=0: the main diagonal (optional)
   #				//i//>0: the //i//-th diagonal __below__ the main diagonal
   #				//i//<0: the //i//-th diagonal __above__ the main diagonal
   # @return Vector<Element>
   user_method diagonal(;$=0) : lvalue_opt : wary : c++;

   # Returns the __anti-diagonal__ of the matrix.
   # @param Int i //i//=0: the main anti_diagonal (optional)
   #				//i//>0: the //i//-th anti_diagonal __below__ the main anti_diagonal
   #				//i//<0: the //i//-th anti_diagonal __above__ the main anti_diagonal
   # @return Vector<Element>
   user_method anti_diagonal(;$=0) : lvalue_opt : wary : c++;

   # Divides every entry by //a// (assuming that every entry is divisible by //a//).
   # @param Int a
   # @return Matrix
   user_method div_exact(*) : lvalue : c++;

   # backward compatibility, see issue #864
   use overload '&{}' => \&deprecated_elem_access;
}

# @category Linear Algebra
# Projects points into the [[null_space|orthogonal complement]] of a subspace given via an orthogonal basis.
# The given points will be overwitten. 
# @param Matrix points will be changed to orthogonal ones 
# @param Matrix orthogonal basis of the subspace
user_function project_to_orthogonal_complement(Matrix:lvalue,Matrix) : void : c++ (include => "polymake/linalg.h");


# @category Algebraic Types
# Use the keyword "All" for all rows or columns, e.g. when constructing a [[minor]].
declare property_type all_rows_or_cols : c++ (name => 'pm::all_selector', include => ["polymake/Matrix.h"], builtin => enum { All });


# @category Data Conversion
# Explicit conversion to a different element type.
# @param Vector v
# @tparam Target
# @return Vector<Target>
# @example > $v = new Vector<Rational>(1/2,2/3,3/4);
# > $vf = convert_to<Float>($v);
# > print $vf;
# | 0.5 0.6666666667 0.75

user_function convert_to<Element>(Vector) : c++ {
   if ($_[0]->type->params->[0]==typeof Element) {
      return $_[0];
   }
}

# @category Data Conversion
# Explicit conversion to a different element type.
# @param Matrix m
# @tparam Target
# @return Matrix<Target>
# @example > $M = new Matrix<Rational>([1/2,2],[3,2/3]);
# > $Mf = convert_to<Float>($M);
# > print $Mf;
# | 0.5 2
# | 3 0.6666666667

user_function convert_to<Element>(Matrix) : c++ {
   if ($_[0]->type->params->[0]==typeof Element) {
      return $_[0];
   }
}

# @category Data Conversion
# Create a [[Matrix]] by repeating the given [[Vector]] as rows.
# @param Vector v
# @param Int i
# @example > $v = new Vector(23,42,666);
# > $M = repeat_row($v,3);
# > print $M;
# | 23 42 666
# | 23 42 666
# | 23 42 666

user_function repeat_row(Vector:anchor,$) : c++ (include => "polymake/Matrix.h");

# @category Data Conversion
# Create a [[Matrix]] by repeating the given [[Vector]] as cols.
# @param Vector v
# @param Int i
# @example > $v = new Vector(23,42,666);
# > $M = repeat_col($v,3);
# > print $M;
# | 23 23 23
# | 42 42 42
# | 666 666 666

user_function repeat_col(Vector:anchor,$) : c++ (include => "polymake/Matrix.h");

# @category Data Conversion
# Convert a [[Vector]] to a [[Matrix]] with a single row.
# @param Vector v
# @return Matrix
# @example This converts a vector into a row and prints it and its type:
# > $v = new Vector([1,2,3,4]);
# > $V = vector2row($v);
# > print $V;
# | 1 2 3 4
# > print $V->type->full_name;
# | Matrix<Rational, NonSymmetric>
user_function vector2row(Vector:anchor) : c++ (include => "polymake/Matrix.h");


# @category Data Conversion
# Convert a [[Vector]] to a [[Matrix]] with a single column.
# @param Vector v
# @return Matrix
# @example This converts a vector into a column and prints it and its type:
# > $v = new Vector([1,2,3,4]);
# > $V = vector2col($v);
# > print $V;
# | 1
# | 2
# | 3
# | 4
# > print $V->type->full_name;
# | Matrix<Rational, NonSymmetric>
user_function vector2col(Vector:anchor) : c++ (include => "polymake/Matrix.h");


# @category Linear Algebra
# Produces a [[SparseMatrix]] from its diagonal.
# @param Vector d the diagonal entries
# @return SparseMatrix
# @example > $v = new Vector(1,2,3,4);
# > $D = diag($v);
# > print $D;
# | (4) (0 1)
# | (4) (1 2)
# | (4) (2 3)
# | (4) (3 4)

user_function diag(Vector:anchor) : c++ (include => "polymake/SparseMatrix.h");


# @category Linear Algebra
# Returns a __block diagonal matrix__ with blocks //m1// and //m2//.
# @param Matrix m1
# @param Matrix m2
# @return SparseMatrix
# @example > $m1 = new Matrix([1,2],[3,4]);
# > $m2 = new Matrix([1,0,2],[3,4,0]);
# > $D = diag($m1,$m2);
# > print $D;
# | (5) (0 1) (1 2)
# | (5) (0 3) (1 4)
# | 0 0 1 0 2
# | 0 0 3 4 0

user_function diag(*:anchor,*:anchor) : c++ (include => "polymake/SparseMatrix.h");


# @category Linear Algebra
# Produces a [[SparseMatrix]] from its anti-diagonal.
# @param Vector d the anti-diagonal entries
# @return SparseMatrix
user_function anti_diag(Vector:anchor) : c++ (include => "polymake/SparseMatrix.h");


# @category Linear Algebra
# Returns a __block anti-diagonal matrix__ with blocks //m1// and //m2//.
# @param Matrix m1
# @param Matrix m2
# @return SparseMatrix
user_function anti_diag(*:anchor,*:anchor) : c++ (include => "polymake/SparseMatrix.h");


# @category Data Conversion
# Returns an array containing the rows of //A//.
# @param Matrix A
# @return ARRAY<Vector>
# @example The following saves the rows of the vertex matrix of a square in
# the variable $w and then prints its contents using a foreach loop and concatenating
# each entry with the string "  ".
# > $w = rows(cube(2)->VERTICES);
# > foreach( @$w ){
# > print @{$_},"  ";
# > }
# | 1-1-1  11-1  1-11  111  
user_function rows(Matrix) : c++;

# @category Data Conversion
# Returns an array containing the columns of //A//.
# @param Matrix A
# @return ARRAY<Vector>
# @example The following saves the columns of the vertex matrix of a square in
# the variable $w and then prints its contents using a foreach loop and concatenating
# each entry with the string "  ".
# > $w = rows(cube(2)->VERTICES);
# > foreach( @$w ){
# > print @{$_},"  ";
# > }
# | 1111  -11-11  -1-111  

user_function cols(Matrix) : c++;

# @category Data Conversion
# Concatenates the rows of //A//.
# @param Matrix A
# @return Vector
# @example Make a vector out of the rows of the vertex matrix of a cube:
# > $v = concat_rows(cube(2)->VERTICES);
# > print $v;
# | 1 -1 -1 1 1 -1 1 -1 1 1 1 1
# @example For a sparse matrix, the resulting vector is sparse, too.
# > $vs = concat_rows(unit_matrix(3));
# > print $vs;
# | (9) (0 1) (4 1) (8 1)

user_function concat_rows(Matrix:lvalue_opt:anchor) : c++;


function permuted_rows(Matrix, *) : c++;

function permuted_inv_rows(Matrix, *) : c++;

function permuted_cols(Matrix, *) : c++;

function permuted_inv_cols(Matrix, *) : c++;


# @category Linear Algebra
# Computes the __transpose__ //A//<sup>T</sup> of a matrix //A//, i.e., (a<sup>T</sup>)<sub>ij</sub> = a<sub>ji</sub>.
# @param Matrix A
# @return Matrix
# @example > $M = new Matrix([1,2,23],[23,22,21]);
# > $Mt = transpose($M);
# > print $Mt;
# | 1 23
# | 2 22
# | 23 21

user_function transpose(Matrix:anchor) : c++ (name => 'T');


# @category Linear Algebra
# Computes the __determinant__ of a matrix using Gauss elimination.
# @param Matrix A
# @return Int det(A)
# @example > print det(unit_matrix(3));
# | 1

user_function det(Matrix:wary) : c++ (include => "polymake/linalg.h");

# @category Linear Algebra
# Computes the __trace__ of a matrix.
# @param Matrix A
# @return Int trace(A)
# @example > $M = new Matrix([1,2,3],[23,24,25],[0,0,1]);
# > print trace($M);
# | 26

user_function trace(Matrix:wary) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Computes the __rank__ of a matrix.
# @param Matrix A
# @return Int
user_function rank(Matrix) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Computes the __inverse__ //A//<sup>-1</sup> of an invertible matrix //A// using Gauss elimination.
# @param Matrix A
# @return Matrix
# @example We save the inverse of a small matrix M in the variable $iM:
# > $M = new Matrix([1,2],[3,4]);
# > $iM = inv($M);
# To print the result, type this:
# > print $iM;
# | -2 1
# | 3/2 -1/2
# As we can see, that is in fact the inverse of M.
# > print $M * $iM;
# | 1 0
# | 0 1

user_function inv(Matrix:wary) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Reduce a vector with a given matrix using Gauss elimination.
# @param Matrix A
# @param Vector b
# @return Vector
user_function reduce(Matrix:wary, Vector:wary) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Normalize a matrix by dividing each row by its length (l2-norm).
# @param Matrix<Float> A
# @return Matrix<Float>
# @example > $A = new Matrix<Float>([1.5,2],[2.5,2.5]);
# > print normalized($A);
# | 0.6 0.8
# | 0.7071067812 0.7071067812

user_function normalized(Matrix) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Computes subsets of the rows and columns of //A// that form a basis for the linear space spanned by //A//.
# @param Matrix A
# @return Pair<Set<Int>, Set<Int>> The first set corresponds to the rows, the second to the columns.
# @example Here we have a nice matrix:
# > $M = new Matrix([[1,0,0,0],[2,0,0,0],[0,1,0,0],[0,0,1,0]]);
# Let's print bases for the row and column space:
# > ($row,$col) = basis($M);
# > print $M->minor($row,All);
# | 1 0 0 0
# | 0 1 0 0
# | 0 0 1 0
# > print $M->minor(All,$col);
# | 1 0 0
# | 2 0 0
# | 0 1 0
# | 0 0 1
user_function basis(Matrix) : returns(@) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Does the same as [[basis]] ignoring the first column of the matrix.
# @param Matrix A
# @return Pair<Set<Int>, Set<Int>> The first set corresponds to the rows, the second to the columns.
user_function basis_affine(Matrix) : returns(@) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Computes a subset of the rows of //A// that form a basis for the linear space spanned by //A//.
# @param Matrix A
# @return Set<Int>
# @example Here we have a nice matrix:
# > $M = new Matrix([[1,0,0,0],[2,0,0,0],[0,1,0,0],[0,0,1,0]]);
# Let's print a basis of its row space:
# > print $M->minor(basis_rows($M),All);
# | 1 0 0 0
# | 0 1 0 0
# | 0 0 1 0
user_function basis_rows(Matrix) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Computes a subset of the columns of //A// that form a basis for the linear space spanned by //A//.
# @param Matrix A
# @return Set<Int>
# @example Here we have a nice matrix:
# > $M = new Matrix([[1,0,0,0],[2,0,0,0],[0,1,0,0],[0,0,1,0]]);
# Let's print a basis of its column space:
# > print $M->minor(All,basis_cols($M));
# | 1 0 0
# | 2 0 0
# | 0 1 0
# | 0 0 1
user_function basis_cols(Matrix) : c++ (include => "polymake/linalg.h");

function numerators(Vector<Rational>) : c++ (include => "polymake/linalg.h");

function numerators(Matrix<Rational>) : c++ (include => "polymake/linalg.h");

function denominators(Vector<Rational>) : c++ (include => "polymake/linalg.h");

function denominators(Matrix<Rational>) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# Creates a unit matrix of given dimension
# @tparam Element default: [[Rational]]
# @param Int d dimension of the matrix
# @return SparseMatrix<Element>
# @example The following stores the 3-dimensional unit matrix (ones on the diagonal and zeros otherwise) in a variable
# and prints it:
# > $M = unit_matrix(3);
# > print $M;
# | (3) (0 1)
# | (3) (1 1)
# | (3) (2 1) 
# @example The following stores the 3-dimensional unit matrix (ones on the diagonal and zeros otherwise) from type Int 
# in a variable and prints it:
# > $M = unit_matrix<Int>(3);
# > print $M->type->full_name;
# | SparseMatrix<Int, Symmetric>
user_function unit_matrix<Element=Rational>($) : c++ (include => "polymake/linalg.h") {
   if ($_[0] < 0) {
      croak( "unit_matrix - invalid dimension" );
   }
}

# @category Linear Algebra
# Creates a zero matrix of given dimensions
# @tparam Element default: [[Rational]]
# @param Int i number of rows
# @param Int j number of columns
# @return SparseMatrix<Element>
# @example The following stores a 2x3 matrix with 0 as entries (from type Rational) in a variable and prints it:
# > $M = zero_matrix(2,3);
# > print $M;
# | 0 0 0 
# | 0 0 0 
# @example The following stores a 2x3 matrix with 0 as entries from type Int in a variable and prints its type:
# > $M = zero_matrix<Int>(2,3);
# > print $M->type->full_name;
# | Matrix<Int, NonSymmetric>
user_function zero_matrix<Element=Rational>($$) : c++ (include => "polymake/linalg.h") {
   if ($_[0] < 0 || $_[1] < 0) {
      croak( "zero_matrix - invalid dimension" );
   }
}

# @category Linear Algebra
# Creates a [[SparseVector]] of given length //d// with a one entry at position //pos// and zeroes elsewhere.
# @tparam Element default: [[Rational]]
# @param Int d the dimension of the vector
# @param Int pos the position of the 1
# @return SparseVector<Element>
# # @example The following stores a vector of dimension 5 with a single 1 (as a Rational) at position 2:
# > $v = unit_vector(5,2);
# > print $v;
# | (5) (2 1)
# @example The following stores a vector of dimension 5 with a single 1 (as a Int) at position 2:
# > $v = unit_vector<Int>(5,2);
# > print $v->type->full_name;
# | SparseVector<Int>
# @example The following concatenates a unit vector of dimension 3 with a 1 at position 2 and a 
# unit vector of dimension 2 with a 1 at position 1:
# > $v = unit_vector(3,2) | unit_vector(2,1);
# > print $v;
# | (5) (2 1) (4 1)
user_function unit_vector<Element=Rational>($$) : c++ (include => "polymake/linalg.h") {
   if ($_[1] < 0 || $_[1] >= $_[0]) {
      croak( "unit_vector - invalid dimension or index out of range" );
   }
}

# @category Linear Algebra
# Creates a vector with all elements equal to zero.
# @param Int d vector dimension.  If omitted, a vector of dimension 0 is created,
#               which can adjust itself when involved in a block matrix operation.
# @tparam Element default: [[Rational]]
# @return Vector<Element>
# @example The following stores a vector of dimension 5 with 0 as entries (from type Rational) in a variable and prints it:
# > $v = zero_vector(5);
# > print $v;
# | 0 0 0 0 0
# @example The following stores a vector of dimension 5 with 0 as entries from type Int in a variable and prints its type:
# > $v = zero_vector<Int>(5);
# > print $v->type->full_name;
# | Vector<Int>
# @example The following concatenates a vector of dimension 2 of ones and a vector of length 2 of zeros:
# > $v = ones_vector(2) | zero_vector(2);
# > print $v;
# | 1 1 0 0
user_function zero_vector<Element=Rational>(;$=0) : c++ (include => "polymake/linalg.h") {
   if ($_[0] < 0) {
      croak( "zero_vector - invalid dimension");
   }
}

# @category Linear Algebra
# Creates a vector with all elements equal to 1.
# @param Int d vector dimension.  If omitted, a vector of dimension 0 is created, which can adjust itself when involved in a block matrix operation.
# @tparam Element default: [[Rational]].
# @return Vector<Element>
# @example To create the all-ones Int vector of dimension 3, do this:
# > $v = ones_vector<Int>(3);
# You can print the result using the print statement:
# > print $v;
# | 1 1 1

user_function ones_vector<Element=Rational>(;$=0) : c++ (include => "polymake/linalg.h") {
   if ($_[0] < 0) {
      croak( "ones_vector - invalid dimension");
   }
}

# @category Linear Algebra
# Compute the __null space__ of a matrix //A//.
# @param Matrix A
# @return Matrix
# @example > $A = new Matrix([1,2,0],[2,0,2]);
# > print null_space($A);
# | -1 1/2 1

user_function null_space(Matrix) : c++ (include => "polymake/linalg.h");

# @category Linear Algebra
# Compute the __null space__ of a vector //b//.
# @param Vector b
# @return Matrix
# @example > $b = new Vector(1,2,3);
# polytope > print null_space($b);
# | -2 1 0
# | -3 0 1

user_function null_space(Vector) : c++ (include => "polymake/linalg.h");

# @category Linear Algebra
# Compute the __lineality space__ of a matrix //A//.
# @param Matrix A
# @return Matrix
# @example > $M = new Matrix([1,1,0,0],[1,0,1,0]);
# > print lineality_space($M);
# | 0 0 0 1

user_function lineality_space(Matrix) : c++ (include => "polymake/linalg.h");

# @category Linear Algebra
# Computes the solution of the system //A//x = //b//
# @param Matrix A must be invertible
# @param Vector b
# @return Vector
# @example from the Wikipedia:
# > $A = new Matrix([3,2,-1],[2,-2,4],[-1,1/2,-1]);
# > $b = new Vector(1,-2,0);
# > print lin_solve($A,$b);
# | 1 -2 -2

user_function lin_solve(Matrix:wary, Vector:wary) : c++ (include => "polymake/linalg.h");

# @category Linear Algebra
# The matrix //A// is totally unimodular if the determinant of each square submatrix equals 0, 1, or -1.
# This is the naive test (exponential in the size of the matrix).
# 
# For a better implementation try Matthias Walter's polymake extension at
#  [[https://github.com/xammy/unimodularity-test/wiki/Polymake-Extension]].
# @param Matrix A
# @return Bool
# @example > $M = new Matrix<Int>([-1,-1,0,0,0,1],[1,0,-1,-1,0,0],[0,1,1,0,-1,0],[0,0,0,1,1,-1]);
# > print totally_unimodular($M);
# | 1

user_function totally_unimodular(Matrix) : c++ (include => "polymake/totally_unimodular.h");

# @category Linear Algebra
# Check whether both matrices are bases of the same linear subspace. 
# Note: It is assumed that they are *bases* of the row space.
# @param Matrix M1
# @param Matrix M2
# @return Bool
# @example > $M1 = new Matrix([1,1,0],[1,0,1],[0,0,1]);
# > $M2 = new Matrix([1,0,0],[0,1,0],[0,0,1]);
# > print equal_bases($M1,$M2);
# | 1

user_function equal_bases(Matrix, Matrix) {
   my ($M1, $M2)=@_;
   # we first need to check for full rank matrices as null_space gives an empty matrix in that case
   return $M1 == $M2 || $M1->rows==$M2->rows && 
          ($M1->rows == $M1->cols || is_zero(null_space($M1) * transpose($M2)));
}

# @category Linear Algebra
# Householder tranformation of Vector b. Only the orthogonal matrix reflection H is returned.
# @param Vector b
# @return Vector
user_function householder_trafo(Vector<Float>) : c++ (include => "polymake/linalg.h");


# @category Linear Algebra
# QR decomposition of a Matrix //M// with rows > cols
# @param Matrix<Float> M
# @return Pair<Matrix,Matrix>
# @example > $M = new Matrix<Float>([23,4],[6,42]);
# > $qr = qr_decomp($M);
# > print $qr->first;
# | 0.9676172724 0.2524218971
# | 0.2524218971 -0.9676172724
# > print $qr->second;
# | 23.76972865 14.47218877
# | 0 -39.63023785
# > print $qr->first * $qr->second ;
# | 23 4
# | 6 42

user_function qr_decomp(Matrix<Float>) : c++ (include => "polymake/linalg.h");

# @topic category property_types/Linear Algebra
# These types are needed as return types of algebraic computations.

# @category Linear Algebra
# Complete result of the __singular value decomposition__ of a matrix //M//,
# such that left_companion * sigma * transpose(right_companion) = //M//
# Contains the following fields:
#	Matrix<Float> sigma: the diagonalized matrix
#	Matrix<Float> left_companion: matrix of left singular vectors
#	Matrix<Float> right_companion: matrix of right singular vectors
declare property_type SingularValueDecomposition : c++ (include => "polymake/linalg.h");

# @category Linear Algebra
# SVD decomposition of a Matrix. Computes the SVD of a matrix into a  diagonal Marix (S), orthogonal square Matrix (U), orthogonal square Matrix (V), such that U*S*V^T=M
# The first element of the output array is S, the second U and the thrid V. 
# @param Matrix<Float> M
# @return SingularValueDecomposition
# @example > $M = new Matrix<Float>([1,2],[23,24]);
# > $SVD = singular_value_decomposition($M);
# The following prints the three matrices, seperated by newline characters.
# > print $SVD->left_companion ,"\n", $SVD->sigma ,"\n", $SVD->right_companion;
# | 0.06414638608 0.9979404998
# | 0.9979404998 -0.06414638608
# | 
# | 33.31011547 0
# | 0 0.6604600341
# | 
# | 0.6909846321 -0.7228694476
# | 0.7228694476 0.6909846321

user_function singular_value_decomposition(Matrix<Float>) : c++ (include => "polymake/linalg.h");

# @category Linear Algebra
# Moore-Penrose Inverse of a Matrix
# @param Matrix M 
# @return Matrix<Float>
user_function moore_penrose_inverse(Matrix<Float>) : c++ (include => "polymake/linalg.h");

# @category Data Conversion
# Return the input vector (which is already in dense form).
# @param Vector v
# @return Vector
user_function dense(Vector) { $_[0] }

# @category Data Conversion
# Return the input matrix (which is already in dense form).
# @param Matrix m
# @return Matrix
user_function dense(Matrix) { $_[0] }


# @category Combinatorics
# Returns the __permuation matrix__ of the permutation given by //p//.
# @tparam Scalar default: [[Int]]
# @param Array<Int> p
# @return Matrix<Scalar>
# @example The following prints the permutation matrix in sparse representation.
# > print permutation_matrix([1,0,3,2]);
# | (4) (1 1)
# | (4) (0 1)
# | (4) (3 1)
# | (4) (2 1)

user_function permutation_matrix<Scalar=Int>(*:anchor) : c++ (include => [ "polymake/permutations.h", "polymake/SparseMatrix.h" ]);

##################################################################################

# @category Algebraic Types
# A SparseVector is an associative container with element indices (coordinates) as keys; elements equal to the default value (ElementType(), which is 0 for most numerical types) are not stored, but implicitly encoded by the gaps in the key set. It is based on an AVL tree.
# 
# The printable representation of a SparseVector looks like a sequence (l) (p<sub>1</sub> v<sub>1</sub>) ... (p<sub>k</sub> v<sub>k</sub>),
# where l is the dimension of the vector and each pair (p<sub>i</sub> v<sub>i</sub>) denotes an entry with value
# v<sub>i</sub> at position p<sub>i</sub>. All other entries are zero.
#
# Use [[dense]] to convert this into its dense form.
#
# You can create a new SparseVector by entering its printable encoding as described above, e.g.:
#	$v = new SparseVector<Int>(<< '.');
#	(6) (1 1) (2 2)
#	.
# @tparam Element
declare property_type SparseVector<Element=Rational> : Vector<Element> : c++ (include => ["polymake/SparseVector.h", "polymake/SparseMatrix.h"]) {

   # The number of non-zero entries.
   # @return Int
   user_method size() : c++;
}


property_type Vector {

   # accept a SparseVector of a different type as a value of Vector property but substitute the desired element type
   type_method coherent_type {
      my ($self, $value)=@_;
      instanceof SparseVector($value) ? typeof SparseVector($self->params->[0]) : undef;
   }
}

# @category Algebraic Types
# A SparseMatrix is a two-dimensional associative array with row and column indices as keys; elements equal to the default value (ElementType(), which is 0 for most numerical types) are not stored, but implicitly encoded by the gaps in the key set. Each row and column is organized as an AVL-tree.
# 
# Use [[dense]] to convert this into its dense form.
#
# You can create a new SparseMatrix by entering its entries row by row, as a list of [[SparseVector|SparseVectors]] e.g.:
#	$A = new SparseMatrix<Int>(<< '.');
#	(5) (1 1)
#	(5) (4 2)
#	(5)
#	(5) (0 3) (1 -1)
#	.
# @tparam Element
# @tparam Sym one of [[Symmetric]] or [[NonSymmetric]], default: [[NonSymmetric]]
declare property_type SparseMatrix<Element=Rational, Sym=NonSymmetric> : Matrix<Element,Sym> : c++ (include => ["polymake/SparseMatrix.h"]) {

   # Removes empty rows and columns.
   # The remaining rows and columns are renumbered without gaps.
   user_method squeeze() : non_const : void : c++;

   # Removes empty rows.
   # The remaining rows are renumbered without gaps.
   user_method squeeze_rows() : non_const : void : c++;

   # Removes empty columns.
   # The remaining columns are renumbered without gaps.
   user_method squeeze_cols() : non_const : void : c++;

   # Resize the matrix
   user_method resize($$) : non_const : void : c++;
}


property_type Matrix {

   # accept a SparseMatrix of a different type as a value of Matrix property but substitute the desired element type
   type_method coherent_type {
      my ($self, $value)=@_;
      instanceof SparseMatrix($value) ? typeof SparseMatrix(@{$self->params}) : undef;
   }
}


function entire(SparseVector:anchor) : c++ : returns(SparseIterator);

# @category Data Conversion
# Convert to an equivalent dense vector of the same element type.
# @tparam Element
# @param SparseVector<Element> v
# @return Vector<Element>
user_function dense<Element>(SparseVector<Element>) { new Vector<Element>(shift) }

# @category Data Conversion
# Convert to an equivalent dense matrix of the same element type.
# @tparam Element
# @param SparseMatrix<Element> m
# @return Matrix<Element>
user_function dense<Element>(SparseMatrix<Element,_>) { new Matrix<Element>(shift) }

# @category Data Conversion
# @param Set S
# @tparam Scalar
user_function toVector<Scalar>(Set:wary:anchor $) : c++ (name => 'same_element_sparse_vector', include => "polymake/SparseVector.h");

# @category Data Conversion
# Convert an IncidenceMatrix to a SparseMatrix.
# @param IncidenceMatrix A
# @tparam Scalar
# @return SparseMatrix<Scalar>
# @example > $M = toMatrix<Int>(cube(2)->VERTICES_IN_FACETS);
# > print $M->type->full_name;
# | SparseMatrix<Int, NonSymmetric>

user_function toMatrix<Scalar>(IncidenceMatrix:anchor) : c++ (name => 'same_element_sparse_matrix', include => "polymake/SparseMatrix.h");

# @category Data Conversion
# Convert to a dense 0/1 matrix.
# @param IncidenceMatrix m
# @return Matrix<Int>
user_function dense(IncidenceMatrix) { dense(toMatrix<Int>(@_)); }

# @category Data Conversion
# Convert to a dense 0/1 vector of a given dimension.
# @param Set s
# @param Int dim
# @return Vector<Int>
user_function dense(Set $) { dense(toVector<Int>(@_)); }

# @category Data Conversion
# Get the positions of non-zero entries of a sparse vector.
# @param SparseVector v
# @return Set<Int>
# @example > $v = new SparseVector(0,1,1,0,0,0,2,0,3);
# > print indices($v);
# | {1 2 6 8}

user_function indices(SparseVector:anchor) : c++ (include => "polymake/Set.h");

# @category Data Conversion
# Get the positions of non-zero entries of a vector.
# @param Vector v
# @return Set<Int>
# @example > print support(new Vector(0,23,0,0,23,0,23,0,0,23));
# | {1 4 6 9}

user_function support(Vector:anchor) : c++ (include => "polymake/linalg.h");

# @category Data Conversion
# Get the positions of non-zero entries of a sparse matrix.
# @param SparseMatrix m
# @return IncidenceMatrix
# @example > $S = new SparseMatrix([1,2,0,0,0,0],[0,0,5,0,0,32]);
# > print index_matrix($S);
# | {0 1}
# | {2 5}

user_function index_matrix(SparseMatrix:anchor) : c++ (include => "polymake/IncidenceMatrix.h");


##################################################################################

# types for tropical addition

# @category Arithmetic
# tropical addition: min
declare property_type Min : c++ (special => 'Min', include => "polymake/TropicalNumber.h") {

    user_method orientation() { return 1; }
    
    user_method apply {my ($o,$x,$y) = @_; return min($x,$y);}
}

# @category Arithmetic
# tropical addition: max
declare property_type Max : c++ (special => 'Max', include => "polymake/TropicalNumber.h") {

    user_method orientation() { return -1; }
    
    user_method apply {my ($o,$x,$y) = @_; return max($x,$y);}
}

# @category Arithmetic
# @tparam Addition
# @tparam Scalar default: [[Rational]]
declare property_type TropicalNumber<Addition, Scalar=Rational> : upgrades( Scalar ) : c++ \
  (operators => '++ + += * *= / /= neg - @compare', include => "polymake/TropicalNumber.h") {
  
    # The orientation of the associated addition, i.e. 
    # +1 if the corresponding 0 is +inf
    # -1 if the corresponding 0 is -inf
    # @return Int 
    user_method orientation() {
      return Addition->orientation();
    }
    
    # The zero element of the tropical semiring of this element.
    # @return Scalar
    user_method zero() : static : c++ (include => "polymake/TropicalNumber.h");
    
  }; 
  
# function is_tropical_addition(Min) {1}
# function is_tropical_addition(Max) {1}
# function is_tropical_addition() {0}
  
##################################################################################

# @category Algebraic Types
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type Ring<Coefficient=Rational, Exponent=Int> \
   : c++ (include => "polymake/Ring.h", operators => '@eq', default_constructor => 'Deserializing') {

   # number of variables and the name stem
   method construct(Int, String) : c++;

   # list of names
   method construct(String+) : c++;

   # ring with polynomial coefficients
   method construct(Ring, Int, String) : c++;

   method construct(Ring, String+) : c++;

   user_method variables() : returns(@) : c++ (include => "polymake/Polynomial.h");

   user_method variable() : c++ (include => "polymake/Polynomial.h");

   user_method n_vars() : c++ (include => "polymake/Polynomial.h");

   user_method names() : c++;

   method id() : c++;

   type_method toXML { &Core::CPlusPlus::serialized_with_id_toXML }
}

# @category Algebraic Types
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type Monomial<Coefficient=Rational, Exponent=Int> \
   : c++ (operators => '@arith_nodiv / ^ ^= @compare', include => "polymake/Polynomial.h", default_constructor => 'Deserializing') {

   method construct(Ring) : c++;

   method construct(Vector,Ring) : c++;
}

# @category Algebraic Types
# A class for __univariate__ monomials.
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type UniMonomial<Coefficient=Rational, Exponent=Int> \
   : c++ (operators => '@arith_nodiv / ^ ^= @compare', include => "polymake/RationalFunction.h") {

   method construct(Ring) : c++;

   method construct(*,Ring) : c++;
}

# @category Algebraic Types
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type Term<Coefficient=Rational, Exponent=Int> \
   : upgrades( Coefficient, Monomial<Coefficient, Exponent> ) \
   : c++ (operators => '@arith @compare', include => "polymake/Polynomial.h", default_constructor => 'Deserializing') {

   method construct(Ring) : c++;

   method construct(Vector,*,Ring) : c++;
}

# @category Algebraic Types
# A class for __univariate__ terms.
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type UniTerm<Coefficient=Rational, Exponent=Int> \
   : upgrades( Coefficient, UniMonomial<Coefficient, Exponent> ) \
   : c++ (operators => '@arith @compare', include => "polymake/RationalFunction.h") {

   method construct(Ring) : c++;

   method construct(*,Ring) : c++;
}

# @category Algebraic Types
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type Polynomial<Coefficient=Rational, Exponent=Int> \
   : upgrades( Term<Coefficient, Exponent> ) \
   : c++ (operators => '@arith @compare', include => "polymake/Polynomial.h", default_constructor => 'Deserializing') {

   method construct(Ring) : c++;

   method construct(*,Ring) : c++;

   method construct(Matrix,*,Ring) : c++;

   # The exponent of the leading monomial.
   # @return Int
   user_method lm_exp() : c++;

   # The __leading coefficient__.
   # @return Int
   user_method lc() : c++;

   user_method print_ordered(Matrix) : c++ : void;

   user_method monomials_as_matrix() : c++;

   user_method get_ring() : c++;

   user_method trivial() : c++;

   user_method coefficients_as_vector() : c++;

}

# @category Algebraic Types
# A class for __univariate__ polynomials.
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type UniPolynomial<Coefficient=Rational, Exponent=Int> \
   : upgrades( UniTerm<Coefficient, Exponent> ) \
   : c++ (operators => '@arith % %= @compare', include => "polymake/RationalFunction.h") {

   method construct(Ring) : c++;

   method construct(*,Ring) : c++;

   method construct(*,*,Ring) : c++;

   user_method deg() : c++;

   user_method lower_deg() : c++;

   # The __leading coefficient__.
   # @return Int
   user_method lc() : c++;

   user_method print_ordered(Matrix) : c++ : void;

   user_method evaluate_float($) : c++;

   user_method evaluate(*;$=1) : c++;
}

# @category Arithmetic
# Returns the __greatest common divisor__ of two univariate polynomials.
# @param UniPolynomial p
# @param UniPolynomial q
# @return UniPolynomial
# @example We first create the polynomial ring with one variable, Rational coefficients and Int exponents.
# $r=new Ring<Rational, Int>(1);
# Then we create two UniPolynomials with said coefficient and exponent type:
# > $p = new UniPolynomial<Rational,Int>([2,2],[3,2],$r);
# > $q = new UniPolynomial<Rational,Int>([6,4],[4,2],$r);
# Printing them reveals what the constructor does:
# > print $p;
# | 2*x^3 + 2*x^2
# > print $q;
# | 4*x^4 + 6*x^2
# Now we can calculate their gcd:
# > print gcd($p,$q);
# | x^2

user_function gcd(UniPolynomial, UniPolynomial) : c++;


# @category Algebraic Types
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Int]]
declare property_type RationalFunction<Coefficient=Rational, Exponent=Int> \
   : upgrades( UniPolynomial<Coefficient, Exponent> ) \
   : c++ (operators => '@arith @eq', include => "polymake/RationalFunction.h") {

   method construct(Ring) : c++;

   method construct(*,*) : c++;

   method construct(*,*,Ring) : c++;
}

# @category Arithmetic
# Returns the __numerator__ of a [[RationalFunction]] //f//.
# @param RationalFunction f
# @return Polynomial
user_function numerator(RationalFunction) : c++;

# @category Arithmetic
# Returns the __denominator__ of a [[RationalFunction]] //f//.
# @param RationalFunction f
# @return Polynomial
user_function denominator(RationalFunction) : c++;

# @category Algebraic Types
# @tparam MinMax type of tropical addition: either [[Min]] or [[Max]]
# @tparam Coefficient default: [[Rational]]
# @tparam Exponent default: [[Rational]]
declare property_type PuiseuxFraction<MinMax, Coefficient=Rational, Exponent=Rational> \
   : upgrades( Coefficient ) \
   : c++ (operators => '@arith @compare', include => "polymake/PuiseuxFraction.h") {

   method construct(type_upgrades_to< RationalFunction<Coefficient,Exponent> >) : c++;

   method construct(type_upgrades_to< UniPolynomial<Coefficient,Exponent> >,type_upgrades_to< UniPolynomial<Coefficient,Exponent> >) : c++;

   # The __valuation__.
   # @return TropicalNumber<MinMax>
   user_method val() : c++;

   # Approximate evaluation at //x//
   # @param Float x
   # @return Float
   user_method evaluate_float($) : c++;

   # Approximate evaluation of a [[Matrix]] at //x//
   # @param Matrix m
   # @param Float x
   # @return Float
   user_method evaluate_float(Matrix,*) : static : c++;

   # Approximate evaluation of a [[Vector]] at //x//
   # @param Vector v
   # @param Float x
   # @return Float
   user_method evaluate_float(Vector,*) : static : c++;

   # Evaluate all [[PuiseuxFraction]]s in a [[Matrix]] at a [[Rational]] number (//x^exp//).
   # Let //explcm// be the lcm of the denominators of all exponents.
   # If there are no denominators or //explcm// divides //exp//, then the evaluation
   # is computed exactly.
   # Otherwise, some rational number close to the root //(x^exp)^-explcm// will be chosen
   # via an intermediate floating point number.
   # @param Matrix m
   # @param Coefficient x
   # @param Int exp (default: 1)
   # @return Matrix<Coefficient>
   user_method evaluate(Matrix,*;$=1) : static : c++;

   # Evaluate all [[PuiseuxFraction]]s in a [[Vector]] at a [[Rational]] number (//x^exp//).
   # Let //explcm// be the lcm of the denominators of all exponents.
   # If there are no denominators or //explcm// divides //exp//, then the evaluation
   # is computed exactly.
   # Otherwise, some rational number close to the root //(x^exp)^-explcm// will be chosen
   # via an intermediate floating point number.
   # @param Vector v
   # @param Coefficient x
   # @param Int exp (default: 1)
   # @return Vector<Coefficient>
   user_method evaluate(Vector,*;$=1) : static : c++;

   # Evaluate a [[PuiseuxFraction]] at a [[Rational]] number (//x^exp//).
   # Let //explcm// be the lcm of the denominators of all exponents.
   # If there are no denominators or //explcm// divides //exp//, then the evaluation
   # is computed exactly.
   # Otherwise, some rational number close to the root //(x^exp)^-explcm// will be chosen
   # via an intermediate floating point number.
   # @param Coefficient x
   # @param Int exp (default: 1)
   # @return Coefficient
   user_method evaluate(*;$=1) : c++;
}

function is_ordered_field_with_unlimited_precision(PuiseuxFraction) { 1 }

# @category Arithmetic
# Returns the __numerator__ of a [[PuiseuxFraction]] //f//.
# @param PuiseuxFraction f
# @return Polynomial
user_function numerator(PuiseuxFraction) : c++;

# @category Arithmetic
# Returns the __denominator__ of a [[PuiseuxFraction]] //f//.
# @param PuiseuxFraction f
# @return Polynomial
user_function denominator(PuiseuxFraction) : c++;


# @topic category property_types/Linear Algebra
# These types are needed as return types of algebraic computations.

# @category Linear Algebra
# Complete result of the __Smith normal form__ computation.
# Contains the following fields:
#	SparseMatrix<Scalar> form: the Smith normal form S of the given matrix //M//
#	List<Pair<Scalar, Int>> torsion:  absolute values of the entries greater than 1 of the diagonal together with their multiplicity
#	Int rank: rank of //M//
#	SparseMatrix<Scalar> left_companion, right_companion: unimodular matrices L and R such that
#                        M = LSR in normal case, or S = LMR in inverted case (as specified in the call to [[smith_normal_form]] function).
declare property_type SmithNormalForm<Scalar> : c++ (include => "polymake/Smith_normal_form.h");

# @category Linear Algebra
# Compute the __Smith normal form__ of a given matrix //M//.
# @param Matrix M must be of integer type
# @param Bool inv optional, if true, compute the inverse of the companion matrices
# @return SmithNormalForm<Integer>
# @example > $M = new Matrix<Integer>([1,2],[23,24]);
# > $SNF = smith_normal_form($M);
# The following line prints the three matrices seperated by newline characters.
# > print $SNF->left_companion ,"\n", $SNF->form ,"\n", $SNF->right_companion;
# | 1 0
# | 23 1
# |
# | 1 0
# | 0 -22
# |
# | 1 2
# | 0 1

user_function smith_normal_form(Matrix; $=0) : c++ (include => "polymake/Smith_normal_form.h");


# Local Variables:
# mode: perl
# cperl-indent-level: 3
# indent-tabs-mode:nil
# End: