Moved parts from EXACUS-AlciX to unify with ConiX-methods

.gitattributes

@@ -258,6 +258,13 @@ Algebraic_kernel_d/doc_tex/Algebraic_kernel_d_ref/intro.tex -text
 Algebraic_kernel_d/dont_submit -text
 Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2.h -text
 Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Algebraic_real_traits.h -text
+Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Approximate_arithmetic_controller.h -text
+Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Best_approximation_cache.h -text
+Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_bfs.h -text
+Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_traits_on_vert_line.h -text
+Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Non_generic_position_exception.h -text
+Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/alg_real_utils.h -text
+Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/enums.h -text
 Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Curve_analysis_2.h -text
 Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Curve_pair_analysis_2.h -text
 Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/LRU_hashed_map.h -text

Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Approximate_arithmetic_controller.h

@@ -0,0 +1,229 @@
+// TODO: Add licence
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL:$
+// $Id: $
+// 
+//
+// Author(s)     : Michael Kerber <mkerber@mpi-inf.mpg.de>
+//                 Ralf Schindlbeck
+//
+// ============================================================================
+
+
+#ifndef CGAL_APPROXIMATE_ARITHMETIC_CONTROLLER
+#define CGAL_APPROXIMATE_ARITHMETIC_CONTROLLER 1
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+template<typename Poly_1_,typename AlgebraicReal>
+class Approximate_arithmetic_controller {
+
+public:
+
+    typedef Poly_1_ Poly_1;
+    typedef AlgebraicReal Algebraic_real;
+ 
+    typedef typename Algebraic_real::Rational Boundary;
+
+    typedef typename Poly_1::NT Coefficient;
+
+    typedef typename boost::numeric::interval<Boundary> Boundary_interval;
+
+    typedef typename CGAL::Coercion_traits<Coefficient, Boundary>::Type 
+    Eval_result;
+
+    typedef typename boost::numeric::interval<Eval_result> Eval_interval;
+
+    //Default constructor
+    Approximate_arithmetic_controller() {}
+
+    // Constrcutor with Polynomial and Algebraic real
+    Approximate_arithmetic_controller(Poly_1 f, Algebraic_real alpha)
+        : _m_f(f), _m_alpha(alpha)
+    {
+#if AcX_USE_DERIVATIVE_OPTION || AcX_USE_BISECTION_OPTION
+        f_x = CGAL::diff(_m_f);
+        flag = false;
+#endif
+
+#if AcX_USE_BISECTION_OPTION
+        flag_r = false;
+#endif
+		iter_refine_steps = 0;
+	}
+
+    //! Obains an interval approximation of f(alpha)
+    Eval_interval interval_approximation() const {
+        Boundary_interval interval(_m_alpha.low(),_m_alpha.high());
+	#if AcX_USE_DERIVATIVE_OPTION || AcX_USE_BISECTION_OPTION
+		//! AcX_USE_DERIVATIVE_OPTION On or AcX_USE_BISECTION_OPTION On;
+		//! flag request, because to save runtime
+		if(!flag)
+		{
+			interval_f_x = evaluate_iv(f_x, interval);
+			CGAL::Sign sign_f_x = CGAL::sign(interval_f_x.lower()) * CGAL::sign(interval_f_x.upper());
+			//! find the sign of the derived function _m_f on their bounds
+			sign_low = f_x.sign_at(_m_alpha.low());
+
+			if(sign_f_x > 0)
+			{
+				//! f_x have no horizontal tangent => interval good ;-)
+				flag = true;
+				//std::cout << "function f is strictly monotone" << std::endl;
+
+			#if AcX_USE_BISECTION_OPTION
+				// help-variables, needed to save the old/new state of the intervals
+				low = _m_alpha.low();
+				high = _m_alpha.high();
+				low_f = _m_f.evaluate(_m_alpha.low());
+				high_f = _m_f.evaluate(_m_alpha.high());
+			#endif
+			}
+			else
+			{
+				//! f_x have a horizontal tangent => interval was overestimated
+				//std::cout << "function f is not monotone" << std::endl;
+				interval_f = evaluate_iv( _m_f, interval );
+				
+			}
+		}
+		if(flag)
+		{
+			//! computing direct, without IA => because the interval isn't
+			//! overestimated
+		#if !AcX_USE_BISECTION_OPTION
+			if(sign_low == CGAL::POSITIVE)
+			{
+				//! slope of _m_f is positive => interval OK
+				interval_f = Eval_interval( _m_f.evaluate(_m_alpha.low()), _m_f.evaluate(_m_alpha.high()) );
+			}
+			else
+			{
+				//! slope of _m_f is negative => switch interval-boundaries
+				interval_f = Eval_interval( _m_f.evaluate(_m_alpha.high()), _m_f.evaluate(_m_alpha.low()) );
+			}
+		#else
+			if(!flag_r)
+			{
+				if(sign_low == CGAL::POSITIVE)
+				{
+					//! slope of _m_f is positive => interval OK
+					interval_f = Eval_interval( low_f, high_f );
+				}
+				else
+				{
+					//! slope of _m_f is negative => switch interval-boundaries
+					interval_f = Eval_interval( high_f, low_f );
+				}
+			}
+			else
+			{
+				//! request if the low boundary hasn't change after the refinement
+				if(low == low_r)
+				{
+					high = high_r;
+					high_f = _m_f.evaluate(_m_alpha.high());
+					if(sign_low == CGAL::POSITIVE)
+					{
+						interval_f = Eval_interval( low_f , high_f );
+					}
+					else
+					{
+						interval_f = Eval_interval( high_f, low_f  );
+					}
+				}
+				//! request if the high boundary hasn't change after the refinement
+				else if(high == high_r)
+				{
+					low = low_r;
+					low_f = _m_f.evaluate(_m_alpha.low());
+					if(sign_low == CGAL::POSITIVE)
+					{
+						interval_f = Eval_interval( low_f , high_f );
+					}
+					else
+					{
+						interval_f = Eval_interval( high_f, low_f  );
+					}
+				}
+				else // for quadratic refinement method
+				{
+					//std::cout << "Refinement method without modified option, step: " << i << std::endl;
+					low = low_r;
+					high = high_r;
+					low_f = _m_f.evaluate(_m_alpha.low());
+					high_f = _m_f.evaluate(_m_alpha.high());
+					if(sign_low == CGAL::POSITIVE)
+					{
+						interval_f = Eval_interval( low_f , high_f );
+					}
+					else
+					{
+						interval_f = Eval_interval( high_f, low_f  );
+					}
+				}
+			}
+		#endif
+		}
+		return interval_f;
+	#else
+		//! AcX_USE_DERIVATIVE_OPTION Off or AcX_USE_BISECTION_OPTION Off;
+		return evaluate_iv( _m_f, interval );
+	#endif
+    }
+
+    //! Refines the algebraic real
+    int refine_value() const {
+        _m_alpha.refine();
+	#if AcX_USE_BISECTION_OPTION
+		// help-variables, needed to save the new intervals after refinement
+		flag_r = true;
+		low_r = _m_alpha.low();
+		high_r = _m_alpha.high();
+	#endif
+		iter_refine_steps++;
+		return iter_refine_steps;
+    }
+
+
+private:
+
+    Poly_1 _m_f;
+    Algebraic_real _m_alpha;
+	//! helpvariables for the "derivate"  and "bisection" option
+#if AcX_USE_DERIVATIVE_OPTION || AcX_USE_BISECTION_OPTION
+	// variable: indicate the first dervative of polynome f
+    Poly_1 f_x;
+	// variable: denotes the sign of the low boundary of the Algebraic_real alpha 
+    mutable CGAL::Sign sign_low;
+	// variable: terms the intervals for the poylnome f and its first derivative
+    mutable Eval_interval interval_f, interval_f_x;
+	// helpvariable: flag for indication, if f is strictly monotone in the boundary
+    mutable bool flag;
+#endif
+	//! helpvariables only for the "bisection" option
+#if AcX_USE_BISECTION_OPTION
+	// variable: denotes the old state of the maximum borders of alpha
+    mutable Boundary low, high;
+	// helpvariable: flag_r for indication, if one refinement step is made
+    mutable bool flag_r;
+	// variable: terms the new state of the maximum borders of alpha
+    mutable Boundary low_r, high_r;
+	// variable: denotes the state of the maximum borders of the polynome f
+    mutable Boundary low_f, high_f;
+#endif
+	// variable, which counts the refinement steps
+    mutable int iter_refine_steps;
+
+}; //class Approximate_arithmetic_controller
+
+} // namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif

Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Best_approximation_cache.h

