diff --git a/DIRECTORY.md b/DIRECTORY.md
index 740a61647..189d6041b 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -15,7 +15,6 @@
   * [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/gaussian_elimination.cpp)
   * [Newton Raphson Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/newton_raphson_method.cpp)
   * [Ordinary Least Squares Regressor](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp)
-  * [Secant Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/secant_method.cpp)
   * [Successive Approximation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/computer_oriented_statistical_methods/successive_approximation.cpp)
 
 ## Data Structure
diff --git a/computer_oriented_statistical_methods/bisection_method.cpp b/computer_oriented_statistical_methods/bisection_method.cpp
index 089828c9a..c93c529d2 100644
--- a/computer_oriented_statistical_methods/bisection_method.cpp
+++ b/computer_oriented_statistical_methods/bisection_method.cpp
@@ -1,25 +1,60 @@
+/**
+ * \file
+ * \brief Solve the equation \f$f(x)=0\f$ using [bisection
+ * method](https://en.wikipedia.org/wiki/Bisection_method)
+ *
+ * Given two points \f$a\f$ and \f$b\f$ such that \f$f(a)<0\f$ and
+ * \f$f(b)>0\f$, then the \f$(i+1)^\text{th}\f$ approximation is given by: \f[
+ * x_{i+1} = \frac{a_i+b_i}{2}
+ * \f]
+ * For the next iteration, the interval is selected
+ * as: \f$[a,x]\f$ if \f$x>0\f$ or \f$[x,b]\f$ if \f$x<0\f$. The Process is
+ * continued till a close enough approximation is achieved.
+ *
+ * \see newton_raphson_method.cpp, false_position.cpp, secant_method.cpp
+ */
 #include <cmath>
 #include <iostream>
+#include <limits>
 
