diff --git a/man.md b/man.md
index 9da1922..00c61ab 100644
--- a/man.md
+++ b/man.md
@@ -3,7 +3,7 @@ title: dc
 section: 1
 header: General Commands Manual
 footer: Marco Cetica
-date: March 25, 2024
+date: March 26, 2024
 ---
 
 
@@ -179,7 +179,7 @@ Pops one value, computes its square root, and pushes that.
 
 **!**
 
-Pops one value, computes its factorial, and pushes that.
+Pops one non-negative integer, computes its factorial, and pushes that.
 
 **pi**
 
@@ -224,6 +224,43 @@ Pops one value, computes its `acos`, and pushes that.
 
 Pops one value, computes its `atan`, and pushes that.
 
+## Statistics Operations
+**dc** supports various common statistics operations, such as permutations, combinations, mean, standard deviation
+summation, sum of squares and linear regression. All statistics functions are limited to **non-negative integers**.
+Accumulating functions use the `X` register.
+
+**gP**
+
+Pops two non-negative integers(that is, >=0) and computes `P_{y, x}`, that is the number of possible different arrangements of
+`y` different items taken in quantities of `x` items at a time. No item shall occur more than once in an arrangmement and different orders of same `x` 
+items are counted separately. The `y` parameter correspond to the second one value popped while the `x` parameter correspond to the first one popped.
+
+**gC**
+
+Pops two non-negative integers(that is, >=0) and computes `C_{y, x}`, that is the number of possible sets of
+`y` different items taken in quantities of `x` items at a time. No item shall occur more than once in a set and different orders of the same `x`
+items are not counted separately. The `y` parameter correspond to the second one value popped while the `x` parameter correspond to the first one popped.
+
+**gN**
+
+Counts the number of accumulated elements of the `X` register's stack and pushes that.
+
+**gs**
+
+Computes Σx of the `X` register's stack and pushes that.
+
+**gS**
+
+Computes Σx^2 of the `X` register's stack and pushes that.
+
+**gM**
+
+Computes x̄(mean) of the `X` register's stack and pushes that.
+
+**gD**
+
+Computes σ(standard deviation) of the `X` register's stack and pushes that.
+
 ## Base Conversion
 
 **pb**
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 23b34bf..66f2016 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,6 +4,7 @@ set(HEADER_FILES
         eval.h
         macro.h
         mathematics.h
+        statistics.h
         bitwise.h
         operation.h
         stack.h
@@ -14,6 +15,7 @@ set(SOURCE_FILES
         eval.cpp
         macro.cpp
         mathematics.cpp
+        statistics.cpp
         bitwise.cpp
         stack.cpp
         adt.cpp
diff --git a/src/adt.cpp b/src/adt.cpp
index c46bd1e..cf447d9 100644
--- a/src/adt.cpp
+++ b/src/adt.cpp
@@ -1,3 +1,4 @@
+#include <numeric>
 #include "adt.h"
 
 #define GET_X this->stack.back()
@@ -121,6 +122,43 @@ namespace dc {
         return this->stack.size();
     }
 
+    /*
+     * Returns the summation of all items
+     */
+    template<typename T>
+    requires is_num_or_str<T>
+    double Stack<T>::summation() {
+        auto sum = std::accumulate(this->stack.begin(), this->stack.end(), 0.0,
+            [](auto accumulator, const T& val) -> double {
+                if constexpr(std::is_same_v<T, std::string>) {
+                    return accumulator + std::stod(val);
+                } else {
+                    return accumulator + val;
+                }
+            });
+
+        return sum;
+    }
+
+    /*
+     * Returns the summation of squares of all items
+     */
+    template<typename T>
+    requires is_num_or_str<T>
+    double Stack<T>::summation_squared() {
+        auto sum = std::accumulate(this->stack.begin(), this->stack.end(), 0.0,
+            [](auto accumulator, const T& val) -> double {
+                if constexpr(std::is_same_v<T, std::string>) {
+                    return accumulator + (std::stod(val) * std::stod(val));
+                } else {
+                    return accumulator + (val * val);
+                }
+            });
+
+        return sum;
+    }
+
+
     /**
      * Returns true if stack is empty
      * false otherwise
diff --git a/src/adt.h b/src/adt.h
index 88ccf7c..9c84134 100644
--- a/src/adt.h
+++ b/src/adt.h
@@ -22,6 +22,8 @@ namespace dc {
         T& at(std::size_t index);
         T& operator[](std::size_t index);
         std::size_t size();
+        double summation();
+        double summation_squared();
         bool empty();
         [[nodiscard]] const std::vector<T>& get_ref() const;
 
diff --git a/src/eval.cpp b/src/eval.cpp
index 8305d8c..256f1cc 100644
--- a/src/eval.cpp
+++ b/src/eval.cpp
@@ -3,6 +3,7 @@
 #include "adt.cpp"
 #include "eval.h"
 #include "mathematics.h"
+#include "statistics.h"
 #include "bitwise.h"
 #include "stack.h"
 #include "macro.h"
@@ -45,7 +46,12 @@ void Evaluate::init_environment() {
     this->op_factory.emplace("@", MAKE_UNIQUE_PTR(Mathematics, OPType::RND));
     this->op_factory.emplace("$", MAKE_UNIQUE_PTR(Mathematics, OPType::INT));
     // Statistical operations
-
+    this->op_factory.emplace("gP", MAKE_UNIQUE_PTR(Statistics, OPType::PERM));
+    this->op_factory.emplace("gC", MAKE_UNIQUE_PTR(Statistics, OPType::COMB));
+    this->op_factory.emplace("gs", MAKE_UNIQUE_PTR(Statistics, OPType::SUMX));
+    this->op_factory.emplace("gS", MAKE_UNIQUE_PTR(Statistics, OPType::SUMXX));
+    this->op_factory.emplace("gM", MAKE_UNIQUE_PTR(Statistics, OPType::MEAN));
+    this->op_factory.emplace("gD", MAKE_UNIQUE_PTR(Statistics, OPType::SDEV));
     // Bitwise operations
     this->op_factory.emplace("{", MAKE_UNIQUE_PTR(Bitwise, OPType::BAND));
     this->op_factory.emplace("}", MAKE_UNIQUE_PTR(Bitwise, OPType::BOR));
diff --git a/src/mathematics.cpp b/src/mathematics.cpp
index d41de16..b402470 100644
--- a/src/mathematics.cpp
+++ b/src/mathematics.cpp
@@ -245,17 +245,17 @@ std::optional<std::string> Mathematics::fn_mod_exp(dc::Stack<std::string> &stack
     auto n = stack[len];
     auto e = stack[len-1];
     auto b = stack[len-2];
-    auto is_n_num = is_num<double>(n);
-    auto is_e_num = is_num<double>(e);
-    auto is_b_num = is_num<double>(b);
+    auto is_n_num = is_num<long long>(n);
+    auto is_e_num = is_num<long long>(e);
+    auto is_b_num = is_num<long long>(b);
 
     // This functions computes
 	// 		c ≡ b^e (mod n)
     if(is_n_num && is_e_num && is_b_num) {
         stack.copy_xyz();
-        auto modulus = std::stol(stack.pop(true));
-        auto exponent = std::stol(stack.pop(true));
-        auto base = std::stol(stack.pop(true));
+        auto modulus = std::stoll(stack.pop(true));
+        auto exponent = std::stoll(stack.pop(true));
+        auto base = std::stoll(stack.pop(true));
 
         if(modulus == 1) {
             stack.push("0");
@@ -502,14 +502,14 @@ std::optional<std::string> Mathematics::fn_fact(dc::Stack<std::string> &stack, c
     // Check whether the entry is a number
     if(is_x_num) {
         stack.copy_xyz();
-        unsigned long factorial = 1;
+        unsigned long long factorial = 1;
         auto val = std::stod(stack.pop(true));
 
-        if(val == 0.0) {
-            factorial = 1;
+        if(val < 0.0) {
+            return "'!' is not defined for negative numbers";
         }
 
-        for(std::size_t i = 2; i <= (std::size_t)val; i++) {
+        for(unsigned long long i = 2; i <= val; i++) {
             factorial *= i;
         }
 
diff --git a/src/statistics.cpp b/src/statistics.cpp
new file mode 100644
index 0000000..261b313
--- /dev/null
+++ b/src/statistics.cpp
@@ -0,0 +1,214 @@
+#include <cmath>
+#include <iomanip>
+
+#include "adt.cpp"
+#include "statistics.h"
+#include "is_num.h"
+
+std::optional<std::string> Statistics::exec(dc::Stack<std::string> &stack, dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs) {
+    std::optional<std::string> err = std::nullopt;
+
+    switch(this->op_type) {
+        case OPType::PERM: err = fn_perm(stack, parameters); break;
+        case OPType::COMB: err = fn_comb(stack, parameters); break;
+        case OPType::SUMX: err = fn_sum(stack, parameters, regs); break;
+        case OPType::SUMXX: err = fn_sum_squared(stack, parameters, regs); break;
+        case OPType::MEAN: err = fn_mean(stack, parameters, regs); break;
+        case OPType::SDEV: err = fn_sdev(stack, parameters, regs); break;
+        default: break;
+    }
+
+    return err;
+}
+
+std::optional<std::string> Statistics::fn_perm(dc::Stack<std::string> &stack, const dc::Parameters &parameters) {
+    // Check if stack has enough elements
+    if(stack.size() < 2) {
+        return "'gP' requires two operands";
+    }
+
+    // Extract two entries from the stack
+    auto len = stack.size()-1;
+    auto head = stack[len];
+    auto second = stack[len-1];
+    auto is_head_num = is_num<long long>(head);
+    auto is_second_num = is_num<long long>(second);
+
+    // Check whether both entries are integers
+    if(is_head_num && is_second_num) {
+        stack.copy_xyz();
+        auto x = std::stoll(stack.pop(true));
+        auto y = std::stoll(stack.pop(true));
+
+        // Compute factorial of numerator and denominator
+        // and check if result is a non-negative integer
+        auto numerator_opt = factorial(y);
+        auto denominator_opt = factorial(y - x);
+        if(numerator_opt == std::nullopt || denominator_opt == std::nullopt) {
+            return "'gP' requires positive integers";
+        }
+
+        unsigned long long permutation = numerator_opt.value() / denominator_opt.value();
+        stack.push(trim_digits(permutation, parameters.precision));
+    } else {
+        return "'gP' requires integers values";
+    }
+
+    return std::nullopt;
+}
+
+std::optional<std::string> Statistics::fn_comb(dc::Stack<std::string> &stack, const dc::Parameters &parameters) {
+    // Check if stack has enough elements
+    if(stack.size() < 2) {
+        return "'gC' requires two operands";
+    }
+
+    // Extract two entries from the stack
+    auto len = stack.size()-1;
+    auto head = stack[len];
+    auto second = stack[len-1];
+    auto is_head_num = is_num<long long>(head);
+    auto is_second_num = is_num<long long>(second);
+
+    // Check whether both entries are integers
+    if(is_head_num && is_second_num) {
+        stack.copy_xyz();
+        auto n = std::stoull(stack.pop(true));
+        auto k = std::stoull(stack.pop(true));
+
+        // Check if combination is non-negative
+        if(n > k) {
+            return "'gC' requires positive integers";
+        }
+
+        // From Knuth's "Seminumerical Algorithms" book
+        unsigned long long combination = 1;
+        for(unsigned long long i = 1; i <= n; i++) {
+            combination *= k--;
+            combination /= i;
+        }
+
+        stack.push(trim_digits(combination, parameters.precision));
+    } else {
+        return "'gC' requires integers values";
+    }
+
+    return std::nullopt;
+}
+
+std::optional<std::string> Statistics::fn_sum(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs) {
+    // Check whether 'x' register exists
+    if(regs.find('X') == regs.end()) {
+        return "Register 'X' is undefined";
+    }
+
+    // Check if register's stack is empty
+    if(regs['X'].stack.empty()) {
+        return "The stack of register 'X' is empty";
+    }
+
+    // Othewise retrieve summation of register's stack
+    auto summation = regs['X'].stack.summation();
+    stack.push(trim_digits(summation, parameters.precision));
+
+    return std::nullopt;
+}
+
+std::optional<std::string> Statistics::fn_sum_squared(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs) {
+    // Check whether 'x' register exists
+    if(regs.find('X') == regs.end()) {
+        return "Register 'X' is undefined";
+    }
+
+    // Check if register's stack is empty
+    if(regs['X'].stack.empty()) {
+        return "The stack of register 'X' is empty";
+    }
+
+    // Othewise retrieve summation of squares of register's stack
+    auto summation_squared = regs['X'].stack.summation_squared();
+    stack.push(trim_digits(summation_squared, parameters.precision));
+
+    return std::nullopt;
+}
+
+std::optional<std::string> Statistics::fn_mean(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs) {
+    // Check whether 'x' register exists
+    if(regs.find('X') == regs.end()) {
+        return "Register 'X' is undefined";
+    }
+
+    // Check if register's stack is empty
+    if(regs['X'].stack.empty()) {
+        return "The stack of register 'X' is empty";
+    }
+
+    // Othewise compute mean
+    auto summation = regs['X'].stack.summation();
+    auto size = regs['X'].stack.size();
+    auto mean = summation / size;
+    stack.push(trim_digits(mean, parameters.precision));
+
+    return std::nullopt;
+}
+
+std::optional<std::string> Statistics::fn_sdev(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs) {
+    // Check whether 'x' register exists
+    if(regs.find('X') == regs.end()) {
+        return "Register 'X' is undefined";
+    }
+
+    // Check if register's stack is empty
+    if(regs['X'].stack.empty()) {
+        return "The stack of register 'X' is empty";
+    }
+
+    // Othewise, compute the mean
+    auto summation = regs['X'].stack.summation();
+    auto count = regs['X'].stack.size();
+    auto mean = summation / count;
+
+    // Then compute the sum of the deviations from the mean and square the result
+    const auto& const_vec = regs['X'].stack.get_ref();
+    double sum_of_deviations = std::accumulate(const_vec.begin(), const_vec.end(), 0.0, 
+        [&](double acc, const std::string& val) {
+            double deviation = std::stod(val) - mean;
+            return acc + std::pow(deviation, 2);
+    });
+    // Then compute the mean of previos values(variance)
+    auto variance = sum_of_deviations / count;
+    
+    // Finally, compute the square root of the variance(standard deviation)
+    auto s_dev = sqrt(variance);
+
+    stack.push(trim_digits(s_dev, parameters.precision));
+    return std::nullopt;
+}
+
+std::optional<unsigned long long> Statistics::factorial(const long long n) {
+    if(n < 0) {
+        return std::nullopt;
+    }
+
+    unsigned long long factorial = 1;
+
+    for(long long i = 2; i <= n; i++) {
+        factorial *= i;
+    }
+
+    return factorial;
+}
+
+std::string Statistics::trim_digits(double number, unsigned int precision) {
+    std::ostringstream oss;
+
+    // Preserve non-zero decimal numbers even when precision is zero
+    if(precision == 0 && std::fmod(number, 1.0) != 0.0) {
+        precision = 2;
+    }
+
+    oss << std::fixed << std::setprecision(static_cast<int>(precision)) << number;
+    std::string s = oss.str();
+
+    return s;
+}
diff --git a/src/statistics.h b/src/statistics.h
new file mode 100644
index 0000000..fd8766e
--- /dev/null
+++ b/src/statistics.h
@@ -0,0 +1,21 @@
+#pragma once
+
+#include "operation.h"
+
+class Statistics : public IOperation {
+public:
+    explicit Statistics(const OPType op_t) : op_type(op_t) {}
+    std::optional<std::string> exec(dc::Stack<std::string> &stack, dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs) override;
+
+private:
+    std::optional<std::string> fn_perm(dc::Stack<std::string> &stack, const dc::Parameters &parameters);
+    std::optional<std::string> fn_comb(dc::Stack<std::string> &stack, const dc::Parameters &parameters);
+    std::optional<std::string> fn_sum(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs);
+    std::optional<std::string> fn_sum_squared(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs);
+    std::optional<std::string> fn_mean(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs);
+    std::optional<std::string> fn_sdev(dc::Stack<std::string> &stack, const dc::Parameters &parameters, std::unordered_map<char, dc::Register> &regs);
+    std::optional<unsigned long long> factorial(const long long n);
+    std::string trim_digits(double number, unsigned int precision);
+
+    OPType op_type;
+};
diff --git a/tests/test_fact b/tests/test_fact
index 9df5fde..a4d290f 100644
--- a/tests/test_fact
+++ b/tests/test_fact
@@ -11,6 +11,11 @@ utest() {
     ACTUAL=$("$PROGRAM" -e '!' 2>&1) || true
     assert_eq "$EXPECTED" "$ACTUAL"
 
+    # Test negative numbers
+    EXPECTED="'!' is not defined for negative numbers"
+    ACTUAL=$("$PROGRAM" -e '-5 !' 2>&1) || true
+    assert_eq "$EXPECTED" "$ACTUAL"
+
     # Test non numerical values
     EXPECTED="'!' requires numeric values"
     ACTUAL=$("$PROGRAM" -e '[ foo ] !' 2>&1) || true