static polynomial div_monomial(const polynomial& divident,
			const polynomial& divisor) {
		int_type degree = divident.degree() - divisor.degree();
		int_type coefficient = divident.data.begin()->second
				/ divisor.data.begin()->second;
		return { {degree, coefficient}};
	}

polynomial operator *(const polynomial &p1, const polynomial &p2)
{
    polynomial p3;
    Iterator<polyterm> poly1, poly2;
    Node<polyterm> newTerm;


    poly1=p1.Begin();
    while(poly1.GetNode()!=NULL)
    {

        poly2 = p2.Begin();
        while(poly2.GetNode()!=NULL)
        {
            newTerm = *poly1.GetNode()* *poly2.GetNode();

            p3.InsertSorted(new Node<polyterm>(newTerm));


            poly2.Next();
        }
        poly1.Next();
    }
    return p3;




}

 friend polynomial derivative(polynomial const& p) {
  int s = p.size();
  if(s == 0) return {};

  polynomial res;
  res.coeffs.resize(s-1);
  for(int n = 0; n < s-1; ++n)
   res.coeffs(n) = p.coeffs(n+1) * CoeffType(n+1);
  return res;
 }

 polynomial operator+(const polynomial<T>& o) const
 {
   int len;
   const polynomial<T> *p;
   size() >= o.size() ? (p=this, len=size()) : (p=&o, len=o.size());
   polynomial<T> out(len);
   std::transform(vals->data(),vals->data()+std::min(size(),o.size()),o.vals->data(),out.vals->data(),[](T a, T b){return a+b;});
   if (size() != o.size()) std::copy_n(&p->element(std::min(size(),o.size())),abs(size()-o.size()),&out(std::min(size(),o.size())));
   return out;
 }

    const polynomial operator+ (polynomial operand) const {
        std::vector<double> result (std::max(degree(), operand.degree()) + 1);

        for (int i = 0; i < result.size(); i++) {
            double x = i <= degree() ? coefficients[i] : 0.0;
            double y = i <= operand.degree() ? operand[i] : 0.0;
            result[i] = x + y;
        }

        return polynomial { result };
    }

 friend polynomial antiderivative(polynomial const& p) {
  int s = p.size();
  if(s == 0) return {};

  polynomial res;
  res.coeffs.resize(s+1);
  res.coeffs(0) = CoeffType{};
  for(int n = 1; n < s+1; ++n)
   res.coeffs(n) = p.coeffs(n-1) / CoeffType(n);
  return res;
 }

void printpoly(ostream &dst, const polynomial<num,Nvar,Ndeg> &p)
{
    int exps[Nvar] = {0};
    for (int i=0; i<p.size; i++)
    {
        p.exponents(i,exps);
        for (int j=0; j<Nvar; j++)
            cout << exps[j] << " ";
        cout << " " << p[i] << endl;
        assert(i == p.term_index(exps));
    }
}

std::vector<polynomial> sturm_chain (polynomial p) {
    std::vector<polynomial> chain (p.degree() + 1);

    chain[0] = p;
    chain[1] = p.derivative();

    for (int i = 2; i < chain.size(); i++) {
        chain[i] = -(chain[i-2] % chain[i-1]);
    }

    return chain;
}

void test_np(const polynomial& test)
{
    unsigned int exponent;

    cout << "Enter exponent: ";
    cin >> exponent;

    cout << "For polynomial: ";
    view(test);
    cout << "next_term(" << exponent << ") = "
	 << test.next_term(exponent) << endl;
    cout << "previous_term(" << exponent << ") = "
	 << test.previous_term(exponent) << endl;
}

    std::pair<polynomial, polynomial> quo_rem (polynomial divisor) const {
        /* Degrees of the remainder and quotient. */
        int kd = divisor.degree();
        int kr = kd - 1;
        int kq = degree() - kd;

        std::vector<double> remainder (kr + 1);
        std::vector<double> quotient (kq + 1);

        auto dbegin = divisor.coefficients.crbegin() + 1;
        auto dend = divisor.coefficients.crend();
        auto qbegin = quotient.end();
        auto qend = quotient.end();

        for (int i = kq; 0 <= i; i--) {
            quotient[i] = (coefficients[kd+i] - ip(dbegin, dend, qbegin, qend)) / divisor[kd];
            --qbegin;
        }

        for (int i = kr; 0 <= i; i--) {
            remainder[i] = coefficients[i] - ip(dbegin, dend, qbegin, qend);
            ++dbegin;
        }

        return { polynomial { quotient }, polynomial { remainder } };
    }

void polynomial::div_and_mod(const polynomial &p, 
		polynomial &result, polynomial &remainder) const {
	remainder = *this;

	while (!remainder.islower(p)) {
		double coe = remainder.beg->coe / p.beg->coe;
		double exp = remainder.beg->exp - p.beg->exp;
		if (coe == 0) break;

		polynomial ans;
		ans.append(coe, exp);
		result.append(coe, exp);

		remainder = remainder.minus(ans.multiply(p));
	}
}

        void sturm_functions(std::vector< polynomial<data__> > &sturm_fun, const polynomial<data__> &f)
        {
          size_t i, n;

          // resize the Sturm function vector
          n=f.degree()+1;
          sturm_fun.resize(n);

          polynomial<data__> *ptemp(f.fp()), pjunk;

          // initialize the Sturm functions
          sturm_fun[0]=f;
          ptemp=f.fp();
          sturm_fun[1]=*ptemp;
          delete ptemp;

          // build the Sturm functions
          sturm_fun[0].divide(std::abs(sturm_fun[0].coefficient(sturm_fun[0].degree())));
          sturm_fun[1].divide(std::abs(sturm_fun[1].coefficient(sturm_fun[1].degree())));
          for (i=2; i<n; ++i)
          {
            if (sturm_fun[i-1].degree()>0)
            {
              data__ small_no, max_c(0);
              typename polynomial<data__>::coefficient_type c;

              sturm_fun[i-1].get_coefficients(c);
              for (typename polynomial<data__>::index_type j=0; j<=sturm_fun[i-1].degree(); ++j)
              {
                if (std::abs(c(j))>max_c)
                  max_c=std::abs(c(j));
              }
              small_no=max_c*std::sqrt(std::numeric_limits<data__>::epsilon());

              pjunk.divide(sturm_fun[i], sturm_fun[i-2], sturm_fun[i-1]);
              sturm_fun[i].negative();
              if (std::abs(sturm_fun[i].coefficient(sturm_fun[i].degree()))>small_no)
                sturm_fun[i].divide(std::abs(sturm_fun[i].coefficient(sturm_fun[i].degree())));
              sturm_fun[i].adjust_zero(small_no);
            }
            else
            {
              sturm_fun[i]=0;
            }
          }
        }

void find_kernels(const polynomial& P, vector<pair<monomial, polynomial>>& kernelmap)
{
    kernelmap.clear();
    kernels(0, P, monomial(), kernelmap);
    if (P.gcd() == monomial())
    {
        kernelmap.push_back(make_pair(monomial(), P));
    }
}

        int descartes_rule(const polynomial<data__> &f, bool positive)
        {
          // catch special case
          if (f.degree()==0)
            return 0;

          // get the coefficients from the polynomial
          std::vector<data__> a(f.degree()+1);
          for (size_t i=0; i<a.size(); ++i)
          {
            a[i]=f.coefficient(i);

            // change the sign of odd powers if want the negative root count
            if (!positive && i%2==1)
              a[i]*=-1;
          }

          return eli::mutil::poly::root::sign_changes(a.begin(), a.end());
        }

polynomial operator + (const polynomial &p, const polynomial &q) {
	polynomial::const_iterator itP, itQ;
	polynomial C;
	itP = p.begin();
	itQ = q.begin();
	while(itP != p.end() && itQ != q.end()) {
		if ((*itP).second == (*itQ).second) {
			term ans(((*itP).first + (*itQ).first), (*itQ).second);
			if (((*itP).first + (*itQ).first) != 0) 
				C.push_back(ans);
			itP++;
			itQ++;
			
		} else if ((*itP).second < (*itQ).second) {
			C.push_back(*itQ);
			itQ++;
		} else {
			C.push_back(*itP);
			itP++;
		}
	}
	while (itP != p.end()) {
		C.push_back(*itP);
		itP++;
	}
	while (itQ != q.end()) {
		C.push_back(*itQ);
		itQ++;
	}


	return C;

}

	polynomial operator -(const polynomial& p1, const polynomial& p2)
	{
		size_t j=p1.degree( );
		size_t q= p2.degree( );
		polynomial answer;

		if(j>=q)
		{
			answer.reserve(j+1);
			size_t i;
			for (i=0; i<=j; i++)
			{
				answer.add_to_coef((p1.coefficient(i)-p2.coefficient(i)),(i));
			}
		}
		else
		{
			answer.reserve(q+1);
			size_t i;
			for (i=0; i<=q; i++)
			{
				answer.add_to_coef((p1.coefficient(i)-p2.coefficient(i)),(i));
			}
		}
		return answer;
	}

 void divide(const polynomial &b)
 {
     int m=b.size(),res[1110];
     for (int i=n;i>=m;i--)
     {
         int t=a[i]/b[m];
         for (int j=i;j>=i-m;j--)
             a[j]-=b[j-i+m]*t;
         res[i-m]=t;
     }
     n-=m;
     for (int i=0;i<=n;i++)
         a[i]=res[i];
 }

void test_assign(polynomial& test)
{
    double coefficient;
    unsigned int exponent;

    cout << "Enter exponent: ";
    cin >> exponent;
    cout << "Enter coefficient: ";
    cin >> coefficient;

    test.assign_coef(coefficient, exponent);
    cout << "After assigning: ";
    view(test);
}

void compute_alternant_error_locator (polynomial&syndrome, gf2m&fld,
                                      uint t, polynomial&out)
{
	if (syndrome.zero()) {
		//ensure no roots
		out.resize (1);
		out[0] = 1;
		return;
	}

	polynomial a, b;

	polynomial x2t; //should be x^2t
	x2t.clear();
	x2t.resize (1, 1);
	x2t.shift (2 * t);

	syndrome.ext_euclid (a, b, x2t, fld, t - 1);
	uint b0inv = fld.inv (b[0]);
	for (uint i = 0; i < b.size(); ++i) b[i] = fld.mult (b[i], b0inv);
	out = b;
	//we don't care about error evaluator
}

        int sturm_count(const polynomial<data__> &f, const data2__ &xmin, const data2__ &xmax)
        {
          // short circuit degenerate cases
          if (xmax<=xmin)
            return 0;

          // catch another degenerate case
          if (f.degree()==0)
            return 0;

          std::vector< polynomial<data__> > sturm_fun;

          // calculate the Sturm functions
          sturm_functions(sturm_fun, f);

          return sturm_count(sturm_fun, xmin, xmax);
        }