using PureFun
using PureFun.Lazy: Stream, @stream

Example: Approximating $\pi$

This example is adapted from Chapter 3 of Structure and Interpretation of Computer Programs

The summation:

\[\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7}\]

gives us a way to approximate $\pi$:

function pi_summands(n, one)
    @stream(Float64, one/n, pi_summands(n+2, -one))
end

function approx_pi(n)
    summands = pi_summands(n, 1)
    sums = accumulate(+, summands, init=0.0)
    map(x -> 4x, sums)
end

π̂ = approx_pi(1)

PureFun.Lazy.NonEmpty{Float64}
4.0
2.666666666666667
3.466666666666667
2.8952380952380956
3.3396825396825403
2.9760461760461765
3.2837384837384844
...

The series converges slowly though. In order to see how close the estimates are:

map(x -> abs(π - x), π̂)

PureFun.Lazy.NonEmpty{Float64}
0.8584073464102069
0.47492598692312615
0.32507401307687367
0.2463545583516975
0.19808988609274714
0.16554647754361662
0.1421458301486913
...

A better approximation

An accelerator is a function that takes a series, and returns a series that converges to the same sum, but more quickly. The Euler transform is an accelerator that works well on series with alternating positive/negative terms, like we have here. It is defined as:

\[S_{n+1} - \frac{(S_{n+1}-S_{n})^{2}}{S_{n-1} - 2S_{n} + S_{n+1}}\]

So:

euler_transform(s₀, s₁, s₂) = s₂ - (s₂ - s₁)^2/(s₀ - 2s₁ + s₂)

function euler_transform(s::Stream)
    @stream(Float64,
          euler_transform(s[1], s[2], s[3]),
          euler_transform(tail(s)))
end

π̂₂ = euler_transform(approx_pi(1))

PureFun.Lazy.NonEmpty{Float64}
3.166666666666667
3.1333333333333337
3.1452380952380956
3.13968253968254
3.1427128427128435
3.1408813408813416
3.142071817071818
...

This converges much more quickly to the true value of $\pi$

map(x -> abs(π - x), π̂₂)

PureFun.Lazy.NonEmpty{Float64}
0.025074013076873847
0.008259320256459368
0.0036454416483024943
0.0019101139072530415
0.0011201891230503414
0.0007113127084514836
0.0004791634820247026
...

and we can reuse the accelerator to keep improving our approximation:

euler_transform(π̂₂)

PureFun.Lazy.NonEmpty{Float64}
3.142105263157895
3.1414502164502167
3.1416433239962656
3.1415712902014277
3.1416028416028423
3.141587320947787
3.141595655236941
...

...

π̂₂ |> euler_transform |> euler_transform |> euler_transform

PureFun.Lazy.NonEmpty{Float64}
3.1415927140337785
3.141592637113005
3.1415926587096235
3.141592651803974
3.141592654277287
3.141592653301986
3.1415926537192127
...

This page was generated using Literate.jl.