Improved math algorithms: update Extended Euclid and binomial exponent implementations

math/binary_exponent.cpp

@@ -5,67 +5,110 @@
  * Calculate \f$a^b\f$ in \f$O(\log(b))\f$ by converting \f$b\f$ to a
  * binary number. Binary exponentiation is also known as exponentiation by
  * squaring.
- * @note This is a far better approach compared to naive method which
- * provide \f$O(b)\f$ operations.
+ *
+ * @note This is a far better approach compared to the naive method which
+ * provides \f$O(b)\f$ operations.
  *
  * Example:
- * </br>10 in base 2 is 1010.
+ *
+ * 10 in base 2 is 1010.
+ *
  * \f{eqnarray*}{
- * 2^{10_d} &=& 2^{1010_b} = 2^8 * 2^2\\
- * 2^1 &=& 2\\
- * 2^2 &=& (2^1)^2 = 2^2 = 4\\
- * 2^4 &=& (2^2)^2 = 4^2 = 16\\
- * 2^8 &=& (2^4)^2 = 16^2 = 256\\
+ * 2^{10_d} &=& 2^{1010_b} = 2^8 \cdot 2^2 \\
+ * 2^1 &=& 2 \\
+ * 2^2 &=& (2^1)^2 = 4 \\
+ * 2^4 &=& (2^2)^2 = 16 \\
+ * 2^8 &=& (2^4)^2 = 256 \\
  * \f}
- * Hence to calculate 2^10 we only need to multiply \f$2^8\f$ and \f$2^2\f$
- * skipping \f$2^1\f$ and \f$2^4\f$.
+ *
+ * Hence to calculate \f$2^{10}\f$, we only need to multiply \f$2^8\f$ and
+ * \f$2^2\f$, skipping \f$2^1\f$ and \f$2^4\f$.
  */
 
 #include <iostream>
+using namespace std;
 
-/// Recursive function to calculate exponent in \f$O(\log(n))\f$ using binary
+/// Recursive function to calculate exponent in O(log(n)) using binary
 /// exponent.
