From: Lee Worden <worden@be...>  20060623 01:23:28

Hi, I haven't been on vxlusers, and I'm having trouble connecting to the archive, so I apologize if this has already been dealt with. My project is generating vnl_real_npolynomial objects dynamically from evolving data structures in a simulation and passing them to vnl_rnpoly_solve, so it is probably stressing the code more than a lot of applications of those objects. I have workarounds for 2 bugs in vnl_real_npolynomial in VXL 1.5.1 and earlier. I'm posting them here in case they're useful to someone or they might be of use in a future release. Feedback is welcome. Thanks for the library  it's great! Lee  first: simplify() doesn't always remove terms with zero coefficient  it fails when there is one or more zero term at the end of the vector. Here is a function that can fix up the object before you use it. The main change is marked with "***". void simplify(vnl_real_npolynomial &poly) { vnl_vector<double> coeffs_ = poly.coefficients(); vnl_matrix<unsigned int> polyn_ = poly.polyn(); unsigned int nvar_ = polyn_.cols(); unsigned int nterms_ = polyn_.rows(); for (unsigned int row1=0; row1<nterms_; ++row1) for (unsigned int row2=row1+1; row2<nterms_; ++row2) { unsigned int col=0; while (col<nvar_ && polyn_(row1,col) == polyn_(row2,col)) ++col; if (col < nvar_) continue; // not all exponents are identical coeffs_(row1) += coeffs_(row2); coeffs_(row2) = 0; } while (nterms_>0 && coeffs_(nterms_1)==0) nterms_; // *** for (unsigned int row=0; row<nterms_; ++row) if (coeffs_(row) == 0) { nterms_; // decrement nterms, and move last element to // vacant place: coeffs_(row) = coeffs_(nterms_); coeffs_(nterms_) = 0; // not really necessary; to keep // coeffs_ consistent for (unsigned int i=0; i<nvar_; ++i) polyn_(row,i) = polyn_(nterms_,i); } if (nterms_ < polyn_.rows()) poly.set(coeffs_.extract(nterms_),polyn_.extract(nterms_,nvar_)); } second: the << operator generates incorrect output for coefficients of 1 or 1 in many cases. The following code works correctly, I think. I also changed it to write, for example, "X0^2 X1" instead of "X0^2X1", though it will still write "X0X1" and so on. // you can use stream << asString(P) instead of stream << P string asString(vnl_real_npolynomial &P) { ostringstream os; vnl_vector<double> coeffs_ = P.coefficients(); vnl_matrix<unsigned int> polyn_ = P.polyn(); unsigned int nvar_ = polyn_.cols(); unsigned int nterms_ = polyn_.rows(); if (nvar_ <= 3) for (unsigned int i=0; i<nterms_; ++i) { if (i>0) { os << ' '; if (coeffs_(i) >= 0) os << "+ "; } if (coeffs_(i) < 0) os << " "; if (fabs(coeffs_(i)) != 1) os << fabs(coeffs_(i)) << ' '; unsigned int totaldeg = 0; if (nvar_ > 0 && polyn_(i,0) > 0) { os << 'X'; totaldeg += polyn_(i,0); } if (nvar_ > 0 && polyn_(i,0) > 1) os << '^' << polyn_(i,0) << ' '; if (nvar_ > 1 && polyn_(i,1) > 0) { os << 'Y'; totaldeg += polyn_(i,1); } if (nvar_ > 1 && polyn_(i,1) > 1) os << '^' << polyn_(i,1) << ' '; if (nvar_ > 2 && polyn_(i,2) > 0) { os << 'Z'; totaldeg += polyn_(i,2); } if (nvar_ > 2 && polyn_(i,2) > 1) os << '^' << polyn_(i,2) << ' '; if (totaldeg == 0 && vcl_fabs(coeffs_(i)) == 1) os << fabs(coeffs_(i)); } else for (unsigned int i=0; i<nterms_; ++i) { if (i>0) { os << ' '; if (coeffs_(i) >= 0) os << "+ "; } if (coeffs_(i) < 0) os << " "; if (fabs(coeffs_(i)) != 1) os << fabs(coeffs_(i)) << ' '; unsigned int totaldeg = 0; for (unsigned int j=0; j<nvar_; ++j) { if (polyn_(i,j) > 0) os << 'X' << j; if (polyn_(i,j) > 1) os << '^' << polyn_(i,j) << ' '; totaldeg += polyn_(i,j); } if (totaldeg == 0 && vcl_fabs(coeffs_(i)) == 1) os << fabs(coeffs_(i)); } //os << endl; return os.str(); } 