julia_bbp.ipynb

Kernel: Julia

Julia version of BBP

Experiments with BBP algorithm in Julia.

First version

"""
For integer \$d\$, \$k\$, and \$j\$ efficiently
compute \$16^{d-k} \\bmod 8k + j ) / (8k + j)\$.
"""
ratio(d, k, j) = powermod(16, d - k, 8k + j) / (8k + j)

"""
Fractional part of real \$x\$.
"""
frac_part(x) = x - floor(x)

"""
For integer \$d\$ and \$j\$ computes the fractional part of
\$\$
    \\sum_{k=0}^d \\frac{16^{d-k} \\bmod 8k+j}{8k+j} +
    \\sum_{k=d+1}^\\infty \\frac{1}{16^{k-d} (8k+j)}.
\$\$
"""
function sigma(d, j)
    r = 0.

    for k in 0:d
        r += ratio(d, k, j)
        if r > 1.0
            r = frac_part(r)
        end
    end

    k = d + 1
    u = 1 / (16^(k - d) * (8k + j)) 
    
    while u > eps(Float64)
        r += u
        k += 1
        u = 1 / ((16^(k-d))*(8k + j))
    end

    return frac_part(r)
end

"""
Compute \$d + 1\$-th hexadecimal digit of the fractional part of \$\\pi\$. 
"""
function digpi(d)
    r = 4*sigma(d, 1) - 2*sigma(d, 4) - sigma(d, 5) - sigma(d, 6)
    return floor(Int, 16*frac_part(r))
end

digpi (generic function with 1 method)

$1\,000\,000$th hexadecimal digit of $\pi$ (after the decimal point).

digpi(999999)

2

A few more digits...

map(digpi, 999999:999999+5)

6-element Array{Int64,1}:
  2
  6
 12
  6
  5
 14

How fast it is? To measure the performance we will use BenchmarkTools.jl.

using BenchmarkTools

Take $d = 10\,000$.

@benchmark digpi(10000)

BenchmarkTools.Trial: 
  samples:          76
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  0.00 bytes
  allocs estimate:  0
  minimum time:     65.35 ms (0.00% GC)
  median time:      65.45 ms (0.00% GC)
  mean time:        65.80 ms (0.00% GC)
  maximum time:     70.41 ms (0.00% GC)

Faster version

It is possible to speed things up a little bit (but not too much).

"""
For integer ``d``, ``k``, and ``j`` efficiently compute
``16^{d-k} \\bmod 8k + j ) / (8k + j)``.
"""
@fastmath function powermod2(a, b, n)
    r = 1
    while b > 0
        if isodd(b)
            r = r * a % n
        end
        b = div(b, 2)
        a = a * a % n
    end
    return r
end

function ratio2(d, k, j)
  @fastmath powermod2(16, d - k, 8k + j) / (8k + j)
end

function sigma2(d::Int64, j::Int64)
    r = 0.
    @fastmath @simd for k in 0:d
        r += ratio2(d, k, j)
        if r > 1.0
            r = frac_part(r)
        end
    end

    kk  = 8. * (d + 1) 
    u   = 1 / 16
    rat = u / kk

    while rat > eps(Float64)
        @fastmath r   += rat
        @fastmath kk  += 8.
        @fastmath u   /= 16.
        @fastmath rat  = u / kk
    end

    return frac_part(r)
end

function digpi2(d::Int64)
    r = 4*sigma2(d, 1) - 2*sigma2(d, 4) - sigma2(d, 5) - sigma2(d, 6)
    return floor(Int, 16*frac_part(r))
end

digpi2 (generic function with 1 method)

digpi2(999999)

2

map(digpi2, 999999:999999+5)

6-element Array{Int64,1}:
  2
  6
 12
  6
  5
 14

Again, speed test for $d = 10\,000$.

@benchmark digpi2(10000)

BenchmarkTools.Trial: 
  samples:          165
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  0.00 bytes
  allocs estimate:  0
  minimum time:     29.49 ms (0.00% GC)
  median time:      29.62 ms (0.00% GC)
  mean time:        30.31 ms (0.00% GC)
  maximum time:     41.13 ms (0.00% GC)

This is twice as fast as the previous Julia version. It is a little bit slower then the Cython version.

Parallel version

"""
Parallel version.
"""
function digpi3(d::Int64)
  if nworkers() < 4
    error("Expects four available workers! Found only $(nworkers()).")
  end

  s1 = remotecall(4, sigma2, d, 1)
  s4 = remotecall(3, sigma2, d, 4)
  s5 = remotecall(2, sigma2, d, 5)
  s6 = remotecall(1, sigma2, d, 6) 
  r = 4*fetch(s1) - 2*fetch(s4) - fetch(s5) - fetch(s6)
  return floor(Int, 16*frac_part(r))
end

digpi3 (generic function with 1 method)

nworkers()

4

addprocs(1)

1-element Array{Int64,1}:
 5

nworkers()

4

@everywhere include("bbp.jl")

digpi3(10000)

8

@benchmark digpi3(10000)

BenchmarkTools.Trial: 
  samples:          248
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  34.16 kb
  allocs estimate:  593
  minimum time:     15.64 ms (0.00% GC)
  median time:      20.72 ms (0.00% GC)
  mean time:        20.22 ms (0.00% GC)
  maximum time:     26.43 ms (0.00% GC)

Test on larger input.

@benchmark digpi3(1000000)

BenchmarkTools.Trial: 
  samples:          3
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  34.88 kb
  allocs estimate:  619
  minimum time:     2.09 s (0.00% GC)
  median time:      2.15 s (0.00% GC)
  mean time:        2.16 s (0.00% GC)
  maximum time:     2.25 s (0.00% GC)

@benchmark digpi2(1000000)

BenchmarkTools.Trial: 
  samples:          2
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  0.00 bytes
  allocs estimate:  0
  minimum time:     4.54 s (0.00% GC)
  median time:      4.60 s (0.00% GC)
  mean time:        4.60 s (0.00% GC)
  maximum time:     4.66 s (0.00% GC)