curves.tcc

CGAL is a collaborative effort of several sites in Europe and Israel. The goal is to make the most i

    VecExpr vE;
    for (int i=0; i<= ydeg; i++) {
        vE.push_back(coeffX[i].eval(x));
    }
    return Polynomial<Expr>(vE);
  }//yPolynomial
    
    But this has many problems.
    Solution below is to have special yExprPolynomial(x).
    *********************************************************** */
  
template <class NT>
Polynomial<NT> BiPoly<NT>::yPolynomial(const NT & x) {
    NT coeffVec[ydeg+1];
    int d=-1;
    for(int i=ydeg; i >= 0 ; i--){
      coeffVec[i] = coeffX[i].eval(x);
      if ((d < 0) && (coeffVec[i] != 0))
	      d = i;
    }
    return Polynomial<NT>(d, coeffVec);
  }

  // Specialized version of yPolynomial for Expressions
template <class NT>
  Polynomial<Expr> BiPoly<NT>::yExprPolynomial(const Expr & x) {
    Expr coeffVec[ydeg+1];
    int d=-1;
    for(int i=ydeg; i >= 0 ; i--){
      coeffVec[i] = coeffX[i].eval(x);
      if ((d < 0) && (coeffVec[i] != 0))
	      d = i;
    }
    return Polynomial<Expr>(d, coeffVec);
  }

  // Specialized version of yPolynomial for BigFloat
  // ASSERTION: BigFloat x is exact
template <class NT>
  Polynomial<BigFloat> BiPoly<NT>::yBFPolynomial(const BigFloat & x) {
    BigFloat coeffVec[ydeg+1];
    int d=-1;
    for(int i=ydeg; i >= 0 ; i--){
      coeffVec[i] = coeffX[i].eval(x);
      if ((d < 0) && (coeffVec[i] != 0))
	      d = i;
    }
    return Polynomial<BigFloat>(d, coeffVec);
  }

  // xPolynomial(y) 
  //   returns the polynomial (in X) when we substitute Y=y
  //   
  //   N.B. May need the
  //   		Polynomial<Expr> xExprPolynomial(Expr y)
  //   version too...
  //

template <class NT>
Polynomial<NT> BiPoly<NT>::xPolynomial(const NT & y) {
    Polynomial<NT> res = coeffX[0];
    NT yPower(y);
    for(int i=1; i <= ydeg ; i++){
      res += coeffX[i].mulScalar(yPower);
      yPower *= y;
    }
    return res;
  }//xPolynomial
  

  // getYdegree()
template <class NT>
int BiPoly<NT>::getYdegree() const {
    return ydeg;
  }
  
  // getXdegree()
template <class NT>
int BiPoly<NT>::getXdegree(){
    int deg=-1;
    for(int i=0; i <=ydeg; i++)
      deg = max(deg, coeffX[i].getTrueDegree());
    return deg;
  }

  // getTrueYdegree
template <class NT>
int BiPoly<NT>::getTrueYdegree(){
    for (int i=ydeg; i>=0; i--){
      coeffX[i].contract();
      if (!zeroP(coeffX[i]))
	return i;
    }
    return -1;	// Zero polynomial
  }//getTrueYdegree


  //eval(x,y)
template <class NT>
Expr BiPoly<NT>::eval(Expr x, Expr y){//Evaluate the polynomial at (x,y)
    Expr e = 0;

    for(int i=0; i <=ydeg; i++){
      e += coeffX[i].eval(x)*pow(y,i);
    }
    return e;
  }//eval

  ////////////////////////////////////////////////////////
  // Polynomial arithmetic (these are all self-modifying)
  ////////////////////////////////////////////////////////
  
  // Expands the nominal y-degree to n;
  //	Returns n if nominal y-degree is changed to n
  //	Else returns -2

template <class NT>
int BiPoly<NT>::expand(int n) {
    if ((n <= ydeg)||(n < 0))
      return -2;
    
    for(int i=ydeg+1; i <=n ;i++)
      coeffX.push_back(Polynomial<NT>::polyZero());
    
    ydeg = n;
    return n;
  }//expand

  // contract() gets rid of leading zero polynomials
  //	and returns the new (true) y-degree;
  //	It returns -2 if this is a no-op

template <class NT>
int BiPoly<NT>::contract() {
    int d = getTrueYdegree();
    if (d == ydeg)
      return (-2);  // nothing to do
    else{
      for (int i = ydeg; i> d; i--)
	coeffX.pop_back();
      ydeg = d;
    }
    return d;
  }//contract

  // Self-assignment
template <class NT>
BiPoly<NT> & BiPoly<NT>::operator=( const BiPoly<NT>& P) {
  if (this == &P)
    return *this;	// self-assignment
  ydeg = P.getYdegree();
  coeffX = P.coeffX;
  return *this;
  }//operator=


  // Self-addition
