poly.tcc

很多二维 三维几何计算算法 C++ 类库

/*
 * Core Library Version 1.7, August 2004
 * Copyright (c) 1995-2004 Exact Computation Project
 * All rights reserved.
 *
 * This file is part of CORE (http://cs.nyu.edu/exact/core/); you may
 * redistribute it under the terms of the Q Public License version 1.0.
 * See the file LICENSE.QPL distributed with CORE.
 *
 * 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.
 *
 *
 * File: Poly.tcc
 * Purpose: 
 *	Template implementations of the functions
 *	of the Polynomial<NT> class (found in Poly.h)
 * 
 * OVERVIEW:  
 *	Each polynomial has a nominal "degree" (this
 *		is an upper bound on the true degree, which
 *		is determined by the first non-zero coefficient).
 *	Coefficients are parametrized by some number type "NT".
 *	Coefficients are stored in the "coeff" array of
 *		length "degree + 1".  
 *	IMPORTANT CONVENTION: the zero polynomial has degree -1
 *		while nonzero constant polynomials have degree 0.
 *	
 * Bugs:
 *	Currently, coefficient number type NT only accept
 *			NT=BigInt and NT=int
 * 
 *	To see where NT=Expr will give trouble,
 *			look for NOTE_EXPR below.
 * 
 * Author: Chee Yap, Sylvain Pion and Vikram Sharma
 * Date:   May 28, 2002
 *
 * WWW URL: http://cs.nyu.edu/exact/
 * Email: [e-mail]
 *
 * $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Core/include/CGAL/CORE/poly/Poly.tcc $
 * $Id: Poly.tcc 36978 2007-03-10 10:37:52Z spion $
 ***************************************************************************/


template <class NT>
int Polynomial<NT>::COEFF_PER_LINE  = 4;           // pretty print parameters
template <class NT>
const char* Polynomial<NT>::INDENT_SPACE ="   ";  // pretty print parameters

// ==================================================
// Polynomial Constructors
// ==================================================

template <class NT>
Polynomial<NT>::Polynomial(void) {
  degree = -1; 	// this is the zero polynomial!
  coeff = NULL;
}

//Creates a polynomial with nominal degree n
template <class NT>
Polynomial<NT>::Polynomial(int n) {
  assert(n>= -1);
  degree = n;
  if (n == -1)
    return;	// return the zero polynomial!
  if (n>=0)
    coeff = new NT[n+1];
  coeff[0]=1;			// otherwise, return the unity polynomial
  for (int i=1; i<=n; i++)
    coeff[i]=0;
}

template <class NT>
Polynomial<NT>::Polynomial(int n, const NT * c) {
  //assert("array c has n+1 elements");
  degree = n;
  if (n >= 0) {
    coeff = new NT[n+1];
    for (int i = 0; i <= n; i++)
      coeff[i] = c[i];
  }
}

//
// Constructor from a vector of NT's
///////////////////////////////////////
template <class NT>
Polynomial<NT>::Polynomial(const VecNT & vN) {
  degree = vN.size()-1;
  if (degree >= 0) {
    coeff = new NT[degree+1];
    for (int i = 0; i <= degree; i++)
      coeff[i] = vN[i];
  }
}

template <class NT>
Polynomial<NT>::Polynomial(const Polynomial<NT> & p) {
  coeff = NULL;//WHY?
  *this = p;	// reduce to assignment operator=
}

// Constructor from a Character string of coefficients
///////////////////////////////////////
template <class NT>
Polynomial<NT>::Polynomial(int n, const char * s[]) {
  //assert("array s has n+1 elements");
  degree = n;
  if (n >= 0) {
    coeff = new NT[n+1];
    for (int i = 0; i <= n; i++)
      coeff[i] = s[i];
  }
}

