While watching the Mathologer masterclass on power sums

I came across a challenge to evaluate the following sum

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 head: Symbol call args: Array{Any}((3,)) 1: Symbol + 2: Expr head: Symbol call args: Array{Any}((3,)) 1: Symbol - 2: Expr head: Symbol call args: Array{Any}((3,)) 1: Symbol + 2: Expr head: Symbol call args: Array{Any}((3,)) 1: Symbol - 2: Expr 3: Expr 3: Expr head: Symbol call 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