Improved efficiency of some contractors

src/core/contractors/codac2_CtcInverse.h

@@ -29,7 +29,10 @@ namespace codac2
         requires IsCtcBaseOrPtr<C,Y>
       CtcInverse(const AnalyticFunction<typename Wrapper<Y>::Domain>& f, const C& ctc_y, bool with_centered_form = true, bool is_not_in = false)
         : _f(f), _ctc_y(ctc_y), _with_centered_form(with_centered_form), _is_not_in(is_not_in)
-      { }
+      {
+        assert_release([&]() { return f.output_size() == size_of(ctc_y); }()
+          && "CtcInverse: invalid dimension of image argument ('y' or 'ctc_y')");
+      }
 
       CtcInverse(const AnalyticFunction<typename Wrapper<Y>::Domain>& f, const Y& y, bool with_centered_form = true, bool is_not_in = false)
         : CtcInverse(f, CtcWrapper_<Y>(y), with_centered_form, is_not_in)

src/core/contractors/codac2_directed_ctc.cpp

@@ -397,8 +397,8 @@ using namespace codac2;
     else*/
     {
       IntervalMatrix Q = gauss_jordan(x1.mid());
-      auto b_tilde = Q*y;
-      auto A_tilde = Q*x1; // should be a tree matrix
+      IntervalVector b_tilde = Q*y;
+      IntervalMatrix A_tilde = Q*x1; // should be a tree matrix
 
       for(int a = 0 ; a < 1 ; a++)
       {
@@ -412,10 +412,7 @@ using namespace codac2;
               if(i != j)
                 u -= x2[j]*A_tilde(k,j);
 
-            if(A_tilde(k,i).contains(0.))
-              continue;
-
-            x2[i] &= u / A_tilde(k,i);
+            MulOp::bwd(u, x2[i], A_tilde(k,i));
           }
         }
       }

src/core/domains/interval/codac2_Interval.h

@@ -794,7 +794,10 @@ namespace codac2
 
     double a = std::max<double>(next_float(-oo),lb());
     double b = std::min<double>(previous_float(oo),ub());
-    return a + (((double)std::rand())/(double)RAND_MAX)*(b-a);
+    double r = a + (((double)std::rand())/(double)RAND_MAX)*(b-a);
+    // The above operation may result in a floating point outside the bounds,
+    // due to floating-point errors. Such possible error is corrected below:
+    return std::max<double>(lb(),std::min<double>(r,ub()));
   }
 
   inline double Interval::rad() const

src/core/functions/analytic/codac2_AnalyticFunction.h

@@ -112,6 +112,26 @@ namespace codac2
         return eval_(x...).da;
       }
 
+      Index output_size() const
+      {
+        if constexpr(std::is_same_v<typename T::Domain,Interval>)
+          return 1;
+
+        else if constexpr(std::is_same_v<typename T::Domain,IntervalVector>)
+        {
+          // A dump evaluation is performed to estimate the dimension
+          // of the image of this function. A natural evaluation is assumed
+          // to be faster.
+          return natural_eval(IntervalVector(this->input_size())).size();
+        }
+
+        else
+        {
+          assert_release(false && "unable to estimate output size");
+          return 0;
+        }
+      }
+
       friend std::ostream& operator<<(std::ostream& os, [[maybe_unused]] const AnalyticFunction<T>& f)
       {
         if constexpr(std::is_same_v<typename T::Domain,Interval>)