//The BNF syntax is the following:-
//    [poly] -> [term]| [term] '+/-' [poly] |
//    		'-' [term] | '-' [term] '+/-' [poly]
//    [term] -> [basic term] | [basic term] [term] | [basic term]*[term]
//    [basic term] -> [number] | 'x' | [basic term] '^' [number]
//                    | '(' [poly] ')' 
//COMMENT: 
//  [number] is assumed to be a BigInt; in the future, we probably
//  want to generalize this to BigFloat, etc.
//
template <class NT>
Polynomial<NT>::Polynomial(const string & s, char myX) {
   string ss(s);
   constructFromString(ss, myX);
}
template <class NT>
Polynomial<NT>::Polynomial(const char * s, char myX) {
   string ss(s);
   constructFromString(ss, myX);
}
template <class NT>
void Polynomial<NT>::constructFromString(string & s, char myX) {
  if(myX != 'x' || myX != 'X'){
    //Replace myX with 'x'.
    unsigned int loc = s.find(myX, 0);
    while(loc != string::npos){
      s.replace(loc,1,1,'x');
      loc = s.find(myX, loc+1);
    }
  }

  coeff = NULL;//Did this to ape the constructor from polynomial above
  *this = getpoly(s);
}

template <class NT>
Polynomial<NT>::~Polynomial() {
  if(degree >= 0)
    delete[] coeff;
}

  ////////////////////////////////////////////////////////
  // METHODS USED BY STRING CONSTRUCTOR ABOVE
  ////////////////////////////////////////////////////////


//Sets the input Polynomial to X^n
template <class NT>
void Polynomial<NT>::constructX(int n, Polynomial<NT>& P){
  Polynomial<NT> q(n);//Nominal degree n
  q.setCoeff(n,NT(1));
  if (n>0) q.setCoeff(0,NT(0));
  P = q;
}

//Returns in P the coeffecient starting from start
template <class NT>
int Polynomial<NT>::getnumber(const char* c, int start, unsigned int len,
			  Polynomial<NT> & P){
  int j=0;
  char *temp = new char[len];
  while(isint(c[j+start])){
    temp[j]=c[j+start];j++;
  }
  temp[j] = '\0';
  NT cf = NT(temp);
  Polynomial<NT> q(0);
  q.setCoeff(0, cf);
  P = q;
  delete[] temp;
  return (j-1+start);//Exactly the length of the number
}

template <class NT>
bool Polynomial<NT>::isint(char c){
  if(c == '0' || c == '1' || c == '2' || c == '3' || c == '4' ||
     c == '5' || c == '6' || c == '7' || c == '8' || c == '9')
    return true;
  else
    return false;
}


//Returns as integer the number starting from start in c
template <class NT>
int Polynomial<NT>::getint(const char* c, int start, unsigned int len,
			  int & n){
  int j=0;
  char *temp = new char[len];
  while(isint(c[j+start])){
    temp[j]=c[j+start];j++;
  }
  temp[j] = '\n';
  n = atoi(temp);
  delete[] temp;
  return (j-1+start);//Exactly the length of the number
}

//Given a string starting with an open parentheses returns the place
// which marks the end of the corresponding closing parentheses.
//Strings of the form (A).
template <class NT>
int Polynomial<NT>::matchparen(const char* cstr, int start){
  int count = 0;
  int j=start;
  
  do{
    if(cstr[j] == '('){
      count++;
    }
    if(cstr[j] == ')'){
      count--;
    }
    j++;      
  }while(count != 0 );//j is one more than the matching ')'
  
  return j-1;
}


