Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

SageMath, Cython and Julia: BBP algorithm

198 views
Tweet
Share
Share
Kernel: SageMath 7.2.rc2

SageMath and Julia

Tomáš Kalvoda, KAM FIT ČVUT, Sage Days 72, May 22, 2016

sage logo julia logo

1. Outline

  1. outline of the talk

  2. how do I use Sage and Julia

  3. comparison of Sage and Julia

  4. practical comparative example (implementation of BBP algorithm)

2. How do I use SageMath/Julia

SageMath:

  • teaching math (basics of calculus and linear algebra)

  • topics for students

Julia:

  • numerical experiments (related to my research)

  • topics for students

3. SageMath and Julia

SageMath

Most important facts were already mentioned yesterday.

  • aims to be an open-source alternative to MMM

  • uses Python to glue together various math-related software

  • contains considerable ammount of original code too

  • its monolithic codebase is hosted on trac server

  • uses Cython for speed (more about this in a moment)

Julia

Motto: "a new approach to technical computing."

  • it attempts to

    • be fast

    • be easy and fun to use

    • replace C and FORTRAN (in numerical computing)

  • relatively small core language library hosted on github

  • abundance of optional user contributed packages (~960)

Key features:

  • JIT (Just In Time) compilation

  • dynamically typed, with optional type annotation

  • multiple-dispatch

We will demonstrate Sage/Cython/Julia on the following...

4. Bailey-Borwin-Plouffe (BBP) formula

The Bailey-Borwin-Plouffe formula was discovered in 1996. It expresses π\pi in the form of a particular infinite series,

π=∑k=0∞116k(48k+1−28k+4−18k+5−18k+6).\pi = \sum_{k=0}^\infty \frac{1}{16^k} \left( \frac{4}{8k+1} - \frac{2}{8k+4} - \frac{1}{8k+5} - \frac{1}{8k+6} \right).

One can use this formula to compute hexadecimal digits of π\pi at an arbitrary position (without the need to compute any of previous digits).

The algorithm

Let dd be a non-negative integer. Multiply the relation by 16d16^d, and split the sum into four parts

16dπ=4∑k=0∞16d−k8k+1−2∑k=0∞16d−k8k+4−∑k=0∞16d−k8k+5−∑k=0∞16d−k8k+6.16^d \pi = 4 \sum_{k=0}^\infty \frac{16^{d-k}}{8k+1} - 2 \sum_{k=0}^\infty \frac{16^{d-k}}{8k+4} - \sum_{k=0}^\infty \frac{16^{d-k}}{8k+5} - \sum_{k=0}^\infty \frac{16^{d-k}}{8k+6}.

We are interested only in the fractional part of the left-hand side,

{16dπ}={4σ(d,1)−2σ(d,4)−σ(d,5)−σ(d,6)},\{ 16^d \pi \} = \big\{ 4 \sigma(d, 1) - 2 \sigma(d, 4) - \sigma(d, 5) - \sigma(d, 6) \big\},

where we have used notation

σ(d,j):={∑k=0∞16d−k8k+j}={∑k=0d16d−k mod (8k+j)8k+j+∑k=d+1∞16d−k8k+j}.\sigma(d, j) := \left\{ \sum_{k=0}^\infty \frac{16^{d - k}}{8k + j} \right\} = \left\{ \sum_{k=0}^{d} \frac{16^{d-k} \bmod (8k + j)}{8k+j} + \sum_{k=d+1}^\infty \frac{16^{d-k}}{8k+j} \right\}.

The denominator of the fraction in the first sum can be computed efficiently using the square and multiply algorithm. The second sum can be evaluated approximately in the given floating point arithmetic since it decays very rapidly.

Fun part: Let us now try to implement this in SageMath and Julia.