8 ms·
20M digits of pi in 1 minute using Julia
- adgjlsfhk1 4y agoThis is pretty cool. One minor gripe is that there were a lot of useless type asserts cluttering things up. every usage of ::BigFloat here wasn't needed.
- juliusgeo 4y agoThis is absolutely true—I was just putting them everywhere possible in an attempt to avoid allocations.
- cmovq 4y agoI don’t know much about Julia, why would it potentially cause allocations?
- adgjlsfhk1 4y agoif the code was type unstable that can cause allocations, but in cases like this it wouldn't. Type instability is the ability of the compiler to predict the types of variables from the types of the inputs. Type unstable code can be found using `@code_warntype` and will cause allocations do to boxing and dynamic dispatch.
- juliusgeo 4y agoThank you for the feedback. I will try to remove as many of the explicit type hints as possible. The main reason I added them everywhere is because I found sometimes there would be instability at super low decimal ranges because of type promotion
- xmcqdpt2 4y agoYeah this is a minor gripe with Julia. It has type inference but it's risky to use.
- andyferris 4y agoTrue, but you can remove a lot of the risk by inspecting the results of type inference with `@code_warntype` in places where you suspect it might be not be working. I believe JET.jl can automate global program analysis, but I haven't used it.
- queuebert 4y agoI've found Julia to be an excellent language for Project Euler [1]. Besides the speed, you can use Unicode identifiers, so the solution can closer follow the math. 1. https://projecteuler.net/ https://projecteuler.net/
- recursive 4y agoThat doesn't seem to be a particularly distinguishing feature. I know that C#, python, and javascript can all use non-ascii identifiers. I believe most widely used languages can also, but I'm less certain of the details.
- martinsmit 4y agoThe difference is that Julia makes it easy to insert these characters (\symbol) so programming with them becomes more natural.
- jb1991 4y agoThere’s also Swift, of course. https://docs.swift.org/swift-book/documentation/the-swift-programming-language/stringsandcharacters#Unicode https://docs.swift.org/swift-book/documentation/the-swift-pr...
- adgjlsfhk1 4y agoI think the distinction is that Julia's parser does a relatively good job assigning code-points to unary operators/binary operators/identifiers as appropriate, combined with the fact that multiple dispatch inherently looks more like math.
- npalli 4y agoJulia will read 2π (fixed:thx kps) as 2 times pi which follows the mathematical notation. In all the other languages you have mentioned, I believe even if you can represent pi in Unicode (π) you will have to write 2*π. This can become very irritating when you get to for ex. python/numpy where you will litter your code with * and ** instead of the conventional mathematical notation. Gets a bit tiring if you are looking at more than a page of code. Aesthetics matter.
- goldenkey 4y agoFor reference: https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrat...
- IncRnd 4y agoHere is a video [1] and associated java source code [2] from pi day in 2018 to generate 1M (or more) digits of pi using the same Gauss-Legendre algorithm. [1] https://www.youtube.com/watch?v=8RONJPOgZjw https://www.youtube.com/watch?v=8RONJPOgZjw [2] https://github.com/R-e-t-u-r-n-N-u-l-l/Java-Algorithms/blob/master/ApproxPi.java https://github.com/R-e-t-u-r-n-N-u-l-l/Java-Algorithms/blob/...
- juliusgeo 4y agoMy algorithm only takes about 1.5 seconds for 1M digits. These are very cool resources though, I will read them.
- IncRnd 4y agoThanks! I haven't run any timings on the link I posted. I suspect julia would be much faster than java, though! I just thought it was neat that your post in March reminded me of an earlier pi day post on the same algorithm. I like the explosion from 1M to 20M.
- juliusgeo 4y agoI would love to see the timings for Java runtimes. The various implementations of BigFloats seem to be very customized to the language, so maybe you could do it faster!
- deleted 4y ago[deleted]
- adgjlsfhk1 4y agoOn my computer the Julia algorithm gives .98 seconds for 1M digits and 40.7 seconds for 20M (zen 2 cpu).
- juliusgeo 4y agoThanks for the data!
- tiffanyh 4y agoIs 1-minute fast?
- juliusgeo 4y agoIn the article I compare my implementation to the MPFR pi implementation which is in pure C. So the fact I’m able to get close to that performance in Julia is impressive(I think).
- Joel_Mckay 4y agoAll variables held the same, Julia often beats C, NumPy, and most other languages. Julia also supports transparent Cluster API and CUDA integration. i.e. the hot-mess you hope to never have to maintain professionally. In my humble opinion, the utility of this high-level language makes it fundamentally different. Thus, it is slowly becoming more mainstream as the ecosystem stabilizes. Remember to have fun =)
- jiggawatts 4y agoI just ran the same thing in Mathematica and I got 20 million digits in just over 7 seconds. You can have "very high level functional programming" and 8x the speed of Julia in a mathematics-oriented language!
- Joel_Mckay 4y agoMathematica does several optimizations before execution, and as such can exploit a priori knowledge of many problems. The danger of course is when Wolfram takes those assumptions into a room of mathematicians. It is controversial for all the wrong reasons. =) It is a common question by the way =) https://stackoverflow.com/questions/60121757/julia-vs-mathematica-numerical-integration-performance https://stackoverflow.com/questions/60121757/julia-vs-mathem...
- Max_Limelihood 4y agoJulia already is a very high-level mathematics-oriented programming language, though. The reason Mathematica is so much faster here is it’s using a different algorithm. When you compare using the same algorithm, Julia is 10-100x faster than Mathematica. https://julialang.org/benchmarks/ https://julialang.org/benchmarks/
- boredemployee 4y agogenuine question: would that amount of digits be useful for any application?
- fisherjeff 4y agohttps://www.jpl.nasa.gov/edu/news/2016/3/16/how-many-decimals-of-pi-do-we-really-need/ https://www.jpl.nasa.gov/edu/news/2016/3/16/how-many-decimal...
- juliusgeo 4y agoThank you for posting this! This is exactly what I was looking for.
- yesenadam 4y agoTLDR: "For JPL's highest accuracy calculations, which are for interplanetary navigation, we use 3.141592653589793"...Using pi to 15 decimal places, you'd get the distance of Voyager 1 wrong by less than 1cm!
- juliusgeo 4y agoNo, it is basically only useful as a programming exercise. You only need like 40 digits to compute the circumference of the universe to the precision of a hydrogen atom.
- edgyquant 4y agoWhat is the circumference of the universe?
- mk_stjames 4y ago2.8×10^39 picometers is circumference of the observable part of our universe, roughly. A sole hydrogen atom has a bohr radius of about 25 picometers. So 40 digits of pi would indeed be in the ballpark for computing anything distance-wise at the circumference of the observable universe down to the hydrogen electron cloud size.
- joseftexas 4y agoIsn't this is more hardware dependent? A 2012 Intel notebook vs a 2022 AMD laptop would yield different results with same code. Also a well tune assembly can do these way faster. A python wrapped with assembly functions and run on a machine 6ghz Intel machine with ramdisk would also outperform this. I always feel this kind of claim is not helpful in day to day usage.
- juliusgeo 4y agoThat is why I compared it to the stdlib implementation.
- Temporary_31337 4y agoJust wanted to mention that SuperPi - the benchmark that is very widely used for at least a decade for measuring PC / CPU single threaded performance uses the same algorithm so not really obscure. That said, multi threaded algorithms make more sense for modern multi core computer architectures, something like this for example: http://www.numberworld.org/y-cruncher/ http://www.numberworld.org/y-cruncher/
- userbinator 4y agoI still remember when the SuperPI 1M <10s record was beaten... with a CPU overclocked to 5GHz+: https://tweakers.net/ext/i/1160952453.gif https://tweakers.net/ext/i/1160952453.gif The world record remains at just under 4s, with 8GHz+ oerclocks: https://hwbot.org/benchmark/superpi_-_1m/halloffame https://hwbot.org/benchmark/superpi_-_1m/halloffame Clock cycle efficiency doesn't seem to have improved much in the last decade.
- xmreg 4y ago[flagged]
- Someone 4y agoNitpick: is that really computing 20M digits correctly? Reading the code gives me the impression that this does all float computations in 20 million digits (or 20 million-ish? If BigFloat internally is binary, it might use a bit more), but I would expect that you need to compute using a few more digits to avoid rounding errors in the last digits. Edit: also, reading https://stackoverflow.com/a/67919533 https://stackoverflow.com/a/67919533, it seems you can do BigFloat(π, precision=20000000) How does that perform?
- SherlockSage 4y agoI ran the following in my Julia REPL julia> @time BigFloat(π, precision=20000000); 11.208211 seconds (60.53 k allocations: 2.932 GiB, 0.16% gc time) Which is pretty good considering I ran this in WSL on my laptop and Mathematica in a different comment took 10 seconds. (Plain Windows took ~26 seconds for some reason?)
- rnestler 4y ago> (Plain Windows took ~26 seconds for some reason?) Maybe the allocations were the reason?
- Someone 4y agoI’m too late to update my own comment, but I now suspect BigFloat(π, precision=20000000) may allocate a Bigfloat with precision of two million digits and then store the constant π which has much lower precision in it. What value does that print?
- juliusgeo 4y agoIt definitely does calculate it otherwise it would diverge significantly at lower digits from the Gauss Legendre and Chudnovsky algorithm that I compared it to. It just defers calculation to MPFR, the C library Julia binds to for BigFloat calculations.
- juliusgeo 4y agoAlso one additional comment is that you are setting precision in number of bits, not in decimal digits. It obviously runs much faster when it is only computing log2(20000000) not 20000000 digits of precision.
- jiggawatts 4y agoIn Wolfram Mathematica on my desktop PC with a 7 year old processor: Timing[N[\[Pi], 20000000];] {10.1094, Null} On a fairly recent laptop CPU I got 7.6 seconds. A fair ways to go still!
- xdavidliu 4y agodoes Mathematica actually compute it or is it just loading the value that it already knows?
- extasia 4y ago7.6 second seems too much to be just loading a value from disk
- jiggawatts 4y agoDefinitely computing. For arbitrary precision evaluation like this, Mathematica has always had cutting-edge algorithms built in. As a random example: large matrix operations have been multi-threaded for ages.
- deleted 4y ago[deleted]
- moelf 4y agoit's not using the same algorithm right?
- juliusgeo 4y agoWhile this is impressive, it is neither the same algorithm, nor even the same type of language (I already know Julia is much slower than C). I don’t think it is much of a relevant comparison for a benchmark when almost nothing is the same except for the expected answer. Perhaps if Gauss Legendre were implemented in Mathematica it would only take 3 seconds! There’s no way for us to know.
- bakul 4y agoThe following implementation of the Chudnovsky algorithm takes about the same time in Gambit-Scheme on M1 Mac for 20M digits (and allocates 16GB): (define (pi digits) (let* ((A 13591409) (B 545140134) (C 640320) (C^3/24 (quotient (expt 640320 3) 24)) (D 12)) (define (-1^n*g n g) (if (odd? n) (- g) g)) (define (split m n) (if (= 1 (- n m)) (let* ((6n (* 6 n)) (g (* (- 6n 5) (- (+ n n) 1) (- 6n 1)))) (list g (* C^3/24 (expt n 3)) (* (-1^n*g n g) (+ (* n B) A)))) (let* ((mid (quotient (+ m n) 2)) (gpq1 (split m mid)) (gpq2 (split mid n)) (g1 (car gpq1)) (p1 (cadr gpq1)) (q1 (caddr gpq1)) (g2 (car gpq2)) (p2 (cadr gpq2)) (q2 (caddr gpq2))) (list (* g1 g2) (* p1 p2) (+ (* q1 p2) (* q2 g1)))))) (let* ((num-terms (inexact->exact (floor (+ 2 (/ digits 14.181647462))))) (sqrt-C (integer-sqrt (* C (expt 100 digits)))) (gpq (split 0 num-terms)) (g (car gpq)) (p (cadr gpq)) (q (caddr gpq))) (quotient (* p C sqrt-C) (* D (+ q (* p A))))))) I am surprised it is much slower in Julia (as per what is noted in the gist).
- juliusgeo 4y agoInteresting! What does Gauss Legendre look like in Gambit-Scheme? I also wonder if perhaps there is a way to get the Julia implementation up to par…
- bakul 4y agoI will give it a try, time permitting. The scheme code above is a variation of something I converted from CommonLisp (orignal code by Robert Smith) in 2013. Its speed depends on bignum multiplication speed. I imagine most serious bignum implementations use the Schönhage–Strassen algorithm for very large numbers.
- juliusgeo 4y agoThanks for looking into this! That is interesting, I wonder if it varies a great deal based on the implementation of the BigFloat rather than the language that it is implemented in.