diff --git a/.vscode/settings.json b/.vscode/settings.json
index 074c4ab03..67fe06477 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -1,6 +1,64 @@
 {
-     "C_Cpp.clang_format_style": "{ BasedOnStyle: Google, UseTab: Never, IndentWidth: 4, TabWidth: 4, AllowShortIfStatementsOnASingleLine: false, IndentCaseLabels: true, ColumnLimit: 80, AccessModifierOffset: -3, AlignConsecutiveMacros: true }",
-     "editor.formatOnSave": true,
-     "editor.formatOnType": true,
-     "editor.formatOnPaste": true
+  "C_Cpp.clang_format_style": "{ BasedOnStyle: Google, UseTab: Never, IndentWidth: 4, TabWidth: 4, AllowShortIfStatementsOnASingleLine: false, IndentCaseLabels: true, ColumnLimit: 80, AccessModifierOffset: -3, AlignConsecutiveMacros: true }",
+  "editor.formatOnSave": true,
+  "editor.formatOnType": true,
+  "editor.formatOnPaste": true,
+  "files.associations": {
+    "array": "cpp",
+    "atomic": "cpp",
+    "*.tcc": "cpp",
+    "bitset": "cpp",
+    "cctype": "cpp",
+    "chrono": "cpp",
+    "cinttypes": "cpp",
+    "clocale": "cpp",
+    "cmath": "cpp",
+    "complex": "cpp",
+    "cstdarg": "cpp",
+    "cstddef": "cpp",
+    "cstdint": "cpp",
+    "cstdio": "cpp",
+    "cstdlib": "cpp",
+    "cstring": "cpp",
+    "ctime": "cpp",
+    "cwchar": "cpp",
+    "cwctype": "cpp",
+    "deque": "cpp",
+    "list": "cpp",
+    "unordered_map": "cpp",
+    "unordered_set": "cpp",
+    "vector": "cpp",
+    "exception": "cpp",
+    "algorithm": "cpp",
+    "functional": "cpp",
+    "iterator": "cpp",
+    "map": "cpp",
+    "memory": "cpp",
+    "memory_resource": "cpp",
+    "numeric": "cpp",
+    "optional": "cpp",
+    "random": "cpp",
+    "ratio": "cpp",
+    "set": "cpp",
+    "string": "cpp",
+    "string_view": "cpp",
+    "system_error": "cpp",
+    "tuple": "cpp",
+    "type_traits": "cpp",
+    "utility": "cpp",
+    "fstream": "cpp",
+    "initializer_list": "cpp",
+    "iomanip": "cpp",
+    "iosfwd": "cpp",
+    "iostream": "cpp",
+    "istream": "cpp",
+    "limits": "cpp",
+    "new": "cpp",
+    "ostream": "cpp",
+    "sstream": "cpp",
+    "stdexcept": "cpp",
+    "streambuf": "cpp",
+    "typeinfo": "cpp",
+    "valarray": "cpp"
+  }
 }
diff --git a/DIRECTORY.md b/DIRECTORY.md
index 781a4cbe0..d2c85ebda 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -223,11 +223,13 @@
 ## Numerical Methods
   * [Bisection Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/bisection_method.cpp)
   * [Brent Method Extrema](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/brent_method_extrema.cpp)
+  * [Composite Simpson Rule](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/composite_simpson_rule.cpp)
   * [Durand Kerner Roots](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/durand_kerner_roots.cpp)
   * [False Position](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/false_position.cpp)
   * [Fast Fourier Transform](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/fast_fourier_transform.cpp)
   * [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/gaussian_elimination.cpp)
   * [Golden Search Extrema](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/golden_search_extrema.cpp)
+  * [Inverse Fast Fourier Transform](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/inverse_fast_fourier_transform.cpp)
   * [Lu Decompose](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/lu_decompose.cpp)
   * [Lu Decomposition](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/lu_decomposition.h)
   * [Midpoint Integral Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/midpoint_integral_method.cpp)
