diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
index 024275b8750..d5f75137eb2 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_polygon_soup.h
@@ -118,43 +118,21 @@ bool simplify_polygon(PointRange& points,
 {
   const std::size_t ini_polygon_size = polygon.size();
 
-  // Start at the last since if two points are identical, the second one gets removed.
-  // By starting at 'last', we ensure that 'to_remove' is ordered from closest to .begin()
-  // to closest to .end()
-  std::size_t last = ini_polygon_size - 1, i = last;
-  bool stop = false;
-  std::vector<std::size_t> to_remove;
-
-  do
+  for(std::size_t i=0; i<polygon.size(); ++i)
   {
-    std::size_t next_i = (i == last) ? 0 : i+1;
-    stop = (next_i == last);
+    const std::size_t s = polygon.size();
+    if(s == 1)
+      break;
 
-    while(polygon[i] == polygon[next_i] || // combinatorial equality
-          traits.equal_3_object()(points[polygon[i]], points[polygon[next_i]])) // geometric equality
+    const std::size_t ni = (i + 1) % s;
+    if(polygon[i] == polygon[ni] ||
+       traits.equal_3_object()(points[polygon[i]], points[polygon[ni]]))
     {
-      to_remove.push_back(next_i);
 #ifdef CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE_PP
-      std::cout << "Duplicate point: polygon[" << next_i << "] = " << polygon[next_i] << std::endl;
+      std::cout << "Duplicate point: polygon[" << ni << "] = " << polygon[ni] << std::endl;
 #endif
-      next_i = (next_i == last) ? 0 : next_i+1;
-
-      // Every duplicate in front of 'last' (circularly-speaking) has already been cleared
-      if(next_i == last)
-      {
-        stop = true;
-        break;
-      }
+      polygon.erase(polygon.begin() + i--);
     }
-
-    i = next_i;
-  }
-  while(!stop);
-
-  while(!to_remove.empty())
-  {
-    polygon.erase(polygon.begin() + to_remove.back());
-    to_remove.pop_back();
   }
 
   const std::size_t removed_points_n = ini_polygon_size - polygon.size();
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_repair_polygon_soup.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_repair_polygon_soup.cpp
index dd1ac38d01a..0a3c01fa52d 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_repair_polygon_soup.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_repair_polygon_soup.cpp
@@ -289,14 +289,24 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
   points.push_back(Point_3(1,3,0)); // #3
   points.push_back(Point_3(0,1,0)); // #4
   points.push_back(Point_3(1,1,0)); // #5
-  points.push_back(Point_3(0,0,0)); // #6
+  points.push_back(Point_3(0,0,0)); // #6 == #0
 
   // ------
   CGAL_polygon polygon;
-  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4); polygon.push_back(0); polygon.push_back(0);
   polygons.push_back(polygon);
 
   std::size_t res = PMP::internal::simplify_polygons_in_polygon_soup<K>(points, polygons);
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
+  assert(res == 1 && polygons.back().size() == 3);
+
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup<K>(points, polygons);
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
   assert(res == 0 && polygons.back().size() == 3);
 
   // ------
@@ -305,6 +315,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
   polygons.push_back(polygon);
 
   res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
   assert(res == 1 && polygons.back().size() == 1);
 
   // ------
@@ -313,6 +324,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
   polygons.push_back(polygon);
 
   res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
   assert(res == 1 && polygons.back().size() == 1);
 
   // ------
@@ -321,8 +333,18 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
   polygons.push_back(polygon);
 
   res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
   assert(res == 1 && polygons.back().size() == 2);
 
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(0); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
+  assert(res == 0 && polygons.back().size() == 4);
+
   // ------
   // Now with the same geometric positions, but different combinatorial information
   polygon.clear();
@@ -330,6 +352,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
   polygons.push_back(polygon);
 
   res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
   assert(res == 1 && polygons.back().size() == 3);
 
   // ------
@@ -338,6 +361,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
   polygons.push_back(polygon);
 
   res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
   assert(res == 1 && polygons.back().size() == 2);
 
   // ------
@@ -346,6 +370,7 @@ void test_simplify_polygons(const bool /*verbose*/ = false)
   polygons.push_back(polygon);
 
   res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  std::cout << "res: " << res << " / size: " << polygons.back().size() << std::endl;
   assert(res == 1 && polygons.back().size() == 3);
 }