C++ Poly::back方法代码示例

Poly divide(Poly const &a, Poly const &b, Poly &r) {
	Poly c;
	r = a; // remainder
	assert(!b.empty());

	const unsigned k = a.degree();
	const unsigned l = b.degree();
	c.resize(k, 0.);

	for(unsigned i = k; i >= l; i--) {
		assert(i >= 0);
		double ci = r.back()/b.back();
		c[i-l] += ci;
		Poly bb = ci*b;
		//std::cout << ci <<"*(" << b.shifted(i-l) << ") = "
		//          << bb.shifted(i-l) << "     r= " << r << std::endl;
		r -= bb.shifted(i-l);
		r.pop_back();
	}
	//std::cout << "r= " << r << std::endl;
	r.normalize();
	c.normalize();

	return c;
}

void reduce(Poly &v)
{
  while (!v.empty() && v.back() % MOD == 0)
  {
    v.pop_back();
  }
}

cdouble laguerre_internal_complex(Poly const & p,
                                  double x0,
                                  double tol,
                                  bool & quad_root) {
    cdouble a = 2*tol;
    cdouble xk = x0;
    double n = p.degree();
    quad_root = false;
    const unsigned shuffle_rate = 10;
//    static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
    unsigned shuffle_counter = 0;
    while(std::norm(a) > (tol*tol)) {
        //std::cout << "xk = " << xk << std::endl;
        cdouble b = p.back();
        cdouble d = 0, f = 0;
        double err = abs(b);
        double abx = abs(xk);
        for(int j = p.size()-2; j >= 0; j--) {
            f = xk*f + d;
            d = xk*d + b;
            b = xk*b + p[j];
            err = abs(b) + abx*err;
        }

        err *= 1e-7; // magic epsilon for convergence, should be computed from tol

        cdouble px = b;
        if(abs(b) < err)
            return xk;
        //if(std::norm(px) < tol*tol)
        //    return xk;
        cdouble G = d / px;
        cdouble H = G*G - f / px;

        //std::cout << "G = " << G << "H = " << H;
        cdouble radicand = (n - 1)*(n*H-G*G);
        //assert(radicand.real() > 0);
        if(radicand.real() < 0)
            quad_root = true;
        //std::cout << "radicand = " << radicand << std::endl;
        if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
            a = - sqrt(radicand);
        else
            a = sqrt(radicand);
        //std::cout << "a = " << a << std::endl;
        a = n / (a + G);
        //std::cout << "a = " << a << std::endl;
        if(shuffle_counter % shuffle_rate == 0)
		{
			//a *= shuffle[shuffle_counter / shuffle_rate];
		}
        xk -= a;
        shuffle_counter++;
        if(shuffle_counter >= 90)
            break;
    }
    //std::cout << "xk = " << xk << std::endl;
    return xk;
}

Poly monic(Poly p)
{
  reduce(p);
  if (is_zero(p))
  {
    return p;
  }
  ll m = p.back();
  ll rev_m = mod_rev(m);
  return rev_m * p;
}

inline bool onEdge(const Point & point, const Poly & poly, const Edges & edges, const S & eps)
{
    for(auto j = edges.begin(), jend = edges.end(); j != jend; ++j)
    {
        auto & p1 = poly[*j];
        auto & p2 = *j == 0 ? poly.back() : poly[*j - 1];
        if(onEdge(point, p1, p2, eps))
            return true;
    }
    return false;
}

inline bool within(const Point & point, const Poly & poly, const Mapping & mapping, const S & ieps)
{
    auto i = std::lower_bound(mapping.begin(), mapping.end(), point.y, detail::CmpY<decltype(*mapping.begin())>());
    if(i == mapping.begin())
    {
        if(ieps < 0 || point.y < i->y - ieps)
            return false;
        return onEdge(point, poly, i->edges, ieps);
    }
    if(i == mapping.end())
    {
        if(ieps < 0)
            return false;
        --i;
        if(point.y > i->y + ieps)
            return false;
        if(!i->edges.empty() && onEdge(point, poly, i->edges, ieps))
            return true;
        --i;
        return onEdge(point, poly, i->edges, ieps);
    }
    S eps = std::abs(ieps);
    if(i->y - point.y < eps)
    {
        if(onEdge(point, poly, i->edges, eps))
            return ieps > 0;
    }
    --i;
    if(onEdge(point, poly, i->edges, eps))
        return ieps > 0;
    bool result = false;
    for(auto j = i->edges.begin(), jend = i->edges.end(); j != jend; ++j)
    {
        auto & p1 = poly[*j];
        auto & p2 = *j == 0 ? poly.back() : poly[*j - 1];
        auto d = p2.y - p1.y;
        if(std::abs(d) > eps)
        {
            auto t = (point.y - p1.y) / d;
            if(point.x + eps > p1.x + (p2.x - p1.x) * t)
                result = !result;
        }
    }
    return result;
}

bool contains(const Poly & poly, const Point & p)
{
    bool result = false;
    const auto * prev = &poly.back();
    for(auto i = poly.begin(), end = poly.end(); i != end; ++i)
    {
        const auto * cur = &*i;
        if(cur->y > prev->y)
        {
            auto dy = cur->y - prev->y;
            if(p.y >= prev->y && p.y < cur->y &&
               (p.x - prev->x) * dy <= (p.y - prev->y) * (cur->x - prev->x))
               result = !result;
        } else {
            auto dy = cur->y - prev->y;
            if(p.y < prev->y && p.y >= cur->y &&
               (p.x - prev->x) * dy >= (p.y - prev->y) * (cur->x - prev->x))
               result = !result;
        }
        prev = cur;
    }
    return result;
}

bool is_monic(const Poly &p)
{
  return (p.back() == 1);
}

// Simplify the polynomial, eliminates the high-degree zero term.
void tidy(Poly& c) {
    while(c.size() > 1 && fabs(c.back()) < EPS)
        c.pop_back();
}

// Simplify the polynomial, eliminates the high-degree zero term.
void tidy(Poly& c) {
    while(!c.empty() && c.back() == 0)
        c.pop_back();
}