From de2e8838aea1bd67a333743e3ff056cf4fb045a3 Mon Sep 17 00:00:00 2001
From: "Yukihiro \"Matz\" Matsumoto" <matz@ruby.or.jp>
Date: Fri, 18 Jul 2025 07:44:37 +0900
Subject: [PATCH] mruby-bigint: replace euclidean gcd with binary gcd algorithm

Replace the traditional Euclidean GCD algorithm with Stein's binary GCD algorithm
for improved performance on large numbers:

- Implement binary GCD (Stein's algorithm) avoiding expensive division operations
- Use bit shifts and subtraction instead of modulo operations
- Handle special cases (zero values) efficiently
- Preserve common factors of 2 for correct results
- Maintain full compatibility with existing rational number functionality

Binary GCD is significantly faster for large numbers as it avoids the costly
division operations used in the Euclidean algorithm, using only bit operations,
addition, and subtraction.

Co-authored-by: Claude <noreply@anthropic.com>
---
 mrbgems/mruby-bigint/core/bigint.c | 67 ++++++++++++++++++++++++++----
 1 file changed, 59 insertions(+), 8 deletions(-)

diff --git a/mrbgems/mruby-bigint/core/bigint.c b/mrbgems/mruby-bigint/core/bigint.c
index 18551703b..68f4d7c8c 100644
--- a/mrbgems/mruby-bigint/core/bigint.c
+++ b/mrbgems/mruby-bigint/core/bigint.c
@@ -1253,21 +1253,72 @@ mpz_abs(mrb_state *mrb, mpz_t *x, mpz_t *y)
     x->sn = 1;
 }
 
+/* Binary GCD algorithm (Stein's algorithm) - faster than Euclidean GCD */
 static void
 mpz_gcd(mrb_state *mrb, mpz_t *gg, mpz_t *aa, mpz_t *bb)
 {
   mpz_t a, b, t;
-  mpz_abs(mrb, &a, aa); mpz_abs(mrb, &b, bb);
+
+  /* Handle special cases */
+  if (zero_p(aa)) {
+    mpz_abs(mrb, gg, bb);
+    return;
+  }
+  if (zero_p(bb)) {
+    mpz_abs(mrb, gg, aa);
+    return;
+  }
+
+  mpz_abs(mrb, &a, aa);
+  mpz_abs(mrb, &b, bb);
   mpz_init(mrb, &t);
 
-  while (b.sn != 0) {
-    mpz_mod(mrb, &t, &a, &b);
-    mpz_set(mrb, &a, &b);
-    mpz_set(mrb, &b, &t);
+  /* Find power of 2 that divides both a and b */
+  size_t shift = 0;
+  while ((a.p[0] & 1) == 0 && (b.p[0] & 1) == 0) {
+    mpz_div_2exp(mrb, &a, &a, 1);
+    mpz_div_2exp(mrb, &b, &b, 1);
+    shift++;
   }
-  trim(&a);
-  mpz_move(mrb, gg, &a);
-  mpz_clear(mrb, &b);
+
+  /* Make a odd */
+  while ((a.p[0] & 1) == 0) {
+    mpz_div_2exp(mrb, &a, &a, 1);
+  }
+
+  /* From here on, a is always odd */
+  do {
+    /* Make b odd */
+    while ((b.p[0] & 1) == 0) {
+      mpz_div_2exp(mrb, &b, &b, 1);
+    }
+
+    /* Now both a and b are odd. Ensure a >= b */
+    if (mpz_cmp(mrb, &a, &b) < 0) {
+      mpz_set(mrb, &t, &a);
+      mpz_set(mrb, &a, &b);
+      mpz_set(mrb, &b, &t);
+    }
+
+    /* Replace a with (a - b) */
+    mpz_sub(mrb, &a, &a, &b);
+
+    /* Remove factors of 2 from the result if it's even */
+    if (a.sz > 0 && (a.p[0] & 1) == 0) {
+      size_t a_trailing = mpz_trailing_zeros(&a);
+      if (a_trailing > 0) {
+        mpz_div_2exp(mrb, &a, &a, a_trailing);
+      }
+    }
+
+  } while (!zero_p(&a));
+
+  /* Restore common factors of 2 */
+  mpz_mul_2exp(mrb, &b, &b, shift);
+
+  trim(&b);
+  mpz_move(mrb, gg, &b);
+  mpz_clear(mrb, &a);
   mpz_clear(mrb, &t);
 }
 #endif