/**************************************************************************** * Core Library Version 1.7, August 2004 * Copyright (c) 1995-2004 Exact Computation Project * All rights reserved. * * This file is part of CGAL (www.cgal.org). * * File: Curves.tcc * * Description: * This file contains the implementations of * functions defined by Curves.h * It is included through Curves.h * Please see Curves.h for more details. * * Author: Vikram Sharma and Chee Yap * Date: April 12, 2004 * * WWW URL: http://cs.nyu.edu/exact/ * Email: exact@cs.nyu.edu * * $URL$ * $Id$ * SPDX-License-Identifier: LGPL-3.0-or-later ***************************************************************************/ #include //CONSTRUCTORS FOR THE BIPOLY CLASS //////////////////////////////////////////////////////// //Constructors //////////////////////////////////////////////////////// template BiPoly::BiPoly(){ // zero polynomial ydeg = -1; } //BiPoly(n) template BiPoly::BiPoly(int n){// creates a BiPoly with nominal y-degree equal to n. // To support this constructor, you need functions // that are equivalent to "setCoeff(i, PolyNT q)" in the Poly Class // You also want "getCoeff(i)" functions. // E.g. BiPoly p(10); // PolyNT q(3, {1,2,3,4}); // p.setCoeff(10, q); // p.setCoeff(0, PolyNT()); // ydeg = n; if (n<0) return; // coeffX not needed for(int i=0; i<= ydeg; i++){ Polynomial temp; coeffX.push_back(temp); } } //BiPoly(vp) template BiPoly::BiPoly(std::vector > vp){ // From vector of Polynomials ydeg = vp.size() - 1; if(ydeg >=0){ coeffX = vp; } } //BiPoly(p, flag): // if true, it converts polynomial p(X) into P(Y) // if false, it creates the polynomial Y-p(X) template BiPoly::BiPoly(Polynomial p, bool flag){ if (flag){ ydeg = p.getTrueDegree(); if(ydeg >=0){ for(int i=0; i<=ydeg; i++){ Polynomial temp(0); temp.setCoeff(0, p.getCoeffi(i)); coeffX.push_back(temp); // does STL make a copy of temp? }//for }//if } else { ydeg = 1; coeffX.push_back(p); coeffX.push_back(Polynomial::polyUnity()); }//else }//BiPoly(p) //BiPoly(deg, d[], C[]): // Takes in a list of list of coefficients. // Each cofficient list represents a polynomial in X // // deg - ydeg of the bipoly // d[] - array containing the degrees of each coefficient (i.e., X poly) // C[] - list of coefficients, we use array d to select the // coefficients template BiPoly::BiPoly(int deg, int *d, NT *C){ ydeg = deg; Polynomial temp; int max=0; for(int i=0; i <=deg; i++) max = core_max(d[i],max); NT *c = new NT[max]; for(int i=0; i<= deg; i++){ for(int j=0; j <=d[i]; j++) c[j] = C[i+j]; temp = Polynomial(d[i],c); coeffX.push_back(temp); } delete[] c; }//BiPoly(deg,d[],C[]) //The BNF syntax is the following:- // [bipoly] -> [term]| [term] +/- [bipoly] // [term] -> [basic term]|[basic term] [term]|[basic term]*[term] // [basic term] -> [number]|'x'|'y'|[basic term]'^'[number] // | '(' [bipoly]')'|'-' //Unary minus is treated as a basic term template BiPoly::BiPoly(const char * s, char myX, char myY){ string ss(s); constructFromString(ss, myX, myY); } template BiPoly::BiPoly(const string & s, char myX, char myY){ string ss(s); constructFromString(ss, myX, myY); } template void BiPoly::constructFromString(string & s, char myX, char myY){ if((myX != 'x' || myX != 'X') && (myY != 'y' || myY != 'Y')){ //Replace myX with 'x' and myY with 'y' in s. unsigned int loc = s.find(myX, 0); while(loc != string::npos){ s.replace(loc,1,1,'x'); loc = s.find(myX, loc+1); } loc = s.find(myY, 0); while(loc != string::npos){ s.replace(loc,1,1,'y'); loc = s.find(myY, loc+1); } } (*this) = (*this).getbipoly(s); } //Constructor from another BiPoly template BiPoly::BiPoly(const BiPoly &P) : ydeg(P.ydeg), coeffX(P.coeffX) { } //Destructor template void BiPoly::deleteCoeffX(){ coeffX.clear(); } template BiPoly::~BiPoly(){ if (ydeg >= 0) deleteCoeffX(); } //////////////////////////////////////////////////////// // METHODS USED BY STRING CONSTRUCTOR ABOVE //////////////////////////////////////////////////////// //Sets the input BiPoly to X^n template void BiPoly::constructX(int n, BiPoly& P){ CGAL_assertion(n>= -1); P.deleteCoeffX();//Clear the present coeffecients Polynomial q(n);//Nominal degree n q.setCoeff(n,NT(1)); if (n>0) q.setCoeff(0,NT(0)); P.coeffX.push_back(q); P.ydeg = 0; } //Sets the input BiPoly to Y^n template void BiPoly::constructY(int n, BiPoly& P){ P.ydeg = 0; P.deleteCoeffX();//Clear the present coeffecients NT cs[] = {1}; Polynomial q(0, cs); P.coeffX.push_back(q);//P is the unit bipoly. Now we just multiply Y^n P.mulYpower(n); } //Returns in P the coeffecient starting from start template int BiPoly::getnumber(const char* c, int start, unsigned int len, BiPoly & 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); P.deleteCoeffX(); Polynomial q(0); q.setCoeff(0, cf); P.coeffX.push_back(q); P.ydeg = 0; delete[] temp; return (j-1+start);//Exactly the length of the number } template bool BiPoly::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 int BiPoly::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 int BiPoly::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 int BiPoly::getbasicterm(string s, BiPoly & P){ const char * cstr = s.c_str(); unsigned int len = s.length(); int i=0; //BiPoly * temp = new BiPoly(); 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] == 'y'||cstr[i] == 'Y'){ constructY(1, P); }else if(cstr[i] =='('){ int oldi = i; i = matchparen(cstr, i); string t = s.substr(oldi+1, i -oldi -1); P = getbipoly(t); }else if(cstr[i] == '-'){//Unary Minus P.ydeg = 0; P.deleteCoeffX();//Clear the present coeffecients NT cs[] = {-1}; Polynomial q(0, cs); P.coeffX.push_back(q); }else{ #ifdef CGAL_CORE_TRACE std::cout <<"ERROR IN PARSING BASIC TERM" << std::endl; #endif } //i+1 points to the beginning of next syntaxtic object in the string. if(cstr[i+1] == '^'){ int n; i = getint(cstr, i+2, len, n); P.pow(n); } return i; } template int BiPoly::getterm(string s, BiPoly & P){ unsigned int len = s.length(); if(len == 0){// Zero BiPoly P = BiPoly(); 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; } BiPoly R; ind = oind + getbasicterm(t, R);//Because the second term is the offset in //t P *= R; } return ind; } template BiPoly BiPoly::getbipoly(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 BiPoly(); } //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; BiPoly 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; NT negone(-1); P.mulScalar(negone); }else{ ind = getterm(s, P); } unsigned int oind =0;//the string between oind and ind is a term while(ind != len -1){ BiPoly 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{ #ifdef CGAL_CORE_TRACE std::cout << "ERROR IN PARSING BIPOLY! " << std::endl; #endif } } return P; } //////////////////////////////////////////////////////// // METHODS //////////////////////////////////////////////////////// // filedump (msg, ofs, com, com2) // where msg and com are strings. // msg is an message and com is the character preceding each line // (e.g., msg="" and com=com2="# ") // This is called by the other dump functions template void BiPoly::dump(std::ostream & os, std::string msg, std::string com, std::string com2) const { if (msg != "") os << msg << std::endl; if(ydeg == -1) { os << com << " Zero Polynomial" << std::endl; return; } bool first = true; os << com; for (int i=0; i <= ydeg; i++){ if (!zeroP(coeffX[i])){ if (i % 3 == 0) os << std::endl << com2 ; // output 3 coefficients per line if (first) first = false; else os << " + "; os << "["; coeffX[i].filedump(os,"","",com2); // Note: first comment string is "" os << "]"; if (i > 0) os <<" * y^"<< i ; } } }//dump // dump(ofs, msg, com, com2) -- dump to file // /* template void BiPoly::dump(std::ofstream & ofs, std::string msg, std::string com, std::string com2) const { dump(ofs, msg, com, com2); } */ // dump(msg, com) -- dump to std output template void BiPoly::dump(std::string msg, std::string com, std::string com2) const { dump(std::cout, msg, com, com2); } /* *********************************************************** We want to substitute X or Y with a general Expression or BigFloat into a BiPoly. E.g., the following produces a Y-polynomial after replacing X by x: BiPoly yPolynomial(const Expr & x) { VecExpr vE; for (int i=0; i<= ydeg; i++) { vE.push_back(coeffX[i].eval(x)); } return Polynomial(vE); }//yPolynomial But this has many problems. Solution below is to have special yExprPolynomial(x). *********************************************************** */ template Polynomial BiPoly::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(d, coeffVec); } // Specialized version of yPolynomial for Expressions template Polynomial BiPoly::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(d, coeffVec); } // Specialized version of yPolynomial for BigFloat // ASSERTION: BigFloat x is exact template Polynomial BiPoly::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(d, coeffVec); } // xPolynomial(y) // returns the polynomial (in X) when we substitute Y=y // // N.B. May need the // Polynomial xExprPolynomial(Expr y) // version too... // template Polynomial BiPoly::xPolynomial(const NT & y) { Polynomial 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 int BiPoly::getYdegree() const { return ydeg; } // getXdegree() template int BiPoly::getXdegree(){ int deg=-1; for(int i=0; i <=ydeg; i++) deg = max(deg, coeffX[i].getTrueDegree()); return deg; } // getTrueYdegree template int BiPoly::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 Expr BiPoly::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 int BiPoly::expand(int n) { if ((n <= ydeg)||(n < 0)) return -2; for(int i=ydeg+1; i <=n ;i++) coeffX.push_back(Polynomial::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 int BiPoly::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 BiPoly & BiPoly::operator=( const BiPoly& P) { if (this == &P) return *this; // self-assignment ydeg = P.getYdegree(); coeffX = P.coeffX; return *this; }//operator= // Self-addition template BiPoly & BiPoly::operator+=( BiPoly& 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 BiPoly & BiPoly::operator-=( BiPoly& 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 BiPoly & BiPoly::operator*=( BiPoly& P) { // *= int d = ydeg + P.getYdegree(); std::vector > vP; Polynomial* c = new Polynomial [d + 1]; for(int i=0; i <=d; i++) c[i] = Polynomial(); 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 BiPoly & BiPoly::mulXpoly( Polynomial & 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 BiPoly & BiPoly::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 BiPoly & BiPoly::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 > vP; if(s > 0){ for(int i=0; i < s; i ++) vP.push_back(Polynomial::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 BiPoly & BiPoly::divXpoly( Polynomial & 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 BiPoly BiPoly::pseudoRemainderY (BiPoly & Q){ contract(); Q.contract(); int qdeg = Q.getYdegree(); Polynomial LC(Q.coeffX[qdeg]), temp1(LC), temp; BiPoly P(Q); std::vector > 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 > vP; for(int i= S.size()-1; i>=0; i--) vP.push_back(S[i]); return BiPoly(vP); }//pseudoRemainder //Partial Differentiation //Partial Differentiation wrt Y template BiPoly & BiPoly::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 BiPoly & BiPoly::differentiateX() { if (ydeg >= 0) for (int i=0; i<=ydeg; i++) coeffX[i].differentiate(); return *this; }// differentiation wrt X template BiPoly & BiPoly::differentiateXY(int m, int n) {//m times wrt X and n times wrt Y CGAL_assertion(m >=0); CGAL_assertion(n >=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 BiPoly & BiPoly::convertXpoly(){ getTrueYdegree(); if(ydeg == -1) return (*this); NT *cs = new NT[ydeg +1]; int xdeg = getXdegree(); std::vector > 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(ydeg, cs)); } delete[] cs; ydeg = xdeg; coeffX = vP; return (*this); } //Set Coeffecient to the polynomial passed as a parameter template bool BiPoly::setCoeff(int i, Polynomial p){ if(i < 0 || i > ydeg) return false; coeffX[i] = p; return true; } template void BiPoly::reverse() { Polynomial 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 Polynomial BiPoly::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(m,cs); }//replaceYwithX template BiPoly& BiPoly::pow(unsigned int n){ if (n == 0) { ydeg = 0; deleteCoeffX(); Polynomial unity(0); coeffX.push_back(unity); }else if(n == 1){ }else{ BiPoly x(*this); while ((n % 2) == 0) { // n is even x *= x; n >>= 1; } BiPoly 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 bool isZeroPinY(BiPoly & 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 template BiPoly gcd( BiPoly& P ,BiPoly& Q){ int m = P.getTrueYdegree(); int n = Q.getTrueYdegree(); //If both Bi-Polys are zero then return zero if( m==-1 && n==-1){ return BiPoly();//Had to do this for the type of //return value o/w causes problem in } //assignment operation if(m < n) { return gcd(Q, P); } //If one of the Bi-Polys are zero then return the other if(n==-1) return P; //When the two BiPolys are just univariate Polynomials in X if(m==0 && n==0){ std::vector > vP; vP.push_back(gcd(P.coeffX[0], Q.coeffX[0])); return BiPoly(vP);//Had to do this for the type of //return value } int delta = m - n; Polynomial a(Q.coeffX[n]); std::vector vN; vN.push_back(pow(BigInt(-1),delta + 1)); Polynomial beta(vN); Polynomial phi(power(a, delta)), t; m = n; BiPoly temp; P.pseudoRemainderY(Q); while(!isZeroPinY(P)){ P.divXpoly(beta); n = P.getTrueYdegree(); delta = m - n; beta = power(phi, delta)*a.mulScalar(pow(BigInt(-1),delta+1)); a = P.coeffX[n]; t = power(phi, delta -1); phi = power(a,delta).pseudoRemainder(t); m = n; temp = Q;//Swapping the pseudo-remainder for the next step Q = P; P = temp; P.pseudoRemainderY(Q); } return Q; }//gcd // resY(P,Q): // Resultant of Bi-Polys P and Q w.r.t. Y. // So the resultant is a polynomial in X template Polynomial resY( BiPoly& P ,BiPoly& Q){ int m = P.getTrueYdegree(); int n = Q.getTrueYdegree(); //If either polynomial is zero, return zero if( m == -1 || n == -1) return Polynomial();//::polyZero(); if(n > m) return resY(Q, P); Polynomial b(Q.coeffX[n]); Polynomial lc(P.coeffX[m]), C, temp; BiPoly r; r = P.pseudoRemainderY(Q); C = b * r.coeffX[r.getTrueYdegree()]; C = C.pseudoRemainder(lc); if(isZeroPinY(P) && n > 0) return Polynomial();//::polyZero(); if(Q.getTrueYdegree() == 0 && isZeroPinY(P)) return power(Q.coeffX[0], m); int l = P.getTrueYdegree(); temp = power(b, m-l).mulScalar(pow(BigInt(-1),m*n))*resY(Q,P); temp = temp.pseudoRemainder(power(C,n)); return temp; }//resY // resX(P,Q): // Resultant of Bi-Polys P and Q w.r.t. X. // So the resultant is a polynomial in Y // We first convert P, Q to polynomials in X. Then // call resY and then turns it back into a polynomial in Y // QUESTION: is this last switch really necessary??? template BiPoly resX( BiPoly& P ,BiPoly& Q){ P.convertXpoly(); Q.convertXpoly(); // Polynomial p(resY(P,Q)); // return BiPoly(p); return(resY(P,Q)); // avoid the switch back }//resX //Equality operator for BiPoly template bool operator==(const BiPoly& P, const BiPoly& Q) { // == BiPoly P1(P); BiPoly Q1(Q); P1.contract(); Q1.contract(); if(P1.getYdegree() != Q1.getYdegree()) return false; else{ for(int i=0; i <= P1.getYdegree() ; i++){ if(P1.coeffX[i] != Q1.coeffX[i]) return false; } } return true; } // Addition P + Q template BiPoly operator+(const BiPoly& P, const BiPoly& Q ) { // + return BiPoly(P) += Q; } // Subtraction P - Q template BiPoly operator-(const BiPoly& P,const BiPoly& Q ) { // + return BiPoly(P) -= Q; } // Multiplication P * Q template BiPoly operator*(const BiPoly& P, const BiPoly& Q ) { // + return BiPoly(P) *= Q; } //////////////////////////////////////////////////////// //Curve Class // extends BiPoly Class //////////////////////////////////////////////////////// template < class NT > Curve::Curve(){} // zero polynomial //Curve(vp): // construct from a vector of polynomials template < class NT > Curve::Curve(std::vector > vp) : BiPoly(vp){ } //Curve(p): // Converts a polynomial p(X) to a BiPoly in one of two ways: // (1) if flag is false, the result is Y-p(X) // (2) if flag is true, the result is p(Y) // The default is (1) because we usually want to plot the // graph of the polynomial p(X) template < class NT > Curve::Curve(Polynomial p, bool flag) : BiPoly(p, flag){ } //Curve(deg, d[], C[]): // Takes in a list of list of coefficients. // Each cofficient list represents a polynomial in X // // deg - ydeg of the bipoly // d[] - array containing the degrees of each coefficient (i.e., X poly) // C[] - list of coefficients, we use array d to select the // coefficients template < class NT > Curve::Curve(int deg, int *d, NT *C) : BiPoly(deg, d, C){ } template < class NT > Curve::Curve(const BiPoly &P) : BiPoly(P){ } //Curve(n) template < class NT > Curve::Curve(int n) : BiPoly(n){// creates a Curve with nominal y-degree equal to n } //Creates a curve from a string (no parentheses, no *) template < class NT > Curve::Curve(const string & s, char myX, char myY) : BiPoly(s, myX, myY){ } template < class NT > Curve::Curve(const char * s, char myX, char myY) : BiPoly(s, myX, myY){ } // verticalIntersections(x, vecI, aprec=0): // // The list vecI is passed an isolating intervals for y's such that (x,y) // lies on the curve. // If aprec is non-zero (!), the intervals have with < 2^{-aprec}. // Return is -2 if curve equation does not depend on Y // -1 if infinitely roots at x, // 0 if no roots at x // 1 otherwise // // ASSERTION: x is an exact BigFloat template < class NT > int Curve::verticalIntersections(const BigFloat & x, BFVecInterval & vI, int aprec) { int d= Curve::getTrueYdegree(); if(d <= 0) return(-2); // This returns a NULL vI, which should be caught by caller Polynomial PY = this->yExprPolynomial(x); // should be replaced // CGAL_assertion(x.isExact()); // Polynomial PY = yBFPolynomial(x); // unstable still d = PY.getTrueDegree(); if(d <= 0) return(d); Sturm Ss(PY); // should be replaced by BigFloat version // Sturm Ss(PY); // unstable still Ss.isolateRoots(vI); int s = vI.size(); if ((aprec != 0) && (s>0)) Ss.newtonRefineAllRoots(vI, aprec); return s; } // plot(eps, x1, y1, x2, y2) // // All parameters have defaults // // Gives the points on the curve at resolution "eps". Currently, // eps is viewed as delta-x step size (but it could change). // The display is done in the rectangale // defined by [(x1, y1), (x2, y2)]. // The output is written into a file in the format specified // by our drawcurve function (see COREPATH/ext/graphics). // // Heuristic: the open polygonal lines end when number of roots // changes... // // TO DO: // (1) need to automatically readjust x- and y-scales // independently // (2) reorder parameters in this order: x1, x2, y1, y2 // [Then y1, y2 can be automatically determined] // (3) should allow the ability to look for interesting // features // // ASSERTION: all input BigFloats are exact // template < class NT > int Curve::plot( BigFloat eps, BigFloat x1, BigFloat y1, BigFloat x2, BigFloat y2, int fileNo){ const char* filename[] = {"data/input", "data/plot", "data/plot2"}; CGAL_assertion(eps.isExact()); // important for plotting... CGAL_assertion(x1.isExact()); CGAL_assertion(y1.isExact()); ofstream outFile; outFile.open(filename[fileNo]); // ought to check if open is successful! outFile << "########################################\n"; outFile << "# Curve equation: \n"; this->dump(outFile,"", "# "); outFile << std::endl; outFile << "# Plot parameters: step size (eps) = " << eps << std::endl; outFile << "# x1 = " << x1 << ",\t y1 = " << y1 << std::endl; outFile << "# x2 = " << x2 << ",\t y2 = " << y2 << std::endl; outFile << "########################################\n"; outFile << "# X-axis " << std::endl; outFile << "o 2" << std::endl; outFile << "0.99 \t 0.99 \t 0.99" << std::endl; outFile << x1 << "\t" << 0 << std::endl; outFile << x2 << "\t" << 0 << std::endl; outFile << "########################################\n"; outFile << "# Y-axis " << std::endl; outFile << "o 2" << std::endl; outFile << "0.99 \t 0.99 \t 0.99" << std::endl; outFile << 0 << "\t" << y1 << std::endl; outFile << 0 << "\t" << y2 << std::endl; // CGAL_assertion(eps>0) int aprec = -(eps.lMSB()).toLong(); // By definition, eps.lMSB() <= lg(eps) BigFloat xCurr = x1; BigFloat xLast = x1; unsigned int numRoots=0; unsigned int numPoints=0; BFVecInterval vI; cout <<"Current value of x " << xCurr << endl; //=================================================================== // FIRST locate the first x-value where there are roots for plotting //=================================================================== do { vI.clear(); if (verticalIntersections(xCurr, vI, aprec) > 0) { numRoots = vI.size(); cout <<"Number of roots at " << xCurr << " are " << numRoots< 0 then there exists a //valid point for plotting return -1; // nothing to plot! } int limit = ((x2 - xCurr + eps)/eps).intValue()+1; //std::cout << "Limit = " << limit << std::endl; machine_double plotCurves[this->getTrueYdegree()][limit];//plot buffer machine_double yval; for (unsigned int i=0; i< numRoots; i++) { yval = (vI[i].first + vI[i].second).doubleValue()/2; if (yval < y1) plotCurves[i][numPoints] = y1.doubleValue(); else if (yval > y2) plotCurves[i][numPoints] = y2.doubleValue(); else plotCurves[i][numPoints] = yval; } vI.clear(); xLast = xCurr - eps; // -eps to compensate for the forward value of xCurr numPoints = 1; //=================================================================== // Get all the curves in a main loop // -- dump the curves when an x-interval is discovered // We define an "x-interval" to be a maximal interval // where the number of curves is constant (this is a heuristic!) // Note that this includes the special case where number is 0. //=================================================================== BigFloat tmp; // used to step from xLast to xCurr in loops while (xCurr < x2) { //main loop //std::cout << "Doing verticalintersec at " << xCurr << std::endl; verticalIntersections(xCurr, vI, aprec); if (vI.size() != numRoots) { // an x-interval discovered! // write previous x-interval to output file outFile << "########################################\n"; outFile << "# New x-interval with " << numRoots << " roots\n"; for (unsigned int i=0; i< numRoots; i++) { outFile << "#=======================================\n"; outFile << "# Curve No. " << i+1 << std::endl; outFile << "o " << numPoints << std::endl; outFile << red_comp(i) << "\t" << green_comp(i) << "\t" << blue_comp(i) << std::endl; tmp = xLast; for (unsigned int j=0; j< numPoints; j++) { outFile << tmp << " \t\t" << plotCurves[i][j] << "\n"; tmp += eps; }//for j }//for i numPoints = 0; // reset numRoots = vI.size(); // reset xLast = xCurr; // reset }//if vI.size() != if (numRoots>0){ // record curr. vertical intersections if numRoots>0 for (unsigned int i=0; i< numRoots; i++) { yval = (vI[i].first + vI[i].second).doubleValue()/2; // HERE SHOULD BE A LOOP TO OUTPUT MORE POINTS IN CASE THE slope IS LARGE // Idea: let previous value of yval be yval-old. // // Two cases: // (1) When i=0: // you need to search backwards and forwards // (yval-old is not defined in this case) // // (2) When i>0: // If |yval-old - yval| > 2eps, we must plot using horizontalIntersection() // We must start from // y = yval-old until // EITHER (y = (vI[i+1].first + vI[i+1].second)/2) // OR (sign of slope changes) // OR (hit the ymax or ymin boundary) // if (yval < y1) plotCurves[i][numPoints] = y1.doubleValue(); else if (yval > y2) plotCurves[i][numPoints] = y2.doubleValue(); else plotCurves[i][numPoints] = yval; }//for i vI.clear(); //Should clear the intersection points numPoints++; }//if xCurr += eps; if (!xCurr.isExact()){ std::ostringstream oss; oss << "xCurr has error! xCurr=" << xCurr; core_error(oss.str(), __FILE__, __LINE__, false); } }//main while loop // Need to flush out the final x-interval: if ((numRoots>0) && (numPoints >0)) { // write to output file outFile << "########################################\n"; outFile << "# New x-interval with " << numRoots << " roots\n"; for (unsigned int i=0; i< numRoots; i++) { outFile << "#=======================================\n"; outFile << "# Curve No. " << i+1 << std::endl; outFile << "o " << numPoints << std::endl; outFile << red_comp(i) << "\t" << green_comp(i) << "\t" << blue_comp(i) << std::endl; tmp = xLast; for (unsigned int j=0; j< numPoints; j++) { outFile << tmp << " \t\t" << plotCurves[i][j] << "\n"; tmp += eps; }//for j }//for i }//if // Put out the final frame (this hides the artificial cut-off of curves outFile << "########################################\n"; outFile << "# FRAME around our figure " << std::endl; outFile << "p 4" << std::endl; outFile << "0.99 \t 0.99 \t 0.99" << std::endl; outFile << x1 << "\t" << y1 << std::endl; outFile << x2 << "\t" << y1 << std::endl; outFile << x2 << "\t" << y2 << std::endl; outFile << x1 << "\t" << y2 << std::endl; outFile << "######### End of File ##################\n"; outFile.close(); return 0; }//plot // selfIntersections(): // this should be another member function that lists // all the self-intersections of a curve // // template // void selfIntersections(BFVecInterval &vI){ // ... // } //////////////////////////////////////////////////////// // Curve helper functions //////////////////////////////////////////////////////// //Xintersections(C, D, vI): // returns the list vI of x-ccordinates of possible intersection points. // Assumes that C & D are quasi-monic.(or generally aligned) template void Xintersections( Curve& P ,Curve& Q, BFVecInterval &vI){ Sturm Ss(resY(P, Q)); Ss.isolateRoots(vI); } //Yintersections(C, D, vI): // similar to Xintersections template void Yintersections( Curve& P ,Curve& Q, BFVecInterval &vI){ Sturm Ss(resX(P, Q)); Ss.isolateRoots(vI); } // Display Intervals // template //DO I NEED THIS OVERHERE AS WELL? void showIntervals(char* s, BFVecInterval &vI) { std::cout << s; for (unsigned int i=0; i< vI.size(); i++) { std::cout << "[ " << vI[i].first << ", " << vI[i].second << " ], " ; } std::cout << std::endl; } /*************************************************************************** */ // END /*************************************************************************** */