From 9ef3362f92a1822172e5cca731b8dd074ae9232f Mon Sep 17 00:00:00 2001 From: "Yukihiro \"Matz\" Matsumoto" Date: Wed, 14 Jan 2026 08:19:23 +0900 Subject: [PATCH] mruby-bigint: consolidate mpn layer for low-level limb operations Rename limb_* functions to follow GMP-style mpn_* naming convention: - limb_zero -> mpn_zero - limb_copy -> mpn_copyi - limb_addmul_1 -> mpn_addmul_1 Add new mpn functions for future optimization work: - mpn_submul_1: multiply-and-subtract (r -= s * limb) - mpn_cmp: compare two same-length limb arrays The mpn layer provides direct limb array operations without mpz memory management overhead, which is essential for implementing advanced algorithms like Newton-Raphson division. Co-authored-by: Claude --- mrbgems/mruby-bigint/core/bigint.c | 66 +++++++++++++++++++++--------- 1 file changed, 46 insertions(+), 20 deletions(-) diff --git a/mrbgems/mruby-bigint/core/bigint.c b/mrbgems/mruby-bigint/core/bigint.c index 6073052b7..1411d49f3 100644 --- a/mrbgems/mruby-bigint/core/bigint.c +++ b/mrbgems/mruby-bigint/core/bigint.c @@ -99,7 +99,7 @@ pool_alloc(mpz_pool_t *pool, size_t limbs) /* Zero n limbs at p */ static inline void -limb_zero(mp_limb *p, size_t n) +mpn_zero(mp_limb *p, size_t n) { memset(p, 0, n * sizeof(mp_limb)); } @@ -124,7 +124,7 @@ mpz_init_heap(mpz_ctx_t *ctx, mpz_t *s, size_t hint) mrb_raise(mrb, E_RUNTIME_ERROR, "bigint size too large"); } s->p = (mp_limb*)mrb_malloc(MPZ_MRB(ctx), hint * sizeof(mp_limb)); - limb_zero(s->p, hint); + mpn_zero(s->p, hint); s->sz = hint; } else { @@ -200,7 +200,7 @@ mpz_realloc(mpz_ctx_t *ctx, mpz_t *x, size_t size) #endif /* Zero-initialize new limbs */ - limb_zero(x->p + old_sz, size - old_sz); + mpn_zero(x->p + old_sz, size - old_sz); x->sz = size; } } @@ -638,7 +638,7 @@ mpz_sub_int(mpz_ctx_t *ctx, mpz_t *x, mrb_int n) /* Multiply-and-add: rp[0..n-1] += s1p[0..n-1] * limb; return carry (high limb) */ static inline mp_limb -limb_addmul_1(mp_limb *rp, const mp_limb *s1p, size_t n, mp_limb limb) +mpn_addmul_1(mp_limb *rp, const mp_limb *s1p, size_t n, mp_limb limb) { #if defined(__SIZEOF_INT128__) && (__SIZEOF_INT128__ == 16) /* Use 128-bit arithmetic with 8x unrolling for maximum efficiency */ @@ -743,6 +743,33 @@ limb_addmul_1(mp_limb *rp, const mp_limb *s1p, size_t n, mp_limb limb) #endif } +/* Multiply-and-subtract: rp[0..n-1] -= s1p[0..n-1] * limb; return borrow */ +static inline mp_limb +mpn_submul_1(mp_limb *rp, const mp_limb *s1p, size_t n, mp_limb limb) +{ + mp_dbl_limb borrow = 0; + for (size_t i = 0; i < n; i++) { + mp_dbl_limb prod = (mp_dbl_limb)s1p[i] * (mp_dbl_limb)limb; + mp_dbl_limb sub = (mp_dbl_limb)rp[i] - LOW(prod) - borrow; + rp[i] = LOW(sub); + /* Borrow is 1 if sub underflowed, plus HIGH(prod) */ + borrow = HIGH(prod) + (sub >> (sizeof(mp_dbl_limb) * 8 - 1)); + } + return (mp_limb)borrow; +} + +/* Compare two same-length limb arrays: returns <0, 0, or >0 */ +static inline int +mpn_cmp(const mp_limb *ap, const mp_limb *bp, size_t n) +{ + while (n-- > 0) { + if (ap[n] != bp[n]) { + return (ap[n] > bp[n]) ? 1 : -1; + } + } + return 0; +} + #define KARATSUBA_THRESHOLD 32 static inline mrb_bool @@ -751,7 +778,7 @@ should_use_karatsuba(size_t x_len, size_t y_len) return x_len >= KARATSUBA_THRESHOLD && y_len >= KARATSUBA_THRESHOLD; } -/* w = u * v (optimized schoolbook using limb_addmul_1) */ +/* w = u * v (optimized schoolbook using mpn_addmul_1) */ static void mpz_mul_basic(mpz_ctx_t *ctx, mpz_t *ww, mpz_t *u, mpz_t *v) { @@ -774,9 +801,9 @@ mpz_mul_basic(mpz_ctx_t *ctx, mpz_t *ww, mpz_t *u, mpz_t *v) mp_limb scalar = b->p[0]; mpz_t w; mpz_init_heap(ctx, &w, a->sz + 1); - limb_zero(w.p, a->sz + 1); + mpn_zero(w.p, a->sz + 1); - mp_limb carry = limb_addmul_1(w.p, a->p, a->sz, scalar); + mp_limb carry = mpn_addmul_1(w.p, a->p, a->sz, scalar); w.p[a->sz] = carry; w.sn = a->sn * b->sn; @@ -787,13 +814,13 @@ mpz_mul_basic(mpz_ctx_t *ctx, mpz_t *ww, mpz_t *u, mpz_t *v) mpz_t w; mpz_init_heap(ctx, &w, a->sz + b->sz); - limb_zero(w.p, a->sz + b->sz); + mpn_zero(w.p, a->sz + b->sz); for (size_t j = 0; j < a->sz; j++) { mp_limb a_limb = a->p[j]; if (a_limb == 0) continue; - mp_limb carry = limb_addmul_1(w.p + j, b->p, b->sz, a_limb); + mp_limb carry = mpn_addmul_1(w.p + j, b->p, b->sz, a_limb); /* Properly handle carry propagation to avoid overflow */ size_t k = j + b->sz; @@ -812,9 +839,9 @@ mpz_mul_basic(mpz_ctx_t *ctx, mpz_t *ww, mpz_t *u, mpz_t *v) /* Allocation-free Karatsuba helper functions */ -/* Copy limbs: dest[0..n-1] = src[0..n-1] */ +/* Copy limbs forward: dest[0..n-1] = src[0..n-1] */ static void -limb_copy(mp_limb *dest, const mp_limb *src, size_t n) +mpn_copyi(mp_limb *dest, const mp_limb *src, size_t n) { if (n > 0) { memcpy(dest, src, n * sizeof(mp_limb)); @@ -857,11 +884,11 @@ static void mpz_mul_basic_limbs(mp_limb *result, const mp_limb *x, size_t x_len, const mp_limb *y, size_t y_len) { - limb_zero(result, x_len + y_len); + mpn_zero(result, x_len + y_len); for (size_t i = 0; i < x_len; i++) { if (x[i] == 0) continue; - mp_limb carry = limb_addmul_1(result + i, y, y_len, x[i]); + mp_limb carry = mpn_addmul_1(result + i, y, y_len, x[i]); if (i + y_len < x_len + y_len) { result[i + y_len] += carry; } @@ -886,7 +913,7 @@ static void mpz_sqr_basic_limbs(mp_limb *result, const mp_limb *x, size_t n) { size_t result_len = 2 * n; - limb_zero(result, result_len); + mpn_zero(result, result_len); /* Step 1: Compute off-diagonal terms xi * xj for i < j */ for (size_t i = 0; i < n; i++) { @@ -895,7 +922,7 @@ mpz_sqr_basic_limbs(mp_limb *result, const mp_limb *x, size_t n) if (i + 1 < n) { /* Compute xi * x[i+1..n-1] and add at position 2*i+1 */ - mp_limb carry = limb_addmul_1(result + 2*i + 1, x + i + 1, n - i - 1, xi); + mp_limb carry = mpn_addmul_1(result + 2*i + 1, x + i + 1, n - i - 1, xi); /* Propagate carry */ size_t k = 2*i + 1 + (n - i - 1); @@ -1064,8 +1091,8 @@ mpz_mul_karatsuba(mpz_ctx_t *ctx, mp_limb *result, /* Step 4: Final assembly: result = z0 + z1*B + z2*B^2 */ size_t result_len = x_len + y_len; - limb_zero(result, result_len); - limb_copy(result, z0, z0_len); + mpn_zero(result, result_len); + mpn_copyi(result, z0, z0_len); limb_add_at(result, result_len, z1, z1_len, half); limb_add_at(result, result_len, z2, z2_len, 2 * half); } @@ -1171,8 +1198,8 @@ mpz_sqr_karatsuba(mpz_ctx_t *ctx, mp_limb *result, const mp_limb *x, size_t n, /* Step 4: Final assembly: result = z0 + z1*B + z2*B^2 */ size_t result_len = 2 * n; - limb_zero(result, result_len); - limb_copy(result, z0, z0_len); + mpn_zero(result, result_len); + mpn_copyi(result, z0, z0_len); limb_add_at(result, result_len, z1, z1_len, half); limb_add_at(result, result_len, z2, z2_len, 2 * half); } @@ -1938,7 +1965,6 @@ cleanup: pool_restore(ctx, pool_state); } - static void udiv(mpz_ctx_t *ctx, mpz_t *qq, mpz_t *rr, mpz_t *xx, mpz_t *yy) {