SPAR/AIModule/BWTA/vendors/CGAL/CGAL/Polynomial/internal/Polynomial_impl.h

Go to the documentation of this file.

// Copyright (c) 2005  Stanford University (USA).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 of the License.
// See the file LICENSE.LGPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Kinetic_data_structures/include/CGAL/Polynomial/internal/Polynomial_impl.h $
// $Id: Polynomial_impl.h 41418 2008-01-03 14:50:13Z spion $
// 
//
// Author(s)     : Daniel Russel <drussel@alumni.princeton.edu>

#ifndef CGAL_POLYNOMIAL_INTERNAL_POLYNOMIAL_CORE_H_
#define CGAL_POLYNOMIAL_INTERNAL_POLYNOMIAL_CORE_H_

#include <CGAL/Polynomial/basic.h>
#include <vector>
#include <iostream>
#include <CGAL/Polynomial/internal/Rational/Evaluate_polynomial.h>
#include <sstream>
#include <algorithm>

#define CGAL_EXCESSIVE(x)

CGAL_POLYNOMIAL_BEGIN_INTERNAL_NAMESPACE

template <class This, class NT_t>
class Polynomial_impl;

template < class T, class NT, class C, class Tr>
inline std::ostream &operator<<(std::basic_ostream<C,Tr> &out,
                                const Polynomial_impl<T, NT> &poly);




template <class This, class NT_t>
class Polynomial_impl
{
  typedef std::vector<NT_t>                  Coefficients;
public:
  typedef NT_t                             NT;
  typedef typename Coefficients::const_iterator  iterator;
  typedef NT                               result_type;
  typedef NT                               argument_type;

  //================
  // CONSTRUCTORS
  //================

  Polynomial_impl() {}

  Polynomial_impl(const NT& c) {
    coefs_.push_back(c);
  }

  template <class Iterator>
  Polynomial_impl(Iterator first, Iterator beyond): coefs_(first, beyond) {}

  //========================
  // ACCESS TO COEFFICIENTS
  //========================


  const NT& operator[](unsigned int i) const
  {
    CGAL_assertion( i < coefs_.size());
    //if (i < coefs_.size()) {
      return coefs_[i];
      /*}
    else {
      return zero_coef();
      }*/
  }

  //=====================================
  //     ITERATORS FOR COEFFICIENTS
  //=====================================

  iterator begin() const
  {
    return coefs_.begin();
  }

  iterator end() const
  {
    return coefs_.end();
  }

  //=========
  // DEGREE
  //=========

  int degree() const
  {
    return static_cast<int>(coefs_.size()) - 1;
  }

  //=============
  // PREDICATES
  //=============

  bool is_constant() const
  {
    return degree() < 1;
  }


  bool is_zero() const
  {
    CGAL_Polynomial_assertion( coefs_.empty() == (coefs_.size()==0) );
    return coefs_.empty();
  }

  //=======================
  //     OPERATORS
  //=======================

  This operator-() const
  {
    This ret;
    ret.coefs_.resize( coefs_.size() );
    for (unsigned int i=0; i < coefs_.size(); ++i) {
      ret.coefs_[i] = -coefs_[i];
    }
    return ret;
  }

  This operator+(const This &o) const
  {
    if (is_zero()) { return o; }
    else if (o.is_zero()) { return This(*this); }
    else {
      This ret;
      unsigned int new_deg = (std::max)(o.degree(), degree());
      ret.coefs_.resize(new_deg + 1);
      unsigned int md= (std::min)(degree(), o.degree());
      for (unsigned int i = 0; i <= md; ++i) {
        ret.coefs_[i] = operator[](i) + o[i];
      }
      for (int i=md+1; i <= degree(); ++i){
        ret.coefs_[i]= operator[](i);
      }
      for (int i=md+1; i <= o.degree(); ++i){
        ret.coefs_[i]= o[i];
      }
      CGAL_EXCESSIVE(std::cout << *this << " - " << o << " + " << ret << std::endl);
      ret.finalize();
      return ret;
    }
  }

