Informal Call for a Julia π - day challenge

[Alan Edelman]

With about a week left, I'd love to find out how many digits of  π  we can get only from Julia.

Perhaps we can coordinate a worldwide distributed computation this week.

[Ivar Nesje]

julia> big(pi)

3.141592653589793238462643383279502884197169399375105820974944592307816406286198e+00 with 256 bits of precision

[Alan Edelman]

with the last two digits being roundoff error

[Luis Benet]

julia> set_bigfloat_precision(10000)
10000

julia> big(pi)
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802759009946576407895126946839835259570982582262052248940772671947826848260147699090264013639443745530506820349625245174939965143142980919065925093722169646151570985838741059788595977297549893016175392846813826868386894277415599185592524595395943104997252468084598727364469584865383673622262609912460805124388439045124413654976278079771569143599770012961608944169486855584840635342207222582848864815845602850601684273945226746767889525213852254995466672782398645659611635488623057745649803559363456817432411251507606947945109659609402522887971089314566913686722874894056010150330861792868092087476091782493858900971490967598526136554978189312978482168299894872265880485756401427047755513237964145152374623436454285844479526586782105114135473573952311342716610213596953623144295248493718711014576540359027993440374200731057853906219838744780847848968332144571386875194350643021845319104848100537061468067491927819119793995206141966342875444064374512371819217999839101591956181467514269123974894090718649423196156794520803e+00 with 10000 bits of precision

[Johan Sigfrids]

Because calling out to mpfr is cheating here is a implementation of Gauss-Legendre in Julia calculating Pi to the desired number of digits (Still relying on BigFloats though).

function gaussLegendrePi(d::Integer)
    prec = get_bigfloat_precision()
    set_bigfloat_precision(iceil(d*log(10)/log(2))+10)
    n = iceil(log2(d))
    a = BigFloat(1)
    b = a / sqrt(BigFloat(2))
    t = a / BigFloat(4)
    x = a
    while n > 0
        y = a
        a = (a + b) / 2
        b = sqrt(b * y)
        t = t - x * (y - a)^2
        x = 2 * x
        n -= 1
    end
    p = ((a + b)^2) / (4 * t)
    println(string(p)[1:d+1])
    set_bigfloat_precision(prec)
    nothing
end

I haven't done any validation beyond checking that it and Mathematica agree on the millionth digit of Pi.

[Simon Byrne]

How about an ijulia notebook of 101 ways to compute pi in Julia, ideally with accompanying graphics? 

[Alan Edelman]

314 ways?

[Simon Byrne]

Here's my nomination: unfortunately, it only works on 0.4 using LLVM >=3.5:

function asmpi()

    Base.llvmcall("""%pi = call double asm "fldpi", "={st},~{dirflag},~{fpsr},~{flags}"() #0

        ret double %pi""", Float64, ())

end

[Hans W Borchers]

I can provide an implementation of the "droplet" algorithm for computing pi, based on work by Borwein, Rabinowitz, and Wagon. It does neither use 'mpfr' nor big floats, only simple arithmetic. It will spit out about 50,000 digits one after the other -- with more digits it gets too slow. A speed-up of the factor 10 or so will be possible with some tricks.
As a test I am abut to reimplement it in Nim. As Nim can also generate C or C++ code, it will be interesting to compare C, Nim, and Julia on this non-trivial, but straightforward algorithm.

[Alan Edelman]

Is this the  one that can be parcelled out to all of us so we can distribute the work?

I realize later digits take longer, so early folks can get more digits

and later digit folks might only do fewer digits

[Hans W Borchers]

Yes, of course.

I have to retrieve the code from a backup drive, so this will have to wait until I'm home and can search for it. Guess it will take until tomorrow evening. I have backuped all my Julia work and removed the files from the main drive. I will return to Julia only when version 0.4 is officially out and running, and most optimization packages and the 'Approxfun' package are working correctly.

At first I was a bit irritated about the length of this Julia "chaos phase", but in the meantime I find it's a nice break and provides an opportunity to go and try out some other new programming languages, like Nim, Rust, or Swift. That is the reason I am reimplementing some of these algorithms in other languages. Well, I'm sure I will return in time.

[Simon Danisch]

Any chance, that you'll write up your findings in a blog post or whatnot!?

[Hans W Borchers]

The following function is an implementation of the famous "droplet" algorithm of Rabinowitz and Wagon for calculating \pi to n (internal) decimal places (no analytic functions, no arbitrary precision, some digits could even be calculated by hand):

    function droplet_pi(n)
        Pi = p = ""
        no_nines = 0
        d = n + 2
        N = ifloor(10*d/3.0 + 1)
        a = fill(2, N+1)
        for l in 1:d
            a = 10 * a
            for i in (N+1):-1:2
                j = 2*i - 1
                q, r = divrem(a[i], j)
                a[i] = r
                a[i-1] += q * (i-1)
            end
            q, r = divrem(a[1], 10)
            a[1] = r
            if q < 9
                Pi *= string(p) * ("9"^no_nines)
                p = q
                no_nines = 0
            elseif q == 9
                no_nines += 1
            elseif q == 10
                p += 1
                Pi *= string(p) * ("0"^no_nines)
                p = 0
                no_nines = 0
            else
                error("droplet_pi: algorithm error!")
            end
        end
        Pi[1:1] * "." * Pi[2:end]
    end

    julia> @time droplet_pi(1000)
    elapsed time: 0.166306075 seconds (27629136 bytes allocated)
    "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348
    2534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555
    9644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610
    4543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600
    1133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807
    4462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190
    7021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577
    8960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403
    4418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908
    3026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534
    9042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019
    89"

This implementation continuously appends digits to a string. It will probably be faster to preallocate a vector of (integer) digits and to combine these into a string in one step at the end of the procedure. Maybe I didn't know how to do this when I wrote this code.