@@ -0,0 +1,122 @@
+// TODO: Add licence
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL:$
+// $Id: $
+// 
+//
+// Author(s)     : Michael Kerber <mkerber@mpi-inf.mpg.de>
+//
+// ============================================================================
+
+#ifndef CGAL_BEST_APPROXIMATION_CACHE
+#define CGAL_BEST_APPROXIMATION_CACHE 1
+
+
+#include <CGAL/basic.h>
+#include <CGAL/Handle_with_policy.h>
+#include <map>
+#include <vector>
+
+CGAL_BEGIN_NAMESPACE
+
+
+  namespace CGALi {
+
+    /*!
+     * \brief Stores computed approximations for coefficients.
+     *
+     * During the execution of the Bitstream Descartes method, approximations
+     * of the coefficients are computed independently by different branches of
+     * the Descartes tree. To avoid computing them several times, the best
+     * konwn approximation for a special coefficient is stored in a \c std::map
+     * structure. The data contains the used precision and the integer that
+     * represents the coefficient for this precision, according to the 
+     * \c CGAL::CGALi::Bitstream_descartes_rndl_tree_traits definition.
+     */
+    template<typename Coefficient_,typename Integer_>
+      class Best_approximation_cache 
+      : public ::CGAL::Handle_with_policy<std::map<Coefficient_,std::pair<long,Integer_> > > {
+
+      public:
+
+      //! The Coefficient type
+      typedef Coefficient_ Coefficient;
+
+      //! The integer type
+      typedef Integer_ Integer;
+
+      typedef std::pair<long,Integer> Data_pair;
+
+      typedef std::map<Coefficient,Data_pair> Map;
+
+      typedef ::CGAL::Handle_with_policy<Map> Base;
+
+      typedef Best_approximation_cache<Coefficient,Integer> Self;
+
+      //! Default constructor
+      Best_approximation_cache() {
+      }
+      
+      //! Copy constructor
+      Best_approximation_cache(const Self& s)
+	: Base(static_cast<const Base&>(s))
+	{}
+      
+      /*! 
+       * \brief Checks whether the coefficient \c c already has already been
+       * approximated.
+       */
+      bool approximation_exists(Coefficient c) const {
+	typename Map::const_iterator it = this->ptr()->find(c);
+	return it!=this->ptr()->end();
+      }
+
+      /*!
+       * \brief The best approximation known for the coefficient <tt>c</tt>
+       * is returned.
+       *
+       * It is necessary to check whether an approximation of <tt>c</tt> exists
+       * using the \c approximation_exists method beforehand!
+       */
+      void get_best_approximation(Coefficient c,long& prec, Integer& val) const{
+	typename Map::const_iterator it = this->ptr()->find(c);
+	CGAL_assertion(it!=this->ptr()->end());
+	prec = it->second.first;
+	val = it->second.second;
+      }
+
+      /*!
+       * \brief Updates the approximation memory
+       *
+       * If an approximation for <tt>c</tt> is known, this method can be used
+       * to store it for later usage. The update is only performed if the
+       * precision is indeed higher than the current precision in the map.
+       */
+      void update_approximation(Coefficient c,long prec, const Integer& val) {
+	typename Map::iterator it = this->ptr()->find(c);
+	if(it!=this->ptr()->end()) {
+	  if(it->second.first >= prec) {
+	    return;
+	  }
+	  else {
+	    Data_pair better_pair = std::make_pair(prec,val);
+	    it->second=better_pair;
+	  }
+	}
+	else {
+	  Data_pair new_pair = std::make_pair(prec,val);
+	  this->ptr()->operator[](c)=new_pair;
+	}
+	
+      }
+    };
+    
+
+  }// namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif // AcX_BEST_APPROXIMATION_CACHE

Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Bitstream_descartes_traits_on_vert_line.h