-static float eq(float i) {
+#define EPSILON \
+    1e-6  // std::numeric_limits<double>::epsilon()  ///< system accuracy limit
+#define MAX_ITERATIONS 50000  ///< Maximum number of iterations to check
+
+/** define \f$f(x)\f$ to find root for
+ */
+static double eq(double i) {
     return (std::pow(i, 3) - (4 * i) - 9);  // original equation
 }
 
-int main() {
-    float a, b, x, z;
+/** get the sign of any given number */
+template <typename T>
+int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
 
-    for (int i = 0; i < 100; i++) {
-        z = eq(i);
-        if (z >= 0) {
-            b = i;
-            a = --i;
+/** main function */
+int main() {
+    double a = -1, b = 1, x, z;
+    int i;
+
+    // loop to find initial intervals a, b
+    for (int i = 0; i < MAX_ITERATIONS; i++) {
+        z = eq(a);
+        x = eq(b);
+        if (sgn(z) == sgn(x)) {  // same signs, increase interval
+            b++;
+            a--;
+        } else {  // if opposite signs, we got our interval
             break;
         }
     }
 
     std::cout << "\nFirst initial: " << a;
     std::cout << "\nSecond initial: " << b;
-    for (int i = 0; i < 100; i++) {
+
+    // start iterations
+    for (i = 0; i < MAX_ITERATIONS; i++) {
         x = (a + b) / 2;
         z = eq(x);
         std::cout << "\n\nz: " << z << "\t[" << a << " , " << b
@@ -31,10 +66,10 @@ int main() {
             b = x;
         }
 
-        if (z > 0 && z < 0.0009)  // stoping criteria
+        if (std::abs(z) < EPSILON)  // stoping criteria
             break;
     }
 
-    std::cout << "\n\nRoot: " << x;
+    std::cout << "\n\nRoot: " << x << "\t\tSteps: " << i << std::endl;
     return 0;
 }
diff --git a/computer_oriented_statistical_methods/false_position.cpp b/computer_oriented_statistical_methods/false_position.cpp
index c5a314508..aebd154bb 100644
--- a/computer_oriented_statistical_methods/false_position.cpp
+++ b/computer_oriented_statistical_methods/false_position.cpp
@@ -1,19 +1,53 @@
+/**
+ * \file
+ * \brief Solve the equation \f$f(x)=0\f$ using [false position
+ * method](https://en.wikipedia.org/wiki/Regula_falsi), also known as the Secant
+ * method
+ *
+ * Given two points \f$a\f$ and \f$b\f$ such that \f$f(a)<0\f$ and
+ * \f$f(b)>0\f$, then the \f$(i+1)^\text{th}\f$ approximation is given by: \f[
+ * x_{i+1} = \frac{a_i\cdot f(b_i) - b_i\cdot f(a_i)}{f(b_i) - f(a_i)}
+ * \f]
+ * For the next iteration, the interval is selected
+ * as: \f$[a,x]\f$ if \f$x>0\f$ or \f$[x,b]\f$ if \f$x<0\f$. The Process is
+ * continued till a close enough approximation is achieved.
+ *
+ * \see newton_raphson_method.cpp, bisection_method.cpp
+ */
 #include <cmath>
 #include <cstdlib>
 #include <iostream>
+#include <limits>
 
-static float eq(float i) {
-    return (pow(i, 3) - (4 * i) - 9);  // origial equation
+#define EPSILON \
+    1e-6  // std::numeric_limits<double>::epsilon()  ///< system accuracy limit
+#define MAX_ITERATIONS 50000  ///< Maximum number of iterations to check
+
+/** define \f$f(x)\f$ to find root for
+ */
+static double eq(double i) {
+    return (std::pow(i, 3) - (4 * i) - 9);  // origial equation
 }
 
+/** get the sign of any given number */
+template <typename T>
+int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+
+/** main function */
 int main() {
-    float a, b, z, c, m, n;
-    system("clear");
-    for (int i = 0; i < 100; i++) {
-        z = eq(i);
-        if (z >= 0) {
-            b = i;
-            a = --i;
+    double a = -1, b = 1, x, z, m, n, c;
+    int i;
+
+    // loop to find initial intervals a, b
+    for (int i = 0; i < MAX_ITERATIONS; i++) {
+        z = eq(a);
+        x = eq(b);
+        if (sgn(z) == sgn(x)) {  // same signs, increase interval
+            b++;
+            a--;
+        } else {  // if opposite signs, we got our interval
             break;
         }
     }
@@ -21,18 +55,20 @@ int main() {
     std::cout << "\nFirst initial: " << a;
     std::cout << "\nSecond initial: " << b;
 
-    for (int i = 0; i < 100; i++) {
-        float h, d;
+    for (i = 0; i < MAX_ITERATIONS; i++) {
         m = eq(a);
         n = eq(b);
+
         c = ((a * n) - (b * m)) / (n - m);
+
         a = c;
         z = eq(c);
-        if (z > 0 && z < 0.09) {  // stoping criteria
+
+        if (std::abs(z) < EPSILON) {  // stoping criteria
             break;
         }
     }
 
-    std::cout << "\n\nRoot: " << c;
+    std::cout << "\n\nRoot: " << c << "\t\tSteps: " << i << std::endl;
     return 0;
 }
diff --git a/computer_oriented_statistical_methods/gaussian_elimination.cpp b/computer_oriented_statistical_methods/gaussian_elimination.cpp
index 0b8bb693d..60b5648ec 100644
--- a/computer_oriented_statistical_methods/gaussian_elimination.cpp
+++ b/computer_oriented_statistical_methods/gaussian_elimination.cpp
@@ -1,17 +1,26 @@
+/**
+ * \file
+ * \brief [Gaussian elimination
+ * method](https://en.wikipedia.org/wiki/Gaussian_elimination)
+ */
 #include <iostream>
 
+/** Main function */
 int main() {
     int mat_size, i, j, step;
 
     std::cout << "Matrix size: ";
     std::cin >> mat_size;
 
+    // create a 2D matrix by dynamic memory allocation
     double **mat = new double *[mat_size + 1], **x = new double *[mat_size];
     for (i = 0; i <= mat_size; i++) {
         mat[i] = new double[mat_size + 1];
-        if (i < mat_size) x[i] = new double[mat_size + 1];
+        if (i < mat_size)
+            x[i] = new double[mat_size + 1];
     }
 
+    // get the matrix elements from user
     std::cout << std::endl << "Enter value of the matrix: " << std::endl;
     for (i = 0; i < mat_size; i++) {
         for (j = 0; j <= mat_size; j++) {
@@ -20,6 +29,7 @@ int main() {
         }
     }
 
+    // perform Gaussian elimination
     for (step = 0; step < mat_size - 1; step++) {
         for (i = step; i < mat_size - 1; i++) {
             double a = (mat[i + 1][step] / mat[step][step]);
@@ -56,7 +66,8 @@ int main() {
 
     for (i = 0; i <= mat_size; i++) {
         delete[] mat[i];
-        if (i < mat_size) delete[] x[i];
+        if (i < mat_size)
+            delete[] x[i];
     }
     delete[] mat;
     delete[] x;
diff --git a/computer_oriented_statistical_methods/newton_raphson_method.cpp b/computer_oriented_statistical_methods/newton_raphson_method.cpp
index 47d276490..318363a39 100644
--- a/computer_oriented_statistical_methods/newton_raphson_method.cpp
+++ b/computer_oriented_statistical_methods/newton_raphson_method.cpp
@@ -1,42 +1,57 @@
+/**
+ * \file
+ * \brief Solve the equation \f$f(x)=0\f$ using [Newton-Raphson
+ * method](https://en.wikipedia.org/wiki/Newton%27s_method)
+ *
+ * The \f$(i+1)^\text{th}\f$ approximation is given by:
+ * \f[
+ * x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}
+ * \f]
+ *
+ * \see bisection_method.cpp, false_position.cpp
+ */
 #include <cmath>
+#include <ctime>
 #include <iostream>
+#include <limits>
 
-static float eq(float i) {
+#define EPSILON \
+    1e-6  // std::numeric_limits<double>::epsilon()  ///< system accuracy limit
+#define MAX_ITERATIONS 50000  ///< Maximum number of iterations to check
+
+/** define \f$f(x)\f$ to find root for
+ */
+static double eq(double i) {
     return (std::pow(i, 3) - (4 * i) - 9);  // original equation
 }
 
-static float eq_der(float i) {
+/** define the derivative function \f$f'(x)\f$
+ */
+static double eq_der(double i) {
     return ((3 * std::pow(i, 2)) - 4);  // derivative of equation
 }
 
+/** Main function */
 int main() {
-    float a, b, z, c, m, n;
+    std::srand(std::time(nullptr));  // initialize randomizer
 
-    for (int i = 0; i < 100; i++) {
-        z = eq(i);
-        if (z >= 0) {
-            b = i;
-            a = --i;
-            break;
-        }
-    }
+    double z, c = std::rand() % 100, m, n;
+    int i;
 
-    std::cout << "\nFirst initial: " << a;
-    std::cout << "\nSecond initial: " << b;
-    c = (a + b) / 2;
+    std::cout << "\nInitial approximation: " << c;
 
-    for (int i = 0; i < 100; i++) {
-        float h;
+    // start iterations
+    for (i = 0; i < MAX_ITERATIONS; i++) {
         m = eq(c);
         n = eq_der(c);
 
         z = c - (m / n);
         c = z;
 
-        if (m > 0 && m < 0.009)  // stoping criteria
+        if (std::abs(m) < EPSILON)  // stoping criteria
             break;
     }
 
-    std::cout << "\n\nRoot: " << z << std::endl;
+    std::cout << "\n\nRoot: " << z << "\t\tSteps: " << i << std::endl;
     return 0;
 }
diff --git a/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp b/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp
index 06bd4ea52..de02b27bb 100644
--- a/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp
+++ b/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp
@@ -1,11 +1,13 @@
 /**
  * @file
+ * \brief Linear regression example using [Ordinary least
+ * squares](https://en.wikipedia.org/wiki/Ordinary_least_squares)
  *
  * Program that gets the number of data samples and number of features per
  * sample along with output per sample. It applies OLS regression to compute
  * the regression output for additional test data samples.
  */
-#include <iomanip>
+#include <iomanip>  // for print formatting
 #include <iostream>
 #include <vector>
 
@@ -52,7 +54,8 @@ inline bool is_square(std::vector<std::vector<T>> const &A) {
     // Assuming A is square matrix
     size_t N = A.size();
     for (size_t i = 0; i < N; i++)
-        if (A[i].size() != N) return false;
+        if (A[i].size() != N)
+            return false;
     return true;
 }
 
@@ -265,7 +268,8 @@ std::vector<std::vector<float>> get_inverse(
         inverse[row] = inverse[row] / divisor;
         // Row transformations
         for (size_t row2 = 0; row2 < N; row2++) {
-            if (row2 == row) continue;
+            if (row2 == row)
+                continue;
             float factor = temp[row2][row];
             temp[row2] = temp[row2] - factor * temp[row];
             inverse[row2] = inverse[row2] - factor * inverse[row];
diff --git a/computer_oriented_statistical_methods/secant_method.cpp b/computer_oriented_statistical_methods/secant_method.cpp
deleted file mode 100644
index c353ef850..000000000
--- a/computer_oriented_statistical_methods/secant_method.cpp
+++ /dev/null
@@ -1,37 +0,0 @@
-#include <cmath>
-#include <iostream>
-
-static float eq(float i) {
-    return (pow(i, 3) - (4 * i) - 9);  // original equation
-}
-
-int main() {
-    float a, b, z, c, m, n;
-    for (int i = 0; i < 100; i++) {
-        z = eq(i);
-        if (z >= 0) {
-            b = i;
-            a = --i;
-            break;
-        }
-    }
-
-    std::cout << "\nFirst initial: " << a;
-    std::cout << "\nSecond initial: " << b;
-    for (int i = 0; i < 100; i++) {
-        float h, d;
-        m = eq(a);
-        n = eq(b);
-
-        c = ((a * n) - (b * m)) / (n - m);
-        a = b;
-        b = c;
-
-        z = eq(c);
-        if (z > 0 && z < 0.09)  // stoping criteria
-            break;
-    }
-
-    std::cout << "\n\nRoot: " << c;
-    return 0;
-}
diff --git a/computer_oriented_statistical_methods/successive_approximation.cpp b/computer_oriented_statistical_methods/successive_approximation.cpp
index efbcc3bbf..351382f24 100644
--- a/computer_oriented_statistical_methods/successive_approximation.cpp
+++ b/computer_oriented_statistical_methods/successive_approximation.cpp
@@ -1,9 +1,22 @@
+/**
+ * \file
+ * \brief Method of successive approximations using [fixed-point
+ * iteration](https://en.wikipedia.org/wiki/Fixed-point_iteration) method
+ */
 #include <cmath>
 #include <iostream>
 
-static float eq(float y) { return ((3 * y) - (cos(y)) - 2); }
-static float eqd(float y) { return ((0.5) * ((cos(y)) + 2)); }
+/** equation 1
+ * \f[f(y) = 3y - \cos y -2\f]
+ */
+static float eq(float y) { return (3 * y) - cos(y) - 2; }
 
+/** equation 2
+ * \f[f(y) = \frac{\cos y+2}{2}\f]
+ */
+static float eqd(float y) { return 0.5 * (cos(y) + 2); }
+
+/** Main function */
 int main() {
     float y, x1, x2, x3, sum, s, a, f1, f2, gd;
     int i, n;