diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_cache.h b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_cache.h
index 62c7950a028..829f00d2a5f 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_cache.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_cache.h
@@ -653,7 +653,7 @@ void _Bezier_cache<NtTraits>::_self_intersection_params
   //   II: Y(t) - Y(s) / (t - s) = 0
   //
   Integer                 *coeffs;
-  int                      i, k;
+  int                      i;
 
   // Consruct the bivariate polynomial that corresponds to Equation I.
   // Note that we represent a bivariate polynomial as a vector of univariate
@@ -667,12 +667,11 @@ void _Bezier_cache<NtTraits>::_self_intersection_params
   coeffs = new Integer [degX];
 
   for (i = 0; i < degX; i++)
-  {
-    for (k = i + 1; k < degX; k++)
-      coeffs[k - i - 1] = nt_traits.get_coefficient (polyX, k);
+    coeffs[i] = nt_traits.get_coefficient(polyX, i + 1);
 
-    coeffsX_st[i] = nt_traits.construct_polynomial (coeffs, degX - i - 1);
-  }
+  for (i = 0; i < degX; i++)
+    coeffsX_st[degX - i - 1] =
+        nt_traits.construct_polynomial(coeffs + i, degX - i - 1);
 
   delete[] coeffs;
 
@@ -685,12 +684,11 @@ void _Bezier_cache<NtTraits>::_self_intersection_params
   coeffs = new Integer [degY];
 
   for (i = 0; i < degY; i++)
-  {
-    for (k = i + 1; k < degY; k++)
-      coeffs[k - i - 1] = nt_traits.get_coefficient (polyY, k);
+    coeffs[i] = nt_traits.get_coefficient(polyY, i + 1);
 
-    coeffsY_st[i] = nt_traits.construct_polynomial (coeffs, degY - i - 1);
-  }
+  for (i = 0; i < degY; i++)
+    coeffsY_st[degY - i - 1] =
+        nt_traits.construct_polynomial(coeffs + i, degY - i - 1);
 
   delete[] coeffs;
 
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_point_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_point_2.h
index 62312301940..feb5def190e 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_point_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_point_2.h
@@ -1641,17 +1641,43 @@ void _Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_make_exact
   const Algebraic      t_min = nt_traits.convert (orig2.point_bound().t_min);
   const Algebraic      t_max = nt_traits.convert (orig2.point_bound().t_max);
 
+  bool self_intersecting = (org_it1->curve().id() == org_it2->curve().id());
+
   for (intr_it = intr_list.begin(); intr_it != intr_list.end(); ++intr_it)
   {
-    if (CGAL::compare (intr_it->s, s_min) != SMALLER &&
-        CGAL::compare (intr_it->s, s_max) != LARGER &&
-        CGAL::compare (intr_it->t, t_min) != SMALLER &&
-        CGAL::compare (intr_it->t, t_max) != LARGER)
+    auto in_bounding_interval =
+      [](const auto& s_, const auto& s_min_, const auto& s_max_) -> bool {
+      return CGAL::compare(s_, s_min_) != SMALLER &&
+             CGAL::compare(s_, s_max_) != LARGER;
+    };
+
+    bool st_in_st_range = in_bounding_interval(intr_it->s, s_min, s_max) &&
+                          in_bounding_interval(intr_it->t, t_min, t_max);
+    bool ts_in_st_range = false;
+
+    if (st_in_st_range)
     {
       // Update the originators.
-      orig1.set_parameter (intr_it->s);
-      orig2.set_parameter (intr_it->t);
+      orig1.set_parameter(intr_it->s);
+      orig2.set_parameter(intr_it->t);
+    }
+    else if (self_intersecting)
+    {
+      // check whether s is in t range, and t is in s range
+      // s and t can be interchanged in case of self intersections
+      ts_in_st_range = in_bounding_interval(intr_it->t, s_min, s_max) &&
+                       in_bounding_interval(intr_it->s, t_min, t_max);
 
+      if (ts_in_st_range)
+      {
+        // Update the originators.
+        orig1.set_parameter(intr_it->t);
+        orig2.set_parameter(intr_it->s);
+      }
+    }
+
+    if (st_in_st_range || ts_in_st_range)
+    {
       // Set the exact point coordinates.
       p_alg_x = new Algebraic (intr_it->x);
       p_alg_y = new Algebraic (intr_it->y);
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_x_monotone_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_x_monotone_2.h
index 5397a42062f..41594edcc2b 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_x_monotone_2.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_geometry_traits/Bezier_x_monotone_2.h
@@ -1139,10 +1139,11 @@ _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::compare_to_left
     Originator_iterator  org = p.get_originator(_curve, _xid);
 
     CGAL_assertion(org != p.originators_end());
-    CGAL_assertion(_inc_to_right != cv._inc_to_right);
 
     if (org->point_bound().type == Bez_point_bound::VERTICAL_TANGENCY_PT)
     {
+      CGAL_assertion(_inc_to_right != cv._inc_to_right);
+
       if (! p.is_exact())
       {
         // Comparison based on the control polygon of the bounded vertical