--- a/include/geometry/kinetic-sort.h	Fri Aug 17 15:08:27 2007 -0400
+++ b/include/geometry/kinetic-sort.h	Fri Aug 24 16:58:25 2007 -0400
@@ -18,7 +18,7 @@
  *  \arg Swap_                is called with an ElementIterator_ when a swap needs to be performed
  */
 template<class ElementIterator_, class TrajectoryExtractor_, 
-		 class Simulator_, class Swap_ = boost::function<void(ElementIterator_ pos)>>
+		 class Simulator_, class Swap_ = boost::function<void(ElementIterator_ pos, Simulator_* simulator)> >
 class KineticSort
 {
 	public:
@@ -50,8 +50,9 @@
 
 		/// \name Core Functionality
 		/// @{
-									KineticSort(Swap swap);
+									KineticSort();
 									KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
+		void						initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator);
 
 		void						insert(iterator pos, ElementIterator f, ElementIterator l, Simulator* simulator);
 		void						erase(iterator pos, Simulator* simulator);
@@ -60,13 +61,11 @@
 		void						swap(iterator pos, Simulator* simulator);
 
 		bool						audit(Simulator* simulator) const;
+		/// @}
 
 		iterator					begin() 									{ return list_.begin(); }
 		iterator					end() 										{ return list_.end(); }
 
-		void						initialize(ElementIterator b, ElementIterator e, Simulator* simulator);
-		/// @}
-
 	private:
 		class SwapEvent;
 		void						schedule_swaps(iterator b, iterator e, Simulator* s);

--- a/include/geometry/kinetic-sort.hpp	Fri Aug 17 15:08:27 2007 -0400
+++ b/include/geometry/kinetic-sort.hpp	Fri Aug 24 16:58:25 2007 -0400
@@ -1,22 +1,21 @@
 template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
-KineticSort(Swap swap):
-	swap_(swap)
+KineticSort()
 {}	
 
 template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
-KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator):
-	swap_(swap)
+KineticSort(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
 {
-	initialize(b, e, simulator);
+	initialize(b, e, swap, simulator);
 }
 
 template<class ElementIterator_, class TrajectoryExtractor_, class Simulator_, class Swap_>
 void
 KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
-initialize(ElementIterator b, ElementIterator e,, Simulator* simulator)
+initialize(ElementIterator b, ElementIterator e, Swap swap, Simulator* simulator)
 {
+	swap_ = swap;
 	for (ElementIterator cur = b; cur != e; ++cur)
 		list_.push_back(Node(cur, simulator->null_key()));
 	schedule_swaps(list_.begin(), list_.end(), simulator);
@@ -80,9 +79,9 @@
 KineticSort<ElementIterator_, TrajectoryExtractor_, Simulator_, Swap_>::
 swap(iterator pos, Simulator* simulator)
 {
-	AssertMsg(boost::next(pos) != list_.end(), "Cannot swap the last element");
+	// TODO: AssertMsg(boost::next(pos) != list_.end(), "Cannot swap the last element");
 
-	swap(pos->element);
+	swap_(pos->element, simulator);
 
 	// Remove events
 	iterator prev = boost::prior(pos);

--- a/include/geometry/polynomial.h	Fri Aug 17 15:08:27 2007 -0400
+++ b/include/geometry/polynomial.h	Fri Aug 24 16:58:25 2007 -0400
@@ -8,6 +8,7 @@
 #include <synaps/usolve/bezier/SlvBzStd.h>
 //#include <synaps/usolve/Sturm.h>
 //#include <synaps/arithm/Infinity.h>
+#include <synaps/arithm/gmp.h>
 
 #include <stack>
 
@@ -30,20 +31,23 @@
 		static RootType				root(const T& r)									{ return SynapsTraits<T>::root(r); }
 		static int					sign_at(const RationalFunction& rf, const RootType& r);
 		static RootType				between(const RootType& r1, const RootType& r2)		{ return SynapsTraits<T>::between(r1,r2); }
+		static bool					sign_at_negative_infinity(const RationalFunction& rf);
 };
 
 template<class T>
 struct SynapsTraits				///< Suitable for double
 {
 		typedef						T													CoefficientType;
+		typedef						T													RootType;
 		typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
-		typedef						SYNAPS::SlvBzStd<CoefficientType>					Solver;
-		typedef						T													RootType;
+		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<>
@@ -51,7 +55,8 @@
 {
 		typedef						ZZ													CoefficientType;
 		typedef						SYNAPS::UPolDse<CoefficientType>					Polynomial;
-		typedef						SYNAPS::Algebraic<CoefficientType>					Solver;
+		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);}
@@ -59,6 +64,24 @@
 		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"

--- a/include/geometry/polynomial.hpp	Fri Aug 17 15:08:27 2007 -0400
+++ b/include/geometry/polynomial.hpp	Fri Aug 24 16:58:25 2007 -0400
@@ -5,8 +5,8 @@
 {
 	typedef SYNAPS::Seq<RootType>		RootSeq;
 
-	RootSeq seq_num = SYNAPS::solve(rf.numerator(), Solver());
-	RootSeq seq_den = SYNAPS::solve(rf.denominator(), Solver());
+	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
@@ -32,3 +32,40 @@
 {
 	return SynapsTraits<T>::sign_at(rf.numerator(), r) * SynapsTraits<T>::sign_at(rf.denominator(), r);
 }
+
+template<class T>
+bool
+UPolynomial<T>::
+sign_at_negative_infinity(const RationalFunction& rf)
+{
+	const Polynomial& num = rf.numerator();
+	const Polynomial& den = rf.denominator();
+	unsigned int ndegree = num.get_degree();
+	unsigned int ddegree = den.get_degree();
+	return (((ndegree + 1) % 2 == 0) ^ (num[ndegree] > 0)) ^
+		   (((ddegree + 1) % 2 == 0) ^ (den[ddegree] > 0));
+}
+
+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;
+}
+

--- a/include/geometry/simulator.h	Fri Aug 17 15:08:27 2007 -0400
+++ b/include/geometry/simulator.h	Fri Aug 24 16:58:25 2007 -0400
@@ -5,8 +5,14 @@
 #include "polynomial.h"
 
 template<class Comparison>  						class IndirectComparison;
-template<class PolyKernel_, class Simulator_>		class Event;
 
+/**
+ * Simulator class. Keeps a queue of events. Infinity is reached if the Event 
+ * at the front of the queue has an empty root stack. Keeps track of current time, 
+ * Event addition, and processes events one by one. Degeneracies are handled by 
+ * assuming that the RationalFunction responsible for the event must be positive
+ * before the Event occurs.
+ */
 template<class PolyKernel_, template<class Event> class EventComparison_ = std::less>
 class Simulator
 {
@@ -56,6 +62,12 @@
 		bool						reached_infinity_;
 };
 
+
+/**
+ * Base class for events. Stores a root stack, subclasses need to define process(). 
+ * Event with an empty root stack compares greater than any other Event, 
+ * pushing those events to the end of the queue.
+ */
 template<class PolyKernel_, template<class Event> class EventComparison_>
 class Simulator<PolyKernel_, EventComparison_>::Event
 {
@@ -80,15 +92,16 @@
 				return root_stack().top() < e.root_stack().top();
 		}
 
-		int							sign_before() const							{ return root_stack().top().sign_low(); }
-		int							sign_after() const							{ return root_stack().top().sign_high(); }
-
 		virtual std::ostream&		print(std::ostream& out) const				{ return out << "Event with " << root_stack_.size() << " roots"; }
 
 	private:
 		RootStack					root_stack_;
 };
 
+/**
+ * Compares elements pointed at by its arguments using the provided Comparison_ 
+ * (which must not take any arguments during construction).
+ */
 template<class Comparison_>
 class IndirectComparison: public std::binary_function<const typename Comparison_::first_argument_type*, 
 													  const typename Comparison_::second_argument_type*, 

--- a/include/geometry/simulator.hpp	Fri Aug 17 15:08:27 2007 -0400
+++ b/include/geometry/simulator.hpp	Fri Aug 24 16:58:25 2007 -0400
@@ -17,8 +17,13 @@
 	Event* ee = new Event_(e);
 	//std::cout << "Solving: " << f << std::endl;
 	PolynomialKernel::solve(f, ee->root_stack());
+	bool sign = PolynomialKernel::sign_at_negative_infinity(f);
 	while (!ee->root_stack().empty() && ee->root_stack().top() < current_time())
+	{
 		ee->root_stack().pop();
+		sign = !sign;
+	}
+	if (sign) ee->root_stack().pop();			// TODO: double-check the logic
 	//std::cout << "Pushing: " << ee->root_stack().top() << std::endl;
 	return queue_.push(ee);
 }