  This operator-(const This &o) const
  {
    if (is_zero()) { return -o; }
    else if (o.is_zero()) { return This(*this); }
    else {
      This ret;
      unsigned int new_deg = (std::max)(o.degree(), degree());
      ret.coefs_.resize( new_deg + 1 );
      unsigned int md= (std::min)(degree(), o.degree());
      for (unsigned int i = 0; i <= md; ++i) {
        ret.coefs_[i] = operator[](i) - o[i];
      }
      for (int i=md+1; i <= degree(); ++i){
        ret.coefs_[i]= operator[](i);
      }
      for (int i=md+1; i <= o.degree(); ++i){
        ret.coefs_[i]= -o[i];
      }
      CGAL_EXCESSIVE(std::cout << *this << " - " << o << " = " << ret << std::endl);
      ret.finalize();
      return ret;
    }
  }

  This operator*(const This &o) const
  {
    if (is_zero()) return This(NT(0));
    else if (o.is_zero()) return o;
    This ret;
    ret.coefs_.resize( degree() + o.degree() + 1 );
    // the following for-loop makes sure that the values on ret.coefs_
    // are properly initialized with zero.
    for (int i = 0; i <= degree() + o.degree(); ++i) {
      ret.coefs_[i] = NT(0);
    }
    for (unsigned int i = 0; i < coefs_.size(); ++i) {
      for (unsigned int j = 0; j < o.coefs_.size(); ++j) {
        NT prev = ret.coefs_[i+j];
        NT result = prev + operator[](i) * o[j];
        ret.coefs_[i+j] = result;
      }
    }
    CGAL_EXCESSIVE(std::cout << *this << " * " << o << " = " << ret << std::endl);
    return ret;
  }

  This operator+(const NT& a) const
  {
    if ( is_zero() ) { return This(a); }

    This res(*this);
    res.coefs_[0] += a;
    CGAL_EXCESSIVE(std::cout << *this << " + " << a << " = " << res << std::endl);
    return res;
  }

  This operator-(const NT& a) const
  {
    return (*this) + (-a);
  }

  This operator*(const NT& a) const
  {
    if ( is_zero() || CGAL_NTS sign(a)==ZERO ) { return This(); }

    This res;
    unsigned int deg = degree();
    res.coefs_.resize(deg + 1, NT(0));
    for (unsigned int i = 0; i <= deg; i++) {
      res.coefs_[i] = coefs_[i] * a;
    }
    return res;
  }

  This operator/(const NT& a) const
  {
    NT inv_a = NT(1) / a;
    return (*this) * inv_a;
  }

  //=======================
  // VALUE OF POLYNOMIAL
  //=======================


  NT operator()(const NT &t) const
  {
    if (is_zero()) return NT(0);
    else if (degree()==0 /* || t==NT(0)*/ ) {
      return operator[](0);
    }
    else {
      return evaluate_polynomial(coefs_, t);
    }
  }



  NT value_at(const NT &t) const
  {
    return operator()(t);
  }


  bool operator!=(const This &o) const
  {
    return !operator==(o);
  }


  bool operator==(const This &o) const
  {
    if (degree() != o.degree()) return false;
    int max_size = (std::max)(o.coefs_.size(), coefs_.size());
    for (int i = 0; i < max_size; ++i) {
      if (o[i] != operator[](i)) return false;
    }
    return true;
  }


  void print()const
  {
    write(std::cout);
  }