diff --git a/numerical_methods/composite_simpson_rule.cpp b/numerical_methods/composite_simpson_rule.cpp
new file mode 100644
index 000000000..9a5e8c180
--- /dev/null
+++ b/numerical_methods/composite_simpson_rule.cpp
@@ -0,0 +1,202 @@
+/**
+ * @file
+ * @brief Implementation of the Composite Simpson Rule for the approximation
+ *
+ * @details The following is an implementation of the Composite Simpson Rule for
+ * the approximation of definite integrals. More info -> wiki:
+ * https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule
+ *
+ * The idea is to split the interval in an EVEN number N of intervals and use as
+ * interpolation points the xi for which it applies that xi = x0 + i*h, where h
+ * is a step defined as h = (b-a)/N where a and b are the first and last points
+ * of the interval of the integration [a, b].
+ *
+ * We create a table of the xi and their corresponding f(xi) values and we
+ * evaluate the integral by the formula: I = h/3 * {f(x0) + 4*f(x1) + 2*f(x2) +
+ * ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)}
+ *
+ * That means that the first and last indexed i f(xi) are multiplied by 1,
+ * the odd indexed f(xi) by 4 and the even by 2.
+ *
+ * In this program there are 4 sample test functions f, g, k, l that are
+ * evaluated in the same interval.
+ *
+ * Arguments can be passed as parameters from the command line argv[1] = N,
+ * argv[2] = a, argv[3] = b
+ *
+ * N must be even number and a<b.
+ *
+ * In the end of the main() i compare the program's result with the one from
+ * mathematical software with 2 decimal points margin.
+ *
+ * Add sample function by replacing one of the f, g, k, l and the assert
+ *
+ * @author [ggkogkou](https://github.com/ggkogkou)
+ *
+ */
+
+#include <cassert>     /// for assert
+#include <cmath>       /// for math functions
+#include <cstdint>     /// for integer allocation
+#include <cstdlib>     /// for std::atof
+#include <functional>  /// for std::function
+#include <iostream>    /// for IO operations
+#include <map>         /// for std::map container
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical algorithms/methods
+ */
+namespace numerical_methods {
+/**
+ * @namespace simpson_method
+ * @brief Contains the Simpson's method implementation
+ */
+namespace simpson_method {
+/**
+ * @fn double evaluate_by_simpson(int N, double h, double a,
+ * std::function<double (double)> func)
+ * @brief Calculate integral or assert if integral is not a number (Nan)
+ * @param N number of intervals
+ * @param h step
+ * @param a x0
+ * @param func: choose the function that will be evaluated
+ * @returns the result of the integration
+ */
+double evaluate_by_simpson(std::int32_t N, double h, double a,
+                           std::function<double(double)> func) {
+    std::map<std::int32_t, double>
+        data_table;  // Contains the data points. key: i, value: f(xi)
+    double xi = a;   // Initialize xi to the starting point x0 = a
+
+    // Create the data table
+    double temp;
+    for (std::int32_t i = 0; i <= N; i++) {
+        temp = func(xi);
+        data_table.insert(
+            std::pair<std::int32_t, double>(i, temp));  // add i and f(xi)
+        xi += h;  // Get the next point xi for the next iteration
+    }
+
+    // Evaluate the integral.
+    // Remember: f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)
+    double evaluate_integral = 0;
+    for (std::int32_t i = 0; i <= N; i++) {
+        if (i == 0 || i == N)
+            evaluate_integral += data_table.at(i);
+        else if (i % 2 == 1)
+            evaluate_integral += 4 * data_table.at(i);
+        else
+            evaluate_integral += 2 * data_table.at(i);
+    }
+
+    // Multiply by the coefficient h/3
+    evaluate_integral *= h / 3;
+
+    // If the result calculated is nan, then the user has given wrong input
+    // interval.
+    assert(!std::isnan(evaluate_integral) &&
+           "The definite integral can't be evaluated. Check the validity of "
+           "your input.\n");
+    // Else return
+    return evaluate_integral;
+}
+
+/**
+ * @fn double f(double x)
+ * @brief A function f(x) that will be used to test the method
+ * @param x The independent variable xi
+ * @returns the value of the dependent variable yi = f(xi)
+ */
+double f(double x) { return std::sqrt(x) + std::log(x); }
+/** @brief Another test function */
+double g(double x) { return std::exp(-x) * (4 - std::pow(x, 2)); }
+/** @brief Another test function */
+double k(double x) { return std::sqrt(2 * std::pow(x, 3) + 3); }
+/** @brief Another test function*/
+double l(double x) { return x + std::log(2 * x + 1); }
+}  // namespace simpson_method
+}  // namespace numerical_methods
+
+/**
+ * \brief Self-test implementations
+ * @param N is the number of intervals
+ * @param h is the step
+ * @param a is x0
+ * @param b is the end of the interval
+ * @param used_argv_parameters is 'true' if argv parameters are given and
+ * 'false' if not
+ */
+static void test(std::int32_t N, double h, double a, double b,
+                 bool used_argv_parameters) {
+    // Call the functions and find the integral of each function
+    double result_f = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::f);
+    assert((used_argv_parameters || (result_f >= 4.09 && result_f <= 4.10)) &&
+           "The result of f(x) is wrong");
+    std::cout << "The result of integral f(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_f << std::endl;
+
+    double result_g = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::g);
+    assert((used_argv_parameters || (result_g >= 0.27 && result_g <= 0.28)) &&
+           "The result of g(x) is wrong");
+    std::cout << "The result of integral g(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_g << std::endl;
+
+    double result_k = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::k);
+    assert((used_argv_parameters || (result_k >= 9.06 && result_k <= 9.07)) &&
+           "The result of k(x) is wrong");
+    std::cout << "The result of integral k(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_k << std::endl;
+
+    double result_l = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::l);
+    assert((used_argv_parameters || (result_l >= 7.16 && result_l <= 7.17)) &&
+           "The result of l(x) is wrong");
+    std::cout << "The result of integral l(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_l << std::endl;
+}
+
+/**
+ * @brief Main function
+ * @param argc commandline argument count (ignored)
+ * @param argv commandline array of arguments (ignored)
+ * @returns 0 on exit
+ */
+int main(int argc, char** argv) {
+    std::int32_t N = 16;  /// Number of intervals to divide the integration
+                          /// interval. MUST BE EVEN
+    double a = 1, b = 3;  /// Starting and ending point of the integration in
+                          /// the real axis
+    double h;             /// Step, calculated by a, b and N
+
+    bool used_argv_parameters =
+        false;  // If argv parameters are used then the assert must be omitted
+                // for the tst cases
+
+    // Get user input (by the command line parameters or the console after
+    // displaying messages)
+    if (argc == 4) {
+        N = std::atoi(argv[1]);
+        a = (double)std::atof(argv[2]);
+        b = (double)std::atof(argv[3]);
+        // Check if a<b else abort
+        assert(a < b && "a has to be less than b");
+        assert(N > 0 && "N has to be > 0");
+        if (N < 16 || a != 1 || b != 3)
+            used_argv_parameters = true;
+        std::cout << "You selected N=" << N << ", a=" << a << ", b=" << b
+                  << std::endl;
+    } else
+        std::cout << "Default N=" << N << ", a=" << a << ", b=" << b
+                  << std::endl;
+
+    // Find the step
+    h = (b - a) / N;
+
+    test(N, h, a, b, used_argv_parameters);  // run self-test implementations
+
+    return 0;
+}
diff --git a/numerical_methods/inverse_fast_fourier_transform.cpp b/numerical_methods/inverse_fast_fourier_transform.cpp
new file mode 100644
index 000000000..d2248be7b
--- /dev/null
+++ b/numerical_methods/inverse_fast_fourier_transform.cpp
@@ -0,0 +1,161 @@
+/**
+ * @file
+ * @brief [An inverse fast Fourier transform
+ * (IFFT)](https://www.geeksforgeeks.org/python-inverse-fast-fourier-transformation/)
+ * is an algorithm that computes the inverse fourier transform.
+ * @details
+ * This algorithm has an application in use case scenario where a user wants find coefficients of
+ * a function in a short time by just using points generated by DFT.
+ * Time complexity
+ * this algorithm computes the IDFT in O(nlogn) time in comparison to traditional O(n^2).
+ * @author [Ameya Chawla](https://github.com/ameyachawlaggsipu)
+ */
+
+#include <cassert>   /// for assert
+#include <cmath>     /// for mathematical-related functions
+#include <complex>   /// for storing points and coefficents
+#include <iostream>  /// for IO operations
+#include <vector>    /// for std::vector
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical algorithms/methods
+ */
+namespace numerical_methods {
+/**
+ * @brief InverseFastFourierTransform is a recursive function which returns list of
+ * complex numbers
+ * @param p List of Coefficents in form of complex numbers
+ * @param n Count of elements in list p
+ * @returns p if n==1
+ * @returns y if n!=1
+ */
+std::complex<double> *InverseFastFourierTransform(std::complex<double> *p, uint8_t n) {
+    if (n == 1) {
+        return p;  /// Base Case To return
+    }
+
+    double pi = 2 * asin(1.0);  /// Declaring value of pi
+
+    std::complex<double> om = std::complex<double>(
+        cos(2 * pi / n), sin(2 * pi / n));  /// Calculating value of omega
+        
+    om.real(om.real()/n); /// One change in comparison with DFT
+    om.imag(om.imag()/n); /// One change in comparison with DFT
+
+    auto *pe = new std::complex<double>[n / 2];  /// Coefficients of even power
+
+    auto *po = new std::complex<double>[n / 2];  /// Coefficients of odd power
+
+    int k1 = 0, k2 = 0;
+    for (int j = 0; j < n; j++) {
+        if (j % 2 == 0) {
+            pe[k1++] = p[j];  /// Assigning values of even Coefficients
+
+        } else
+            po[k2++] = p[j];  /// Assigning value of odd Coefficients
+    }
+
+    std::complex<double> *ye =
+        InverseFastFourierTransform(pe, n / 2);  /// Recursive Call
+
+    std::complex<double> *yo =
+        InverseFastFourierTransform(po, n / 2);  /// Recursive Call
+
+    auto *y = new std::complex<double>[n];  /// Final value representation list
+
+    k1 = 0, k2 = 0;
+
+    for (int i = 0; i < n / 2; i++) {
+        y[i] =
+            ye[k1] + pow(om, i) * yo[k2];  /// Updating the first n/2 elements
+        y[i + n / 2] =
+            ye[k1] - pow(om, i) * yo[k2];  /// Updating the last n/2 elements
+
+        k1++;
+        k2++;
+    }
+ 
+    if(n!=2){
+     
+        delete[] pe;
+        delete[] po;
+        
+    }
+
+    delete[] ye;  /// Deleting dynamic array ye
+    delete[] yo;  /// Deleting dynamic array yo
+    return y;
+}
+
+}  // namespace numerical_methods
+
+/**
+ * @brief Self-test implementations
+ * @details
+ * Declaring two test cases and checking for the error
+ * in predicted and true value is less than 0.000000000001.
+ * @returns void
+ */
+static void test() {
+    /* descriptions of the following test */
+
+    auto *t1 = new std::complex<double>[2];  /// Test case 1
+    auto *t2 = new std::complex<double>[4];  /// Test case 2
+
+    t1[0] = {3, 0};
+    t1[1] = {-1, 0};
+    t2[0] = {10, 0};
+    t2[1] = {-2, -2};
+    t2[2] = {-2, 0};
+    t2[3] = {-2, 2};
+
+    uint8_t n1 = 2;
+    uint8_t n2 = 4;
+    std::vector<std::complex<double>> r1 = {
+        {1, 0}, {2, 0}};  /// True Answer for test case 1
+
+    std::vector<std::complex<double>> r2 = {
+        {1, 0}, {2, 0}, {3, 0}, {4, 0}};  /// True Answer for test case 2
+
+    std::complex<double> *o1 = numerical_methods::InverseFastFourierTransform(t1, n1);
+    
+    std::complex<double> *o2 = numerical_methods::InverseFastFourierTransform(t2, n2);
+
+    for (uint8_t i = 0; i < n1; i++) {
+        assert((r1[i].real() - o1[i].real() < 0.000000000001) &&
+               (r1[i].imag() - o1[i].imag() <
+                0.000000000001));  /// Comparing for both real and imaginary
+                                   /// values for test case 1
+        
+    }
+
+    for (uint8_t i = 0; i < n2; i++) {
+        assert((r2[i].real() - o2[i].real() < 0.000000000001) &&
+               (r2[i].imag() - o2[i].imag() <
+                0.000000000001));  /// Comparing for both real and imaginary
+                                   /// values for test case 2
+        
+    }
+    
+    
+    delete[] t1;
+    delete[] t2;
+    delete[] o1;
+    delete[] o2;
+    std::cout << "All tests have successfully passed!\n";
+}
+
+/**
+ * @brief Main function
+ * @param argc commandline argument count (ignored)
+ * @param argv commandline array of arguments (ignored)
+ * calls automated test function to test the working of fast fourier transform.
+ * @returns 0 on exit
+ */
+
+int main(int argc, char const *argv[]) {
+    test();  //  run self-test implementations
+             //  with 2 defined test cases
+    return 0;
+}