--- a/tests/geometry/test-kinetic-sort.cpp	Fri Aug 17 15:08:27 2007 -0400
+++ b/tests/geometry/test-kinetic-sort.cpp	Fri Aug 24 16:58:25 2007 -0400
@@ -4,10 +4,11 @@
 #include <iostream>
 
 #include <boost/utility.hpp>
+#include <boost/bind.hpp>
 
-//typedef		double							FieldType;
+typedef		double							FieldType;
 //typedef		ZZ								FieldType;
-typedef		QQ								FieldType;
+//typedef		QQ								FieldType;
 typedef 	UPolynomial<FieldType>			PolyKernel;
 typedef		PolyKernel::Polynomial			Polynomial;
 typedef 	std::list<Polynomial>			SortDS;
@@ -30,7 +31,7 @@
 	FieldType					vv;
 };
 
-void swap(SortDS* s, SortDSIterator i)
+void swap(SortDS* s, SortDSIterator i, SimulatorFT* simulator)
 {
 	std::cout << "Swapping " << *i << " " << *boost::next(i) << std::endl;
 	s->splice(i, *s, boost::next(i));
@@ -45,6 +46,9 @@
 	list.push_back(Polynomial("x^3 - 3"));
 	list.push_back(Polynomial("x^2 - 2*x - 2"));
 	list.push_back(Polynomial("2*x - 4"));
+	list.push_back(Polynomial("x"));
+	list.push_back(Polynomial("-x + 4"));
+	list.push_back(Polynomial("2"));
 	//list.sort(EvaluatedComparison(simulator.current_time()));
 	list.sort(EvaluatedComparison(0));
 
@@ -53,15 +57,14 @@
 		std::cout << *cur << std::endl;
 
 	// Setup kinetic sort
-	KineticSortDS	ks(list.begin(), list.end(), std::bind1st(&swap, &list), &simulator);
+	KineticSortDS	ks(list.begin(), list.end(), boost::bind(swap, &list, _1, _2), &simulator);
 
-	while(!simulator.reached_infinity() && simulator.current_time() < 1)
+	while(!simulator.reached_infinity() && simulator.current_time() < 4)
 	{
-		//std::cout << "Current time before: " << simulator.current_time() << std::endl;
-		if (!ks.audit(&simulator)) return 1;
+		std::cout << "Current time before: " << simulator.current_time() << std::endl;
+		//if (!ks.audit(&simulator)) return 1;
 		//simulator.print(std::cout << "Auditing ");
 		simulator.process();
-		//std::cout << "Current time after: " << simulator.current_time() << std::endl;
+		std::cout << "Current time after: " << simulator.current_time() << std::endl;
 	}
-
 }