/*******************************************************************************
+
+  LEDA 3.5
+
+  _rat_vector.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
#include <LEDA/rat_vector.h>
#include <LEDA/integer_matrix.h>

rat_vector::rat_vector(int d) 
{ PTR= new geo_rep(d);
  for(int i=0; i<d; i++) ptr()->v[d] = 0;
  ptr()->v[d] = 1;
}

rat_vector::rat_vector(integer a,integer b,integer D)
{ if (D==0)
    error_handler(1,"rat_vector::constructor: D must be nonzero.");
  if (D<0)
    PTR= new geo_rep(-a,-b,-D);
  else
    PTR= new geo_rep(a,b,D);
}

rat_vector::rat_vector(integer a,integer b,integer c,integer D)
{ if (D==0)
    error_handler(1,"rat_vector::constructor: D must be nonzero.");
  if (D<0)
    PTR= new geo_rep(-a,-b,-c,-D);
  else
    PTR= new geo_rep(a,b,c,D);
}

rat_vector::rat_vector(const integer_vector&c,integer D)
{ if (D==0)
    error_handler(1,"rat_vector::constructor: D must be nonzero.");
  if (D<0)
    PTR= new geo_rep(-c,-D);
  else
    PTR= new geo_rep(c,D);
}

rat_vector::rat_vector(const integer_vector&c)
{ int d = c.dim();
  integer D = c[d];
  if (D==0)
    error_handler(1,"rat_point::constructor: D must be nonzero");
  if (D<0)
    PTR = new geo_rep(-c);
  else
    PTR= new geo_rep(c);
}




// sn (02/1997)
rational  rat_vector::sqr_length() const
{ int      d = ptr()->dim;
  integer* V = ptr()->v;
  integer nom = 0;
  for(int i=0; i<d; i++) nom += V[i]*V[i];
  integer den = V[d]*V[d];
  return rational(nom,den);
}

// sn (02/1997)
vector rat_vector::to_vector() const
{ int      d = ptr()->dim;
  integer* V = ptr()->v;
  double   w = V[d].to_double();
  vector vec(d);
  for(int i=0; i<d; i++) vec[i] = V[i].to_double()/w;
  return vec;
}


rat_vector rat_vector::rotate90() const
{ if (ptr()->dim != 2)
      error_handler(1,"rat_vector::rotate90: dimension must be two. ");
  return rat_vector(-hcoord(1),hcoord(0),hcoord(2)); 
 }



rat_vector rat_vector::d2(integer a, integer b, integer D)
{
  if (D == 0)
    error_handler(1, "rat_vector::d2: denominator must not be zero.");
  return rat_vector(a, b, D);
}

rat_vector rat_vector::d3(integer a, integer b, integer c, integer D)
{
  rat_vector d3v(3);
  if (D == 0)
    error_handler(1, "rat_vector::d3: denominator must not be zero.");
  if (D < 0)
    d3v.ptr()->init4(-a, -b, -c, -D);
  else
    d3v.ptr()->init4(a, b, c, D);
  return d3v;
}

rat_vector rat_vector::unit(int i, int d)
{
  if (i < 0 || i >= d)
    error_handler(1, "rat_vector::unit: i out of range.");
  rat_vector uv(d);
  uv.ptr()->v[i] = 1;
  return uv;
}

rat_vector rat_vector::zero(int d)
{
  return rat_vector(d);
}





void rat_vector::print(ostream & out) const
{
  out << ptr();
}

void rat_vector::read(istream & in)
{
  int d = dim();
  if (refs() > 1)
    operator = (rat_vector(d));
  in >> ptr();
  if (hcoord(d) == 0)
    error_handler(1, "operator>>: denominator of vector must be nonzero.");
  if (hcoord(d) < 0)
    ptr()->negate(d + 1);
}






rat_vector rat_vector::scale(integer m, integer n) const
{
  int d = dim();
  rat_vector result(d);
  result.ptr()->v[d] = ptr()->v[d] * n;
  integer g = gcd(ptr()->v[d], m);
  result.ptr()->v[d] /= g;
  m /= g;
  for (int i = 0; i < d; i++)
    result.ptr()->v[i] = ptr()->v[i] * m;
  return result;
}

void rat_vector::self_scale(integer m, integer n)
{
  int d = dim();
  ptr()->v[d] *= n;
  integer g = gcd(ptr()->v[d], m);
  ptr()->v[d] /= g;
  m /= g;
  for (int i = 0; i < d; i++)
    ptr()->v[i] *= m;
}

rat_vector operator *(int n, const rat_vector & p) {
  return p.scale(n, 1);
}

rat_vector operator *(integer n, const rat_vector & p) {
  return p.scale(n, 1);
}

rat_vector operator *(rational r, const rat_vector & p) {
  return p.scale(r.numerator(), r.denominator());
}


rat_vector operator / (const rat_vector & p, int n) {
  return p.scale(1, n);
}

rat_vector operator / (const rat_vector & p, integer n) {
  return p.scale(1, n);
}

rat_vector operator / (const rat_vector & p, rational r) {
  return p.scale(r.denominator(), r.numerator());
}
rat_vector & rat_vector::operator *= (integer n) {
  self_scale(n, 1);
  return *this;
}
rat_vector & rat_vector::operator *= (int n) {
  self_scale(n, 1);
  return *this;
}
rat_vector & rat_vector::operator *= (rational r) {
  self_scale(r.numerator(), r.denominator());
  return *this;
}
rat_vector & rat_vector::operator /= (integer n) {
  self_scale(1, n);
  return *this;
}
rat_vector & rat_vector::operator /= (int n) {
  self_scale(1, n);
  return *this;
}
rat_vector & rat_vector::operator /= (rational r) {
  self_scale(r.denominator(), r.numerator());
  return *this;
}


rational operator *(const rat_vector v, const rat_vector & w) {
  int d = v.dim();
  if (d != w.dim())
    error_handler(1, "inner product: dimensions disagree.");
  integer nom = 0;
  for (int i = 0; i < d; i++)
    nom += v.hcoord(i) * w.hcoord(i);
  integer denom = v.hcoord(d) * w.hcoord(d);
  return rational(nom, denom);
}




rat_vector operator + (const rat_vector & v, const rat_vector & w) {
  rat_vector res(v.dim());
  c_add(res.ptr(), v.ptr(), w.ptr());
  return res;
}
rat_vector & rat_vector::operator += (const rat_vector & w) {
  int d = dim();
  rat_vector old(*this);
  if (ptr()->count > 2)
    *this = rat_vector(d);
  c_add(ptr(), old.ptr(), w.ptr());
  return *this;
}

rat_vector operator - (const rat_vector & v, const rat_vector & w) {
  rat_vector res(v.dim());
  c_sub(res.ptr(), v.ptr(), w.ptr());
  return res;
}
rat_vector & rat_vector::operator -= (const rat_vector & w) {
  int d = dim();
  rat_vector old(*this);
  if (ptr()->count > 2)
    *this = rat_vector(d);
  c_sub(ptr(), old.ptr(), w.ptr());
  return *this;
}

rat_vector rat_vector::operator - () const
{
  int d = dim();
  rat_vector result(d);
  result.ptr()->copy(ptr());
  result.ptr()->negate(d);
  return result;
}





bool contained_in_linear_hull(const array <rat_vector > &A,
   const rat_vector & x)
{
  int al = A.low();
  int k = A.high() - al + 1;	/* |A| contains |k| vectors */
  int d = A[al].dim();
  integer_matrix M(d, k);
  integer_vector b(d);
  for (int i = 0; i < d; i++) {
    b[i] = x.hcoord(i);
    for (int j = 0; j < k; j++)
      M(i, j) = A[al + j].hcoord(i);
  }
  return is_solvable(M, b);
}



int linear_rank(const array <rat_vector > &A)
{
  int al = A.low();
  int k = A.high() - al + 1;	/* |A| contains |k| vectors */
  int d = A[al].dim();
  integer_matrix M(d, k);
  for (int i = 0; i < d; i++)
    for (int j = 0; j < k; j++)
      M(i, j) = A[al + j].hcoord(i);
  return rank(M);
}




bool linearly_independent(const array <rat_vector > &A)
{
  return (linear_rank(A) == A.high() - A.low() + 1);
}




array <rat_vector > linear_base(const array <rat_vector > &A)
{
  int al = A.low();
  int k = A.high() - al + 1;	/* |A| contains |k| vectors */
  int d = A[al].dim();
  integer denom;
  integer_matrix M(d, k);
  for (int j = 0; j < k; j++)
  { for (int i = 0; i < d; i++)
      M(i, j) = A[al + j].hcoord(i);
   }

  array <int >indcols;
  independent_columns(M, indcols);

  int indcolsdim = indcols.high() + 1;
  array <rat_vector > L(indcolsdim);
  for (int i = 0; i < indcolsdim; i++)
    L[i] = rat_vector(M.col(indcols[i]), 1);
  return L;
}