@@ -0,0 +1,687 @@
+// TODO: Add licence
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL:$
+// $Id: $
+// 
+//
+// Author(s)     : Michael Kerber <mkerber@mpi-inf.mpg.de>
+//
+// ============================================================================
+
+#ifndef CGAL_BITSTREAM_DESCARTES_TRAITS_ON_VERT_LINE
+#define CGAL_BITSTREAM_DESCARTES_TRAITS_ON_VERT_LINE 1
+
+#include <CGAL/Algebraic_kernel_1.h>
+#include <CGAL/Arithmetic_kernel.h>
+#include <CGAL/Algebraic_kernel_d/Real_embeddable_extension.h>
+#include <CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/alg_real_utils.h>
+#include <CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Best_approximation_cache.h>
+
+CGAL_BEGIN_NAMESPACE 
+
+    // Workaround for min/max. This is just temporary!
+    // Necessary to prevent compiler confusion about std::min/max 
+    // and CGAL::min/max. The better solution is to put the CGAL functions into
+    // some seperate namespace and include it into CGAL-namepsace with the
+    // using command. The workaround works only for Sqrt_extensions
+template< class NT, class Root >
+inline Sqrt_extension<NT,Root> min BOOST_PREVENT_MACRO_SUBSTITUTION
+( const Sqrt_extension<NT,Root>& x , 
+  const Sqrt_extension<NT,Root>& y) {
+    return Min<Sqrt_extension<NT,Root> >() (x,y);
+}
+
+template< class NT, class Root >
+inline Sqrt_extension<NT,Root> max BOOST_PREVENT_MACRO_SUBSTITUTION
+( const Sqrt_extension<NT,Root>& x , 
+  const Sqrt_extension<NT,Root>& y) {
+    return Min<Sqrt_extension<NT,Root> >() (x,y);
+}
+
+namespace CGALi {
+
+
+// Representation of the Bitstream Descartes traits class
+template<typename Coefficient_,typename AlgebraicReal, typename Integer_>
+class Bitstream_descartes_traits_on_vert_line_rep {
+   
+public:
+
+    typedef AlgebraicReal Algebraic_real;
+    typedef Integer_ Integer;
+    typedef Coefficient_ Coefficient;
+
+    // Coefficient is a polynomial type
+    typedef typename Coefficient::NT Coeff_NT;
+
+    typedef typename Algebraic_real::Rational Boundary;
+
+    typedef typename Algebraic_real::Rational Rational;
+
+    typedef typename CGAL::Coercion_traits<Coeff_NT,Rational>::Type
+    Coercion;
+
+    typedef typename CGAL::Coercion_traits<Coeff_NT,Rational>::Cast
+    Rational_to_coercion_cast;
+
+    typedef typename 
+    CGAL::Get_arithmetic_kernel<Rational>::Arithmetic_kernel::
+    Bigfloat_interval BFI;
+
+    typedef typename 
+    CGAL::Get_arithmetic_kernel<Rational>::Arithmetic_kernel::
+    Bigfloat BF;
+
+    typedef CGAL::CGALi::Best_approximation_cache<Coefficient,Integer> 
+    Best_approximation_cache;
+
+    typedef Bitstream_descartes_traits_on_vert_line_rep
+    <Coefficient,Algebraic_real,Integer> Self;
+
+    Bitstream_descartes_traits_on_vert_line_rep
+    (Algebraic_real alpha,bool use_approx_memory) 
+	: alpha_(alpha), use_approx_memory(use_approx_memory)
+    {	
+    } 
+
+
+    /*
+      Bitstream_descartes_traits_on_vert_line_rep
+      (const Self& s) 
+      : alpha_(s.alpha)
+      {	
+      } 
+    */
+
+    bool uses_approx() const {
+	return use_approx_memory;
+    }
+
+    Best_approximation_cache approx_mem() const {
+	return approximations_;
+    }
+     
+
+    Algebraic_real alpha() const {
+	return alpha_;
+    }
+
+    class Approximator {
+
+    private:
+	
+	Algebraic_real alpha_;
+
+	Best_approximation_cache approximations_;
+
+	bool use_approx_mem;
+
+        long curr_prec;
+
+    public:
+	Approximator(Algebraic_real alpha) : alpha_(alpha),use_approx_mem(true), curr_prec(4) {};
+
+	Approximator(Algebraic_real alpha,
+		     Best_approximation_cache approximations) 
+            : alpha_(alpha), approximations_(approximations), use_approx_mem(true), curr_prec(4) {};
+
+	Approximator() {};
+
+        template<typename NT>
+        class Coeff_to_bfi_functor {
+
+        public:
+
+            typedef NT argument_type;
+            typedef BFI result_type;
+
+            BFI operator() (NT x) const {
+                
+                return CGAL::convert_to_bfi(x);
+
+            }
+
+        };
+        
+
+	typedef boost::numeric::interval<Coercion> Interval;
+
+        long log_width(Algebraic_real alpha) const {
+            CGAL_assertion(! alpha.is_rational());
+            double d = CGAL::to_double(alpha.high()- alpha.low());
+            if(d != 0.0) {
+                return CGAL::CGALi::ceil_log2_abs(d);
+            }
+            return 0;
+        }
+    
+#if AcX_USE_NO_BFI_APPROX_IN_BITSTREAM_TRAITS
+#warning uses no bfi-approx!
+
+Integer operator() (Coefficient f, long p) {
+
+            //AcX_DSTREAM("Called approximator with " << f << " and " << p << std::endl);
+            //AcX_DSTREAM("Alpha: " << (CGAL::to_double(alpha_.high()-alpha_.low())) << " " << rational_to_double(alpha_.high()) << std::endl);
+
+
+      
+            if(use_approx_mem) {
+                if(approximations_.approximation_exists(f)) {
+                    long prec;
+                    Integer approx;
+                    approximations_.get_best_approximation(f,prec,approx);
+                    if(prec>=p) {
+                        //AcX_DSTREAM("use older approx.." << std::flush);
+                        Integer divisor = CGAL::ipower((Integer)2,prec-p);
+                        Integer return_value = CGAL::div(approx,divisor);
+                        //AcX_DSTREAM("done" << std::endl);
+                        return return_value;
+                    }
+                }
+            }
+            
+            //AcX_DSTREAM("Have to approximate..." << std::endl);
+	 
+            Coercion bound;
+            
+            Rational_to_coercion_cast rational_to_coercion_cast;
+
+            if(p<0) {
+                bound = rational_to_coercion_cast
+                    (Rational(CGAL::ipower((Integer)2,-p)));
+            }
+            else {
+                typename CGAL::Fraction_traits<Rational>::Compose compose;
+                bound =
+                    rational_to_coercion_cast
+                    (compose(1,CGAL::ipower(Integer(2),p)));
+            }
+            Interval f_alpha_iv;
+
+
+
+            if(alpha_.is_rational()) {
+                Coercion f_alpha=f.evaluate(alpha_.low());
+                f_alpha_iv=Interval(f_alpha,f_alpha);
+                
+            }
+            else {
+                typedef Approximate_arithmetic_controller
+                    <Coefficient,Algebraic_real>
+                    Approximate_controller;
+
+                Approximate_controller approx_controller(f,alpha_);
+
+                while(true) {
+
+                    f_alpha_iv = approx_controller.interval_approximation();
+	      
+                    //AcX_DSTREAM("done" << std::endl);
+                    if(CGAL::compare
+                       (f_alpha_iv.upper()-f_alpha_iv.lower(),bound) ==
+                       CGAL::SMALLER) {
+                        break;
+                    }
+                    //AcX_DSTREAM("refine..." << std::flush);
+	      
+                    approx_controller.refine_value();
+	      
+                    //AcX_DSTREAM("done" << std::endl);
+                }
+            }
+       
+
+            Integer return_value;
+            typename CGAL::Fraction_traits<Coercion>::Decompose decompose;
+            typedef typename 
+                CGAL::Fraction_traits<Coercion>::Numerator_type Numerator;
+            typedef typename 
+                CGAL::Fraction_traits<Coercion>::Denominator_type Denominator;
+            if(sign(f_alpha_iv.lower())==CGAL::ZERO) {
+                return_value = Integer(0);
+            }
+            else if(CGAL::sign(f_alpha_iv.lower()) != 
+                    CGAL::sign(f_alpha_iv.upper())) {
+                return_value = Integer(0);
+            }
+            else {
+                Numerator num;
+                Denominator denom;
+
+                decompose(f_alpha_iv.upper(),num,denom);                
+
+
+                if(p<0) {
+                    denom *= CGAL::ipower(Integer(2),-p);
+                }
+                else {
+                    num *= CGAL::ipower(Integer(2),p);
+                }
+                return_value = CGAL::div(num,denom);
+
+            }
+
+            if(use_approx_mem) {
+                approximations_.update_approximation(f,p,return_value);
+            }
+            //AcX_DSTREAM(" returning " << " " <<  return_value << std::endl);
+
+            return return_value;
+	}
+    
+
+#else
+
+        CGAL::Polynomial<BFI> _convert_polynomial_to_bfi(Coefficient f) const {
+            std::vector<BFI> coeffs;
+            for(int i = 0; i <= f.degree(); i++) {
+                coeffs.push_back(CGAL::convert_to_bfi(f[i]));
+            }
+            return CGAL::Polynomial<BFI>(coeffs.begin(), coeffs.end());   
+        }
+        
+        
+        Integer operator() (Coefficient f, long p) {
+            
+            typename CGAL::CGALi::Float_traits<BF>::Get_exponent get_exp;
+            typename CGAL::CGALi::Float_traits<BF>::Get_mantissa get_m;
+
+            long old_prec = CGAL::get_precision(BFI());
+            
+            CGAL::Polynomial<BFI> f_bfi;
+            BFI alpha_bfi, f_alpha_bfi;
+            
+            long prec = 16;
+            
+            Integer return_value;
+            
+            long wbit = 0;
+
+            while(true) {
+                CGAL::set_precision(BFI(),prec);
+                
+                f_bfi = _convert_polynomial_to_bfi(f);
+                alpha_bfi = CGAL::convert_to_bfi(alpha_);
+                
+                f_alpha_bfi = f_bfi.evaluate(alpha_bfi);
+                
+                if(!CGAL::singleton(f_alpha_bfi)) {
+                    long ceil = CGAL::CGALi::ceil_log2_abs(f_alpha_bfi);
+                    long signi = CGAL::get_significant_bits(f_alpha_bfi);
+                    wbit   = ceil - signi + p;
+                     
+                } 
+                
+                if(wbit<-5 || CGAL::singleton(f_alpha_bfi)) {
+                    break;
+                } else {
+                    prec*=2;
+                }
+            }
+            BF lower = CGAL::lower(f_alpha_bfi);
+            
+            long shift = - (p + get_exp(lower)); 
+            Integer bfi_m(get_m(lower)); 
+             if( shift > 0 ){
+                while(shift>63) {
+                    bfi_m = (bfi_m >> 63);
+                    shift-=63;
+                }
+                bfi_m = (bfi_m >> shift);
+            }else{
+                // add 0 bits 
+                bfi_m = (bfi_m << -shift);   
+            }   
+            return_value = bfi_m;
+            CGAL::set_precision(BFI(),old_prec);
+            
+            return return_value;
+        }
+#endif
+
+    };
+
+
+    class Lower_bound_log2_abs {
+
+    private:
+	Algebraic_real alpha_;
+
+    public:
+	Lower_bound_log2_abs(Algebraic_real alpha) : alpha_(alpha) {}
+	Lower_bound_log2_abs() {};
+
+	typedef boost::numeric::interval<Coercion> Interval;
+
+	long operator() (Coefficient f) {
+            //std::cout << "Called lower_bound_log2_abs with " 
+            //          << f << std::flush;
+          
+            CGAL_assertion(! alpha_.is_root_of(f));
+
+            Interval f_alpha_iv;
+            Coercion abs_lower,abs_upper;
+
+            typedef CGAL::CGALi::Approximate_arithmetic_controller
+                <Coefficient,Algebraic_real> Approximate_controller;
+          
+            Approximate_controller approx_controller(f,alpha_);
+
+            while(true) {
+                f_alpha_iv = approx_controller.interval_approximation();
+                CGAL::Sign lower_sign = CGAL::sign(f_alpha_iv.lower());
+                if(CGAL::sign(f_alpha_iv.upper())==lower_sign) {
+                    if(lower_sign==CGAL::POSITIVE) {
+                        abs_lower=f_alpha_iv.lower();
+                        abs_upper=f_alpha_iv.upper();
+                    }
+                    else {
+                        abs_lower=CGAL::abs(f_alpha_iv.upper());
+                        abs_upper=CGAL::abs(f_alpha_iv.lower());
+                    }
+    
+                    BFI bfi_low = CGAL::convert_to_bfi(abs_lower),
+                        bfi_high = CGAL::convert_to_bfi(abs_upper);
+                    long lower_bound = CGAL::CGALi::floor_log2_abs(bfi_low),
+                        upper_bound = CGAL::CGALi::ceil_log2_abs(bfi_high);
+                    
+                    CGAL_assertion(upper_bound>=lower_bound);
+                    if(upper_bound-lower_bound <=2) {
+                        //std::cout << "returning " << lower_bound << std::endl;
+                        return lower_bound;
+                    }
+                }
+                approx_controller.refine_value();
+            }
+	}
+
+    };
+
+    class Upper_bound_log2_abs_approximator {
+
+    private:
+	Algebraic_real alpha_;
+
+        std::vector<long> coefficients_for_alpha;
+
+	int refinements_of_alpha;
+
+	bool zero_test_enabled;
+
+	int refinement_limit;
+
+	// Stores id of polynomials which are known to vanish (or not to 
+	// vanish) at alpha
+	std::vector<long> zeroes,non_zeroes;
+
+    public:
+	Upper_bound_log2_abs_approximator(Algebraic_real alpha) 
+            : alpha_(alpha),refinements_of_alpha(0),zero_test_enabled(false)
+        {}
+
+        Upper_bound_log2_abs_approximator() {};
+
+        typedef boost::numeric::interval<Boundary> Boundary_interval;
+	typedef boost::numeric::interval<Coercion> Interval;
+
+	bool initial_upper_bound
+        (Coefficient f, long& ub_log2_abs,bool& is_certainly_zero) {
+            return improve_upper_bound(f,ub_log2_abs,is_certainly_zero);
+	}
+
+	bool improve_upper_bound
+        (const Coefficient f, long& ub_log2_abs,bool& is_certainly_zero) {
+            //AcX_DSTREAM("improve upper bound.." << f << std::flush);
+            if(std::find(coefficients_for_alpha.begin(),
+                         coefficients_for_alpha.end(),
+                         f.id())!=coefficients_for_alpha.end()) {
+                coefficients_for_alpha.clear();
+                //AcX_DSTREAM("refine.." << std::flush);
+                alpha_.bisect();
+                ++refinements_of_alpha;
+                if(refinements_of_alpha>=refinement_limit) {
+                    zero_test_enabled=true;
+                }
+                //AcX_DSTREAM("done.." << std::flush);
+            }
+            coefficients_for_alpha.push_back(f.id());
+            if(zero_test_enabled) {
+                if(std::find(zeroes.begin(),
+                             zeroes.end(),
+                             f.id())!=zeroes.end()) {
+                    is_certainly_zero=true;
+                    return true;
+                }
+                else if(std::find(non_zeroes.begin(),
+                                  non_zeroes.end(),
+                                  f.id())!=non_zeroes.end()) {
+                    is_certainly_zero=false;
+                }
+                else {
+                    bool zero = CGAL::CGALi::is_root_of(alpha_,f);
+                    if(zero) {
+                        zeroes.push_back(f.id());
+                        is_certainly_zero=true;
+                        return true;
+                    }
+                    else {
+                        non_zeroes.push_back(f.id());
+                        is_certainly_zero=false;
+                    }
+                }
+            }
+            else {
+                is_certainly_zero=false;
+            }
+            Boundary_interval alpha_iv(alpha_.low(),alpha_.high());
+            Interval f_alpha_iv = evaluate_iv(f,alpha_iv);
+            CGAL::Sign lower_sign = CGAL::sign(f_alpha_iv.lower());
+            bool iv_contains_zero = lower_sign 
+                != CGAL::sign(f_alpha_iv.upper());
+            Coercion abs_upper
+                =(CGAL::abs(f_alpha_iv.lower())<CGAL::abs(f_alpha_iv.upper())) 
+                ? CGAL::abs(f_alpha_iv.upper()) 
+                : CGAL::abs(f_alpha_iv.lower());
+            Coercion abs_lower(0);
+            if(! iv_contains_zero) {
+                abs_lower = (CGAL::abs(f_alpha_iv.lower())
+                             < CGAL::abs(f_alpha_iv.upper())) 
+                    ? CGAL::abs(f_alpha_iv.lower()) 
+                    : CGAL::abs(f_alpha_iv.upper());
+            }
+            //Numerator upper_num;
+            //Denominator upper_denom;
+            //decompose(abs_upper,upper_num,upper_denom);
+            if(CGAL::sign(abs_upper)==CGAL::ZERO) {
+                is_certainly_zero=true;
+                return true;
+            }
+            ub_log2_abs = CGAL::CGALi::ceil_log2_abs
+                (CGAL::convert_to_bfi(abs_upper));
+            //ub_log2_abs = num_ceil_log2_abs(upper_num) - 
+            //    denom_floor_log2_abs(upper_denom);
+
+            if(! iv_contains_zero) {
+                //Numerator lower_num;
+                //Denominator lower_denom;
+                //decompose(abs_lower,lower_num,lower_denom);
+                //long lb_log2_abs = num_floor_log2_abs(lower_num) 
+                //    - denom_ceil_log2_abs(lower_denom);
+                long lb_log2_abs 
+                    = CGAL::CGALi::ceil_log2_abs
+                    (CGAL::convert_to_bfi(abs_lower));
+                CGAL_assertion(ub_log2_abs >= lb_log2_abs);
+                //AcX_DSTREAM("Upper: " << ub_log2_abs << " Lower: " << lb_log2_abs << std::endl); 
+                //AcX_DSTREAM(((ub_log2_abs - lb_log2_abs) <= 2) << std::endl);
+	    
+                return ((ub_log2_abs - lb_log2_abs) <= 2);
+            }
+            else {
+                //AcX_DSTREAM("Upper: " << ub_log2_abs << std::endl);
+                return false;
+            }
+	}
+    };
+
+private:
+      
+    Algebraic_real alpha_;
+
+    Best_approximation_cache approximations_;
+
+    bool use_approx_memory;
+
+};
+
+/*! 
+ * \brief Traits class for the Bitstream Descartes method over algebraic 
+ * extensions.
+ *
+ * Approximates coefficients of the polynomial 
+ * <tt>g=f|<sub>x=a</sub></tt> 
+ * using interval arithmetic. The \c AlgebraicReal object \c a is fixed
+ * for the traits class. Each coefficient is represented by an univariate
+ * integer polynomial. For a more detailled specification of its member 
+ * functions, see the documentation of 
+ * \c CGAL::CGALi::Bitstream_descartes_rndl_tree
+ *
+ * Using this traits class, bivariate polynomials can be isolated with 
+ * AcX::Bitstream_descartes_bfs, the inner variable \c x is interpreted as
+ * the algebraic number \c a.
+ */
+template<typename Coefficient_,typename AlgebraicReal,typename Integer_>
+class Bitstream_descartes_traits_on_vert_line
+    : public ::CGAL::Handle_with_policy
+    <CGAL::CGALi::Bitstream_descartes_traits_on_vert_line_rep
+      <Coefficient_,AlgebraicReal,Integer_> > {
+      
+public:
+      
+    //! Coefficient type
+    typedef Coefficient_ Coefficient;
+
+    //! Algebraic real type
+    typedef AlgebraicReal Algebraic_real;
+
+    //! Integer type
+    typedef Integer_ Integer;
+      
+    typedef Bitstream_descartes_traits_on_vert_line
+    <Coefficient,Algebraic_real,Integer> Self;
+      
+    typedef CGAL::CGALi::Bitstream_descartes_traits_on_vert_line_rep
+    <Coefficient,AlgebraicReal,Integer> Rep;
+      
+    typedef  ::CGAL::Handle_with_policy<Rep> Handle;
+
+    typedef typename Rep::Boundary Boundary;
+      
+    //! Functor type for the approximation of coefficients
+    typedef typename Rep::Approximator Approximator;
+      
+    //! Functor type for lower bounds on coefficients
+    typedef typename Rep::Lower_bound_log2_abs Lower_bound_log2_abs;
+      
+    //! Functor type for upper bounds on coefficients
+    typedef typename Rep::Upper_bound_log2_abs_approximator
+    Upper_bound_log2_abs_approximator;
+      
+    typedef typename Rep::Rational Rational;
+      
+    // We can just take these from NT_traits
+    typedef typename CGAL::Real_embeddable_traits<Integer>::Sign
+    Sign;
+    typedef typename CGAL::CGALi::Real_embeddable_extension<Integer>
+        ::Ceil_log2_abs Ceil_log2_abs_Integer;
+    typedef typename CGAL::CGALi::Real_embeddable_extension<long>
+        ::Ceil_log2_abs Ceil_log2_abs_long;
+
+      
+    /*!
+     * \brief Constructor for the polynomial <tt>f|<sub>x=a</sub></tt>
+     *
+     * Initialises an instance such that all bivariate polynomials are
+     * interpreted as univariate polynomials in the algebraic extension
+     * <tt>Z[a]</tt>.
+     * If the \c use_approx_mem flag is set, the best known approximation
+     * of each coefficient is stored in a 
+     * \c CGAL::CGALi::Intern::Best_approximation_cache instance.
+     */
+    Bitstream_descartes_traits_on_vert_line
+    (AlgebraicReal alpha=AlgebraicReal(),bool use_approx_mem = true) 
+	: Handle(alpha,use_approx_mem)
+    {	
+    } 
+      
+    /*
+      Bitstream_descartes_traits_on_vert_line
+      (const Self& s) 
+      : Handle(s)
+      {	
+      } 
+    */
+      
+    Algebraic_real alpha() const {
+        return this->ptr()->alpha();
+    }
+
+    //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree.      
+    Approximator approximator_object() const {
+	if(this->ptr()->uses_approx()) {
+            return Approximator(this->ptr()->alpha(),this->ptr()->approx_mem());
+	}
+	else {
+            return Approximator(this->ptr()->alpha());
+	}
+    }
+      
+    //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree.            
+    Upper_bound_log2_abs_approximator 
+    upper_bound_log2_abs_approximator_object() const {
+	return Upper_bound_log2_abs_approximator(this->ptr()->alpha());
+    }
+      
+    //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree.            
+    Lower_bound_log2_abs lower_bound_log2_abs_object() const {
+	return Lower_bound_log2_abs(this->ptr()->alpha());
+    }
+      
+      
+    //! See the documentation of CGAL::CGALi::Bitstream_descartes_rndl_tree. 
+    class Boundary_creator {
+	
+    public:
+	
+	Boundary_creator() {}
+	
+	Boundary operator() (Integer x,long p) {
+            Integer num=x, denom,two(2),q,r;
+            if(p <0) {
+                CGAL::div_mod(num,two,q,r);
+                while(r==Integer(0) && p<0) {
+                    num=q;
+                    p++;
+                    CGAL::div_mod(num,two,q,r);
+                }
+                denom = CGAL::ipower(Integer(2),-p);
+            }
+            else {
+                num*=CGAL::ipower(Integer(2),p);
+                denom=1;
+            }
+            Rational b(num,denom);
+            CGAL::simplify(b);
+            return b;
+	} 
+    };
+      
+};
+    
+} // namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif

Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Non_generic_position_exception.h