template <class NT>
BiPoly<NT> & BiPoly<NT>::operator+=( BiPoly<NT>& P) { // +=

    int d = P.getYdegree();
    if (d > ydeg)
      expand(d);
    for (int i = 0; i<=d; i++)
      coeffX[i] += P.coeffX[i];

  return *this;
  }//operator+=
   
  // Self-subtraction
template <class NT>
BiPoly<NT> & BiPoly<NT>::operator-=( BiPoly<NT>& P) { // -=
    int d = P.getYdegree();
    if (d > ydeg)
      expand(d);
  for (int i = 0; i<=d; i++)
    coeffX[i] -= P.coeffX[i];

  return *this;
  }//operator-=

  // Self-multiplication
template <class NT>
BiPoly<NT> & BiPoly<NT>::operator*=( BiPoly<NT>& P) { // *=
    int d = ydeg + P.getYdegree();

    std::vector<Polynomial<NT> > vP;

    Polynomial<NT>* c = new Polynomial<NT> [d + 1];
    for(int i=0; i <=d; i++)
      c[i] = Polynomial<NT>();
    
    for (int i = 0; i<=P.getYdegree(); i++)
      for (int j = 0; j<=ydeg; j++) {
	if(!zeroP(P.coeffX[i]) && !zeroP(coeffX[j]))
	  c[i+j] += P.coeffX[i] * coeffX[j];
      }

    for(int i=0; i <= d; i++)
      vP.push_back(c[i]);

    delete[] c;
    coeffX.clear();
    coeffX = vP;
    ydeg = d;
    return *this;
  }//operator*=
  
  // Multiply by a polynomial in X
template <class NT>
BiPoly<NT> & BiPoly<NT>::mulXpoly( Polynomial<NT> & p) {
    contract();
    if (ydeg == -1)
      return *this;

    for (int i = 0; i<=ydeg ; i++)
      coeffX[i] *= p;
    return *this;
  }//mulXpoly

  //Multiply by a constant
template <class NT>
BiPoly<NT> & BiPoly<NT>::mulScalar( NT & c) {
    for (int i = 0; i<=ydeg ; i++)
      coeffX[i].mulScalar(c);
    return *this;
  }//mulScalar

  // mulYpower: Multiply by Y^i (COULD be a divide if i<0)
template <class NT>
BiPoly<NT> & BiPoly<NT>::mulYpower(int s) {
  // if s >= 0, then this is equivalent to
  // multiplying by Y^s;  if s < 0, to dividing by Y^s
  if (s==0)
    return *this;
  int d = s+getTrueYdegree();
  if (d < 0) {
    ydeg = -1;
    deleteCoeffX();
    return *this;
  }

  std::vector<Polynomial<NT> > vP;
  if(s > 0){
    for(int i=0; i < s; i ++)
      vP.push_back(Polynomial<NT>::polyZero());
    for(int i=s; i<=d; i++)
      vP.push_back(coeffX[i-s]);
  }

  if (s < 0) {
    for (int j=-1*s; j <= d-s; j++)
      vP.push_back(coeffX[j]);
  }

  coeffX = vP;
  ydeg= d;
  return *this;
  }//mulYpower


  
  // Divide by a polynomial in X.
  // We replace the coeffX[i] by the pseudoQuotient(coeffX[i], P)
template <class NT>
BiPoly<NT> & BiPoly<NT>::divXpoly( Polynomial<NT> & p) {
    contract();
    if (ydeg == -1)
      return *this;

    for (int i = 0; i<=ydeg ; i++)
      coeffX[i] = coeffX[i].pseudoRemainder(p);

    return *this;
  }// divXpoly


  
  //Using the standard definition of pseudRemainder operation.
  //	--No optimization!
template <class NT>
BiPoly<NT> BiPoly<NT>::pseudoRemainderY (BiPoly<NT> & Q){
    contract();
    Q.contract();
    int qdeg = Q.getYdegree();
    Polynomial<NT> LC(Q.coeffX[qdeg]), temp1(LC), temp;
    BiPoly<NT> P(Q);

    std::vector<Polynomial<NT> > S;//The quotient in reverse order
    if(ydeg < qdeg){
      return (*this);
    }
    (*this).mulXpoly(temp1.power(ydeg - qdeg +1));

    while(ydeg >= qdeg){
      temp1 = coeffX[ydeg];
      temp = temp1.pseudoRemainder(LC);
      S.push_back(temp);
      P.mulXpoly(temp);
      P.mulYpower(ydeg - qdeg);
      *this -= P;
      contract();
      P = Q;//P is used as a temporary since mulXpower and mulYpower are
            // self modifying
    }

    //Correct the order of S
    std::vector<Polynomial<NT> > vP;

    for(int i= S.size()-1; i>=0; i--)
      vP.push_back(S[i]);


    return BiPoly<NT>(vP);
    
  }//pseudoRemainder

  //Partial Differentiation
  //Partial Differentiation wrt Y
