20M digits of pi in 1 minute using Julia

https://gist.github.com/juliusgeo/41811563811a6e523086e514ef2bec4a

by juliusgeo • 3 years ago

109 77 3 years ago

gist.github.com

adgjlsfhk17• 3 years ago

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

juliusgeo18• 3 years ago

This is absolutely true—I was just putting them everywhere possible in an attempt to avoid allocations.

cmovq• 3 years ago

I don’t know much about Julia, why would it potentially cause allocations?

xmcqdpt2• 3 years ago

Yeah this is a minor gripe with Julia. It has type inference but it's risky to use.

queuebert• 3 years ago

I'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/

recursive• 3 years ago

That 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• 3 years ago

The difference is that Julia makes it easy to insert these characters (\symbol) so programming with them becomes more natural.

adgjlsfhk17• 3 years ago

I 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• 3 years ago

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

moelf2• 3 years ago

the Julia REPL allows you to type unicode with \blah, also if you use `?` help mode, you can give it a unicode sequence and it will tell you how to type it, so it's nice quality of life

    help?> Ψ₂ˣ
    "Ψ₂ˣ" can be typed by \Psi<tab>\_2<tab>\^x<tab>

andyferris2• 3 years ago

While true, one distinguishing features is the sheer number of unicode binary and unary operators that can be used and overloaded. Take for example ∪ for set union and ∩ for set intersection (eg `set1 ∪ set2`).

(You can also type these easily via latex-style completions, e.g. with \cup<TAB> for ∪ or \cap<TAB> for ∩)

goldenkey• 3 years ago

For reference: https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrat...

IncRnd2• 3 years ago

Here 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

[2] https://github.com/R-e-t-u-r-n-N-u-l-l/Java-Algorithms/blob/...

juliusgeo18• 3 years ago

My algorithm only takes about 1.5 seconds for 1M digits. These are very cool resources though, I will read them.

IncRnd2• 3 years ago

Thanks! 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.

adgjlsfhk17• 3 years ago

On my computer the Julia algorithm gives .98 seconds for 1M digits and 40.7 seconds for 20M (zen 2 cpu).

juliusgeo18• 3 years ago

Thanks for the data!

tiffanyh• 3 years ago

Is 1-minute fast?

juliusgeo18• 3 years ago

In 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_Mckay3• 3 years ago

All 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 =)

jiggawatts3• 3 years ago

I 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!

bboygravity• 3 years ago

But last I checked there's Julia for Python? So you could get the best of both world (nice syntax of Python and the speed of Julia).

ftxbro• 3 years ago

When I googled digits of pi I found http://www.numberworld.org/y-cruncher/ where they get 25M digits in about one second (rather than one minute) depending on what kind of computer you use.

juliusgeo18• 3 years ago

y-cruncher is the state of the art but it is also not written in a language like julia, and doubtlessly occupies many more lines of code

ewams• 3 years ago

I've done 2.5 billion digits in 88 seconds with ycruncher. That's about 28.4 million digits a second?

http://ewams.net/?date=2020/04/17&view=sockets_vs_cores_for_...

Though that is with a broadwell CPU and about 3 years ago. Might see if I can get my hands on a 6348 ice lake CPU and see how she burns.

boredemployee• 3 years ago

genuine question: would that amount of digits be useful for any application?

fisherjeff• 3 years ago

https://www.jpl.nasa.gov/edu/news/2016/3/16/how-many-decimal...

juliusgeo18• 3 years ago

Thank you for posting this! This is exactly what I was looking for.

yesenadam• 3 years ago

TLDR: "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!

juliusgeo18• 3 years ago

No, 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• 3 years ago

What is the circumference of the universe?

melagonster• 3 years ago

https://github.com/philipl/pifs

but this is a joke hahaha

joseftexas• 3 years ago

Isn'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.

juliusgeo18• 3 years ago

That is why I compared it to the stdlib implementation.

Temporary_31337• 3 years ago

Just 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/

userbinator• 3 years ago

I still remember when the SuperPI 1M <10s record was beaten... with a CPU overclocked to 5GHz+: 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

Clock cycle efficiency doesn't seem to have improved much in the last decade.

Someone2• 3 years ago

Nitpick: 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, it seems you can do

  BigFloat(π, precision=20000000)

How does that perform?

SherlockSage• 3 years ago

I 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• 3 years ago

(Plain Windows took ~26 seconds for some reason?)

Maybe the allocations were the reason?

Someone2• 3 years ago

I’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?

juliusgeo18• 3 years ago

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

juliusgeo18• 3 years ago

BigFloat internally is not binary, it uses an array because it is made for arbitrary precision. There will absolutely be small errors at the end but it’s easier to say 20M than 1.999999M digits of pi.

jiggawatts3• 3 years ago

In 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• 3 years ago

does Mathematica actually compute it or is it just loading the value that it already knows?

extasia• 3 years ago

7.6 second seems too much to be just loading a value from disk

jiggawatts3• 3 years ago

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

moelf2• 3 years ago

it's not using the same algorithm right?

juliusgeo18• 3 years ago

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

bakul2• 3 years ago

The 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).

juliusgeo18• 3 years ago

Interesting! 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…

bakul2• 3 years ago

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

20M digits of pi in 1 minute using Julia

Related Stories