void cut(poly &p){
    bool flag = 1;
    while(flag && p.s() > eps){
        flag = 0;
        rep(i, n)
            if(cross(l[i], p)){
                poly *L = new poly();
                int st = 0, sz = p.ver.size();
                for(; !cross(line(p.ver[st], p.ver[(st + 1) % sz]), l[i]); st = (st + 1) % sz);
                L->ver.pb(crossPoint(line(p.ver[st], p.ver[(st + 1) % sz]), l[i]));
                for(st = (st + 1) % sz; !cross(line(p.ver[st], p.ver[(st + 1) % sz]), l[i]); st = (st + 1) % sz)
                    L->ver.pb(p.ver[st]);
                L->ver.pb(p.ver[st]);
                L->ver.pb(crossPoint(line(p.ver[st], p.ver[(st + 1) % sz]), l[i]));
                cut(*L);

                L->ver.clear();
                L->ver.pb(crossPoint(line(p.ver[st], p.ver[(st + 1) % sz]), l[i]));
                for(st = (st + 1) % sz; !cross(line(p.ver[st], p.ver[(st + 1) % sz]), l[i]); st = (st + 1) % sz)
                    L->ver.pb(p.ver[st]);
                L->ver.pb(p.ver[st]);
                L->ver.pb(crossPoint(line(p.ver[st], p.ver[(st + 1) % sz]), l[i]));
                p.ver = L->ver;

                delete(L);
                flag = 1;
            }
    }
    double s = p.s();
    if(s > eps && s < inf)
        ans.pb(s);
}

poly mult(poly& p1,poly& p2){
	poly ret(p1.size()+p2.size()-1,0);
	for(int i=0;i<p1.size();i++){
		for(int j=0;j<p2.size();j++){
			ret[i+j]+=p1[i]*p2[j];
		}
	}
	return ret;
}

double bisect(const poly &p, double lower, double upper) {
  if (upper - lower < 1e-5) return lower;
  double mid = (lower + upper) / 2;
  if (sgn(p.eval(mid)) == sgn(p.eval(lower))) {
    return bisect(p, mid, upper);
  } else {
    return bisect(p, lower, mid);
  }
}

// Test whether point a is inside poly p
bool is_inside(const point & a, const poly & p) {
	unsigned i;
	int n = p.size();
	bool tmp, b = false;
	for (i = 0; i < p.size(); i++) {
		tmp = intersect(a, e, p[i], p[(i + 1) % n]);
		b = (b != tmp);
	}
	return b;
}

poly operator*(const poly &a, const poly &b)
{
    poly ret(max(0, (int)a.size() + (int)b.size() - 1), 0.0);
    for (int i = 0; i < a.size(); i++) {
        for (int j = 0; j < b.size(); j++) {
            ret[i + j] += a[i] * b[j];
        }
    }
    return ret;
}

void								MeshData::GenerateNormals(
	const bool&							InvertNormals)
{
	//	Sources
	const AttribData<vec3> & Vertices	= m_VertexArray;
	const AttribData<poly> & Faces		= m_Polygons;

	vec3* pNormals = new vec3[Vertices.Size()];
	memset(pNormals, 0, Vertices.Size() * sizeof(vec3));

	//	Go through the polylists
	const uint NrPolyLists = Faces.Size();
	for(uint p=0; p<NrPolyLists; p++) {

		const poly Polylist = Faces.GetElem(p);
		const uint NrFaces = Polylist.NrFaces();
		for(uint f=0; f<NrFaces; f++)
		{
			//	Get vertices
			const uint i1 = Polylist.GetFaces().GetElem(f).v1;
			const uint i2 = Polylist.GetFaces().GetElem(f).v2;
			const uint i3 = Polylist.GetFaces().GetElem(f).v3;

			const vec3 v1 = Vertices.GetElem(i1);
			const vec3 v2 = Vertices.GetElem(i2);
			const vec3 v3 = Vertices.GetElem(i3);

			//	Get edges
			const vec3 edge1 = v2-v1; const vec3 edge2 = v3-v1;

			//	Set Crossproduct
			vec3 facenormal;
			if(InvertNormals) {
				facenormal=normalize(glm::cross(edge1, edge2));
				facenormal = facenormal * -1.0f;
			} else
				facenormal=normalize(glm::cross(edge1, edge2));

			//	Add vertexnormals
			pNormals[i1]+=facenormal;
			pNormals[i2]+=facenormal;
			pNormals[i3]+=facenormal;
		}
	}

	//	Normalize
	vec3* pN = pNormals;
	for(uint n=0; n<Vertices.Size(); n++, pN++)
		*pN = normalize(*pN);

	MeshData::SetNormalArray().Load(pNormals, Vertices.Size());

	delete [] pNormals;
}

