# Evaluating Arbitrary Precision Integer Expressions in Julia using Metaprogramming

While watching the Mathologer masterclass on power sums

I came across a challenge to evaluate the following sum
$$1^{10}+2^{10}+\cdots+1000^{10}$$

This can be easily evaluated via brute force in Julia to yield

function power_cal(n)
a=big(0)
for i=1:n
a+=big(i)^10
end
a
end

julia> power_cal(1000)
91409924241424243424241924242500

Note the I had to use big to makes sure we are using BigInt in the computation. Without that, we would be quickly running in into an overflow issue and we will get a very wrong number.

In the comment section of the video I found a very elegant solution to the above sum, expressed as

(1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 – 1000^7 + 1000^5-(1/2) * 1000^3 + (5/66) * 1000 = 91409924241424243424241924242500

If I try to plug this into the Julia, I get

julia> (1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 - 1000^7 + 1000^5-(1/2) * 1000^3 + (5/66) * 1000
-6.740310541071357e18

This negative answer is not surprising at all, because we obviously ran into an overflow. We can, of course, go through that expression and modify all instances of Int64 with BigInt by wrapping it in the big function. But that would be cumbersome to do by hand.

## The power of Metaprogramming

In Julia, metaprogramming allows you to write code that creates code, the idea here to manipulate the abstract syntax tree (AST) of a Julia expression. We start to by “quoting” our original mathematical expressing into a Julia expression. In the at form it is not evaluated yet, however we can always evaluate it via eval.

julia> ex1=:((1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 - 1000^7 + 1000^5-(1/2) * 1000^3 + (5/66) * 1000)

:((((((1 / 11) * 1000 ^ 11 + (1 / 2) * 1000 ^ 10 + (5 / 6) * 1000 ^ 9) - 1000 ^ 7) + 1000 ^ 5) - (1 / 2) * 1000 ^ 3) + (5 / 66) * 1000)

julia> dump(ex1)
Expr
args: Array{Any}((3,))
1: Symbol +
2: Expr
args: Array{Any}((3,))
1: Symbol -
2: Expr
args: Array{Any}((3,))
1: Symbol +
2: Expr
args: Array{Any}((3,))
1: Symbol -
2: Expr
3: Expr
3: Expr
args: Array{Any}((3,))
1: Symbol ^
2: Int64 1000
3: Int64 5

julia> eval(ex1)
-6.740310541071357e18

The output of dump show the full AST in all its glory (…well almost the depth is a bit truncated). Notice that here all our numbers are interpreted as Int64.
Now we walk through the AST and change all occurrences of Int64 with BigInt by using the big function.

function makeIntBig!(ex::Expr)
args=ex.args
for i in eachindex(args)
if args[i] isa Int64
args[i]=big(args[i])
end
if args[i] isa Expr
makeIntBig!(args[i])
end
end
end

julia> ex2=copy(ex1)
:((((((1 / 11) * 1000 ^ 11 + (1 / 2) * 1000 ^ 10 + (5 / 6) * 1000 ^ 9) - 1000 ^ 7) + 1000 ^ 5) - (1 / 2) * 1000 ^ 3) + (5 / 66) * 1000)

julia> makeIntBig!(ex2)

julia> eval(ex2)
9.14099242414242434242419242425000000000000000000000000000000000000000000000014e+31

We see an improvement here, but the results are not very satisfactory. The divisions yield BigFloat results, which had a tiny bit of floating point errors. Can we do better?
Julia has support for Rational expressions baked in. We can use that improve the results. We just need to search for call expressions the / symbol and replace it by the // symbol. For safety we just have to makes sure the operands are as subtype of Integer.

function makeIntBig!(ex::Expr)
args=ex.args
if ex.head == :call && args[1]==:/ &&
length(args)==3 &&
all(x->typeof(x) <: Integer,args[2:end])
args[1]=://
args[2]=big(args[2])
args[3]=big(args[3])
else
for i in eachindex(args)
if args[i] isa Int64
args[i]=big(args[i])
end
if args[i] isa Expr
makeIntBig!(args[i])
end
end
end
end

julia> ex2=copy(ex1);

julia> makeIntBig!(ex2)

julia> eval(ex2)
91409924241424243424241924242500//1

Now that is much better! We have not lost any precision and we ended us with a Rational expression.

Finally, we can build a macro so the if we run into such expressions in the future and we want to evaluate them, we could just conveniently call it.

macro eval_bigInt(ex)
makeIntBig!(ex)
quote # Removing the denominator if it is redundant
local x=\$ex
(x isa Rational && x.den==1) ? x.num : x
end
end

and we can now simply evaluate our original expression as

julia> @eval_bigInt (1/11) * 1000^11 + (1/2) * 1000^10 + (5/6) * 1000^9 - 1000^7 + 1000^5-(1/2) * 1000^3 + (5/66) * 1000
91409924241424243424241924242500

This site uses Akismet to reduce spam. Learn how your comment data is processed.