template <class NT>
BiPoly<NT> & BiPoly<NT>::differentiateY() {	
  if (ydeg >= 0) {
    for (int i=1; i<=ydeg; i++)
      coeffX[i-1] = coeffX[i].mulScalar(i);
    ydeg--;
  }
  return *this;
  }// differentiation wrt Y

template <class NT>
BiPoly<NT> & BiPoly<NT>::differentiateX() {	
    if (ydeg >= 0)
      for (int i=0; i<=ydeg; i++)
	coeffX[i].differentiate();
    
    return *this;
  }// differentiation wrt X

template <class NT>
BiPoly<NT> & BiPoly<NT>::differentiateXY(int m, int n) {//m times wrt X and n times wrt Y
    assert(i >=0); assert(j >=0);
    for(int i=1; i <=m; i++)
      (*this).differentiateX();
    for(int i=1; i <=n; i++)
      (*this).differentiateY();
    
    return *this;
  }

  //Represents the bivariate polynomial in (R[X])[Y] as a member
  //of (R[Y])[X].
  //But since our polynomials in X can only have NT coeff's thus
  // to represent the above polynomial we switch X and Y once
  // the conversion has been done.
  //NOTE: This is different from replacing X by Y which was
  //      done in the case of the constructor from a polynomial in X
  //Need to calculate resultant wrt X.
template <class NT>
BiPoly<NT> & BiPoly<NT>::convertXpoly(){
    getTrueYdegree();
    if(ydeg == -1) return (*this);
    NT *cs = new NT[ydeg +1];
    int xdeg = getXdegree();
    std::vector<Polynomial<NT> > vP;

    for(int i=0; i<=xdeg; i++){
      for(int j=0; j<=ydeg; j++){
	cs[j] = coeffX[j].getCoeffi(i);
      }
      
      vP.push_back(Polynomial<NT>(ydeg, cs));
    }
    delete[] cs;
      
    ydeg = xdeg;
    coeffX = vP;
    return (*this);
  }

  //Set Coeffecient to the polynomial passed as a parameter
template <class NT>
bool BiPoly<NT>::setCoeff(int i, Polynomial<NT> p){
    if(i < 0 || i > ydeg)
      return false;
    coeffX[i] = p;
    return true;
  }
  
template <class NT>
void BiPoly<NT>::reverse() {
    Polynomial<NT> tmp;
    for (int i=0; i<= ydeg/2; i++) {
      tmp = coeffX[i];
      coeffX[i] =   coeffX[ydeg-i];
      coeffX[ydeg-i] = tmp;
    }
  }//reverse

  //replaceYwithX()
  //   used used when the coeffX in BiPoly are constants,
  //   to get the corresponding univariate poly in X
  //   E.g., Y^2+2*Y+9 will be converted to X^2+2*X+9
template <class NT>
Polynomial<NT>  BiPoly<NT>::replaceYwithX(){
    int m = getTrueYdegree();
    NT *cs = new NT[m+1];
    for(int i=0; i <= m ; i++)
      cs[i]=coeffX[i].getCoeffi(0);
    delete[] cs;

    return Polynomial<NT>(m,cs);
  }//replaceYwithX
  
template <class NT>
BiPoly<NT>&  BiPoly<NT>::pow(unsigned int n){

  if (n == 0) {
    ydeg = 0;
    deleteCoeffX();
    Polynomial<NT> unity(0);
    coeffX.push_back(unity);
  }else if(n == 1){

  }else{
    BiPoly<NT> x(*this);

    while ((n % 2) == 0) { // n is even
      x *= x;
      n >>= 1;
    }
    BiPoly<NT> u(x);
    while (true) {
      n >>= 1;
      if (n == 0){
	(*this) = u;
	return (*this);
      }
      x *= x;
      if ((n % 2) == 1) // n is odd
        u *= x;
    }
    (*this) = u;
  }
  return (*this);
}//pow



  ////////////////////////////////////////////////////////
  // Helper Functions
  ////////////////////////////////////////////////////////

// isZeroPinY(P)
//	checks whether a Bi-polynomial is a zero Polynomial
template <class NT>
bool isZeroPinY(BiPoly<NT> & P){
    if(P.getTrueYdegree() == -1)
      return true;
    return false;
}

// gcd(P,Q)
//   This gcd is based upon the subresultant PRS to avoid
//   exponential coeffecient growth and gcd computations, both of which 
//   are expensive since the coefficients are polynomials