[David van Leeuwen]

Intriguing... 

searches seem to refer to this as a Spigot algorithm, but I find no references to droplet. 

---david

[Steven G. Johnson]

If we're doing a 31 ways to to compute pi, I would suggest that most of them should be more aimed at fun than trying to compute a huge number of digits.  Like Monte Carlo integration of the area of a unit circle, or summing a slowly converging series, using quadgk(x -> 1/(1+x^2), 0,1)[1]*4, or other methods with no concern for efficiency.

[Alan Edelman]

fun and/or digits (wherever your fun lies!)

buffon needle could be one 

[Simon Byrne]

On Monday, 9 March 2015 13:06:46 UTC, Alan Edelman wrote:

buffon needle could be one 

That has the nice benefit of leading to a discussion of how "lucky" Lazzarini was in getting his approximation:

http://en.wikipedia.org/wiki/Buffon's_needle#Estimating_.CF.80

[Hans W Borchers]

True, 'spigot' is another name for this class of algorithms. Borwein and Bailey, in their book "Mathematics by Experiment", say:

    A spigot method for a numerical constant is one that can produce digits one by one ("drop by drop").

so perhaps the term 'droplet' stems from this "image in our mind's eye."

@Alan
I may have misunderstood you. This is *not* the Borwein & Borwein approach using their intricate `arctan` series (a BBP-like formula). That approach generates 16 bits per loop, but strangely enough reconstructing the decimal representation takes more time than computing the binary representation itself.

@Steven
Do you think it was no fun to implement the droplet/spigot algorithm? Then you may be completely wrong.

[Steven G. Johnson]

On Monday, March 9, 2015 at 11:24:06 AM UTC-4, Hans W Borchers wrote:

@Steven
Do you think it was no fun to implement the droplet/spigot algorithm? Then you may be completely wrong.

I'm sure it was fun.  The point I meant to make was that fun is not restricted to efficient algorithms, and to focus exclusively on the latter would be unfortunate. 

[Hans W Borchers]

Steven mentioned he would be interested in computing pi through a slowly converging series. Well, I think the Leibniz series might be a good example:

    pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 -+ ...

Let us define the first hundred terms of this series and see what their sum is:

    julia> k = [1:100];
    julia> leibniz = (-1).^(k-1) ./ (2*k-1);

    julia> 4*sum(leibniz)
    3.1315929035585537

This is not even correct in the second digit after the decimal point.

What can we do to accelerate the convergence of an alternating sum? A very powerful method has been proposed by Henri Cohen, F. Rodriguez Villegas, and Don Zagier (2000). The following function is my implementation of 'Algorithm 1' from their paper:

    function sumalt(s_alt, n)
        b = 2^(2*n-1)
        c = b
        s = 0.0
        for k in (n-1):-1:0
            t = s_alt[k+1]
            s = s + c*t
            b = b * (2*k+1) * (k+1) / (2 * (n-k) * (n+k))
            c = c + b
        end
        s = s / c
        return(s)
    end

Let us apply this to the first 20(!) terms of the Leibniz series:

    julia> 4*sumalt(leibniz, 20)
    3.1415926535897936

The result is pi with an accuracy of 2*eps() = 4.44e-16 or 15 digits.

Using Big Floats with 256 bit, we could easily reach 75 correct digits. Of course, `sumalt` is general and can be applied in other contexts as well, for example to compute the `[z]eta` function for complex values very accurately.

[Simon Byrne]

On a pi-related note: one algorithm that would be useful in Base is a remainder modulo-pi. At the moment we call out to a routine in openspecfun, but this requires allocating a 2-element array on each call (and also has some issues on certain compilers (see, for example https://github.com/JuliaLang/julia/issues/8799).

Having this would also make it possible to implement a faster "sincos" function (https://github.com/JuliaLang/julia/issues/10442).

I've made a start on translating it here:

https://github.com/JuliaLang/julia/compare/master...simonbyrne:sincos

However it still requires an implementation of the Payne-Hanek range reduction for large (>1e6) arguments. For reference, the openspecfun code is available here:

https://github.com/JuliaLang/openspecfun/blob/master/rem_pio2/k_rem_pio2.c

Anyone interested?

-Simon

[David P. Sanders]

El sábado, 7 de marzo de 2015, 11:59:47 (UTC-6), Alan Edelman escribió:

With about a week left, I'd love to find out how many digits of  π  we can get only from Julia.

Perhaps we can coordinate a worldwide distributed computation this week.

In a different direction, it seems not to be well known that we can calculate *rigorous bounds*

on π using just floating-point aritmetic. I have made an example notebook that may be viewed here:

http://nbviewer.ipython.org/gist/dpsanders/3e69e4a2d4f5b5727258

[Comments of course welcome!]

Best,

David.

[Shashi Gowda]

Late to the party but here's a simple monte carlo simulation visualized

using Reactive, Interact, Compose

randpoint() = (rand(), rand())

isin(x, y) = 0.25 > (x-0.5)^2 + (y-0.5)^2

Compose.set_default_graphic_size(8inch, 8inch)

@manipulate for switch=false,

    t=fpswhen(switch, 10),

    point=lift(x -> randpoint(), t),

    points=foldl(push!, Any[randpoint()], point),

    incount=foldl((x, y) -> x + isin(y...), 1, point),

    count = foldl((x, y) -> x+1, 1, point)

    compose(context(),

        (context(), text(0, 0.1, "π ≊ $(4*incount/count)")),

        [(context(), fill("red"), circle(x, y, 0.005)) for (x, y) in points]...,

        (context(), fill("lightblue"), circle()))

end

Also http://arxiv.org/pdf/1404.1499v2.pdf ;)