-int binExpo(int a, int b) {
-    if (b == 0) {
+// Use long long for avoiding overflow.
+long long binary_exponent_recursive(long long base, long long exponent) {
+    if (exponent == 0) {
         return 1;
     }
-    int res = binExpo(a, b / 2);
-    if (b % 2) {
-        return res * res * a;
+    long long half = binary_exponent_recursive(base, exponent / 2);
+    if (exponent % 2) {
+        return half * half * base;
     } else {
-        return res * res;
+        return half * half;
     }
 }
 
 /// Iterative function to calculate exponent in \f$O(\log(n))\f$ using binary
 /// exponent.
-int binExpo_alt(int a, int b) {
-    int res = 1;
-    while (b > 0) {
-        if (b % 2) {
-            res = res * a;
+/// Use long long for avoiding Overflow
+long long binary_exponent_iterative(long long base, long long exponent) {
+    long long res = 1;
+    while (exponent > 0) {
+        if (exponent % 2) {
+            res *= base;
         }
-        a = a * a;
-        b /= 2;
+        base = base * base;
+        exponent >>= 1;
     }
     return res;
 }
 
-/// Main function
-int main() {
-    int a, b;
-    /// Give two numbers a, b
-    std::cin >> a >> b;
-    if (a == 0 && b == 0) {
-        std::cout << "Math error" << std::endl;
-    } else if (b < 0) {
-        std::cout << "Exponent must be positive !!" << std::endl;
-    } else {
-        int resRecurse = binExpo(a, b);
-        /// int resIterate = binExpo_alt(a, b);
+//@brief Self-test examples for verifying
+void test() {
+    // Test 0
+    long long expected0 = 1024;
+    cout << "Test 0" << endl;
+    cout << "Input: base = 2 and exponent = 10" << endl;
+    cout << "Expected: " << expected0 << endl;
+    cout << "Recursive: " << binary_exponent_recursive(2, 10) << endl;
+    cout << "Iterative: " << binary_exponent_iterative(2, 10) << endl;
+    cout << endl;
 
-        /// Result of a^b (where '^' denotes exponentiation)
-        std::cout << resRecurse << std::endl;
-        /// std::cout << resIterate << std::endl;
-    }
+    // Test 1
+    long long expected1 = 2187;
+    cout << "Test 1" << endl;
+    cout << "Input: base = 3 and exponent = 7" << endl;
+    cout << "Expected: " << expected1 << endl;
+    cout << "Recursive: " << binary_exponent_recursive(3, 7) << endl;
+    cout << "Iterative: " << binary_exponent_iterative(3, 7) << endl;
+    cout << endl;
+
+    // Test 2
+    long long expected2 = 16777216;
+    cout << "Test 2" << endl;
+    cout << "Input: base = 4 and exponent = 12" << endl;
+    cout << "Expected: " << expected2 << endl;
+    cout << "Recursive: " << binary_exponent_recursive(4, 12) << endl;
+    cout << "Iterative: " << binary_exponent_iterative(4, 12) << endl;
+    cout << endl;
+
+    // Test 3
+    long long expected3 = 30517578125;
+    cout << "Test 3" << endl;
+    cout << "Input: base = 5 and exponent = 15" << endl;
+    cout << "Expected: " << expected3 << endl;
+    cout << "Recursive: " << binary_exponent_recursive(5, 15) << endl;
+    cout << "Iterative: " << binary_exponent_iterative(5, 15) << endl;
+    cout << endl;
+
+    // Test 4
+    long long expected4 = 3656158440062976;
+    cout << "Test 4" << endl;
+    cout << "Input: base = 6 and exponent = 20" << endl;
+    cout << "Expected: " << expected4 << endl;
+    cout << "Recursive: " << binary_exponent_recursive(6, 20) << endl;
+    cout << "Iterative: " << binary_exponent_iterative(6, 20) << endl;
+    cout << endl;
 }
+
+/// Main function
+// returns 0 in case of no error.
+int main() {
+    test();  // run self-test examples
+    return 0;
+}

math/extended_euclid_algorithm.cpp

@@ -1,97 +1,161 @@
 /**
  * @file
- * @brief GCD using [extended Euclid's algorithm]
+ * @brief GCD using extended Euclid's algorithm
  * (https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm)
  *
- * Finding coefficients of a and b ie x and y in  Bézout's identity
- * \f[\text{gcd}(a, b) = a \times x + b \times y \f]
- * This is also used in finding Modular
- * multiplicative inverse of a number. (A * B)%M == 1 Here B is the MMI of A for
- * given M, so extendedEuclid (A, M) gives B.
+ * Finding coefficients x and y for Bézout's identity:
+ * \f[
+ * \gcd(a, b) = a \times x + b \times y
+ * \f]
+ *
+ * This algorithm is also used to compute the Modular Multiplicative Inverse
+ * (MMI). If (A × B) % M == 1, then B is MMI(A, M), and
+ * extended_euclid_recursive(A, M) provides B.
  */
-#include <algorithm>  // for swap function
+
 #include <iostream>
-#include <cstdint>
+using namespace std;
 
 /**
- * function to update the coefficients per iteration
- * \f[r_0,\,r = r,\, r_0 - \text{quotient}\times r\f]
+ * @brief Recursive Extended Euclid Algorithm
+ *
+ * This method recursively applies the idea that:
+ * \f[
+ * \gcd(a, b) = \gcd(b,\, a \bmod b)
+ * \f]
+ *
+ * Each recursive step reduces the problem size by replacing (a, b)
+ * with (b, a % b). Once b becomes 0, the algorithm has reached the base
+ * case, where gcd = a and the Bézout coefficients (x, y) are known.
  *
- * @param[in,out] r signed or unsigned
- * @param[in,out] r0 signed or unsigned
- * @param[in] quotient  unsigned
+ * On returning from recursion, the coefficients are updated in reverse
+ * order using:
+ * \f[
+ * x = y_1,\qquad
+ * y = x_1 - \left\lfloor \frac{a}{b} \right\rfloor \cdot y_1
+ * \f]
+ *
+ * This back-substitution step reconstructs the Bézout identity
+ * for the original pair (a, b).
+ *
+ * @param[in] a first number
+ * @param[in] b second number
+ * @param[out] gcd greatest common divisor
+ * @param[out] x coefficient of a
+ * @param[out] y coefficient of b
  */
-template <typename T, typename T2>
-inline void update_step(T *r, T *r0, const T2 quotient) {
-    T temp = *r;
-    *r = *r0 - (quotient * temp);
-    *r0 = temp;
+void extended_euclid_recursive(long long a, long long b, long long& gcd,
+                               long long& x, long long& y) {
+    if (b == 0) {
+        gcd = a;
+        x = 1;
+        y = 0;
+        return;
+    }
+
+    long long x1, y1;
+    extended_euclid_recursive(b, a % b, gcd, x1, y1);
+
+    x = y1;
+    y = x1 - (a / b) * y1;
 }
 
 /**
- * Implementation using iterative algorithm from
- * [Wikipedia](https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Pseudocode)
+ * @brief Iterative Extended Euclid Algorithm
+ *
+ * This version performs the same transformations as the recursive method
+ * but unfolds the process into a loop. It repeatedly updates (a, b)
+ * using the identity:
+ * \f[
+ * (a, b) \rightarrow (b,\; a - \lfloor a/b \rfloor \cdot b)
+ * \f]
+ *
+ * Alongside (a, b), it maintains two pairs of coefficients (x0, y0) and
+ * (x1, y1) which evolve according to the same quotient used during each
+ * division step. Once b becomes zero, the surviving values (x0, y0)
+ * form the Bézout coefficients:
+ * \f[
+ * a \times x_0 + b_{\text{original}} \times y_0 = \gcd(a, b)
+ * \f]
  *
- * @param[in] A unsigned
- * @param[in] B unsigned
- * @param[out] GCD unsigned
- * @param[out] x signed
- * @param[out] y signed
+ * This approach avoids recursion and makes the sequence of updates
+ * easier to trace step-by-step.
+ *
+ * @param[in] a first number
+ * @param[in] b second number
+ * @param[out] gcd greatest common divisor
+ * @param[out] x coefficient of a
+ * @param[out] y coefficient of b
  */
-template <typename T1, typename T2>
-void extendedEuclid_1(T1 A, T1 B, T1 *GCD, T2 *x, T2 *y) {
-    if (B > A)
-        std::swap(A, B);  // Ensure that A >= B
-
-    T2 s = 0, s0 = 1;
-    T2 t = 1, t0 = 0;
-    T1 r = B, r0 = A;
-
-    while (r != 0) {
-        T1 quotient = r0 / r;
-        update_step(&r, &r0, quotient);
-        update_step(&s, &s0, quotient);
-        update_step(&t, &t0, quotient);
+void extended_euclid_iterative(long long a, long long b, long long& gcd,
+                               long long& x, long long& y) {
+    long long x0 = 1, y0 = 0, x1 = 0, y1 = 1;
+
+    while (b != 0) {
+        long long q = a / b, temp = b;
+        b = a - q * b;
+        a = temp;
+
+        temp = x1;
+        x1 = x0 - q * x1;
+        x0 = temp;
+
+        temp = y1;
+        y1 = y0 - q * y1;
+        y0 = temp;
     }
-    *GCD = r0;
-    *x = s0;
-    *y = t0;
+
+    gcd = a;
+    x = x0;
+    y = y0;
 }
 
 /**
- * Implementation using recursive algorithm
- *
- * @param[in] A unsigned
- * @param[in] B unsigned
- * @param[out] GCD unsigned
- * @param[in,out] x signed
- * @param[in,out] y signed
+ * @brief Self-test examples in your preferred style
  */
-template <typename T, typename T2>
-void extendedEuclid(T A, T B, T *GCD, T2 *x, T2 *y) {
-    if (B > A)
-        std::swap(A, B);  // Ensure that A >= B
-
-    if (B == 0) {
-        *GCD = A;
-        *x = 1;
-        *y = 0;
-    } else {
-        extendedEuclid(B, A % B, GCD, x, y);
-        T2 temp = *x;
-        *x = *y;
-        *y = temp - (A / B) * (*y);
-    }
+void test() {
+    // Test 0
+    long long gcd0, x0, y0;
+    cout << "Test 0" << endl;
+    cout << "Input: a = 30 and b = 20" << endl;
+
+    extended_euclid_recursive(30, 20, gcd0, x0, y0);
+    cout << "Recursive => gcd: " << gcd0 << "  x: " << x0 << "  y: " << y0
+         << endl;
+
+    extended_euclid_iterative(30, 20, gcd0, x0, y0);
+    cout << "Iterative => gcd: " << gcd0 << "  x: " << x0 << "  y: " << y0
+         << endl;
+
+    // Test 1
+    long long gcd1, x1, y1;
+    cout << "Test 1" << endl;
+    cout << "Input: a = 101 and b = 23" << endl;
+
+    extended_euclid_recursive(101, 23, gcd1, x1, y1);
+    cout << "Recursive => gcd: " << gcd1 << "  x: " << x1 << "  y: " << y1
+         << endl;
+
+    extended_euclid_iterative(101, 23, gcd1, x1, y1);
+    cout << "Iterative => gcd: " << gcd1 << "  x: " << x1 << "  y: " << y1
+         << endl;
+
+    // Test 2
+    long long gcd2, x2, y2;
+    cout << "Test 2" << endl;
+    cout << "Input: a = 55 and b = 34" << endl;
+
+    extended_euclid_recursive(55, 34, gcd2, x2, y2);
+    cout << "Recursive => gcd: " << gcd2 << "  x: " << x2 << "  y: " << y2
+         << endl;
+
+    extended_euclid_iterative(55, 34, gcd2, x2, y2);
+    cout << "Iterative => gcd: " << gcd2 << "  x: " << x2 << "  y: " << y2
+         << endl;
 }
 
 /// Main function
 int main() {
-    uint32_t a, b, gcd;
-    int32_t x, y;
-    std::cin >> a >> b;
-    extendedEuclid(a, b, &gcd, &x, &y);
-    std::cout << gcd << " " << x << " " << y << std::endl;
-    extendedEuclid_1(a, b, &gcd, &x, &y);
-    std::cout << gcd << " " << x << " " << y << std::endl;
+    test();  // run self-test examples
     return 0;
-}
+}