C++ Polynomial::eval方法代码示例

bool Quadric::intersect(const Point3D &eye, const Point3D &_ray, HitReporter &hr) const
{
  Polynomial<2> x, y, z;
  const Vector3D ray = _ray - eye;
  x[0] = eye[0];
  x[1] = ray[0];
  y[0] = eye[1];
  y[1] = ray[1];
  z[0] = eye[2];
  z[1] = ray[2];

  const Polynomial<2> eqn = A * x * x + B * x * y + C * x * z
    + D * y * y + E * y * z + F * z * z + G * x + H * y + J * z + K;

  double ts[2];
  auto nhits = eqn.solve(ts);

  if(nhits > 0)
  {
    const Polynomial<2> ddx = 2 * A * x + B * y + C * z + G,
			ddy = 2 * D * y + B * x + E * z + H,
			ddz = 2 * F * z + C * x + E * y + J;

    for(int i = 0; i < nhits; i++)
    {
      // Figure out the normal.
      const double t = ts[i];
      const Point3D pt = eye + t * ray;
      if(!predicate(pt))
	continue;

      Vector3D normal(ddx.eval(t), ddy.eval(t), ddz.eval(t));

      // Figure out the uv.
      Point2D uv;
      Vector3D u, v;
      get_uv(pt, normal, uv, u, v);

      if(!hr.report(ts[i], normal, uv, u, v))
	return false;
    }
  }

  return true;
}

int main()
{

Polynomial a ({2,1,0,1,0}) ;
Polynomial b ;
Polynomial c (a);

cout<<a<<endl;
cout<<a.deriv()<<endl;

cout<<endl;
cout<<endl;
cout<<endl;


cout<<a<<endl;
b.setcoefs({1,1,1,1});

cout<<b.getcoef(0)<<endl;
cout<<b<<endl;
cout<< b.eval(2)<<endl;
cout<<endl;
cout<<endl;
cout<<endl;

cout<<b<<endl;
cout<<b*3<<endl;
cout<<endl;
cout<<endl;
cout<<endl;

cout<<b<<endl;
cout<<b+3<<endl;


cout<<endl;

return 0;

}

void isolate_root_with_cluster(double& low, double& hgh, Polynomial *p)
{
  dprint(low);
  dprint(hgh);
  // ensure starting conditions
  //
  PointVal pv_low(low,p);
  // double p_low = eval(low);
  if(pv_low.v() < 0.)
  {
    Polynomial p_minus = -(*p);
    p_minus.isolate_root(low,hgh);
    return;
  }

  // initialize algorithm
  //
  // initialize root bracket.
  //
  PointVal pv_hgh(hgh,p);
  RootBracketClusterState rb(pv_low,pv_hgh);
  dprint(rb.low().x());
  dprint(rb.low().v());
  dprint(rb.hgh().x());
  dprint(rb.hgh().v());
  //
  Polynomial Dp = p->get_derivative();

  // requirements of algorithm:
  // - second-order convergence in all cases
  // - high performance in generic cases
  // - handles repeated polynomial roots
  // - converges to a root between low and hgh
  // - residual of answer is less than machine precision

  // conditions maintained by algorithm:
  // * low < hgh
  // * p(low) > EPSILON
  // * p(hgh) < -EPSILON
  // using_hgh = fabs(p(hgh)) > fabs(p(low))

  // state of algorithm is specified by:
  // * low
  // * hgh
  // * using_hgh
  // * num_roots_in_cluster = assumed odd number of roots in cluster

  // iteration of algorithm:
  // use newton iteration on smaller side to look for root
  // - check that guess is between low and hgh
  //   and that new value has smaller residual
  //   else switch to secant method
  // - if newton does not change the sign
  //   increment num_roots_in_cluster until you leave interval,
  //   the sign changes, or the residual increases.
  //  - if the sign does change
  //    decrement num_roots_in_cluster until num_roots_in_cluster=1
  //    or the sign does not change (if residual does not shrink
  //    sufficiently modify num_roots_in_cluster).
  //  - do one iteration of bisection if residual did not
  //    shrink by factor of e.g. .5

  // assumed number of roots in root cluster we are converging towards
  goto start_newton_iteration;

 start_newton_iteration:
  // initialize newton iteration
  // choose smaller endpoint as starting value.
  goto newton_iteration;

 done:
   dprint(rb.num_evals());
   low = rb.low().x();
   hgh = rb.hgh().x();
   //return rb.x();
   return;

 // newton's method begins with a point which is one end of a root bracket.
 // it determines a displacement based on the value of the function
 // and its derivative at the point.  To handle root clusters
 // it assumes that the function is approximated by a shifted
 // and scaled monomial.
 newton_iteration:
  {
    if(rb.done()) goto done;
    dprintf("using newton: x0=%24.16e, p(x0)=%24.16e",rb.pv().x(),rb.pv().v());

    const PointVal start_pv0(rb.pv());
    //
    const double Dpx0 = Dp.eval(rb.pv().x());
    // bail if the derivative is too small
    if(fabs(Dpx0) < 1e-6)
    {
      dprintf("bad derivative: %f",Dpx0);
      goto newton_bailed;
    }
    double displacement = -rb.pv().v()/Dpx0;
   compute_x0:
    double x0 = rb.pv().x() + rb.num_roots_in_cluster()*displacement;
    // if Newton prediction is outside interval
//.........这里部分代码省略.........