# 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.