diff --git a/mrbgems/mruby-bigint/core/bigint.c b/mrbgems/mruby-bigint/core/bigint.c index 0e5397752..246af4c06 100644 --- a/mrbgems/mruby-bigint/core/bigint.c +++ b/mrbgems/mruby-bigint/core/bigint.c @@ -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)