From 512fffdac82280f99307d8bf97cbf99e24126615 Mon Sep 17 00:00:00 2001 From: "Yukihiro \"Matz\" Matsumoto" Date: Fri, 9 Jan 2026 22:09:51 +0900 Subject: [PATCH] mruby-bigint: fix division bug with non-standard qhat refinement The udiv function had two buggy modifications to Knuth's Algorithm D: 1. A "three-limb pre-adjustment" that only decremented qhat once 2. A "3-limb refinement" loop with incorrect carry handling These caused incorrect quotients for certain decimal divisions like 10^52 / 10^26. Restored standard Knuth Algorithm D which uses only 2-limb qhat refinement with correction via subtract and add-back. Co-authored-by: Claude --- mrbgems/mruby-bigint/core/bigint.c | 52 ++---------------------------- 1 file changed, 3 insertions(+), 49 deletions(-) diff --git a/mrbgems/mruby-bigint/core/bigint.c b/mrbgems/mruby-bigint/core/bigint.c index ae665ae24..c9cacf9ea 100644 --- a/mrbgems/mruby-bigint/core/bigint.c +++ b/mrbgems/mruby-bigint/core/bigint.c @@ -1965,60 +1965,14 @@ udiv(mpz_ctx_t *ctx, mpz_t *qq, mpz_t *rr, mpz_t *xx, mpz_t *yy) rhat = dividend_val % z; } else { - /* Two limbs available - use enhanced estimation */ + /* Two limbs available - standard Knuth estimation */ mp_dbl_limb dividend_val = ((mp_dbl_limb)x.p[j+yd] << DIG_SIZE) + x.p[j+yd-1]; qhat = dividend_val / z; rhat = dividend_val % z; - - /* Three-limb pre-adjustment when available */ - if (yd >= 2 && j+yd-2 < x.sz && y.p[yd-2] != 0) { - mp_dbl_limb y_second = y.p[yd-2]; - mp_dbl_limb x_third = x.p[j+yd-2]; - - if (qhat > 0) { - mp_dbl_limb left = qhat * y_second; - mp_dbl_limb right = (rhat << DIG_SIZE) + x_third; - - if (qhat >= ((mp_dbl_limb)1 << DIG_SIZE) || left > right) { - qhat--; - rhat += z; - } - } - } } - /* Enhanced qhat refinement step */ - if (yd > 2) { // Now considering at least 3 limbs of divisor - mp_dbl_limb y_second = y.p[yd-2]; - mp_dbl_limb y_third = y.p[yd-3]; // New: third limb of divisor - mp_dbl_limb x_third = (j+yd-2 < x.sz) ? x.p[j+yd-2] : 0; - mp_dbl_limb x_fourth = (j+yd-3 < x.sz) ? x.p[j+yd-3] : 0; // New: fourth limb of dividend - - // Initial check with 2 limbs - mp_dbl_limb left_side = qhat * y_second; - mp_dbl_limb right_side = (rhat << DIG_SIZE) + x_third; - - while (qhat >= ((mp_dbl_limb)1 << DIG_SIZE) || (left_side > right_side)) { - qhat--; - rhat += z; - if (rhat >= ((mp_dbl_limb)1 << DIG_SIZE)) break; - left_side -= y_second; - right_side = (rhat << DIG_SIZE) + x_third; - } - - // Additional check with 3 limbs (new refinement) - left_side = qhat * y_third; - right_side = (rhat << DIG_SIZE) + x_fourth; - - while (qhat >= ((mp_dbl_limb)1 << DIG_SIZE) || (left_side > right_side)) { - qhat--; - rhat += z; - if (rhat >= ((mp_dbl_limb)1 << DIG_SIZE)) break; - left_side -= y_third; - right_side = (rhat << DIG_SIZE) + x_fourth; - } - } - else if (yd == 2) { // Original 2-limb check + /* Standard Knuth Algorithm D qhat refinement (2-limb check) */ + if (yd >= 2) { mp_dbl_limb y_second = y.p[yd-2]; mp_dbl_limb x_third = (j+yd-2 < x.sz) ? x.p[j+yd-2] : 0; mp_dbl_limb left_side = qhat * y_second;