From 715f6af296ee0ed576299e28eac06318debed77b Mon Sep 17 00:00:00 2001
From: Imbus <>
Date: Sun, 1 Jun 2025 14:34:58 +0200
Subject: [PATCH] Some math related code for calculatin binomial coef,
 nth-fibonacci and sin

---
 bincoef.c | 38 ++++++++++++++++++++++++++++++++++
 fibmat.c  | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 sin.c     | 36 ++++++++++++++++++++++++++++++++
 3 files changed, 136 insertions(+)
 create mode 100644 bincoef.c
 create mode 100644 fibmat.c
 create mode 100644 sin.c

diff --git a/bincoef.c b/bincoef.c
new file mode 100644
index 0000000..7d945e2
--- /dev/null
+++ b/bincoef.c
@@ -0,0 +1,38 @@
+#include <stdio.h>
+
+/**
+ * @brief Computes the binomial coefficient "n choose k" (nCk).
+ *
+ * This function calculates the number of ways to choose k elements from a set
+ * of n elements without repetition and without order. It uses an efficient
+ * multiplicative approach to avoid large intermediate factorials.
+ *
+ * @param n The total number of elements.
+ * @param k The number of elements to choose.
+ * @return The computed binomial coefficient (n choose k), or 0 if k > n.
+ */
+unsigned long long binomial_coefficient(unsigned int n, unsigned int k);
+
+unsigned long long binomial_coefficient(unsigned int n, unsigned int k) {
+    if (k > n)
+        return 0;
+    if (k == 0 || k == n)
+        return 1;
+
+    if (k > n - k)
+        k = n - k;
+
+    unsigned long long result = 1;
+    for (unsigned int i = 1; i <= k; ++i) {
+        result *= n - (k - i);
+        result /= i;
+    }
+
+    return result;
+}
+
+int main() {
+    unsigned int n = 10, k = 3;
+    printf("C(%u, %u) = %llu\n", n, k, binomial_coefficient(n, k));
+    return 0;
+}
diff --git a/fibmat.c b/fibmat.c
new file mode 100644
index 0000000..2937a31
--- /dev/null
+++ b/fibmat.c
@@ -0,0 +1,62 @@
+#include <stdio.h>
+#include <stdint.h>
+
+/**
+ * @brief Computes the n-th Fibonacci number using matrix exponentiation.
+ *
+ * This implementation uses the identity:
+ *   [F(n+1) F(n)  ] = [1 1]^n
+ *   [F(n)   F(n-1)]   [1 0]
+ *
+ * The matrix is exponentiated in O(log n) time using exponentiation by squaring.
+ *
+ * @param n The index of the Fibonacci number to compute.
+ * @return The n-th Fibonacci number.
+ */
+uint64_t fibonacci_matrix(uint32_t n);
+
+// 2x2 matrix structure for Fibonacci computation
+typedef struct {
+    uint64_t a, b;
+    uint64_t c, d;
+} FibMatrix;
+
+/**
+ * @brief Multiplies two 2x2 matrices.
+ */
+static FibMatrix matrix_multiply(FibMatrix x, FibMatrix y) {
+    FibMatrix result;
+    result.a = x.a * y.a + x.b * y.c;
+    result.b = x.a * y.b + x.b * y.d;
+    result.c = x.c * y.a + x.d * y.c;
+    result.d = x.c * y.b + x.d * y.d;
+    return result;
+}
+
+/**
+ * @brief Raises a 2x2 matrix to the power of n using exponentiation by squaring.
+ */
+static FibMatrix matrix_power(FibMatrix base, uint32_t n) {
+    FibMatrix result = {1, 0, 0, 1}; // Identity matrix
+    while (n > 0) {
+        if (n % 2 == 1)
+            result = matrix_multiply(result, base);
+        base = matrix_multiply(base, base);
+        n /= 2;
+    }
+    return result;
+}
+
+uint64_t fibonacci_matrix(uint32_t n) {
+    if (n == 0) return 0;
+    FibMatrix base = {1, 1, 1, 0};
+    FibMatrix result = matrix_power(base, n - 1);
+    return result.a;
+}
+
+int main() {
+    for (uint32_t i = 0; i <= 20; ++i) {
+        printf("F(%u) = %lu\n", i, fibonacci_matrix(i));
+    }
+    return 0;
+}
diff --git a/sin.c b/sin.c
new file mode 100644
index 0000000..3b67808
--- /dev/null
+++ b/sin.c
@@ -0,0 +1,36 @@
+#include <stdio.h>
+
+#define HALFPI 1.5707963268
+
+// Compute factorial iteratively
+double factorial(int n) {
+    double result = 1.0;
+    for (int i = 2; i <= n; ++i) {
+        result *= i;
+    }
+    return result;
+}
+
+double abs_double(double x) { return x < 0 ? -x : x; }
+
+// SICP-style iterative approximation for sin(x)
+double sin_iter(double x) {
+    double term = x;   // First term of the series
+    double sum = term; // Initial sum
+    // double prev_sum;
+    int n = 1; // Starting from x^3/3!
+
+    do {
+        // prev_sum = sum;
+        term *= -x * x / ((2 * n) * (2 * n + 1)); // Next term in series
+        sum += term;
+        ++n;
+    } while (abs_double(term) > 1e-10); // Stop when term is sufficiently small
+
+    return sum;
+}
+
+int main() {
+    printf("Approximated sin(pi/2) = %.10f\n", sin_iter(HALFPI));
+    return 0;
+}