poly operator+(const poly &a, const poly &b)
{
    poly ret = a;
    if (ret.size() < b.size()) {
        ret.resize(b.size(), 0.0);
    }
    for (int i = 0; i < b.size(); i++) {
        ret[i] += b[i];
    }
    return ret;
}

void normalize(poly &p)
{
    int minx=1000 , miny=1000 ;
    for(int i=0;i<p.size();i++)
    {
        minx=min(minx,p[i].x) ;
        miny=min(miny,p[i].y) ;
    }
    poly q ;
    for(int i=0;i<p.size();i++)
        q.push_back((P){p[i].x-minx,p[i].y-miny}) ;
    p=q ;
}

void read( poly& pp , string ss , LL mod ){
  pp.clear();
  if( ss[ 0 ] != '-' ) ss = '+' + ss;
  size_t st = 0;
  while( st < ss.length() ){
    size_t lst = st;
    st ++;
    while( st < ss.length() and
           ss[ st ] != '+' and
           ss[ st ] != '-' ) st ++;
    pp.push_back( parse( ss , lst , st , mod ) );
  }
}

poly add(poly p1,poly p2){
	if(p1.size()>=p2.size()){
		for(int i=0;i<p2.size();i++){
			p1[i]+=p2[i];
		}
		return p1;
	}
	else{
		for(int i=0;i<p1.size();i++){
			p2[i]+=p1[i];
		}
		return p2;
	}
}

poly graham(poly p){
	int i,j,n = p.size();
	poly g;
	pivo = *min_element(p.begin(), p.end(), cmp_pivo);
	sort( p.begin(), p.end(), cmp_radial);
	// n previous declared
	for( i=j=0; i<n;i++){
		while( j>=2 && ccw( g[j-2] , g[j-1], p[i]) >=0){
			g.pop_back(); j--;
		}
		g.push_back(p[i]); j++;
	}
	return g;
}

inline poly subtract(poly &a, poly &b){
    int N = (int)a.size(), M = (int)b.size();
    int K = min(N, M);

    poly ans;
    if (N == 0 && M == 0)return ans;
    ans.resize(max(N, M));

    for (int i = 0; i < K; i++)ans[i] = a[i] - b[i];
    for (int i = K; i < N; i++)ans[i] = a[i];
    for (int i = K; i < M; i++)ans[i] = -b[i];

    reduce(ans);
    return ans;
}

std::pair<poly, cdouble> horner(const poly & a, const cdouble & x) {
  int n = a.size();
  poly b = poly(std::max(1, n - 1));
  for (int i = n - 1; i > 0; i--)
    b[i - 1] = a[i] + (i < n - 1 ? b[i] * x : 0);
  return std::make_pair(b, a[0] + b[0] * x);
}

inline poly multiply(poly &a, poly &b){
    int N = (int)a.size(), M = (int)b.size();

    poly ans;
    if (N == 0 || M == 0)return ans;
    ans.assign(N + M - 1, 0);

    for (int i = 0; i < N; i++){
        for (int j = 0; j < M; j++){
            ans[i + j] += a[i] * b[j];
        }
    }
    reduce(ans);

    return ans;
}

poly hull( poly p ) {
	int k = 0, n = p.size();
	sort( p.begin(), p.end());
	poly h( 2 * n );
	for(int i = 0; i < n; i++){
		while( k >= 2 && ccw( h[k-2], h[k-1], p[i] ) <= 0 ) k--;
		h[ k++ ] = p[ i ];
	}
	int tmp = k + 1;
	for(int i = n - 2; i >= 0; i--){
		while( k >= tmp && ccw( h[k-2], h[k-1], p[i] ) <= 0 ) k--;
		h[ k++ ] = p[ i ];
	}
	h.resize( k );
	return h;
}

poly derivative(const poly & p) {
  int n = p.size();
  poly r(std::max(1, n - 1));
  for(int i = 1; i < n; i++)
    r[i - 1] = p[i] * cdouble(i);
  return r;
}

double dist( poly p ){
	double res = 0;
	for(int i = 0; i < p.size() - 1; i++){
		res += dist( p[i], p[i + 1]);
	}
	return res;
}

ld integrate(poly p){
	ld res=0;
	for(int i=0;i<p.size();i++){
		res+=(p[i]/(i+1));
	}
	return res;
}

double eval(const poly &p, double x)
{
    double ret = 0;
    for (int i = (int)p.size() - 1; i >= 0; i--) {
        ret = ret * x + p[i];
    }
    return ret;
}

void rotat(poly &p)
{
    poly q ;
    for(int i=0;i<p.size();i++)
        q.push_back((P){-p[i].y,p[i].x}) ;
    normalize(q) ;
    p=q ;
}