Give Eigen a scoring strategy for pivot selection with intervals.

Number_types/include/CGAL/Interval_nt.h

@@ -52,6 +52,7 @@
 #include <CGAL/FPU.h>
 #include <CGAL/IO/io.h>
 #include <iostream>
+#include <boost/operators.hpp>
 
 namespace CGAL {
 
@@ -1274,6 +1275,46 @@ namespace Eigen {
     template<bool b>
       struct significant_decimals_impl<CGAL::Interval_nt<b> >
       : significant_decimals_impl<typename CGAL::Interval_nt<b>::value_type> { };
+
+    // Without this, when computing some decompositions for a matrix of
+    // intervals, Eigen looks for the largest element in a column (for
+    // instance). There may easily be 2 equal, slightly imprecise numbers that
+    // could equally well be used as pivots, but Eigen ends up spuriously
+    // throwing in the comparison between them. So we provide a different
+    // strategy for picking the pivot.
+    template<typename> struct scalar_score_coeff_op;
+    template<bool b> struct scalar_score_coeff_op<CGAL::Interval_nt<b> > {
+      // If all coeffs can be 0, it is essential to designate as the best one
+      // that can be non-zero and has a non-zero score, if there is one.
+      struct result_type : boost::totally_ordered1<result_type> {
+        CGAL::Interval_nt<b> i;
+        result_type():i(){}
+        result_type(CGAL::Interval_nt<b> j):i(j){}
+        friend bool operator<(result_type x, result_type y){
+          if(x.i.inf()==0){
+            if(y.i.inf()==0)return x.i.sup()<y.i.sup(); // [0,0]<[0,1]
+            else return true; // [0,*]<[1,*]
+          }
+#if 0
+          // The following is already handled by the general formula below
+          if(y.i.inf()==0)return false; // [0,*]<[1,*]
+#endif
+          // Both numbers are guaranteed non-zero. With double people usually
+          // pick the biggest number. Here we choose the tightest interval.
+          // This is purely heuristic, it doesn't matter much if overflow makes
+          // us do random choices.
+          // Best is largest inf/sup (ideally 1)
+          // Risk of {over,under}flow
+          return x.i.inf()*y.i.sup() < y.i.inf()*x.i.sup();
+        }
+        // Only used as: if(max==Score(0))
+        friend bool operator==(result_type x, result_type y){
+          // Throw if we don't know if the max coeff is 0
+          return x.i == y.i;
+        }
+      };
+      result_type operator()(CGAL::Interval_nt<b> const&x)const{return abs(x);}
+    };
   }
 }