admin/mruby-mruby
1
0
Fork 0
mirror of https://github.com/mruby/mruby synced 2026-06-08 16:11:16 +00:00
Code Issues Packages Projects Releases Wiki Activity

mruby-bigint: fix division bug with non-standard qhat refinement

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
Yukihiro "Matz" Matsumoto
2026-01-09 22:09:51 +09:00
parent 1f590521b0
commit 512fffdac8
1 changed files with 3 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mrbgems/mruby-bigint/core/bigint.c
+3 -49
View File
@@ -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;