template <class NT>
int Polynomial<NT>::getbasicterm(string & s, Polynomial<NT> & P){
  const char * cstr = s.c_str();
  unsigned int len = s.length();
  int i=0;
  //Polynomial<NT> * temp = new Polynomial<NT>();

  if(isint(cstr[i])){
    i = getnumber(cstr, i, len, P);
  }else if(cstr[i] == 'x'||cstr[i] == 'X'){
    constructX(1, P);
  }else if(cstr[i] =='('){
    int oldi = i;
    i = matchparen(cstr, i);
    string t = s.substr(oldi+1, i -oldi -1);
    P = getpoly(t);
  }else{
    std::cout <<"ERROR IN PARSING BASIC TERM" << std::endl;
  }
  //i+1 points to the beginning of next syntactic object in the string.
 if(cstr[i+1] == '^'){
    int n;
    i = getint(cstr, i+2, len, n);
    P.power(n);
  }
  return i;
}


template <class NT>
int Polynomial<NT>::getterm(string & s, Polynomial<NT> & P){
  unsigned int len = s.length();
  if(len == 0){// Zero Polynomial
    P=Polynomial<NT>();
    return 0;
  }
  unsigned int ind, oind;
  const char* cstr =s.c_str();
  string t;
  //P will be used to accumulate the product of basic terms.
  ind = getbasicterm(s, P);
  while(ind != len-1 && cstr[ind + 1]!='+' && cstr[ind + 1]!='-' ){
    //Recursively get the basic terms till we reach the end or see
    // a '+' or '-' sign.
    if(cstr[ind + 1] == '*'){
      t = s.substr(ind + 2, len - ind -2);
      oind = ind + 2;
    }else{
      t = s.substr(ind + 1, len -ind -1);
      oind = ind + 1;
    }

    Polynomial<NT> R;
    ind = oind + getbasicterm(t, R);//Because the second term is the offset in
                                     //t
    P *= R;
  }

  return ind;
}

template <class NT>
Polynomial<NT> Polynomial<NT>::getpoly(string & s){

    //Remove white spaces from the string
    unsigned int cnt=s.find(' ',0);
    while(cnt != string::npos){
      s.erase(cnt, 1);
      cnt = s.find(' ', cnt);
    }

    unsigned int len = s.length();
    if(len <= 0){//Zero Polynomial
      return Polynomial<NT>();
    }

    //To handle the case when there is one '=' sign
    //Suppose s is of the form s1 = s2. Then we assign s to
    //s1 + (-1)(s2) and reset len
    unsigned int loc;
    if((loc=s.find('=',0)) != string::npos){
      s.replace(loc,1,1,'+');
      string s3 = "(-1)(";
      s.insert(loc+1, s3);
      len = s.length();
      s.insert(len, 1, ')');
    }
    len = s.length();

    const char *cstr = s.c_str();
    string t;
    Polynomial<NT> P;
    // P will be the polynomial in which we accumulate the
    //sum and difference of the different terms.
    unsigned int ind;
    if(cstr[0] == '-'){
      t = s.substr(1, len);
      ind = getterm(t,P) + 1;
      P.negate();
    }else{
      ind = getterm(s, P);
    }
    unsigned int oind =0;//the string between oind and ind is a term
    while(ind != len -1){
      Polynomial<NT> R;
      t = s.substr(ind + 2, len -ind -2);
      oind = ind;
      ind = oind + 2 + getterm(t, R);
      if(cstr[oind + 1] == '+')
		P += R;
      else if(cstr[oind + 1] == '-')
		P -= R;
      else
	std::cout << "ERROR IN PARSING POLY! " << std::endl;
    }

    return (P);
}

  ////////////////////////////////////////////////////////
  // METHODS
  ////////////////////////////////////////////////////////

// assignment:
template <class NT>
Polynomial<NT> & Polynomial<NT>::operator=(const Polynomial<NT>& p) {
  if (this == &p)
    return *this;	// self-assignment
  degree = p.getDegree();
  delete[] coeff;
  coeff = new NT[degree +1];
  for (int i = 0; i <= degree; i++)
    coeff[i] = p.coeff[i];
  return *this;
}

// getTrueDegree
template <class NT>
int Polynomial<NT>::getTrueDegree() const {
  for (int i=degree; i>=0; i--) {
    if (sign(coeff[i]) != 0)
      return i;
  }
  return -1;	// Zero polynomial
}

