where the first is the original script, while the second is the modified one (the full code is below). However, I ended up still needing to use some low-level functions in the final script; maybe there is a better way to do it with MutableArithmetics? (cc @blegat).
It’s my first attempt at using the package, although I’ve wanted to for awhile, since it looks very cool.
Full code
using MutableArithmetics, BenchmarkTools
@eval MutableArithmetics begin
mutability(::Type{Rational{BigInt}}) = IsMutable()
mutable_copy(x::Rational{BigInt}) = deepcopy(x)
# zero
promote_operation(::typeof(zero), ::Type{Rational{BigInt}}) = Rational{BigInt}
mutable_operate!(::typeof(zero), x::Rational{BigInt}) = (mutable_operate!(zero, x.num); mutable_operate!(one, x.den); x)
# one
promote_operation(::typeof(one), ::Type{Rational{BigInt}}) = Rational{BigInt}
mutable_operate!(::typeof(one), x::Rational{BigInt}) = (mutable_operate!(one, x.num); mutable_operate!(one, x.den); x)
function mutable_operate_to!(output::BigInt, ::typeof(gcd), x::BigInt, y::BigInt)
Base.GMP.MPZ.gcd!(output, x, y)
end
function mutable_operate!(::typeof(gcd), x::BigInt, y::BigInt)
Base.GMP.MPZ.gcd!(x, y)
end
function mutable_operate_to!(output::BigInt, ::typeof(div), x::BigInt, y::BigInt)
Base.GMP.MPZ.tdiv_q!(output, x, y)
end
function mutable_operate!(::typeof(div), x::BigInt, y::BigInt)
Base.GMP.MPZ.tdiv_q!(x, y)
end
# +, -
promote_operation(::Union{typeof(+), typeof(-)}, ::Rational{BigInt}, ::Rational{BigInt}) = Rational{BigInt}
buffer_for(::Union{typeof(+), typeof(-)}, ::Rational{BigInt}, ::Rational{BigInt}) = BigInt()
function mutable_operate_to!(output::Rational{BigInt}, op::Function, x::Rational{BigInt}, y)
return mutable_buffered_operate_to!(buffer_for(op, x, y), output, op, x, y)
end
function mutable_buffered_operate_to!(buffer::BigInt, output::Rational{BigInt}, op::Union{typeof(+), typeof(-)}, x::Rational{BigInt}, y::Rational{BigInt})
# non-mutating algorithm:
# g = gcd(x.den, y.den)
# xd = div(x.den, g)
# yd = div(y.den, g)
# op(x.num * yd, y.num * xd) // (x.den * yd)
mutable_operate_to!(buffer, gcd, x.den, y.den) # buffer = g = gcd(x.den, y.den)
mutable_operate_to!(output.num, div, x.den, buffer) # output.num = xd = div(x.den, g)
mutable_operate_to!(output.den, div, y.den, buffer) # output.den = yd = div(y.den, g)
# `buffer` no longer needs to store `g`
# Now we need to store `xd`
Base.GMP.MPZ.set!(buffer, output.num)
mutable_operate_to!(output.num, *, x.num, output.den) # output.num = x.num * yd
mutable_operate!(*, buffer, y.num) # buffer = xd * y.num
mutable_operate!(op, output.num, buffer) # output.num = op(x.num * yd, xd * y.num)
mutable_operate!(*, output.den, x.den)
# Now `output` is complete, but isn't in lowest terms yet
mutable_operate_to!(buffer, gcd, output.num, output.den)
mutable_operate!(div, output.num, buffer)
mutable_operate!(div, output.den, buffer)
output
end
function mutable_operate!(op::Function, x::Rational{BigInt}, args::Vararg{Any, N}) where N
mutable_operate_to!(x, op, mutable_copy(x), args...)
end
buffer_for(::typeof(*), ::Rational{BigInt}, ::Rational{BigInt}) = (BigInt(), BigInt())
function mutable_buffered_operate_to!(buffer::Tuple{BigInt, BigInt}, output::Rational{BigInt}, ::typeof(*), x::Rational{BigInt}, y::Rational{BigInt})
# non-mutating algorithm:
# g1 = gcd(x.num, y.den)
# xn = div(x.num, g1)
# yd = div(y.den, g1)
# g2 = gcd(x.den, y.num)
# xd = div(x.den, g2)
# yn = div(y.num, g2)
# (xn * yn) // (xd * yd)
g1 = buffer[1]
mutable_operate_to!(g1, gcd, x.num, y.den) # g1 = gcd(x.num, y.den)
mutable_operate_to!(output.num, div, x.num, g1) # output.num = xn = div(x.num, g1)
mutable_operate_to!(output.den, div, y.den, g1) # output.den = yd = div(y.den, g1)
# Now g1 is no longer needed
g2 = buffer[1]
mutable_operate_to!(g2, gcd, x.den, y.num) # g2 = gcd(x.den, y.num)
yn = buffer[2]
mutable_operate_to!(yn, div, y.num, g2) # yn = div(y.num, g2)
mutable_operate!(*, output.num, yn) # output.num is finished; equals xn*yn
xd = buffer[2]
mutable_operate_to!(xd, div, x.den, g2) # xd = div(x.den, g2)
mutable_operate!(*, output.den, xd) # output.den = xd * yd
# Now `output` is complete, but isn't in lowest terms yet
g = buffer[1]
mutable_operate_to!(g, gcd, output.num, output.den)
mutable_operate!(div, output.num, g)
mutable_operate!(div, output.den, g)
output
end
buffer_for(::typeof(*), ::Rational{BigInt}, ::BigInt) = BigInt()
function mutable_buffered_operate_to!(buffer::BigInt, output::Rational{BigInt}, ::typeof(*), x::Rational{BigInt}, y::BigInt)
g = buffer
mutable_operate_to!(output.num, *, x.num, y)
Base.GMP.MPZ.set!(output.den, x.den)
# Now `output` is complete, but isn't in lowest terms yet
mutable_operate_to!(g, gcd, output.num, output.den)
mutable_operate!(div, output.num, g)
mutable_operate!(div, output.den, g)
output
end
end
using Test
@testset begin
a = big(5)//big(2)
b = big(3)
c = big(1)//big(1)
MutableArithmetics.mutable_operate_to!(c, *, a, b)
@test c == a*b
a = big(2) // big(3)
b = big(5) // big(7)
MutableArithmetics.mutable_operate_to!(c, -, a, b)
@test c == a - b
MutableArithmetics.mutable_operate_to!(c, +, a, b)
@test c == a + b
MutableArithmetics.mutable_operate_to!(c, *, a, b)
@test c == a * b
a = big(2)//big(5)
MutableArithmetics.one!(a)
@test a == big(1)//big(1)
MutableArithmetics.zero!(a)
@test a == big(0)//big(1)
end
function bernoulli(n,A)
for m = 0:n
A[m+1] = 1 // (m+1)
for j = m:-1:1
A[j] = j * (A[j] - A[j+1])
end
end
return A[1]
end
function bernoulli2(n,A)
buffer_int = big(1)
buffer_rat = big(1)//big(1)
for m = 0:n
A[m+1] = 1 // (m+1)
for j = m:-1:1
Base.GMP.MPZ.set!(buffer_rat.num, A[j].num)
Base.GMP.MPZ.set!(buffer_rat.den, A[j].den)
MutableArithmetics.mutable_buffered_operate_to!(buffer_int, A[j], -, buffer_rat, A[j+1])
Base.GMP.MPZ.set!(buffer_rat.num, A[j].num)
Base.GMP.MPZ.set!(buffer_rat.den, A[j].den)
MutableArithmetics.mutable_buffered_operate_to!(buffer_int, A[j], *, buffer_rat, big(j))
end
end
return A[1]
end
function bernoulli_display(n, b; do_print = false)
A = Vector{Rational{BigInt}}(undef, n + 1)
B=map(x->b(x,A), 0:n)
for (i, val) in enumerate(B)
(val.num != 0) && do_print && println("B($(i-1))", lpad(val.num,45+Int(i<10)), " / ", val.den)
end
end
@btime bernoulli_display(60, bernoulli)
@btime bernoulli_display(60, bernoulli2)
## set `do_print = true` to confirm they agree
# bernoulli_display(60, bernoulli; do_print = true)
# bernoulli_display(60, bernoulli2; do_print = true)
. Everything ok. It was a mistake of mine with using 
