From 8e91554c6d5102c9ac84bf0051011b6773ffb35d Mon Sep 17 00:00:00 2001 From: "Yukihiro \"Matz\" Matsumoto" Date: Thu, 23 Apr 2026 22:52:05 +0900 Subject: [PATCH] mruby-bigint: rewrite mpz_gcd main loop as binary Stein algorithm Replace the classical Euclidean main loop (mpz_mod per iteration) with Stein's binary GCD: subtract + factor out trailing 2s on odd-maintained operands. Keep an mpz_mod fallback for heavily unbalanced pairs (the smaller operand has at least two fewer limbs) where one long division replaces many Stein subtracts. The old loop allocated a temporary mpz_t every iteration to hold the mod result; the Stein loop is allocation-free thanks to the in-place paths in mpz_sub and mpz_div_2exp. This matches the "memory first" priority and also happens to be faster on typical inputs because mpz_mod's setup cost dominates when the quotient is small (the classical Fibonacci-neighbor worst case). Also add a small mpz_swap helper used by the new loop. Benchmark (bin/mruby benchmark/bm_bigint_gcd.rb, median of 3): case before after ratio single-limb 42 ms 50 ms 1.19 (fast-path unchanged; noise) fib(200) vs fib(201) 82 ms 23 ms 0.28 balanced ~700-bit shared 46 ms 32 ms 0.70 unbalanced big vs small 35 ms 28 ms 0.80 power-of-2 path 38 ms 42 ms 1.11 (fast-path unchanged; noise) balanced ~2800-bit shared 19 ms 15 ms 0.79 Co-authored-by: Claude --- mrbgems/mruby-bigint/core/bigint.c | 45 +++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 10 deletions(-) diff --git a/mrbgems/mruby-bigint/core/bigint.c b/mrbgems/mruby-bigint/core/bigint.c index 668611dd9..15c2031b1 100644 --- a/mrbgems/mruby-bigint/core/bigint.c +++ b/mrbgems/mruby-bigint/core/bigint.c @@ -348,6 +348,14 @@ mpz_move(mpz_ctx_t *ctx, mpz_t *y, mpz_t *x) x->sz = 0; } +static inline void +mpz_swap(mpz_t *a, mpz_t *b) +{ + mpz_t tmp = *a; + *a = *b; + *b = tmp; +} + static size_t digits(mpz_t *x) { @@ -4821,9 +4829,11 @@ mpz_power_of_2_p(mpz_t *x) return (limb != 0) && ((limb & (limb - 1)) == 0); } -/* Hybrid GCD: a binary (Stein) prelude factors out common powers of 2 - and handles single-limb / power-of-2 fast paths; the multi-limb main - loop is classical Euclidean using mpz_mod. */ +/* Binary GCD (Stein's algorithm): factor out common powers of 2, + then iterate on odd operands with subtract + trailing-zero shift. + For heavily unbalanced pairs (one operand has at least two more + limbs than the other) a single Euclidean step via mpz_mod replaces + many Stein subtracts. */ static void mpz_gcd(mpz_ctx_t *ctx, mpz_t *gg, mpz_t *aa, mpz_t *bb) { @@ -4895,14 +4905,29 @@ mpz_gcd(mpz_ctx_t *ctx, mpz_t *gg, mpz_t *aa, mpz_t *bb) mpz_div_2exp(ctx, &a, &a, a_zeros); mpz_div_2exp(ctx, &b, &b, b_zeros); - /* Euclidean algorithm for multi-limb numbers */ + /* Stein main loop. Invariant: a and b are positive and odd. + Euclidean fallback when b has >=2 more limbs than a. */ while (!zero_p(&b)) { - mpz_t temp; - mpz_init_temp(ctx, &temp, a.sz); - mpz_mod(ctx, &temp, &a, &b); - mpz_move(ctx, &a, &b); - mpz_move(ctx, &b, &temp); - mpz_clear(ctx, &temp); + if (mpz_cmp(ctx, &a, &b) > 0) { + mpz_swap(&a, &b); + } + if (b.sz >= a.sz + 2) { + mpz_t temp; + mpz_init_temp(ctx, &temp, a.sz); + mpz_mod(ctx, &temp, &b, &a); + mpz_move(ctx, &b, &temp); + mpz_clear(ctx, &temp); + if (zero_p(&b)) break; + size_t bz = mpz_trailing_zeros(&b); + if (bz > 0) + mpz_div_2exp(ctx, &b, &b, bz); + } + else { + mpz_sub(ctx, &b, &b, &a); + if (zero_p(&b)) break; + size_t bz = mpz_trailing_zeros(&b); + mpz_div_2exp(ctx, &b, &b, bz); + } } mpz_mul_2exp(ctx, gg, &a, shift); mpz_clear(ctx, &a);