//get i'th Coeff. We check whether i is not greater than the
// true degree and if not then return coeff[i] o/w 0
template <class NT>
NT Polynomial<NT>::getCoeffi(int i) const {
  int deg = getTrueDegree();
  if(i > deg)
      return NT(0);
  return coeff[i];
}

// ==================================================
// polynomial arithmetic
// ==================================================

// Expands the nominal degree to n;
//	Returns n if nominal degree is changed to n
//	Else returns -2
template <class NT>
int Polynomial<NT>::expand(int n) {
  if ((n <= degree)||(n < 0))
    return -2;
  int i;
  NT * c = coeff;
  coeff = new NT[n+1];
  for (i = 0; i<= degree; i++)
    coeff[i] = c[i];
  for (i = degree+1; i<=n; i++)
    coeff[i] = 0;
  delete[] c;
  degree = n;
  return n;
}

// contract() gets rid of leading zero coefficients
//	and returns the new (true) degree;
//	It returns -2 if this is a no-op
template <class NT>
int Polynomial<NT>::contract() {
  int d = getTrueDegree();
  if (d == degree)
    return (-2);  // nothing to do
  else
    degree = d;
  NT * c = coeff;
  coeff = new NT[d+1];
  for (int i = 0; i<= d; i++)
    coeff[i] = c[i];
  delete[] c;
  return d;
}

template <class NT>
Polynomial<NT> & Polynomial<NT>::operator+=(const Polynomial<NT>& p) { // +=
  int d = p.getDegree();
  if (d > degree)
    expand(d);
  for (int i = 0; i<=d; i++)
    coeff[i] += p.coeff[i];

  return *this;
}
template <class NT>
Polynomial<NT> & Polynomial<NT>::operator-=(const Polynomial<NT>& p) { // -=
  int d = p.getDegree();
  if (d > degree)
    expand(d);
  for (int i = 0; i<=d; i++)
    coeff[i] -= p.coeff[i];
  return *this;
}

// SELF-MULTIPLICATION
// This is quadratic time multiplication!
template <class NT>
Polynomial<NT> & Polynomial<NT>::operator*=(const Polynomial<NT>& p) { // *=
  int d = degree + p.getDegree();
  NT * c = new NT[d+1];
  for (int i = 0; i<=d; i++)
    c[i] = 0;

  for (int i = 0; i<=p.getDegree(); i++)
    for (int j = 0; j<=degree; j++) {
      c[i+j] += p.coeff[i] * coeff[j];
    }
  degree = d;
  delete[] coeff;
  coeff = c;
  return *this;
}

// Multiply by a scalar
template <class NT>
Polynomial<NT> & Polynomial<NT>::mulScalar( const NT & c) {
  for (int i = 0; i<=degree ; i++)
    coeff[i] *= c;
  return *this;
}

// mulXpower: Multiply by X^i (COULD be a divide if i<0)
template <class NT>
Polynomial<NT> & Polynomial<NT>::mulXpower(int s) {
  // if s >= 0, then this is equivalent to
  // multiplying by X^s;  if s < 0, to dividing by X^s
  if (s==0)
    return *this;
  int d = s+getTrueDegree();
  if (d < 0) {
    degree = -1;
    delete[] coeff;
    coeff = NULL;
    return *this;
  }
  NT * c = new NT[d+1];
  if (s>0)
    for (int j=0;  j <= d; j++) {
      if (j <= degree)
        c[d-j] = coeff[d-s-j];
      else
        c[d-j] = 0;
    }
  if (s<0) {
    for (int j=0; j <= d; j++)
      c[d-j] = coeff[d-s-j];  // since s<0, (d-s-j) > (d-j) !!
  }
  delete[] coeff;
  coeff = c;
  degree = d;
  return *this;
}//mulXpower