diff --git a/Sample/HardParticle/Pyramid2.cpp b/Sample/HardParticle/Pyramid2.cpp
index 83718a588d846826f1aaf6d984c3aabacaf1a0fa..3b42525201af0dbb75dd88f4a3f2dc386b10d734 100644
--- a/Sample/HardParticle/Pyramid2.cpp
+++ b/Sample/HardParticle/Pyramid2.cpp
@@ -64,8 +64,8 @@ std::string Pyramid2::validate() const
     double W = m_width / 2;
     double w = m_width / 2 * (1 - s);
 
-    double zcom =
-        m_height * (.5 - (r + s) / 3 + r * s / 4) / (1 - (r + s) / 2 + r * s / 3); // center of mass
+    // center of mass
+    double zcom = m_height * 0.5 * (6 - 4 * (r + s) + 3 * r * s) / (6 - 3 * (r + s) + 2 * r * s);
 
     setPolyhedron(topology, -zcom,
                   {// base:
diff --git a/Sample/HardParticle/Pyramid3.cpp b/Sample/HardParticle/Pyramid3.cpp
index aec2e912ca388e3faade4b949afb7ca32ab70dd4..a43c46f3a02a8d995758a79803cb6cbda879b604 100644
--- a/Sample/HardParticle/Pyramid3.cpp
+++ b/Sample/HardParticle/Pyramid3.cpp
@@ -55,14 +55,15 @@ std::string Pyramid3::validate() const
 
     double a = m_base_edge;
     double as = a / 2;
-    double ac = a / sqrt(3) / 2;
-    double ah = a / sqrt(3);
+    double ac = a * (1. / sqrt(3) / 2);
+    double ah = a * (1. / sqrt(3));
     double b = a * (1 - r);
     double bs = b / 2;
-    double bc = b / sqrt(3) / 2;
-    double bh = b / sqrt(3);
+    double bc = b * (1. / sqrt(3) / 2);
+    double bh = b * (1. / sqrt(3));
 
-    double zcom = m_height * (.5 - 2 * r / 3 + r * r / 4) / (1 - r + r * r / 3); // center of mass
+    // center of mass
+    double zcom = m_height * (1. / 4) * (6 - 8 * r + 3 * r * r) / (3 - 3 * r + r * r);
 
     setPolyhedron(topology, -zcom,
                   {// base:
diff --git a/Sample/HardParticle/Pyramid4.cpp b/Sample/HardParticle/Pyramid4.cpp
index 81437173a5aebb3ee101904a20495629ea138f22..43497b046332dbb4b69452d78010fc4ca2154e32 100644
--- a/Sample/HardParticle/Pyramid4.cpp
+++ b/Sample/HardParticle/Pyramid4.cpp
@@ -58,7 +58,8 @@ std::string Pyramid4::validate() const
     double a = m_base_edge / 2;
     double b = a * (1 - r);
 
-    double zcom = m_height * (.5 - 2 * r / 3 + r * r / 4) / (1 - r + r * r / 3); // center of mass
+    // center of mass
+    double zcom = m_height * (1. / 4) * (6 - 8 * r + 3 * r * r) / (3 - 3 * r + r * r);
 
     setPolyhedron(topology, -zcom,
                   {// base: