mruby-bigint: eliminate mpz_sqrt_pool duplication

Removed mpz_sqrt_pool function (203 lines) and its forward declaration
to eliminate code duplication. mpz_sqrt 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 22:33:13 +09:00
parent 7902368d86
commit 715696623e
-204
View File
@@ -37,7 +37,6 @@ typedef struct mpz_pool {
/* Forward declarations */
static int mpz_mul_sliding_window(mrb_state *mrb, mpz_t *result, mpz_t *first, mpz_t *second);
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);
/* Memory allocation tracking for benchmarking */
typedef struct allocation_stats {
@@ -2708,12 +2707,6 @@ mpz_barrett_reduce(mrb_state *mrb, mpz_t *r, mpz_t *x, mpz_t *m, mpz_t *mu)
static void
mpz_sqrt(mrb_state *mrb, mpz_t *z, mpz_t *x)
{
/* Try pool-based square root first for eligible operands */
if (mpz_sqrt_pool(mrb, z, x)) {
return; /* Success with pool-based square root */
}
/* Fallback to traditional algorithm */
mrb_assert(x->sn >= 0);
if (x->sz == 0) {
@@ -2751,203 +2744,6 @@ mpz_sqrt(mrb_state *mrb, mpz_t *z, mpz_t *x)
mpz_clear(mrb, &t);
}
/* Pool-based square root using Newton-Raphson method with stack memory */
static int
mpz_sqrt_pool(mrb_state *mrb, mpz_t *z, mpz_t *x)
{
mrb_assert(x->sn >= 0);
if (x->sz == 0) {
// sqrt(0) = 0
z->sn = 0;
z->sz = 0;
return 1; /* Success - trivial case */
}
/* Only use pools for medium-sized operands that will benefit from stack allocation */
size_t x_limbs = x->sz;
if (x_limbs < 4 || x_limbs > 64) {
return 0; /* Use traditional square root */
}
/* Estimate space needed: s (~x_limbs/2 + 1), t (~x_limbs/2 + 1), quotient (~x_limbs), remainder (~x_limbs) */
size_t estimated_s_size = (x_limbs + 1) / 2 + 2; /* sqrt result size + margin */
size_t temp_space = estimated_s_size * 3 + x_limbs * 2; /* s, t, division temps */
if (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 s, t, quotient, remainder;
int pool_success = 0;
/* Initialize temporary variables in pool */
mpz_init_pool(mrb, &s, pool, estimated_s_size);
mpz_init_pool(mrb, &t, pool, estimated_s_size);
mpz_init_pool(mrb, &quotient, pool, x_limbs);
mpz_init_pool(mrb, &remainder, pool, x_limbs);
/* Verify all pool allocations succeeded */
if (MPZ_POOL_VERIFY_4(s, t, quotient, remainder, pool)) {
/* Estimate initial value: 1 << (bit_length(x) / 2) */
size_t xbits = mpz_bits(x);
size_t sbit = (xbits + 1) / 2;
/* Initialize s = 1 << sbit using pool memory */
if (sbit == 0) {
s.p[0] = 1;
s.sz = 1;
s.sn = 1;
}
else {
size_t limb_shift = sbit / DIG_SIZE;
size_t bit_shift = sbit % DIG_SIZE;
/* Clear s array */
for (size_t i = 0; i < s.sz; i++) {
s.p[i] = 0;
}
if (limb_shift < s.sz) {
s.p[limb_shift] = ((mp_limb)1) << bit_shift;
s.sz = limb_shift + 1;
s.sn = 1;
}
else {
/* Fallback if shift too large */
s.p[0] = 1;
s.sz = 1;
s.sn = 1;
}
}
/* Newton-Raphson iteration: s = (s + x / s) / 2 */
int max_iterations = 100; /* Safety limit */
for (int iter = 0; iter < max_iterations; iter++) {
/* t = x / s using pool-based division */
if (udiv_pool(mrb, &quotient, &remainder, x, &s)) {
/* Pool division succeeded - copy quotient to t */
for (size_t i = 0; i < quotient.sz && i < t.sz; i++) {
t.p[i] = quotient.p[i];
}
t.sz = (quotient.sz < t.sz) ? quotient.sz : t.sz;
t.sn = quotient.sn;
}
else {
/* Pool division failed - fallback to traditional */
pool_success = 0;
break;
}
/* t = t + s using pool-based addition */
mpz_t temp_sum;
mpz_init_pool(mrb, &temp_sum, pool, estimated_s_size + 1);
if (is_pool_memory(&temp_sum, pool)) {
mpz_add(mrb, &temp_sum, &t, &s);
/* Copy sum back to t */
for (size_t i = 0; i < temp_sum.sz && i < t.sz; i++) {
t.p[i] = temp_sum.p[i];
}
t.sz = (temp_sum.sz < t.sz) ? temp_sum.sz : t.sz;
t.sn = temp_sum.sn;
MPZ_POOL_CLEANUP(mrb, temp_sum, pool);
}
else {
/* Pool addition failed - fallback to traditional */
MPZ_POOL_CLEANUP(mrb, temp_sum, pool);
pool_success = 0;
break;
}
/* t = t / 2 (right shift by 1 bit) */
mp_limb carry = 0;
for (size_t i = t.sz; i > 0; i--) {
size_t idx = i - 1;
mp_limb current = t.p[idx];
t.p[idx] = (current >> 1) | carry;
carry = (current & 1) << (DIG_SIZE - 1);
}
/* Trim leading zeros */
while (t.sz > 1 && t.p[t.sz - 1] == 0) {
t.sz--;
}
if (t.sz == 0) {
t.sz = 1;
t.p[0] = 0;
t.sn = 0;
}
/* Check convergence: if t >= s, we're done */
int cmp_result = 0;
if (t.sz > s.sz) {
cmp_result = 1;
}
else if (t.sz < s.sz) {
cmp_result = -1;
}
else {
for (size_t i = t.sz; i > 0; i--) {
size_t idx = i - 1;
if (t.p[idx] > s.p[idx]) {
cmp_result = 1;
break;
}
else if (t.p[idx] < s.p[idx]) {
cmp_result = -1;
break;
}
}
}
if (cmp_result >= 0) {
/* Converged: t >= s */
pool_success = 1;
break;
}
/* s = t for next iteration */
for (size_t i = 0; i < t.sz && i < s.sz; i++) {
s.p[i] = t.p[i];
}
s.sz = (t.sz < s.sz) ? t.sz : s.sz;
s.sn = t.sn;
}
if (pool_success) {
/* Copy final result from pool to heap-allocated output */
mpz_realloc(mrb, z, s.sz);
for (size_t i = 0; i < s.sz; i++) {
z->p[i] = s.p[i];
}
z->sz = s.sz;
z->sn = s.sn;
}
}
/* Pool cleanup is automatic */
MPZ_POOL_CLEANUP(mrb, s, pool);
MPZ_POOL_CLEANUP(mrb, t, pool);
MPZ_POOL_CLEANUP(mrb, quotient, pool);
MPZ_POOL_CLEANUP(mrb, remainder, pool);
pool_storage.active = 0;
if (pool_success) {
return 1; /* Success - used pool memory for square root! */
}
else {
return 0; /* Pool allocation failed, fallback */
}
} while(0);
return 0; /* Should not reach here */
}
/* Barrett reduction for efficient modular arithmetic with repeated operations */