mruby-bigint: eliminate mpz_gcd_pool duplication

Removed mpz_gcd_pool function (299 lines) and its forward declaration
to eliminate code duplication. mpz_gcd now uses heap allocation only.
Pool support should be restored in future using unified approach.

Co-authored-by: Claude <noreply@anthropic.com>
This commit is contained in:
Yukihiro "Matz" Matsumoto
2025-07-25 20:07:01 +09:00
parent 3d1db70052
commit 69541a156d
-306
View File
@@ -39,7 +39,6 @@ static int mpz_mul_sliding_window(mrb_state *mrb, mpz_t *result, mpz_t *first, m
static int udiv_pool(mrb_state *mrb, mpz_t *qq, mpz_t *rr, mpz_t *xx, mpz_t *yy);
static int mpz_sqrt_pool(mrb_state *mrb, mpz_t *z, mpz_t *x);
static int mpz_powm_pool(mrb_state *mrb, mpz_t *zz, mpz_t *x, mpz_t *ex, mpz_t *n);
static int mpz_gcd_pool(mrb_state *mrb, mpz_t *gg, mpz_t *aa, mpz_t *bb);
/* Memory allocation tracking for benchmarking */
typedef struct allocation_stats {
@@ -2670,12 +2669,6 @@ mpz_power_of_2_p(mpz_t *x)
static void
mpz_gcd(mrb_state *mrb, mpz_t *gg, mpz_t *aa, mpz_t *bb)
{
/* Try pool-based GCD first for eligible operands */
if (mpz_gcd_pool(mrb, gg, aa, bb)) {
return; /* Success with pool-based GCD */
}
/* Fallback to traditional algorithm */
mpz_t a, b;
/* Handle special cases */
@@ -3066,305 +3059,6 @@ mpz_sub_pool(mrb_state *mrb, mpz_t *result, mpz_t *a, mpz_t *b, mpz_pool_t *pool
}
/* Pool-based GCD using binary GCD algorithm with Lehmer acceleration */
static int
mpz_gcd_pool(mrb_state *mrb, mpz_t *gg, mpz_t *aa, mpz_t *bb)
{
/* Handle special cases first - no pool needed */
if (zero_p(aa)) {
mpz_abs(mrb, gg, bb);
return 1; /* Success - trivial case */
}
if (zero_p(bb)) {
mpz_abs(mrb, gg, aa);
return 1; /* Success - trivial case */
}
/* Fast path for single-limb numbers - no pool needed */
if (aa->sz <= 1 && bb->sz <= 1) {
mp_limb a_limb = (aa->sz == 0) ? 0 : aa->p[0];
mp_limb b_limb = (bb->sz == 0) ? 0 : bb->p[0];
mp_limb result = limb_gcd(a_limb, b_limb);
mpz_init(mrb, gg);
if (result == 0) {
gg->sn = 0;
gg->sz = 0;
}
else {
mpz_realloc(mrb, gg, 1);
gg->p[0] = result;
gg->sn = 1;
}
return 1; /* Success */
}
/* Only use pools for medium-sized operands that will benefit from stack allocation */
size_t max_limbs = (aa->sz > bb->sz) ? aa->sz : bb->sz;
if (max_limbs < 4 || max_limbs > 32) {
return 0; /* Use traditional GCD */
}
/* Fast paths for powers of 2 - delegate to traditional algorithm for clean memory management */
if (mpz_power_of_2_p(aa) || mpz_power_of_2_p(bb)) {
return 0; /* Use traditional GCD to avoid mixing pool/heap memory */
}
/* Estimate space needed for all GCD temporaries */
size_t estimated_temp_size = max_limbs + 2; /* Working copy size estimation */
size_t lehmer_temp_space = estimated_temp_size * 6; /* 6 temporaries for Lehmer */
size_t total_temp_space = estimated_temp_size * 2 + lehmer_temp_space; /* a, b + Lehmer temps */
if (total_temp_space > BIGINT_POOL_DEFAULT_SIZE / 2) {
return 0; /* Pool too small for all temporaries */
}
do {
mpz_pool_t pool_storage = {0};
pool_storage.capacity = BIGINT_POOL_DEFAULT_SIZE;
pool_storage.active = 1;
mpz_pool_t *pool = &pool_storage;
mpz_t a, b;
int pool_success = 0;
/* Initialize working copies in pool */
mpz_init_pool(mrb, &a, pool, estimated_temp_size);
mpz_init_pool(mrb, &b, pool, estimated_temp_size);
/* Verify pool allocation succeeded */
if (!is_pool_memory(&a, pool) || !is_pool_memory(&b, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
/* Copy absolute values to working variables using pool-based operations */
if (mpz_abs_copy(mrb, &a, aa) && mpz_abs_copy(mrb, &b, bb)) {
/* Find power of 2 that divides both a and b */
size_t a_zeros = mpz_trailing_zeros(&a);
size_t b_zeros = mpz_trailing_zeros(&b);
size_t shift = (a_zeros < b_zeros) ? a_zeros : b_zeros;
/* Remove common factors of 2 using manual bit shifting */
if (shift > 0) {
if (!mpz_div_2exp_pool(mrb, &a, &a, shift, pool) ||
!mpz_div_2exp_pool(mrb, &b, &b, shift, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
}
/* Remove remaining factors of 2 from a */
if (a_zeros > shift) {
if (!mpz_div_2exp_pool(mrb, &a, &a, a_zeros - shift, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
}
/* Remove remaining factors of 2 from b */
if (b_zeros > shift) {
if (!mpz_div_2exp_pool(mrb, &b, &b, b_zeros - shift, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
}
/* Use Lehmer's algorithm for large multi-limb numbers (> 3 limbs) */
if (a.sz > 3 && b.sz > 3) {
/* Extract the two most significant limbs for approximation */
mp_limb a_high = a.p[a.sz - 1];
mp_limb a_low = a.p[a.sz - 2];
mp_limb b_high = b.p[b.sz - 1];
mp_limb b_low = b.p[b.sz - 2];
/* Perform Lehmer reduction on double-precision approximations */
mp_limb u0 = 1, u1 = 0, v0 = 0, v1 = 1;
while (b_high > 0) {
/* Calculate quotient using double-precision approximation */
mp_limb q;
if (a_high == b_high) {
q = (a_low >= b_low) ? 1 : 0;
}
else {
/* Approximate quotient from most significant limbs */
q = a_high / (b_high + 1);
}
if (q == 0) break;
/* Check if applying this quotient would cause overflow */
mp_limb max_limb = (mp_limb)(-1);
if (u1 > 0 && q > max_limb / u1) break;
if (v1 > 0 && q > max_limb / v1) break;
/* Update transformation matrix */
mp_limb t;
t = u0 - q * u1; u0 = u1; u1 = t;
t = v0 - q * v1; v0 = v1; v1 = t;
t = a_high - q * b_high; a_high = b_high; b_high = t;
/* Stop if coefficients get too large */
if (u1 == 0 && v1 == 0) break;
}
/* Apply the transformation if it's non-trivial using pool memory */
if (u1 != 0 || v1 != 0) {
mpz_t temp_a, temp_b, u0_a, v0_b, u1_a, v1_b;
/* Initialize all Lehmer temporaries in pool */
mpz_init_pool(mrb, &temp_a, pool, estimated_temp_size);
mpz_init_pool(mrb, &temp_b, pool, estimated_temp_size);
mpz_init_pool(mrb, &u0_a, pool, estimated_temp_size);
mpz_init_pool(mrb, &v0_b, pool, estimated_temp_size);
mpz_init_pool(mrb, &u1_a, pool, estimated_temp_size);
mpz_init_pool(mrb, &v1_b, pool, estimated_temp_size);
/* Verify all Lehmer pool allocations succeeded */
int lehmer_pool_ok = MPZ_POOL_VERIFY_6(temp_a, temp_b, u0_a, v0_b, u1_a, v1_b, pool);
if (lehmer_pool_ok) {
/* Set initial values using pool-based copy operations */
if (mpz_set_pool(mrb, &u0_a, &a, pool) && mpz_set_pool(mrb, &v0_b, &b, pool) &&
mpz_set_pool(mrb, &u1_a, &a, pool) && mpz_set_pool(mrb, &v1_b, &b, pool)) {
/* Compute u0*a, v0*b, u1*a, v1*b using pool-based multiplication */
if (mpz_mul_int_pool(mrb, &u0_a, u0, pool) && mpz_mul_int_pool(mrb, &v0_b, v0, pool) &&
mpz_mul_int_pool(mrb, &u1_a, u1, pool) && mpz_mul_int_pool(mrb, &v1_b, v1, pool)) {
/* temp_a = u0*a + v0*b using addition */
mpz_add(mrb, &temp_a, &u0_a, &v0_b);
mpz_add(mrb, &temp_b, &u1_a, &v1_b);
if (1) {
/* Update a and b using pool-based set operations */
if (mpz_set_pool(mrb, &a, &temp_a, pool) && mpz_set_pool(mrb, &b, &temp_b, pool)) {
/* Lehmer transformation successful */
}
else {
pool_success = 0;
}
}
else {
pool_success = 0;
}
}
else {
pool_success = 0;
}
}
else {
pool_success = 0;
}
}
else {
/* Lehmer pool allocation failed */
pool_success = 0;
}
/* Cleanup Lehmer temporaries - ALWAYS clean up regardless of success/failure */
MPZ_POOL_CLEANUP(mrb, temp_a, pool);
MPZ_POOL_CLEANUP(mrb, temp_b, pool);
MPZ_POOL_CLEANUP(mrb, u0_a, pool);
MPZ_POOL_CLEANUP(mrb, v0_b, pool);
MPZ_POOL_CLEANUP(mrb, u1_a, pool);
MPZ_POOL_CLEANUP(mrb, v1_b, pool);
if (pool_success == 0) {
goto cleanup_gcd;
}
/* Ensure a >= b after transformation */
if (mpz_cmp(mrb, &a, &b) < 0) {
/* In-place swap - just swap the mpz_t structures */
mpz_t temp_holder = a;
a = b;
b = temp_holder;
}
}
}
/* Main binary GCD loop using pool-based operations */
do {
/* Make b odd efficiently using pool-based division */
if (b.sz > 0 && (b.p[0] & 1) == 0) {
size_t b_trailing = mpz_trailing_zeros(&b);
if (b_trailing > 0) {
if (!mpz_div_2exp_pool(mrb, &b, &b, b_trailing, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
}
}
/* Now both a and b are odd. Ensure a >= b */
if (mpz_cmp(mrb, &a, &b) < 0) {
/* In-place swap without temporary variable */
mpz_t temp_holder = a;
a = b;
b = temp_holder;
}
/* Replace a with (a - b) using pool-based subtraction */
if (!mpz_sub_pool(mrb, &a, &a, &b, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
/* 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) {
if (!mpz_div_2exp_pool(mrb, &a, &a, a_trailing, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
}
}
} while (!zero_p(&a));
/* Restore common factors of 2 using pool-based multiplication */
if (shift > 0) {
if (!mpz_mul_2exp_pool(mrb, &b, &b, shift, pool)) {
pool_success = 0;
goto cleanup_gcd;
}
}
/* Copy final result from pool to heap-allocated output */
trim(&b);
mpz_realloc(mrb, gg, b.sz);
for (size_t i = 0; i < b.sz; i++) {
gg->p[i] = b.p[i];
}
gg->sz = b.sz;
gg->sn = b.sn;
pool_success = 1;
}
else {
/* Pool absolute value operations failed */
pool_success = 0;
}
cleanup_gcd:
/* Pool cleanup is automatic */
MPZ_POOL_CLEANUP(mrb, a, pool);
MPZ_POOL_CLEANUP(mrb, b, pool);
pool_storage.active = 0;
if (pool_success == 1) {
return 1; /* Success - used pool memory for GCD! */
}
else {
return 0; /* Pool allocation failed, fallback */
}
} while(0);
return 0; /* Should not reach here */
}
static size_t
mpz_bits(const mpz_t *x)