@@ -0,0 +1,44 @@
+// TODO: Add licence
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL:$
+// $Id: $
+// 
+//
+// Author(s)     : Michael Kerber <mkerber@mpi-inf.mpg.de>
+//
+// ============================================================================
+
+
+
+#ifndef CGAL_NON_GENERIC_POSITION_EXCEPTION
+#define CGAL_NON_GENERIC_POSITION_EXCEPTION
+
+
+CGAL_BEGIN_NAMESPACE
+
+  namespace CGALi {
+
+    /*! 
+     * \brief Exception class for not sufficiently generic positions.
+     *
+     * Must be thrown whenever a curve cannot be analysed because its position
+     * is not "good enough".
+     */
+    class Non_generic_position_exception {
+      
+    public:
+
+      //! Default constructible
+      Non_generic_position_exception() {}
+      
+    };
+
+  } // namespace CGALi
+
+CGAL_END_NAMESPACE
+
+
+#endif

Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/alg_real_utils.h

@@ -0,0 +1,691 @@
+// TODO: Add licence
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL:$
+// $Id: $
+// 
+//
+// Author(s)     : Michael Kerber <mkerber@mpi-inf.mpg.de>
+//
+// ============================================================================
+
+
+#ifndef CGAL_ALG_REAL_UTILS
+#define CGAL_ALG_REAL_UTILS 1
+
+#include <iterator>
+
+#include <CGAL/basic.h>
+#include <CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/enums.h>
+#include <CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/Approximate_arithmetic_controller.h>
+#include <CGAL/Arithmetic_kernel.h>
+#include <CGAL/interval_support.h>
+#include <CGAL/Algebraic_kernel_d/Float_traits.h>
+
+#include <CGAL/Algebraic_kernel_d/Real_embeddable_extension.h>
+
+#include <CGAL/Coercion_traits.h>
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+/* 
+ * \brief Function for merging two sets
+ *
+ * This function is similar to the \c std::union_set operation.
+ * Additionally, it provides a sequence of AcX::Three_valued
+ * providing information to which input set the corresponding root
+ * in the merged sequence belonged
+ *
+ * The BinaryFunction must have the Result type CGAL::Comparison_result.
+ */
+template<typename BinaryFunction,
+      typename InputIterator1,typename InputIterator2,
+      typename OutputIterator1,typename OutputIterator2>
+std::pair<OutputIterator1,OutputIterator2>
+set_union_with_source(InputIterator1 first_begin,
+                      InputIterator1 first_end,
+                      InputIterator2 second_begin,
+                      InputIterator2 second_end,
+                      OutputIterator1 merged_values,
+                      OutputIterator2 merged_values_info,
+                      BinaryFunction compare) {
+      
+    InputIterator1 first_it=first_begin;
+    InputIterator2 second_it=second_begin;
+
+    while((first_it != first_end) || (second_it!=second_end)) {
+	if(first_it == first_end) {
+            *merged_values=*second_it++;
+            ++merged_values;
+            *merged_values_info++ = CGAL::CGALi::ROOT_OF_SECOND_SET;
+            continue;
+	}
+	if(second_it == second_end) {
+            *merged_values++=*first_it++;
+            *merged_values_info++ = CGAL::CGALi::ROOT_OF_FIRST_SET;
+            continue;
+	}
+
+	CGAL::Comparison_result c = compare(*first_it,*second_it);
+
+	if(c==CGAL::EQUAL) {
+            *merged_values++=*first_it++;
+            ++second_it;
+            *merged_values_info++ = CGAL::CGALi::ROOT_OF_BOTH_SETS;
+            continue;
+	}
+	if(c==CGAL::SMALLER) {
+            *merged_values++=*first_it++;
+            *merged_values_info++ = CGAL::CGALi::ROOT_OF_FIRST_SET;
+            continue;
+	}
+	if(c==CGAL::LARGER) {
+            *merged_values++=*second_it++;
+            *merged_values_info++ = CGAL::CGALi::ROOT_OF_SECOND_SET;
+            continue;
+	}
+    }
+    return std::make_pair(merged_values,merged_values_info);
+}
+
+
+
+/*! \brief Tries to find a SIMPLE rational q with a<q<b.
+ *
+ * In this context, simple means that the denominator of <tt>q</tt>
+ * is a power of two, and is not too big. There is no guarantee to find
+ * the rational value between <tt>a</tt> and <tt>b</tt> of minimal
+ * bit size.
+ */
+template<typename Algebraic_real>
+typename Algebraic_real::Rational
+simple_rational_between(const Algebraic_real& a,
+                        const Algebraic_real&b) {
+
+    //srb.start();
+    typedef typename Algebraic_real::Rational Rational;
+    typename CGAL::Fraction_traits<Rational>::Compose compose;
+    typedef typename 
+        CGAL::Get_arithmetic_kernel<Rational>::Arithmetic_kernel AK;
+    typedef typename AK::Bigfloat Bigfloat;
+    typedef typename AK::Bigfloat_interval Bigfloat_interval;
+    typedef typename AK::Integer Integer;
+
+    long old_prec = CGAL::get_precision(Bigfloat_interval());
+
+    CGAL_assertion(a<b);
+    //std::cout << "Intermediate1: " << CGAL::to_double(a) << " " << CGAL::to_double(b) << std::endl;
+    /*
+     * First, refine a and b until their isolating intervals are disjoint
+     * Therefore, the bigger interval is refined in each substep
+     */
+    //srb_a.start();
+    if(a.high() >= b.low()) {
+        Rational size_a=a.high()-a.low(),
+            size_b=b.high() - b.low();
+        while(a.high() >= b.low()) {
+            if(size_a < size_b) {
+                b.refine();
+                size_b=b.high() - b.low();
+            } else {
+                a.refine();
+                size_a=a.high()-a.low();
+            }
+        }
+    }
+    //srb_a.stop();
+
+    //srb_b.start();
+    Bigfloat x=CGAL::upper(CGAL::convert_to_bfi(a.high())),
+        y=CGAL::lower(CGAL::convert_to_bfi(b.low()));
+    if(x>=y) {
+        Rational size_a=a.high() - a.low(),
+            size_b=b.high() - b.low(),
+            size_max = size_a>size_b ? size_a : size_b,
+            size_int = b.low()-a.high();
+        while(x>=y) {
+            //std::cout << "x and y: " <<  x << " and " << y << std::endl;
+            //std::cout << "sizes: " << CGAL::to_double(size_int) << " " << CGAL::to_double(size_max) << std::endl;
+            if(size_int>size_max) {
+                CGAL::set_precision(Bigfloat_interval(),
+                                   2*CGAL::get_precision(Bigfloat_interval()));
+                x=CGAL::upper(CGAL::convert_to_bfi(a.high()));
+                y=CGAL::lower(CGAL::convert_to_bfi(b.low()));
+            } else {
+                if(size_a < size_b) {
+                    b.refine();
+                    size_b=b.high() - b.low();
+                    y=CGAL::lower(CGAL::convert_to_bfi(b.low()));
+                } else {
+                    a.refine();
+                    size_a=a.high()-a.low();    
+                    x=CGAL::upper(CGAL::convert_to_bfi(a.high()));
+                }
+                size_max = size_a>size_b ? size_a : size_b;
+                size_int = b.low()-a.high();
+            }
+        }
+    }
+    //srb_b.stop();
+    //std::cout << "Intermediate2: " << x << " " << y << std::endl;
+    typename CGAL::CGALi::Float_traits<Bigfloat>::Get_mantissa mantissa;
+    typename CGAL::CGALi::Float_traits<Bigfloat>::Get_exponent exponent;
+
+    Integer x_m = mantissa(x),
+        y_m=mantissa(y);
+    long x_e = exponent(x),  y_e = exponent(y);
+    //std::cout << "Floats1: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
+    if((x_m > 0 &&  y_m < 0) || x_m < 0 && y_m > 0) {
+        //srb.stop();
+        return Rational(0);
+    }
+    bool negative=false;
+    if(x_m<=0 && y_m <=0) {
+        x_m=-x_m;
+        y_m=-y_m;
+        std::swap(x_m,y_m);
+        std::swap(x_e,y_e);
+        negative=true;
+    }
+    // Now, we have that (x_m,x_e) represents a number smaller than (y_m,y_e)
+    //srb_c.start();
+    //std::cout << "Floats2: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
+
+    // As long as the mantissa is even, simplify
+    while(x_m != 0 && (x_m & 1)==0 ) {
+        x_m=x_m >> 1;
+        x_e++;
+    }
+    while(y_m != 0 && (y_m & 1)==0 ) {
+        y_m=y_m >> 1;
+        y_e++;
+    }
+    //srb_c.stop();
+    //std::cout << "Floats3: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
+
+    // Bring both numbers to a common exponent
+    //srb_d.start();
+    long min_e = x_e < y_e ? x_e : y_e;
+    while(x_e > min_e) {
+        x_m=x_m << 1;
+        x_e--;
+    }
+    while(y_e > min_e) {
+        y_m=y_m << 1;
+        y_e--;
+    }
+    //srb_d.stop();
+    CGAL_assertion(y_e==x_e && x_e==min_e);
+
+    //std::cout << "Floats4: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
+
+    // Avoid mantissas to have difference one
+    if(y_m-x_m==Integer(1)) {
+        x_m=x_m << 1;
+        y_m=y_m << 1;
+        x_e--;
+        y_e--;
+        min_e--;
+    }
+    //std::cout << "Floats5: " << x_m << " " << x_e << " and " << y_m << " " << y_e << std::endl;
+    Integer final_mantissa(0);
+    //srb_e.start();
+    long x_log = x_m==Integer(0) ? -1 : CGAL::CGALi::floor_log2_abs(x_m),
+        y_log = y_m==Integer(0) ? -1 : CGAL::CGALi::floor_log2_abs(y_m),
+        old_log = y_log;
+    while(x_log==y_log) {
+        //std::cout << "here" << std::endl;
+        while(old_log > y_log) {
+            final_mantissa = final_mantissa << 1;
+            old_log--;
+        }
+        CGAL_assertion((x_m & ((Integer(1) << x_log) - 1)) == x_m - CGAL::ipower(Integer(2),x_log));
+        x_m = x_m & ((Integer(1) << x_log) - 1); // x_m - CGAL::ipower(Integer(2),x_log);
+        y_m = y_m & ((Integer(1) << y_log) - 1); // y_m - CGAL::ipower(Integer(2),y_log);
+
+        final_mantissa++;
+        old_log=y_log;
+        x_log = x_m==0 ? -1 : CGAL::CGALi::floor_log2_abs(x_m);
+        y_log = y_m==0 ? -1 : CGAL::CGALi::floor_log2_abs(y_m);
+    }
+    //srb_e.stop();
+    // Now, x_log != y_log, in fact, y_log is greater
+    CGAL_assertion(x_log<y_log);
+    //srb_f.start();
+    while(old_log > y_log) {
+        final_mantissa = final_mantissa << 1;
+        old_log--;
+    }
+    if((y_m & ((Integer(1) << y_log) - 1 ))==0) { // y_m - CGAL::ipower(Integer(2),y_log)==0) {
+        // Now, the constructed value would be equal to
+        while(y_log!=0 && x_log==y_log-1) {
+            final_mantissa = final_mantissa << 1;
+            final_mantissa++;
+            y_log--;
+            x_m = x_m==0 ? 0 : x_m & ((Integer(1) << x_log) - 1); //x_m - CGAL::ipower(Integer(2),x_log);
+            x_log = x_m==0 ? -1 : CGAL::CGALi::floor_log2_abs(x_m);
+        }
+        final_mantissa = final_mantissa << 1;
+        final_mantissa++;
+        y_log--;
+    } else {
+        final_mantissa++;
+    }
+    //srb_f.stop();
+    min_e += y_log;
+    Rational rat_between;
+    //std::cout << "Min_e: " << min_e << std::endl;
+    if(min_e > 0) {
+        rat_between = compose(final_mantissa << min_e,
+                              Integer(1));
+    } else {
+        rat_between = compose(final_mantissa, Integer(1) << -min_e);
+    }
+    if(negative) {
+        rat_between = -rat_between;
+    }
+    //std::cout << "Result: " << a.high() << " " << rat_between << " " << b.low() << std::endl;
+    CGAL_assertion(a.high() < rat_between);
+    CGAL_assertion(b.low() > rat_between);
+    CGAL::set_precision(Bigfloat_interval(),old_prec);
+    //srb.stop();
+    return rat_between;
+}
+
+
+/*!
+ * \brief Produces intermediate rational values for a list of 
+ * algebraic reals.
+ *
+ * For a list of Algebraic real values with \c n elements, a list with 
+ * <tt>n+1</tt> elements of rational values is given such that the 
+ * <tt>i</tt>th element is
+ * between the <tt>i</tt>th and the <tt>(i+1)</tt>th element of the input list
+ *
+ * The input list must be in increasing order
+ */
+template<typename InputIterator, typename OutputIterator> 
+OutputIterator find_intermediate_values(InputIterator start,
+                                        InputIterator end,
+                                        OutputIterator output) {
+    typedef typename std::iterator_traits<InputIterator>::value_type 
+	Alg_real;
+    typedef typename Alg_real::Rational Rational;
+    if(start==end) {
+	// Empty vector, create one element
+	*output=Rational(0);
+	++output;
+	return output;
+    }
+    switch( CGAL::sign( start->low() ) ) {
+    case(CGAL::ZERO): {
+        *output = Rational(-1);
+        break;
+    }
+    case(CGAL::POSITIVE): {
+        *output = Rational(0);
+        break;
+    }
+    case(CGAL::NEGATIVE): {
+        Alg_real small_value = Alg_real(Rational(2)*start->low());
+        *output 
+            = simple_rational_between(small_value,*start);
+    }
+    }
+    output++;
+      
+    ++output;
+    InputIterator it_1(start),it_2(start);
+    ++it_2;
+    while(it_2 != end) {
+	CGAL_assertion(it_1->compare(*it_2)==CGAL::SMALLER);
+	Rational beta 
+            = simple_rational_between(*it_1,*it_2);
+	*output=beta;
+	++output;
+	++it_1;
+	++it_2;
+    }
+    switch( CGAL::sign( it_1->high() ) ) {
+    case(CGAL::ZERO): {
+        *output = Rational(1);
+        break;
+    }
+    case(CGAL::NEGATIVE): {
+        *output = Rational(0);
+        break;
+    }
+    case(CGAL::POSITIVE): {
+        Alg_real big_value = Alg_real(Rational(2)*it_1->high());
+        *output 
+            = simple_rational_between(*it_1,big_value);
+    }
+    }
+    output++;
+    return output;
+}
+
+/*!
+ * \brief finds a Rational value left of an Algebraic real alpha
+ */
+template<typename AlgebraicReal> typename AlgebraicReal::Rational 
+simple_rational_left_of(AlgebraicReal ar) {
+
+    typedef AlgebraicReal Algebraic_real;
+    typedef typename Algebraic_real::Rational Rational;
+
+    switch( CGAL::sign( ar ) ) {
+    case(CGAL::ZERO): {
+        return Rational(-1);
+        break;
+    }
+    case(CGAL::POSITIVE): {
+        return Rational(0);
+        break;
+    }
+    case(CGAL::NEGATIVE): {
+        Algebraic_real small_value 
+            = Algebraic_real(Rational(2)*ar.low());
+    
+        return simple_rational_between(small_value,ar);
+        // = small_value.rational_between(ar);
+        //= ar.low()-1;
+    }
+    }
+    // never reached
+    return Rational(0);
+}
+
+/*!
+ * \brief finds a Rational value rightt of an Algebraic real alpha
+ */
+template<typename AlgebraicReal> typename AlgebraicReal::Rational 
+simple_rational_right_of(AlgebraicReal ar) {
+    
+    typedef AlgebraicReal Algebraic_real;
+    typedef typename Algebraic_real::Rational Rational;
+
+    return -simple_rational_left_of(-ar);
+
+}
+  
+/*!
+ * \brief Tries to get the (non-zero) sign of <tt>f(alpha)</tt> with interval
+ * arithmetic.
+ *
+ * Performs refinements of alpha and evaluates \c f at \c alpha. If
+ * zero is still in the interval when the isolating interval of \c alpha
+ * is smaller then \c eps, CGAL::ZERO
+ * is returned, denoting that the sign could not be computed. Otherwise,
+ * the correct sign of \c f(alpha) is returned
+ *
+ * If \c n is set to zero (default), \c alpha is refined 
+ * until the sign can be determined correctly. This may lead to an 
+ * infinite loop, if <tt>f(alpha)=0</tt>
+ *
+ */
+template<typename AlgebraicReal,typename Poly_1_,typename Boundary_>
+CGAL::Sign estimate_sign_at(AlgebraicReal alpha,
+                            const Poly_1_& f,
+                            const Boundary_ eps
+                            =Boundary_(0)) {
+    typedef Boundary_ Boundary;
+    typedef Poly_1_ Poly_1;
+    typedef AlgebraicReal Algebraic_real;
+    typedef typename CGAL::Coercion_traits<
+        typename Poly_1::NT,
+        Boundary>::Type Interval_boundary;
+    typename CGAL::Coercion_traits<typename Poly_1::NT,
+                                   Boundary>::Cast boundary_cast;
+    typedef boost::numeric::interval<Interval_boundary> Interval;
+    Interval ev;
+    CGAL::Sign sign=CGAL::ZERO;
+
+    typedef CGAL::CGALi::Approximate_arithmetic_controller
+        <Poly_1,Algebraic_real> Approximate_controller;
+
+    Approximate_controller approx_controller(f,alpha);
+    
+    while(sign==CGAL::ZERO) {
+        ev = approx_controller.interval_approximation();
+
+        if(! boost::numeric::in_zero(ev)) {
+            sign=CGAL::compare(ev.lower(),0);
+            CGAL_assertion(sign!=CGAL::ZERO);
+        }
+        else {
+            if(CGAL::compare(alpha.high()-alpha.low(),
+                             boundary_cast(eps))==CGAL::SMALLER) {
+                break;
+            }
+            approx_controller.refine_value();
+
+        }
+    }
+    return sign;
+}
+     
+
+/*!
+ * \brief Evaluates poynomial at intervals with the Horner scheme
+ */
+template<typename Poly1_,typename Boundary_>
+boost::numeric::interval
+<typename CGAL::Coercion_traits<typename Poly1_::NT, Boundary_>::Type >
+evaluate_iv(Poly1_ f,boost::numeric::interval<Boundary_> iv) {
+
+    typedef Boundary_ Boundary;
+    typedef Poly1_ Poly1;
+    typedef typename 
+        CGAL::Coercion_traits<typename Poly1::NT, Boundary_>::Type Coercion;
+	//typename CGAL::Coercion_traits<typename Poly1::NT, Boundary_>::Cast cast;
+    CGAL_assertion(f.degree()>=0);
+#if !AcX_USE_AFFINE_ARITHMETIC
+	typename CGAL::Coercion_traits<typename Poly1::NT, Boundary_>::Cast cast;
+    typedef boost::numeric::interval<Coercion> Coercion_interval;
+    Coercion_interval iv_cast(cast(iv.lower()),cast(iv.upper()));
+    int n=f.degree();
+    Coercion_interval ret(cast(f[n]),cast(f[n]));
+    for(int i=n-1;i>=0;i--) {
+        ret *= iv_cast;
+        ret += Coercion_interval(cast(f[i]),cast(f[i]));
+    }
+    return ret;
+#else
+	typedef boost::numeric::interval<Coercion> Coercion_interval;
+	int n=f.degree();
+	Boundary x0 = (iv.lower() + iv.upper()) / 2;
+	Boundary x1 = iv.upper() - x0;
+
+	//! Now here four test-implementations of Affine Arithmetic
+	// AF1 is the first form, which should avoid the fact, that affine forms
+	// grow by doing non-affine operations
+	if(AcX_USE_AFFINE_ARITHMETIC == 1)
+	{
+		//! AF 1
+		//std::cout << "Using Affine Arithmetic -> AF1" << std::endl;
+		Boundary y0 = f[n];
+		Boundary y1 = 0, y2 = 0, spread = 0;
+		for(int i = n-1; i >= 0; i--)
+		{
+			y2 = abs(x0)*y2 + abs(x1)*(abs(y1)+abs(y2));
+			y1 = x0*y1 + x1*y0;
+			y0 = x0*y0 + f[i];
+		}
+		spread = abs(y1) + abs(y2);
+		return  Coercion_interval(y0-spread, y0+spread);
+	}
+	// AF2 is the second form, which should avoid the fact described in AF1
+	// and to decrease the difficulty of even powers
+	else if(AcX_USE_AFFINE_ARITHMETIC == 2)
+	{
+		//! AF 2
+		//std::cout << "Using Affine Arithmetic -> AF2" << std::endl;
+		Boundary y0 = f[n];
+		Boundary y1 = 0, k2 = 0, k3 = 0, e1 = 0, e2 = 0, e3 = 0, low = 0, high = 0, tmp = 0;
+		for(int i = n-1; i >= 0; i--)
+		{
+			if(x0 > 0)
+			{
+				k2 = x0 * e2;
+				k3 = x0 * e3;
+			}
+			else
+			{
+				k2 = -x0 * e3;
+				k3 = -x0 * e2;
+			}
+			tmp = x1 * y1;
+			if(tmp > 0)
+			{
+				k2 = k2 + tmp;
+			}
+			else
+			{
+				k3 = k3 + abs(tmp);
+			}
+			e1 = abs(x0)*e1 + abs(x1)*(e1 + e2 + e3);
+			e2 = k2;
+			e3 = k3;
+			y1 = x0*y1 + x1*y0;
+			y0 = x0*y0 + f[i];
+		}
+		low = abs(y1) + e1 + e3;
+		high = abs(y1) + e1 + e2;
+		return  Coercion_interval(y0-low, y0+high);
+	}
+	// QF (quadratic form) is the third form, which should avoid the facts of
+	// AF1 and decreases better the problem of even powers than AF2
+	else if(AcX_USE_AFFINE_ARITHMETIC == 3)
+	{
+		//! QF
+		//std::cout << "Using Affine Arithmetic -> QF" << std::endl;
+		Boundary y0 = f[n];
+		Boundary y1 = 0, z1 = 0, e1 = 0, spread = 0, low = 0, high = 0;
+		for(int i = n-1; i >= 0; i--)
+		{
+			e1 = abs(x0)*e1 + abs(x1)*e1 + abs(x1*z1);
+			z1 = x0*z1 + x1*y1;
+			y1 = x0*y1 + x1*y0;
+			y0 = x0*y0 + f[i];
+		}
+		spread = abs(y1) + e1;
+		if(z1 > 0)
+		{
+			high = spread + z1;
+		}
+		else
+		{
+			low = spread - z1;
+		}
+		return  Coercion_interval(y0-low, y0+high);
+	}
+	// Here we create a own version of the implementation for Affine Arithmetic
+	// We sort the powers to null, even and odd and compute the intervals
+	else if(AcX_USE_AFFINE_ARITHMETIC == 4)
+	{
+		//! AF own creation
+		//std::cout << "Using Affine Arithmetic -> own AF" << std::endl;
+		Boundary mid_point = (iv.lower() + iv.upper()) / 2;
+		Boundary width = iv.upper() - mid_point;
+		typedef CGAL::Polynomial<Boundary> Rat_poly_1;
+		typedef CGAL::Polynomial<Coercion> Coer_poly_1;
+		Rat_poly_1 affine_poly(mid_point, width);
+		Coer_poly_1 poly_aa = f.evaluate(affine_poly);
+		int n = poly_aa.degree();
+		CGAL_assertion(n>=0);
+		Coercion low = poly_aa[0];
+		Coercion high = poly_aa[0];
+		for(int i = 1; i <= n ; i++)
+		{
+			CGAL::Sign sign = CGAL::sign(poly_aa[i]);
+			if(sign == CGAL::ZERO)
+			{
+				continue;
+			}
+			if(i%2 == 0)
+			{
+				if(sign == CGAL::POSITIVE)
+				{
+					high = high + poly_aa[i];
+				}
+				else
+				{
+					low = low + poly_aa[i];
+				}
+			}
+			else
+			{
+				if(sign == CGAL::POSITIVE)
+				{
+					low = low - poly_aa[i];
+					high = high + poly_aa[i];
+				}
+				else
+				{
+					low = low + poly_aa[i];
+					high = high - poly_aa[i];
+				}
+			}
+		}
+		return  boost::numeric::interval<Coercion>(low, high);
+	}
+#endif
+}
+
+//! \brief replaces the \c NiX::is_root_of function, using \c AcX::ntl_gcd
+template<typename Polynomial,typename AlgebraicReal>
+bool is_root_of(const AlgebraicReal & x,Polynomial p) {
+    if(x.is_rational()) {
+        typedef typename AlgebraicReal::Rational Rational;
+        Rational exact_value = x.rational();
+        typedef typename CGAL::Coercion_traits<Rational, 
+            typename Polynomial::NT>::Type
+            Type;
+        return CGAL::sign(p.evaluate(exact_value))==CGAL::ZERO; 
+    }
+    else {
+        CGAL_precondition(p.degree()>=0);
+      
+        if(p.degree()==0) return p.is_zero();
+        Polynomial g=CGAL::CGALi::gcd_utcf(p,x.polynomial());
+        return g.sign_at(x.low())!=g.sign_at(x.high());
+      
+    }
+    return false; 
+}
+
+/*
+ * \brief Removes the leading term of the polynomial \c f as long as it
+ * vanishes at \c alpha
+ *
+ */
+template<typename Poly_2, typename Algebraic_real>
+Poly_2 poly_non_vanish_leading_term(const Poly_2& pol,Algebraic_real alpha) {
+    Poly_2 f(pol);
+    while(true) {
+	if(CGAL::CGALi::is_root_of(alpha,f.lcoeff())) {
+            typename Poly_2::const_iterator poly_end = f.end();
+            if(f.begin()==poly_end) {
+                break;
+            }
+            poly_end--;
+            f=Poly_2(f.begin(),poly_end);
+	}
+	else {
+            break;
+	}
+    }
+    return f;
+}
+	  
+} // namespace CGALi
+
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_ALG_REAL_UTILS

Algebraic_kernel_d/include/CGAL/Algebraic_curve_kernel_2/Bitstream_descartes_at_x/enums.h

@@ -0,0 +1,41 @@
+// TODO: Add licence
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL:$
+// $Id: $
+// 
+//
+// Author(s)     : Michael Kerber <mkerber@mpi-inf.mpg.de>
+//
+// ============================================================================
+
+
+#ifndef AcX_ENUMS_H
+#define AcX_ENUMS_H 1
+
+CGAL_BEGIN_NAMESPACE
+
+namespace CGALi {
+
+  enum Three_valued {
+
+    ROOT_OF_FIRST_SET = 1,
+    ROOT_OF_BOTH_SETS = 0,
+    ROOT_OF_SECOND_SET=-1,
+
+  };
+
+  enum Degeneracy_strategy {
+
+    SHEAR_STRATEGY = 0,
+    EXCEPTION_STRATEGY = 1,
+
+  };
+
+} // namespace CGALi
+
+CGAL_END_NAMESPACE
+
+#endif