diff --git a/Math/Functions/Functions.php b/Math/Functions/Functions.php
index 70289ad69..ee1780b08 100644
--- a/Math/Functions/Functions.php
+++ b/Math/Functions/Functions.php
@@ -24,6 +24,14 @@ namespace phpOMS\Math\Functions;
  */
 final class Functions
 {
+    /**
+     * Epsilon for float comparison.
+     *
+     * @var float
+     * @since 1.0.0
+     */
+    public const EPSILON = 4.88e-04;
+
     /**
      * Constructor.
      *
@@ -268,7 +276,7 @@ final class Functions
             $sum  += $term / (2 * $i + 1);
 
             ++$i;
-        } while ($sum !== 0.0 && \abs($term / $sum) > 0.0000001);
+        } while ($sum !== 0.0 && \abs($term / $sum) > self::EPSILON);
 
         return 2 / \sqrt(\M_PI) * $sum;
     }
@@ -311,7 +319,7 @@ final class Functions
             $n += 0.5;
             $q1 = $q2;
             $q2 = $b / $d;
-        } while (\abs($q1 - $q2) / $q2 > 0.0000001);
+        } while (\abs($q1 - $q2) / $q2 > self::EPSILON);
 
         return 1 / \sqrt(\M_PI) * \exp(-$value * $value) * $q2;
     }
diff --git a/Math/Geometry/Shape/D2/Polygon.php b/Math/Geometry/Shape/D2/Polygon.php
index 76baef5d1..669802e21 100644
--- a/Math/Geometry/Shape/D2/Polygon.php
+++ b/Math/Geometry/Shape/D2/Polygon.php
@@ -30,7 +30,7 @@ final class Polygon implements D2ShapeInterface
      * @var float
      * @since 1.0.0
      */
-    public const EPSILON = 0.00001;
+    public const EPSILON = 4.88e-04;
 
     /**
      * Coordinates.
diff --git a/Math/Matrix/EigenvalueDecomposition.php b/Math/Matrix/EigenvalueDecomposition.php
index 1a24139f7..8d10eed72 100644
--- a/Math/Matrix/EigenvalueDecomposition.php
+++ b/Math/Matrix/EigenvalueDecomposition.php
@@ -30,6 +30,14 @@ use phpOMS\Math\Geometry\Shape\D2\Triangle;
  */
 final class EigenvalueDecomposition
 {
+    /**
+     * Epsilon for float comparison.
+     *
+     * @var float
+     * @since 1.0.0
+     */
+    public const EPSILON = 4.88e-04;
+
     /**
      * Dimension m
      *
@@ -279,14 +287,13 @@ final class EigenvalueDecomposition
 
         $f    = 0.0;
         $tst1 = 0.0;
-        $eps  = 0.00001;
 
         for ($l = 0; $l < $this->m; ++$l) {
             $tst1 = \max($tst1, \abs($this->D[$l]) + \abs($this->E[$l]));
             $m    = $l;
 
             while ($m < $this->m) {
-                if (\abs($this->E[$m]) <= $eps * $tst1) {
+                if (\abs($this->E[$m]) <= self::EPSILON * $tst1) {
                     break;
                 }
 
@@ -344,7 +351,7 @@ final class EigenvalueDecomposition
                     $p           = -$s * $s2 * $c3 * $el1 * $this->E[$l] / $dl1;
                     $this->E[$l] = $s * $p;
                     $this->D[$l] = $c * $p;
-                } while (\abs($this->E[$l]) > $eps * $tst1);
+                } while (\abs($this->E[$l]) > self::EPSILON * $tst1);
             }
 
             $this->D[$l] += $f;
@@ -506,7 +513,6 @@ final class EigenvalueDecomposition
         $n       = $nn - 1;
         $low     = 0;
         $high    = $nn - 1;
-        $eps     = 0.00001;
         $exshift = 0.0;
         $p       = 0;
         $q       = 0;
@@ -535,7 +541,7 @@ final class EigenvalueDecomposition
                     $s = $norm;
                 }
 
-                if (\abs($this->H[$l][$l - 1]) < $eps * $s) {
+                if (\abs($this->H[$l][$l - 1]) < self::EPSILON * $s) {
                     break;
                 }
 
@@ -656,7 +662,7 @@ final class EigenvalueDecomposition
                     $r /= $s;
 
                     if ($m === $l
-                        || \abs($this->H[$m][$m - 1]) * (\abs($q) + \abs($r)) < $eps * (\abs($p) * (\abs($this->H[$m - 1][$m - 1]) + \abs($z) + \abs($this->H[$m + 1][$m + 1])))
+                        || \abs($this->H[$m][$m - 1]) * (\abs($q) + \abs($r)) < self::EPSILON * (\abs($p) * (\abs($this->H[$m - 1][$m - 1]) + \abs($z) + \abs($this->H[$m + 1][$m + 1])))
                     ) {
                         break;
                     }
@@ -774,7 +780,7 @@ final class EigenvalueDecomposition
                         $l = $i;
 
                         if ($this->E[$i] == 0) {
-                            $this->H[$i][$n] = $w != 0 ? -$r / $w : -$r / ($eps * $norm);
+                            $this->H[$i][$n] = $w != 0 ? -$r / $w : -$r / (self::EPSILON * $norm);
                         } else {
                             $x                   = $this->H[$i][$i + 1];
                             $y                   = $this->H[$i + 1][$i];
@@ -785,7 +791,7 @@ final class EigenvalueDecomposition
                         }
 
                         $t = \abs($this->H[$i][$n]);
-                        if (($eps * $t) * $t > 1) {
+                        if ((self::EPSILON * $t) * $t > 1) {
                             for ($j = $i; $j <= $n; ++$j) {
                                 $this->H[$j][$n] = $this->H[$j][$n] / $t;
                             }
@@ -836,7 +842,7 @@ final class EigenvalueDecomposition
                             $vi = ($this->D[$i] - $p) * 2.0 * $q;
 
                             if ($vr == 0 & $vi == 0) {
-                                $vr = $eps * $norm * (\abs($w) + \abs($q) + \abs($x) + \abs($y) + \abs($z));
+                                $vr = self::EPSILON * $norm * (\abs($w) + \abs($q) + \abs($x) + \abs($y) + \abs($z));
                             }
 
                             $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
@@ -855,7 +861,7 @@ final class EigenvalueDecomposition
                         }
 
                         $t = \max(\abs($this->H[$i][$n - 1]), \abs($this->H[$i][$n]));
-                        if (($eps * $t) * $t > 1) {
+                        if ((self::EPSILON * $t) * $t > 1) {
                             for ($j = $i; $j <= $n; ++$j) {
                                 $this->H[$j][$n - 1] = $this->H[$j][$n - 1] / $t;
                                 $this->H[$j][$n]     = $this->H[$j][$n] / $t;
diff --git a/Math/Matrix/Matrix.php b/Math/Matrix/Matrix.php
index 7081e2ae3..0fd51886b 100644
--- a/Math/Matrix/Matrix.php
+++ b/Math/Matrix/Matrix.php
@@ -29,6 +29,14 @@ use phpOMS\Math\Matrix\Exception\InvalidDimensionException;
  */
 class Matrix implements \ArrayAccess, \Iterator
 {
+    /**
+     * Epsilon for float comparison.
+     *
+     * @var float
+     * @since 1.0.0
+     */
+    public const EPSILON = 4.88e-04;
+
     /**
      * Matrix.
      *
@@ -326,7 +334,7 @@ class Matrix implements \ArrayAccess, \Iterator
 
         for ($i = 0; $i < $nDim; ++$i) {
             for ($j = 0; $j < $mDim; ++$j) {
-                if (!$selected[$j] && \abs($matrix[$j][$i]) > 0.0001) {
+                if (!$selected[$j] && \abs($matrix[$j][$i]) > self::EPSILON) {
                     break;
                 }
             }
@@ -342,7 +350,7 @@ class Matrix implements \ArrayAccess, \Iterator
             }
 
             for ($k = 0; $k < $mDim; ++$k) {
-                if ($k !== $j && \abs($matrix[$k][$i]) > 0.0001) {
+                if ($k !== $j && \abs($matrix[$k][$i]) > self::EPSILON) {
                     for ($p = $i + 1; $p < $nDim; ++$p) {
                         $matrix[$k][$p] -= $matrix[$j][$p] * $matrix[$k][$i];
                     }
diff --git a/Math/Number/Numbers.php b/Math/Number/Numbers.php
index e7ba65b66..6b702613c 100644
--- a/Math/Number/Numbers.php
+++ b/Math/Number/Numbers.php
@@ -24,6 +24,10 @@ namespace phpOMS\Math\Number;
  */
 final class Numbers
 {
+    public const SFLOAT = 1.175494351E-38;
+
+    public const EPSILON = 4.88e-04;
+
     /**
      * Constructor.
      *
@@ -94,7 +98,7 @@ final class Numbers
      */
     public static function isSquare(int $n) : bool
     {
-        return \abs(((int) \sqrt($n)) * ((int) \sqrt($n)) - $n) < 0.001;
+        return \abs(((int) \sqrt($n)) * ((int) \sqrt($n)) - $n) < self::EPSILON;
     }
 
     /**
diff --git a/Math/Number/Prime.php b/Math/Number/Prime.php
index 64695bed4..95679ddff 100644
--- a/Math/Number/Prime.php
+++ b/Math/Number/Prime.php
@@ -24,6 +24,14 @@ namespace phpOMS\Math\Number;
  */
 final class Prime
 {
+    /**
+     * Epsilon for float comparison.
+     *
+     * @var float
+     * @since 1.0.0
+     */
+    public const EPSILON = 4.88e-04;
+
     /**
      * Constructor.
      *
@@ -47,7 +55,7 @@ final class Prime
     {
         $mersenne = \log($n + 1, 2);
 
-        return $mersenne - (int) $mersenne < 0.00001;
+        return $mersenne - (int) $mersenne < self::EPSILON;
     }
 
     /**
diff --git a/Math/Stochastic/Distribution/LogDistribution.php b/Math/Stochastic/Distribution/LogDistribution.php
index 29145bebb..3f7ce1907 100644
--- a/Math/Stochastic/Distribution/LogDistribution.php
+++ b/Math/Stochastic/Distribution/LogDistribution.php
@@ -25,6 +25,14 @@ use phpOMS\Math\Functions\Beta;
  */
 final class LogDistribution
 {
+    /**
+     * Epsilon for float comparison.
+     *
+     * @var float
+     * @since 1.0.0
+     */
+    public const EPSILON = 4.88e-04;
+
     /**
      * Get probability mass function.
      *
@@ -53,9 +61,9 @@ final class LogDistribution
     public static function getCdf(float $p, int $k) : float
     {
         // This is a workaround!
-        // Actually 0 should be used instead of 0.0001.
+        // Actually 0 should be used instead of self::EPSILON.
         // This is only used because the incomplete beta function doesn't work for p or q = 0
-        return 1 + Beta::incompleteBeta($p, $k + 1, 0.0001) / \log(1 - $p);
+        return 1 + Beta::incompleteBeta($p, $k + 1, self::EPSILON) / \log(1 - $p);
     }
 
     /**