include/geometry/polynomial.hpp
 author Dmitriy Morozov Tue, 26 Feb 2008 16:45:25 -0500 branch ar changeset 72 35500d7f9fca parent 71 a6c5cb5a17cc child 87 2c2e2f3b5d15 permissions -rw-r--r--
KineticSort can deal with equal trajectories: UPolynomial::sign_at_negative_infinity() can handle 0 and so can Simulator::add()

#include <utilities/log.h>

template<class T>
void
UPolynomial<T>::
solve(const RationalFunction& rf, RootStack& stack)
{
typedef SYNAPS::Seq<RootType>		RootSeq;

AssertMsg(stack.empty(), "Stack must be empty before solve");

RootSeq seq_num = SYNAPS::solve(SynapsTraits<T>::convert(rf.numerator()), Solver());
RootSeq seq_den = SYNAPS::solve(SynapsTraits<T>::convert(rf.denominator()), Solver());

// TODO: assert that all roots in seq_den have positive multiplicity
// TODO: deal with multiplicities for the numerator
for (typename RootSeq::const_reverse_iterator cur = seq_num.rbegin();
cur != seq_num.rend();
++cur)
{
if (SynapsTraits<T>::multiplicity(*cur) % 2 != 0)
{
if (!stack.empty() && stack.top() == *cur)			// get rid of even multiplicities
// TODO: add logging information for this
stack.pop();
else
stack.push(*cur);
}
}
}

template<class T>
int
UPolynomial<T>::
sign_at(const RationalFunction& rf, const RootType& r)
{
return SynapsTraits<T>::sign_at(rf.numerator(), r) * SynapsTraits<T>::sign_at(rf.denominator(), r);
}

template<class T>
int
UPolynomial<T>::
sign_at_negative_infinity(const RationalFunction& rf)
{
const Polynomial& num = rf.numerator();
const Polynomial& den = rf.denominator();
int ndegree = num.get_degree();
int ddegree = den.get_degree();
if (ndegree == -1) return 0;          // ndegree == -1 => num == 0, and 0 is 0 at -infinity
return !((((ndegree + 1) % 2 == 0) ^ (num[ndegree] > 0)) ^
(((ddegree + 1) % 2 == 0) ^ (den[ddegree] > 0))) ? 1:-1;
}

SynapsTraits<QQ>::RootType
SynapsTraits<QQ>::
root(const CoefficientType& r)
{
ZZ p[2] = { -SYNAPS::numerator(r), SYNAPS::denominator(r) };
return SYNAPS::solve(SolverPolynomial(2, p), Solver(), 0);
}

SynapsTraits<QQ>::SolverPolynomial
SynapsTraits<QQ>::
convert(const Polynomial& p)
{
SolverPolynomial result(1, p.size() - 1);
assert(result.size() == p.size());
ZZ denominator_product = 1;
for (unsigned int i = 0; i < p.size(); ++i)
denominator_product *= SYNAPS::denominator(p[i]);
for (unsigned int i = 0; i < p.size(); ++i)
result[i] = SYNAPS::numerator(p[i]) * denominator_product /
SYNAPS::denominator(p[i]);
return result;
}