{"id":638,"date":"2019-11-09T14:47:28","date_gmt":"2019-11-09T14:47:28","guid":{"rendered":"http:\/\/perfectionatic.org\/?p=638"},"modified":"2020-06-09T23:49:23","modified_gmt":"2020-06-09T23:49:23","slug":"evaluating-bigint-expressions-in-julia-using-metaprogramming","status":"publish","type":"post","link":"https:\/\/perfectionatic.org\/?p=638","title":{"rendered":"Evaluating Arbitrary Precision Integer Expressions in Julia using Metaprogramming"},"content":{"rendered":"\n<p>While watching the <a href=\"https:\/\/www.youtube.com\/channel\/UC1_uAIS3r8Vu6JjXWvastJg\">Mathologer<\/a>  masterclass  on power sums <\/p>\n\n\n\n<figure class=\"wp-block-embed-youtube wp-block-embed is-type-video is-provider-youtube wp-embed-aspect-16-9 wp-has-aspect-ratio\"><div class=\"wp-block-embed__wrapper\">\n<iframe loading=\"lazy\" title=\"Power sum MASTER CLASS: How to sum quadrillions of powers ... by hand! (Euler-Maclaurin formula)\" width=\"740\" height=\"416\" src=\"https:\/\/www.youtube.com\/embed\/fw1kRz83Fj0?feature=oembed\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share\" referrerpolicy=\"strict-origin-when-cross-origin\" allowfullscreen><\/iframe>\n<\/div><\/figure>\n\n\n\n<p>I came across a challenge to evaluate the following sum  <br>\\(1^{10}+2^{10}+\\cdots+1000^{10}\\)<\/p>\n\n\n\n<p>This can be easily evaluated via brute force in Julia to yield <\/p>\n\n\n\n<pre lang=\"julia\">function power_cal(n)\n  a=big(0)\n  for i=1:n\n    a+=big(i)^10\n  end\n  a \nend\n\njulia> power_cal(1000)\n91409924241424243424241924242500\n<\/pre>\n\n\n\n<p>Note the I had to use <code>big<\/code> to makes sure we are using <code>BigInt<\/code> in the computation. Without that, we would be quickly running in into an overflow issue and we will get a very wrong number. <\/p>\n\n\n\n<p>In the <a href=\"https:\/\/www.youtube.com\/watch?v=fw1kRz83Fj0&amp;lc=Ugy4AM5t06yx8hQF1md4AaABAg\">comment<\/a> section of the video I found a very elegant solution to the above sum, expressed as <\/p>\n\n\n\n<p>(1\/11) * 1000^11 + (1\/2) * 1000^10 + (5\/6) * 1000^9 &#8211; 1000^7 + 1000^5-(1\/2) * 1000^3 + (5\/66) * 1000 = 91409924241424243424241924242500<\/p>\n\n\n\n<p> If I try to plug this into the Julia, I get <\/p>\n\n\n\n<pre lang=\"julia\">julia> (1\/11) * 1000^11 + (1\/2) * 1000^10 + (5\/6) * 1000^9 - 1000^7 + 1000^5-(1\/2) * 1000^3 + (5\/66) * 1000\n-6.740310541071357e18\n<\/pre>\n\n\n\n<p>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 <code>Int64<\/code> with <code>BigInt<\/code> by wrapping it in the <code>big<\/code> function. But that would be cumbersome to do by hand.<\/p>\n\n\n\n<h2 class=\"wp-block-heading\">The power of Metaprogramming <\/h2>\n\n\n\n<p>In Julia, <a href=\"https:\/\/docs.julialang.org\/en\/v1\/manual\/metaprogramming\/\">metaprogramming<\/a> allows you to write <a href=\"https:\/\/nbviewer.jupyter.org\/github\/dpsanders\/intermediate_julia_2019\/blob\/master\/4.%20Metaprogramming.ipynb\">code that creates code<\/a>, the idea here to manipulate the abstract syntax tree (AST) of a Julia expression. We start to by &#8220;quoting&#8221; our original mathematical expressing into a Julia expression. In the at form it is not evaluated yet, however we can always evaluate it via <code>eval<\/code>. <\/p>\n\n\n\n<pre lang=\"julia\">\njulia> ex1=:((1\/11) * 1000^11 + (1\/2) * 1000^10 + (5\/6) * 1000^9 - 1000^7 + 1000^5-(1\/2) * 1000^3 + (5\/66) * 1000)\n\n:((((((1 \/ 11) * 1000 ^ 11 + (1 \/ 2) * 1000 ^ 10 + (5 \/ 6) * 1000 ^ 9) - 1000 ^ 7) + 1000 ^ 5) - (1 \/ 2) * 1000 ^ 3) + (5 \/ 66) * 1000)\n\njulia> dump(ex1)\nExpr\n  head: Symbol call\n  args: Array{Any}((3,))\n    1: Symbol +\n    2: Expr\n      head: Symbol call\n      args: Array{Any}((3,))\n        1: Symbol -\n        2: Expr\n          head: Symbol call\n          args: Array{Any}((3,))\n            1: Symbol +\n            2: Expr\n              head: Symbol call\n              args: Array{Any}((3,))\n                1: Symbol -\n                2: Expr\n                3: Expr\n            3: Expr\n              head: Symbol call\n              args: Array{Any}((3,))\n                1: Symbol ^\n                2: Int64 1000\n                3: Int64 5\n     \n\njulia> eval(ex1)\n-6.740310541071357e18\n<\/pre>\n\n\n\n<p>The output of <code>dump<\/code> show the full AST in all its glory (&#8230;well almost the depth is a bit truncated). Notice that here all our numbers are interpreted as <code>Int64<\/code>. <br>Now we walk through the AST and change all occurrences of <code>Int64<\/code> with <code>BigInt<\/code> by using the <code>big<\/code> function.<\/p>\n\n\n\n<pre lang=\"julia\">function makeIntBig!(ex::Expr)\n  args=ex.args\n  for i in eachindex(args)\n       if args[i] isa Int64\n          args[i]=big(args[i])\n       end\n       if args[i] isa Expr\n          makeIntBig!(args[i])\n       end\n   end\nend\n\njulia> ex2=copy(ex1)\n:((((((1 \/ 11) * 1000 ^ 11 + (1 \/ 2) * 1000 ^ 10 + (5 \/ 6) * 1000 ^ 9) - 1000 ^ 7) + 1000 ^ 5) - (1 \/ 2) * 1000 ^ 3) + (5 \/ 66) * 1000)\n\njulia> makeIntBig!(ex2)\n\njulia> eval(ex2)\n9.14099242414242434242419242425000000000000000000000000000000000000000000000014e+31\n<\/pre>\n\n\n\n<p>We see an improvement here, but the results are not very satisfactory. The divisions yield <code>BigFloat<\/code> results, which had a tiny bit of floating point errors. Can we do better? <br>Julia has support for <code><a href=\"https:\/\/docs.julialang.org\/en\/v1\/manual\/complex-and-rational-numbers\/#Rational-Numbers-1\">Rational<\/a><\/code> expressions baked in. We can use that improve the results. We just need to search for <code>call<\/code> expressions the <code>\/<\/code> symbol and replace it by the <code>\/\/<\/code> symbol. For safety we just have to makes sure the operands are as subtype of <code>Integer<\/code>. <\/p>\n\n\n\n<pre lang=\"julia\">\nfunction makeIntBig!(ex::Expr)\n    args=ex.args\n    if ex.head == :call && args[1]==:\/ && \n              length(args)==3 && \n              all(x->typeof(x) <: Integer,args[2:end]) \n        args[1]=:\/\/\n        args[2]=big(args[2])\n        args[3]=big(args[3])\n    else \n        for i in eachindex(args)\n           if args[i] isa Int64\n              args[i]=big(args[i])\n           end\n           if args[i] isa Expr\n              makeIntBig!(args[i])\n           end\n        end\n    end\nend\n\njulia> ex2=copy(ex1);\n\njulia> makeIntBig!(ex2)\n\njulia> eval(ex2)\n91409924241424243424241924242500\/\/1\n<\/pre>\n\n\n\n<p>Now that is much better! We have not lost any precision and we ended us with a <code>Rational<\/code> expression. <\/p>\n\n\n\n<p>Finally, we can build a  <code>macro<\/code> so the if we run into such expressions in the future and we want to evaluate them, we could just conveniently call it. <\/p>\n\n\n\n<pre lang=\"julia\"> \nmacro eval_bigInt(ex)\n   makeIntBig!(ex)\n   quote # Removing the denominator if it is redundant \n      local x=$ex\n      (x isa Rational && x.den==1) ? x.num : x \n   end\nend\n<\/pre>\n\n\n\n<p>and we can now simply evaluate our original expression as <\/p>\n\n\n\n<pre lang=\"julia\"> \njulia> @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\n91409924241424243424241924242500\n<\/pre>\n\n\n\n<p><\/p>\n","protected":false},"excerpt":{"rendered":"<p>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 [&hellip;]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[73,62,71],"tags":[92],"class_list":["post-638","post","type-post","status-publish","format-standard","hentry","category-julialang","category-math","category-programming","tag-julia"],"_links":{"self":[{"href":"https:\/\/perfectionatic.org\/index.php?rest_route=\/wp\/v2\/posts\/638","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/perfectionatic.org\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/perfectionatic.org\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/perfectionatic.org\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/perfectionatic.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=638"}],"version-history":[{"count":28,"href":"https:\/\/perfectionatic.org\/index.php?rest_route=\/wp\/v2\/posts\/638\/revisions"}],"predecessor-version":[{"id":694,"href":"https:\/\/perfectionatic.org\/index.php?rest_route=\/wp\/v2\/posts\/638\/revisions\/694"}],"wp:attachment":[{"href":"https:\/\/perfectionatic.org\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=638"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/perfectionatic.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=638"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/perfectionatic.org\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=638"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}