Barrett reduction requires x < 2^(2*bits(m)); the gate condition only
required x->sz >= y->sz + 2, which let x.sz reach 25 limbs against a
4-limb modulus. When the precondition is violated, mpz_barrett_reduce
silently truncates high limbs and returns garbage.
Integer#remainder, which routes through mpz_mod, was affected:
(3**500).remainder((2**100)+3) returned the wrong value. Integer#%
took the udiv path via mpz_mmod and was unaffected.
Add x->sz <= 2 * y->sz to the gate so out-of-range inputs fall
through to the general udiv path.
Co-authored-by: Claude <noreply@anthropic.com>
REDC's input T must satisfy T < R*N; mpz_powm_montgomery() computed
T = base * R^2 without first reducing base modulo n. When base >= n,
T exceeds R*N and REDC silently truncates upper limbs, producing a
wrong base_mont and ultimately a wrong result.
Use mpz_mmod (general division path) to pre-reduce base mod n before
the multiplication by R^2.
Visible effect: (2**160).pow(2, (2**40)+1) returned 0 instead of 1.
A separate pre-existing issue in mpz_mod's Barrett path (broken
precondition check) means that path could return a wrong reduction
for large operands; mmod sidesteps that path entirely.
Co-authored-by: Claude <noreply@anthropic.com>
mpz_t treats sn (sign) as the canonical "is zero" flag (zero_p(x) :=
(x)->sn == 0). When an arithmetic operation produces a value whose
limbs trim to zero size, sn must be reset to 0 to preserve the
invariant. Several call sites already enforced this locally (e.g.
mpz_sub line 616); make trim() responsible so every caller benefits.
Without this, an inconsistent zero bignum (sn!=0, sz=0) can flow into
mpz_sqr, miss the zero_p guard, and reach mpz_init_heap with hint=0
where mpn_zero(NULL, 0) invokes UB (memset() declares its first
argument nonnull). The trip survives at runtime on glibc but is
formally undefined behavior, flagged by UBSan via Integer#pow(b,e,m)
with specific operands.
close#6849
Co-authored-by: Claude <noreply@anthropic.com>
Cover zero operands, power-of-2 fast path, negative operands,
balanced multi-limb pairs with a shared Fibonacci factor, highly
unbalanced pairs (to exercise the Euclidean fallback), and
Fibonacci neighbors (always coprime). Declare a test dependency
on mruby-numeric-ext since Integer#gcd is defined there.
Co-authored-by: Claude <noreply@anthropic.com>
Replace the classical Euclidean main loop (mpz_mod per iteration)
with Stein's binary GCD: subtract + factor out trailing 2s on
odd-maintained operands. Keep an mpz_mod fallback for heavily
unbalanced pairs (the smaller operand has at least two fewer limbs)
where one long division replaces many Stein subtracts.
The old loop allocated a temporary mpz_t every iteration to hold the
mod result; the Stein loop is allocation-free thanks to the in-place
paths in mpz_sub and mpz_div_2exp. This matches the "memory first"
priority and also happens to be faster on typical inputs because
mpz_mod's setup cost dominates when the quotient is small (the
classical Fibonacci-neighbor worst case).
Also add a small mpz_swap helper used by the new loop.
Benchmark (bin/mruby benchmark/bm_bigint_gcd.rb, median of 3):
case before after ratio
single-limb 42 ms 50 ms 1.19 (fast-path
unchanged;
noise)
fib(200) vs fib(201) 82 ms 23 ms 0.28
balanced ~700-bit shared 46 ms 32 ms 0.70
unbalanced big vs small 35 ms 28 ms 0.80
power-of-2 path 38 ms 42 ms 1.11 (fast-path
unchanged;
noise)
balanced ~2800-bit shared 19 ms 15 ms 0.79
Co-authored-by: Claude <noreply@anthropic.com>
The preceding comment claimed "Binary GCD (Stein's algorithm)",
but the multi-limb main loop is classical Euclidean using mpz_mod.
Only the prelude (factoring out common 2s, single-limb and
power-of-2 fast paths) is Stein-flavored. Describe what the code
actually does so readers are not misled.
Co-authored-by: Claude <noreply@anthropic.com>
Per gemini-code-assist review on #6791: replace member-by-member
assignment with C89/C++98-style positional aggregate initialization,
keeping the macro a single declaration and matching the existing
pool_storage initializer style in the same macro.
Co-authored-by: Claude <noreply@anthropic.com>
MPZ_CTX_INIT used a compound literal with designated initializers:
mpz_ctx_t ctx##_struct = ((mpz_ctx_t){.mrb = ..., .pool = ...});
Both features are C99-only, and are not accepted by legacy C++
compilers (notably gcc 4.x) when mruby is pulled into a C++
translation unit via the -cxx.cxx wrapper.
Replace with plain member assignment so the macro expands to code
that is valid under C89/C++98 as well.
Co-authored-by: Claude <noreply@anthropic.com>
mpz_abs() internally allocates via mpz_init_heap(), so
pre-allocating abs_x/abs_y with mpz_init_temp() leaked the
original allocations. let mpz_abs() handle allocation directly.
also use divide-first formula (abs_x/gcd)*abs_y to reduce
intermediate product size, and add bint_norm() for the result.
reported by OSS-Fuzz (clusterfuzz-testcase-6501272051318784).
Co-authored-by: Claude <noreply@anthropic.com>
Barrett and Montgomery reduction compute 2^(2k) internally where k
is the modulus bit length. When this exceeds MRB_BIGINT_BIT_LIMIT,
mrb_raise() via longjmp skips cleanup of allocated temporaries.
Add early modulus size check before any heap allocation.
Co-authored-by: Claude <noreply@anthropic.com>
The work buffer size in mpz_montgomery_reduce() was calculated as
x_len + k + 2, which assumed x_len >= k. However, when R^2 mod n
produces a small result, x_len can be much smaller than k.
The Montgomery reduction loop writes k limbs at work[i] for each
iteration i=0..k-1, so the maximum index accessed is work[2k-1].
This requires at least 2k limbs in the work buffer.
Fixed by ensuring work_size is at least 2*k+2 limbs when x_len < k.
Also initialize b->as.heap before mpz_move in bint_set() to ensure
the destination mpz_t has valid initial state.
Co-authored-by: Claude <noreply@anthropic.com>
mpz_mul_sparse allocated temporary mpz_t variables (shifted, temp) that
were leaked when an exception was raised (e.g., RangeError from shift
width too large). bint_mul had the same issue with its output mpz_t z.
Wrap both functions with MRB_ENSURE to guarantee cleanup runs regardless
of exceptions, following the existing pattern used by mpz_mul_all_ones.
Co-authored-by: Claude <noreply@anthropic.com>
Function was only called once and contained just 2 lines of code.
Inlining directly reduces code size and improves clarity.
Co-authored-by: Claude <noreply@anthropic.com>
Add CPython-style parsing for base-10 string to integer conversion:
- Parse 9 digits at a time into decimal-base array
- Convert decimal-base to binary in single pass
- Use memory pool for temporary decimal buffer
- Use realloc for result buffer to reduce allocations
Also add digit_pairs lookup table for faster to_s output.
Performance: 2-5x faster for to_i, 60% fewer allocations.
Co-authored-by: Claude <noreply@anthropic.com>
This improves to_s performance for medium-sized bigints (40-50 limbs,
~800-1000 digits) by approximately 5x by enabling the divide-and-conquer
algorithm earlier.
Benchmark results:
40 limbs (772 digits): 88 us -> 18 us (5x faster)
50 limbs (964 digits): 134 us -> 25 us (5.4x faster)
Co-authored-by: Claude <noreply@anthropic.com>
Benchmarks show the previous threshold of 50 was too low, causing
Toom-3's setup overhead to outweigh its asymptotic benefits for
medium-sized numbers. Raising to 100 limbs provides:
- 2x faster at 300 limbs (192 -> 96 us)
- 2.5x faster at 120 limbs (42 -> 17 us)
- 2.6x faster at 80 limbs (31 -> 12 us)
Co-authored-by: Claude <noreply@anthropic.com>
When multiplying numbers where one is significantly larger than the other
(at least 2x size difference), split the larger number into chunks matching
the smaller number's size, multiply each chunk, and combine results. This
avoids pathological performance when Toom-3 pads asymmetric operands with
zeros.
Benchmarks show 6-16x speedup for size ratios from 10:1 to 40:1, with no
regression for symmetric cases.
Co-authored-by: Claude <noreply@anthropic.com>
mpz_init_heap() already allocates the requested size, so immediately
calling mpz_realloc() with the same size is a no-op. Remove these
redundant calls from mpz_and, mpz_or, mpz_xor, mpz_mod_2exp, and
mpz_abs.
Co-authored-by: Claude <noreply@anthropic.com>
When the output and input are the same variable, avoid unnecessary
heap allocations by modifying in place:
- mpz_neg: just flip the sign
- mpz_abs: just make sign positive
- ulshift: use mpn_lshift in-place (safe since it processes high-to-low)
Co-authored-by: Claude <noreply@anthropic.com>
Add mpz_sqr_toom3() that performs Toom-3 squaring with reduced memory
and computation:
- Only evaluates x (not y), reducing evaluation buffer from 6 to 3
- Uses recursive squaring instead of multiplication for all 5 products
- Simplifies interpolation since squared values are always positive
The specialized squaring uses the same Toom-3 structure but avoids
redundant computation when both operands are the same number.
Co-authored-by: Claude <noreply@anthropic.com>
The mpn_divexact_3 function used ap[i] in the borrow computation after
writing to rp[i]. When rp == ap (in-place operation), this read the
modified value instead of the original input, causing incorrect borrow
propagation.
This bug caused Toom-3 multiplication to produce wrong results for
certain input patterns where t6 - t5 had non-zero values followed by
zeros. The corrupted r3 coefficient then propagated errors to the final
result.
Fix by saving the original ap[i] value before writing rp[i].
Co-authored-by: Claude <noreply@anthropic.com>
move in-place optimization directly into mpz_sub() so callers
just use mpz_sub(ctx, x, x, y) and get automatic optimization.
remove separate mpz_sub_inplace() function.
Co-authored-by: Claude <noreply@anthropic.com>
when destination equals source in mpz_div_2exp(), use memmove
and mpn_rshift in-place instead of allocating a temporary.
reduces sqrt allocations by 49% since Newton iteration uses
in-place division by 2 on each iteration.
Co-authored-by: Claude <noreply@anthropic.com>
add usub_inplace() and mpz_sub_inplace() for allocation-free
subtraction when the minuend is larger than the subtrahend.
apply to Mersenne multiplication which reduces allocations by
17% and improves performance by 7-9% for small numbers.
Co-authored-by: Claude <noreply@anthropic.com>
skip two's complement conversion in mpz_and, mpz_or, mpz_xor when both
operands are positive. this avoids the per-limb make_2comp overhead and
provides up to 1.6x speedup for large bigints.
Co-authored-by: Claude <noreply@anthropic.com>
convert mpz_init_heap to mpz_init_temp for temporary quotient and
remainder variables in div_limb. these variables are now allocated
from the memory pool when possible, reducing heap allocation overhead.
the div_limb function already uses pool_save/pool_restore, so these
temporary variables are proper candidates for pool allocation.
Co-authored-by: Claude <noreply@anthropic.com>
Replace Karatsuba multiplication (O(n^1.585)) with Toom-3 (O(n^1.465))
for large number multiplication. Toom-3 splits numbers into thirds and
evaluates at 5 points, providing better asymptotic performance.
Threshold is 50 limbs (~1600 bits). For operands below threshold or
highly asymmetric sizes, schoolbook multiplication is used.
Co-authored-by: Claude <noreply@anthropic.com>
Remove automatic downgrade to 16-bit limbs on 32-bit Windows.
Modern compilers (including MSVC) have supported uint64_t for decades.
MRB_NO_MPZ64BIT remains available for constrained platforms.
Adjust BATCH_DIVISOR and BATCH_DIGITS for 16-bit limb compatibility:
- 32-bit limbs: 10^9 (9 digits per batch)
- 16-bit limbs: 10^4 (4 digits per batch)
Co-authored-by: Claude <noreply@anthropic.com>
Add trim() after mpz_set in early return paths to prevent propagation
of inflated sz values. When an mpz_t has sz larger than actual allocated
limbs, copying it without trim causes subsequent operations to read
beyond allocated memory.
Fixed functions:
- mpz_add: when one operand is zero
- mpz_neg: when copying operand
- mpz_mod_2exp: when x < 2^e
Co-authored-by: Claude <noreply@anthropic.com>
Use mrb_protect_error API instead of direct MRB_TRY/MRB_CATCH to handle
exceptions in mpz_mul_all_ones and mpz_to_s_dc. This maintains C++
compatibility (issue #6702) while ensuring temporary mpz_t allocations
are properly freed even when exceptions occur.
Co-authored-by: Claude <noreply@anthropic.com>
Remove MRB_TRY/MRB_CATCH exception handling from bigint.c to fix
compilation errors when using mruby-bigint in C++ projects with
MRB_USE_CXX_EXCEPTION enabled.
The exception handling was added for cleanup on error, but it requires
throw.h which doesn't work when a C file is compiled in a C++ context
with C++ exceptions enabled. Accepting potential memory leaks on
exception (rare) is preferable to breaking C++ builds.
Fixes#6702
Co-authored-by: Claude <noreply@anthropic.com>
Replace the two-buffer swap pattern in D&C to_s base case with
in-place division using new mpn_div10_9 function. This eliminates
the q_base scratch buffer and reduces per-iteration overhead.
Compilers optimize the constant division by 10^9 to multiplication
and shift operations for better performance.
Co-authored-by: Claude <noreply@anthropic.com>
Replace per-call lo allocation with depth-indexed lo_stack buffers
that are reused across recursion levels. Each buffer is allocated
on first use at that depth with appropriate size.
This reduces 493 malloc/free calls (5%) and 580KB of memory (2%)
for large number to_s conversions while maintaining performance.
Co-authored-by: Claude <noreply@anthropic.com>
Instead of allocating a separate hi buffer for the upper part of the
split, reuse q5 by shifting it in place after extracting the lower
bits to q5_low.
This eliminates one mpz_t allocation per recursive call:
- Extract q5_low = q5 mod 2^k first (copy lower bits)
- Shift q5 right in place (memmove + mpn_rshift) to get hi
- q5 now serves as hi for the recursive call
Benchmark results (2.4M bit number):
- Memory: -3.5% (5.58MB -> 5.38MB peak heap)
- Instructions: -16.5%
- Speed: unchanged (within measurement noise)
Co-authored-by: Claude <noreply@anthropic.com>
Replace array indexing x.p[i+j] with pointer arithmetic *xp++ in the
hot inner loop of Knuth Algorithm D division. This avoids recalculating
the index i+j on every iteration.
Profiling showed the inner loop accounts for ~80% of udiv execution time,
with the array indexing contributing significant overhead.
Benchmark improvement: ~4% faster (7.80s -> 7.48s for 2.4M bit to_s).
Co-authored-by: Claude <noreply@anthropic.com>
When converting large numbers to strings using D&C algorithm, the
base case extracts digits in batches of 9 (for 32-bit limbs). The
extraction logic: 1 leading digit + 4 pairs (8 digits) = 9 digits.
The pair extraction loop condition `pos >= 2` exits when pos < 2,
but when the remaining batch still has value and pos == 1, that
final digit was being lost and replaced with '0' by the padding loop.
This caused roundtrip failures (x.to_s.to_i != x) for numbers just
above the D&C threshold (1000 digits), where the split boundary
produced a lo part requiring exactly the right number of digits to
trigger this edge case.
Co-authored-by: Claude <noreply@anthropic.com>
Yukihiro "Matz" Matsumoto