  template <class charT, class Traits>
  void read(std::basic_istream<charT, Traits> &in) {
    bool pos=(in.peek()!='-');
    if (in.peek() == '+' || in.peek() == '-') {
      char c;
      in >> c;
    }
    char v='\0';
    while (true) {
      char vc, t;
      NT coef;
      // eat
      in >> coef;
      if (in.fail()) {
        //std::cout << "Read ";
        //write(std::cout);
        //std::cout << std::endl;
        return;
      }
      //std::cout << "Read " << coef << std::endl;

      //NT cp1= coef+NT(1);
      unsigned int pow;
      char p= in.peek();
      if (in.peek() =='*') {
        in >> t >> vc;
        if (t != '*') {
          in.setstate(std::ios_base::failbit);
          return;
          //return in;
        }
        if ( !(v=='\0' || v== vc)) {
          in.setstate(std::ios_base::failbit);
          return;
          //return in;
        }
        v=vc;
        p=in.peek();
        if (in.peek() =='^') {
          char c;
          in  >> c >> pow;
        }
        else {
          pow=1;
        }
      }
      else {
        pow=0;
      }

      if (coefs_.size() <=pow) {
        coefs_.resize(pow+1, NT(0));
      }

      if (!pos) coef=-coef;
      //std::cout << "Read 2 " << coef << "*t^" << pow << std::endl;
      coefs_[pow]=coef;
      //std::cout << "Read " << coefs_[pow] << std::endl;

      char n= in.peek();
      if (n=='+' || n=='-') {
        pos= (n!='-');
        char e;
        in >> e;
      } else {
        /*bool f= in.fail();
          bool g= in.good();
          bool e= in.eof();
          bool b= in.bad();*/
        // This is necessary since peek could have failed if we hit the end of the buffer
        // it would better to do without peek, but that is too messy
        in.clear();
        //std::cout << f << g << e << b<<std::endl;
        break;
      }
    };
  }


  template <class C, class T>
  void write(std::basic_ostream<C,T> &out) const
  {
    std::basic_ostringstream<C,T> s;
    s.flags(out.flags());
    s.imbue(out.getloc());
    s.precision(12);
    if (coefs_.size()==0) {
      s << "0";
    }
    else {
      for (unsigned int i=0; i< coefs_.size(); ++i) {
        if (i==0) {
          if (coefs_[i] != 0) s << coefs_[i];
        }
        else {
          if ( coefs_[i] != 0 ) {
            if (coefs_[i] >0) s << "+";
            s << coefs_[i] << "*t";
            if (i != 1) {
              s << "^" << i;
            }
          }
        }
      }
    }

    out << s.str();
  }


  /*
    typedef typename Coefficients::const_iterator Coefficients_iterator;
    Coefficients_iterator coefficients_begin() const {
    return coefs_.begin();

    }

    Coefficients_iterator coefficients_end() const {
    return coefs_.end();
    }*/

protected:

  std::size_t size() const {return coefs_.size();}

  static const NT & zero_coef() {
    static NT z(0);
    return z;
  }


  Coefficients coefs_;
};

template < class T, class NT,  class C, class Tr>
inline std::ostream &operator<<(std::basic_ostream<C, Tr> &out,
                                const Polynomial_impl<T, NT> &poly)
{
  poly.write(out);
  return out;
}


template < class T, class NT, class C, class Tr>
inline std::istream &operator>>(std::basic_istream<C, Tr> &in,
                                Polynomial_impl<T, NT> &poly)
{
  poly.read(in);
  //std::cout << "Read " << poly << std::endl;
  return in;
}


template <class T, class NT>
inline T operator*(const NT &a, const Polynomial_impl<T, NT> &poly)
{
  return (poly * a);
}


template <class T, class NT>
inline T operator+(const NT &a, const Polynomial_impl<T, NT> &poly)
{
  return (poly + a);
}


template < class T, class NT>
inline T  operator+(int a, const Polynomial_impl<T, NT> &poly)
{
  return (poly + NT(a));
}


template <class T, class NT>
inline T operator-(const NT &a, const Polynomial_impl<T, NT> &poly)
{
  return -(poly - a);
}


template <class T, class NT>
inline T operator-(int a, const Polynomial_impl<T, NT> &poly)
{
  return -(poly - NT(a));
}


#undef CGAL_EXCESSIVE

/*template <class NT>
class Input_rep<Polynomial<NT> , CGAL::Maple_format_tag > {
  Input_rep(Polynomial<NT>& p): p_(p){}
  std::istream & operator()(std::istream &in) const {
    
  }
  };*/

CGAL_POLYNOMIAL_END_INTERNAL_NAMESPACE
#endif