curves.tcc

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

    if(m==0 && n==0){
      std::vector<Polynomial<NT> > vP;
      vP.push_back(gcd(P.coeffX[0], Q.coeffX[0]));
      return BiPoly<NT>(vP);//Had to do this for the type of
                            //return value
      }
    
    int delta = m - n;
    Polynomial<NT> a(Q.coeffX[n]);
    std::vector<NT> vN;
    vN.push_back(pow(BigInt(-1),delta + 1));
    Polynomial<NT> beta(vN);
    Polynomial<NT> phi(power(a, delta)), t;
    m = n;
    BiPoly<NT> 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 <class NT>
Polynomial<NT>  resY( BiPoly<NT>& P ,BiPoly<NT>& Q){

  int m = P.getTrueYdegree();
  int n = Q.getTrueYdegree();

  //If either polynomial is zero, return zero
  if( m == -1 || n == -1) return Polynomial<NT>();//::polyZero();

  if(n > m)
    return resY(Q, P);

  Polynomial<NT> b(Q.coeffX[n]);
  Polynomial<NT> lc(P.coeffX[m]), C, temp;
  BiPoly<NT> r;
  
  r = P.pseudoRemainderY(Q);
  C = b * r.coeffX[r.getTrueYdegree()];
  C = C.pseudoRemainder(lc);

  if(isZeroPinY(P) && n > 0)
    return Polynomial<NT>();//::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 <class NT>
BiPoly<NT>  resX( BiPoly<NT>& P ,BiPoly<NT>& Q){
  P.convertXpoly();
  Q.convertXpoly();

  // Polynomial<NT> p(resY(P,Q));
  // return BiPoly<NT>(p);
  return(resY(P,Q)); // avoid the switch back
}//resX


//Equality operator for BiPoly
template <class NT>
bool operator==(const BiPoly<NT>& P, const BiPoly<NT>& Q) {	// ==
  BiPoly<NT> P1(P);
  BiPoly<NT> 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 <class NT>
BiPoly<NT> operator+(const BiPoly<NT>& P, const BiPoly<NT>& Q ) { // +
  return BiPoly<NT>(P) += Q;
}
  // Subtraction P - Q
template <class NT>
BiPoly<NT> operator-(const  BiPoly<NT>& P,const  BiPoly<NT>& Q ) { // +
  return BiPoly<NT>(P) -= Q;
}

  // Multiplication P * Q
template <class NT>
BiPoly<NT> operator*(const  BiPoly<NT>& P, const BiPoly<NT>& Q ) { // +
  return BiPoly<NT>(P) *= Q;
}


  ////////////////////////////////////////////////////////
  //Curve Class
  //  	extends BiPoly Class
  ////////////////////////////////////////////////////////

template < class NT >
Curve<NT>::Curve(){} // zero polynomial
  
  //Curve(vp):
  //    construct from a vector of polynomials
template < class NT >
Curve<NT>::Curve(std::vector<Polynomial<NT> > vp)
	  : BiPoly<NT>(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<NT>::Curve(Polynomial<NT> p, bool flag)
	  : BiPoly<NT>(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<NT>::Curve(int deg, int *d, NT *C)
	  : BiPoly<NT>(deg, d, C){
  }

template < class NT >
Curve<NT>::Curve(const BiPoly<NT> &P)
	  : BiPoly<NT>(P){
  }

  //Curve(n)
template < class NT >
Curve<NT>::Curve(int n)
	  : BiPoly<NT>(n){// creates a Curve with nominal y-degree equal to n
  }

  //Creates a curve from a string (no parentheses, no *)
template < class NT >
Curve<NT>::Curve(const string & s, char myX, char myY)
	  : BiPoly<NT>(s, myX, myY){
  }
template < class NT >
Curve<NT>::Curve(const char * s, char myX, char myY)
	  : BiPoly<NT>(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<NT>::verticalIntersections(const BigFloat & x, BFVecInterval & vI,
	int aprec) {
    int d= Curve<NT>::getTrueYdegree();
    if(d <= 0) return(-2);
    	   // This returns a NULL vI, which should be caught by caller
	
    Polynomial<Expr> PY = this->yExprPolynomial(x); // should be replaced
    // assert(x.isExact());
    // Polynomial<BigFloat> PY = yBFPolynomial(x); // unstable still

    d = PY.getTrueDegree();
    if(d <= 0) return(d);

    Sturm<Expr> SS(PY); // should be replaced by BigFloat version
    // Sturm<BigFloat> 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<NT>::plot( BigFloat eps, BigFloat x1,
	BigFloat y1, BigFloat x2, BigFloat y2, int fileNo){

  const char* filename[] = {"data/input", "data/plot", "data/plot2"};

  assert(eps.isExact()); // important for plotting...
  assert(x1.isExact());
  assert(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;
  // assert(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<<endl;
    }
    xCurr += eps;

  } while ((numRoots == 0) && (xCurr < x2));//numRoots <= ydeg

  if (numRoots == 0 && x2 <= xCurr) {//if 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::cout<<"xCurr has error! xCurr=" << xCurr << std::endl;
   }//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 <class NT>
//  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 <class NT>
void  Xintersections( Curve<NT>& P ,Curve<NT>& Q, BFVecInterval &vI){
  Sturm<NT> SS(resY(P, Q));
  SS.isolateRoots(vI);
}

//Yintersections(C, D, vI):
//	similar to Xintersections
template <class NT>
void  Yintersections( Curve<NT>& P ,Curve<NT>& Q, BFVecInterval &vI){
  Sturm<NT> SS(resX(P, Q));
  SS.isolateRoots(vI);
}

// Display Intervals
// 
template <class NT>//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
/*************************************************************************** */