include/geometry/polynomial.h
 author Dmitriy Morozov Tue, 02 Feb 2010 23:43:46 -0800 branch dev changeset 193 7198c9bf279d parent 87 2c2e2f3b5d15 permissions -rw-r--r--
StaticPersistence::pair_simplices() has an option to not display progress (LSVineyard uses it during consistency checking)

#ifndef __POLYNOMIAL_H__
#define __POLYNOMIAL_H__

#include <synaps/upol.h>
#include <synaps/upol/gcd.h>
#include <synaps/usolve/Algebraic.h>
#include <synaps/usolve/algebraic/sign_fct.h>
#include <synaps/usolve/bezier/SlvBzStd.h>
//#include <synaps/usolve/Sturm.h>
//#include <synaps/arithm/Infinity.h>
#include <synaps/arithm/gmp.h>

#include <stack>

#include "rational-function.h"

template<class T> class SynapsTraits;

template<class T>
class UPolynomial
{
public:
typedef						typename SynapsTraits<T>::Polynomial				Polynomial;
typedef						RationalFunction<Polynomial>						RationalFunctionP;
typedef                     RationalFunctionP                                   Function;

typedef						typename SynapsTraits<T>::Solver					Solver;
typedef						typename SynapsTraits<T>::RootType					RootType;
typedef						std::stack<RootType>								RootStack;

static void					solve(const RationalFunctionP& rf, RootStack& stack);
static RootType				root(const T& r)									{ return SynapsTraits<T>::root(r); }
static int					sign_at(const RationalFunctionP& rf, const RootType& r);
static RootType				between(const RootType& r1, const RootType& r2)		{ return SynapsTraits<T>::between(r1,r2); }
static int					sign_at_negative_infinity(const RationalFunctionP& rf);
};

template<class T>
struct SynapsTraits				///< Suitable for double
{
typedef						T													CoefficientType;
typedef						T													RootType;
typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
typedef						SYNAPS::SlvBzStd<RootType, CoefficientType>			Solver;
typedef						Polynomial											SolverPolynomial;

static RootType				root(CoefficientType r)								{ return r; }
static unsigned int			multiplicity(RootType r)							{ return 1; }
static int					sign_at(const Polynomial& p, RootType r)			{ return SYNAPS::UPOLDAR::sign_at(p, r); }
static RootType				between(RootType r1, RootType r2)					{ return (r1 + r2)/2; }
static SolverPolynomial		convert(const Polynomial& p)						{ return p; }
};

template<>
struct SynapsTraits<ZZ>
{
typedef						ZZ													CoefficientType;
typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
typedef						Polynomial											SolverPolynomial;
typedef						SYNAPS::Algebraic<ZZ>								Solver;
typedef						Solver::root_t										RootType;

static RootType				root(const CoefficientType& r)						{ CoefficientType p[2] = {-r, 1}; return SYNAPS::solve(Polynomial(2, p), Solver(), 0);}
static unsigned int			multiplicity(const RootType& r)						{ return r.multiplicity(); }
static int					sign_at(const Polynomial& p,
const RootType& r)							{ return SYNAPS::ALGEBRAIC::sign_at(p, r); }
//static RootType				between(const RootType& r1, const RootType& r2) 	{ RootType r = r1; r += r2; r /= root(2); return r; }
static SolverPolynomial		convert(const Polynomial& p)						{ return p; }
};

template<>
struct SynapsTraits<QQ>
{
typedef						QQ													CoefficientType;
typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
typedef						SYNAPS::Algebraic<ZZ>								Solver;
typedef						SYNAPS::UPolDse<ZZ>									SolverPolynomial;
typedef						Solver::root_t										RootType;

static RootType				root(const CoefficientType& r);
static unsigned int			multiplicity(const RootType& r)						{ return r.multiplicity(); }
static int					sign_at(const Polynomial& p,
const RootType& r)							{ return SYNAPS::ALGEBRAIC::sign_at(convert(p), r); }
//static RootType				between(const RootType& r1, const RootType& r2) 	{ RootType r = r1; r += r2; r /= root(2); return r; }
static SolverPolynomial		convert(const Polynomial& p);
};

#include "polynomial.hpp"

#endif