diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h
index 73dc2983aef..458d23466e3 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h
@@ -1039,19 +1039,21 @@ bool polygon_soup_snap_rounding(PointRange &points,
 *     \cgalParamDefault{false}
 *   \cgalParamNEnd
 *   \cgalParamNBegin{snap_grid_size}
-*     \cgalParamDescription{A value `gs` used to scale the points to `[-2^gs, 2^gs]` before rounding them on integers. Used only if `apply_iterative_snap_rounding()` is set to `true`.}
+*     \cgalParamDescription{A value `gs` used to scale the points to `[-2^gs, 2^gs]` before rounding them on integers. Used only if `apply_iterative_snap_rounding()` is set to `true`}
 *     \cgalParamType{unsigned int}
 *     \cgalParamDefault{23}
 *     \cgalParamExtra{Must be lower than 52.}
 *   \cgalParamNEnd
 *   \cgalParamNBegin{number_of_iterations}
-*     \cgalParamDescription{Maximum number of iteration performed by the snap algorithm. Use only if `apply_iterative_snap_rounding` is true.}
+*     \cgalParamDescription{Maximum number of iterations performed by the snap algorithm. Used only if `apply_iterative_snap_rounding` is true.}
 *     \cgalParamType{unsigned int}
 *     \cgalParamDefault{15}
 *   \cgalParamNEnd
 * \cgalNamedParamsEnd
 *
-* \return `true` if the modified triangle soup is free from self-intersection, and `false` otherwise. The return value is only meaningful when the option `apply_snap_rounding()` is set to true.
+* \return `true` if `apply_iterative_snap_rounding` is set to `false`; otherwise, return `true` if the modified triangle soup is free from
+* self-intersection, and `false` if the algorithm was unable to provide such a triangle soup within the number of iterations.
+*
 */
 template <class PointRange, class TriangleRange, class NamedParameters = parameters::Default_named_parameters>
 bool autorefine_triangle_soup(PointRange& soup_points,
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h
index e7a2d5c6821..c1fa810db6c 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/triangle_soup_snap_rounding.h
@@ -22,7 +22,6 @@
 
 #include <CGAL/Fraction_traits.h>
 #include <CGAL/Lazy_exact_nt.h>
-#include <CGAL/Gmpq.h>
 
 namespace CGAL
 {
@@ -33,11 +32,11 @@ namespace Polygon_mesh_processing
 namespace internal
 {
 
-template <class NT> double ceil(Lazy_exact_nt< NT > x);
-template <class NT> double ceil(NT x);
+template <class NT> double double_ceil(Lazy_exact_nt< NT > x);
+template <class NT> double double_ceil(NT x);
 
 template <class NT>
-double ceil(Lazy_exact_nt< NT > x){
+double double_ceil(Lazy_exact_nt< NT > x){
   // If both sides are in the same ceil, return this ceil
   double ceil_left=std::ceil(to_interval(x).first);
   if(ceil_left==std::ceil(to_interval(x).second))
@@ -48,11 +47,11 @@ double ceil(Lazy_exact_nt< NT > x){
   if(ceil_left==std::ceil(to_interval(x).second))
     return ceil_left;
   // If not return the ceil of the exact value
-  return ceil( x.exact());
+  return double_ceil( x.exact());
 };
 
 template <class NT>
-double ceil(NT x){
+double double_ceil(NT x){
   using FT = Fraction_traits<NT>;
   if constexpr(std::is_same<typename FT::Is_fraction, Tag_true>::value){
     // If NT is a fraction, the ceil value is the result of the euclidian division of the numerator and the denominator.
@@ -60,8 +59,8 @@ double ceil(NT x){
     typename FT::Denominator_type denom;
     typename FT::Decompose()(x,num,denom);
     div_mod(num, denom, r, e);
-    if((r<0) && e!=0) //If the result is negative, the ceil value is one below
-      return to_double(r-1);
+    if((r>=0) && e!=0) //If the result is positive, the ceil value is one above
+      return to_double(r+1);
     return to_double(r);
   } else {
     // Return the ceil of the approximation
@@ -184,7 +183,7 @@ bool polygon_soup_snap_rounding(PointRange &points,
   {
     // Scale the coordinate, round to nearest integer and scale back
     // TODO replace this ceil by the one of Algebraic_fondation when it will be add to master
-    return internal::ceil((x * scale) - 0.5) / scale;
+    return internal::double_ceil((x * scale) - 0.5) / scale;
     // return ceil((x * scale) - 0.5) / scale;
   };
   auto snap_p = [scale, snap](const Point_3 &p)