From 51da16df99bfe975b82e64c6e870c64022096b01 Mon Sep 17 00:00:00 2001 From: Lilith Orion Hafner Date: Mon, 1 Sep 2025 15:27:31 -0500 Subject: [PATCH 1/5] New `prod(::AbstractArray{BigInt})` implementation --- base/gmp.jl | 49 +++++++++++++++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 18 deletions(-) diff --git a/base/gmp.jl b/base/gmp.jl index e4d9294766aaa..7d3db6ead812c 100644 --- a/base/gmp.jl +++ b/base/gmp.jl @@ -605,10 +605,14 @@ Number of ones in the binary representation of abs(x). """ count_ones_abs(x::BigInt) = iszero(x) ? 0 : MPZ.mpn_popcount(x) +"Assumes `x != 0`. Computes top_set_bit(abs(x))." +function _top_set_bit(x::BigInt) + x.size * BITS_PER_LIMB - GC.@preserve x leading_zeros(unsafe_load(x.d, abs(x.size))) +end function top_set_bit(x::BigInt) isnegative(x) && throw(DomainError(x, "top_set_bit only supports negative arguments when they have type BitSigned.")) iszero(x) && return 0 - x.size * sizeof(Limb) << 3 - leading_zeros(GC.@preserve x unsafe_load(x.d, x.size)) + _top_set_bit(x) end divrem(x::BigInt, y::BigInt, ::typeof(RoundToZero) = RoundToZero) = MPZ.tdiv_qr(x, y) @@ -680,24 +684,33 @@ sum(arr::Union{AbstractArray{BigInt}, Tuple{BigInt, Vararg{BigInt}}}) = foldl(MPZ.add!, arr; init=BigInt(0)) function prod(arr::AbstractArray{BigInt}) - # compute first the needed number of bits for the result, - # to avoid re-allocations; - # GMP will always request n+m limbs for the result in MPZ.mul!, - # if the arguments have n and m limbs; so we add all the bits - # taken by the array elements, and add BITS_PER_LIMB to that, - # to account for the rounding to limbs in MPZ.mul! - # (BITS_PER_LIMB-1 would typically be enough, to which we add - # 1 for the initial multiplication by init=1 in foldl) - nbits = BITS_PER_LIMB - for x in arr - iszero(x) && return zero(BigInt) - xsize = abs(x.size) - lz = GC.@preserve x leading_zeros(unsafe_load(x.d, xsize)) - nbits += xsize * BITS_PER_LIMB - lz + any(iszero, arr) && return zero(BigInt) + _prod(arr, firstindex(arr), lastindex(arr)) +end +function _prod(arr::AbstractArray{BigInt}, lo, hi) + if hi - lo + 1 <= 16 + # compute first the needed number of bits for the result, + # to avoid re-allocations; + # GMP will always request n+m limbs for the result in MPZ.mul!, + # if the arguments have n and m limbs; so we add all the bits + # taken by the array elements, and add BITS_PER_LIMB to that, + # to account for the rounding to limbs in MPZ.mul! + # (BITS_PER_LIMB-1 would typically be enough, to which we add + # 1 for the initial multiplication by init) + nbits = BITS_PER_LIMB + for i in lo:hi + nbits += _top_set_bit(arr[i]) + end + init = BigInt(; nbits) + MPZ.set_si!(init, 1) + for i in lo:hi + MPZ.mul!(init, arr[i]) + end + init + else + mid = (lo + hi) รท 2 + MPZ.mul!(_prod(arr, lo, mid), _prod(arr, mid+1, hi)) end - init = BigInt(; nbits) - MPZ.set_si!(init, 1) - foldl(MPZ.mul!, arr; init) end factorial(n::BigInt) = !isnegative(n) ? MPZ.fac_ui(n) : throw(DomainError(n, "`n` must not be negative.")) From 8a2387606af212bf347a64acf138c6358316ffac Mon Sep 17 00:00:00 2001 From: Lilith Orion Hafner Date: Mon, 1 Sep 2025 19:00:39 -0500 Subject: [PATCH 2/5] bugfix --- base/gmp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/gmp.jl b/base/gmp.jl index 7d3db6ead812c..af982ca1257d4 100644 --- a/base/gmp.jl +++ b/base/gmp.jl @@ -607,7 +607,7 @@ count_ones_abs(x::BigInt) = iszero(x) ? 0 : MPZ.mpn_popcount(x) "Assumes `x != 0`. Computes top_set_bit(abs(x))." function _top_set_bit(x::BigInt) - x.size * BITS_PER_LIMB - GC.@preserve x leading_zeros(unsafe_load(x.d, abs(x.size))) + abs(x.size) * BITS_PER_LIMB - GC.@preserve x leading_zeros(unsafe_load(x.d, abs(x.size))) end function top_set_bit(x::BigInt) isnegative(x) && throw(DomainError(x, "top_set_bit only supports negative arguments when they have type BitSigned.")) From a5511b3eccce6642ad03dd05b9566b85d6fad565 Mon Sep 17 00:00:00 2001 From: Lilith Orion Hafner Date: Tue, 2 Sep 2025 11:42:12 -0500 Subject: [PATCH 3/5] Simplify by using a looser approximation for size --- base/gmp.jl | 20 +++++--------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/base/gmp.jl b/base/gmp.jl index af982ca1257d4..7aa5a3764ad92 100644 --- a/base/gmp.jl +++ b/base/gmp.jl @@ -605,14 +605,10 @@ Number of ones in the binary representation of abs(x). """ count_ones_abs(x::BigInt) = iszero(x) ? 0 : MPZ.mpn_popcount(x) -"Assumes `x != 0`. Computes top_set_bit(abs(x))." -function _top_set_bit(x::BigInt) - abs(x.size) * BITS_PER_LIMB - GC.@preserve x leading_zeros(unsafe_load(x.d, abs(x.size))) -end function top_set_bit(x::BigInt) isnegative(x) && throw(DomainError(x, "top_set_bit only supports negative arguments when they have type BitSigned.")) iszero(x) && return 0 - _top_set_bit(x) + x.size * BITS_PER_LIMB - GC.@preserve x leading_zeros(unsafe_load(x.d, abs(x.size))) end divrem(x::BigInt, y::BigInt, ::typeof(RoundToZero) = RoundToZero) = MPZ.tdiv_qr(x, y) @@ -690,18 +686,12 @@ end function _prod(arr::AbstractArray{BigInt}, lo, hi) if hi - lo + 1 <= 16 # compute first the needed number of bits for the result, - # to avoid re-allocations; - # GMP will always request n+m limbs for the result in MPZ.mul!, - # if the arguments have n and m limbs; so we add all the bits - # taken by the array elements, and add BITS_PER_LIMB to that, - # to account for the rounding to limbs in MPZ.mul! - # (BITS_PER_LIMB-1 would typically be enough, to which we add - # 1 for the initial multiplication by init) - nbits = BITS_PER_LIMB + # to avoid re-allocations + nlimbs = 0 for i in lo:hi - nbits += _top_set_bit(arr[i]) + nlimbs += arr[i].size end - init = BigInt(; nbits) + init = BigInt(; nlimbs*BITS_PER_LIMB) MPZ.set_si!(init, 1) for i in lo:hi MPZ.mul!(init, arr[i]) From 1e27175a05d85bbb3e0df606b409af6cd16d8de1 Mon Sep 17 00:00:00 2001 From: Lilith Orion Hafner Date: Wed, 3 Sep 2025 12:29:57 -0500 Subject: [PATCH 4/5] Fix typo Co-authored-by: Alex Arslan --- base/gmp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/gmp.jl b/base/gmp.jl index 7aa5a3764ad92..ed82f300f3fcb 100644 --- a/base/gmp.jl +++ b/base/gmp.jl @@ -691,7 +691,7 @@ function _prod(arr::AbstractArray{BigInt}, lo, hi) for i in lo:hi nlimbs += arr[i].size end - init = BigInt(; nlimbs*BITS_PER_LIMB) + init = BigInt(; nbits=nlimbs*BITS_PER_LIMB) MPZ.set_si!(init, 1) for i in lo:hi MPZ.mul!(init, arr[i]) From bec1ac147ccb85069ef6aa9ee6119a2f30cd9f01 Mon Sep 17 00:00:00 2001 From: Lilith Orion Hafner Date: Sun, 7 Sep 2025 12:20:41 -0500 Subject: [PATCH 5/5] Back out unrelated changes --- base/gmp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/gmp.jl b/base/gmp.jl index ed82f300f3fcb..b2cac7d1c4760 100644 --- a/base/gmp.jl +++ b/base/gmp.jl @@ -608,7 +608,7 @@ count_ones_abs(x::BigInt) = iszero(x) ? 0 : MPZ.mpn_popcount(x) function top_set_bit(x::BigInt) isnegative(x) && throw(DomainError(x, "top_set_bit only supports negative arguments when they have type BitSigned.")) iszero(x) && return 0 - x.size * BITS_PER_LIMB - GC.@preserve x leading_zeros(unsafe_load(x.d, abs(x.size))) + x.size * sizeof(Limb) << 3 - leading_zeros(GC.@preserve x unsafe_load(x.d, x.size)) end divrem(x::BigInt, y::BigInt, ::typeof(RoundToZero) = RoundToZero) = MPZ.tdiv_qr(x, y)