CIS-4230/5230 Homework #4 (Computing Pi)

Due: Thursday, March 14, 2024 (not really, it's actually Friday, April 19, 2024)

In this assignment we will compute π to high precision. We will attempt to compute one million decimal digits of π, but we may have to tone down that goal depending on how things go. In addition to using an advanced algorithm, we will also parallelize the computation using a combination of MPI and OpenMP.

This assignment is ambitious. It is also not well-specified. Regard the information below more as a collection of notes than as specific instructions to follow.

The Algorithms

There are several advanced methods for computing π. For this assignment we will focus on possibilities: using Ramanujan's π formula, which entails evaluating an infinite series.

The terms in Ramanujan's formula are independent and can be computed in parallel in a relatively straightforward way. Also, each term adds about eight digits of precision to the final answer. Thus, about 125,000 terms are needed to one million decimal digits of precision.

In general, and this is typical of these kinds of computations, the computation gets more difficult as the desired precision increases. In this case, the terms in Ramanujan's formula get more difficult to evaluate. Here "difficult" means more time-consuming and more memory intensive.

As an aside, the best known algorithm for computing π is the Chudnovsky algorithm. It is an iterative method and more complex to parallelize because each iteration depends on the previous iteration. However, the current world record of 100 trillion digits was computed using Chudnovsky. That computation finished on March 21, 2022, and was done by Emma Haruka Iwao, a computer scientist at Google. She used a program called y-cruncher that implements the algorithm, running the program on a machine with 128 PEs, 864 GB of RAM and 663 TB of storage. The computation required 158 days.

The Approach

We will need to use arbitrary precision integers and arbitrary precision floating point values. Here "arbitrary precision" means "as many digits as required." Implementing this functionality is itself a non-trivial exercise. We will use two well-established and highly optimized libraries: GMP (GNU Multi-Precision Arithmetic Library) and MPFR (GNU Multi-Precision Floating Point Reliable Library).

The "reliable" part of MPFR's name refers to the fact that it faithfully extends IEEE 754 floating point semantics, including proper rounding behavior, to extended precision values. It will produce the same results on all systems on which it runs. GMP contains support for extended precision floating point values as well, but it is less reliable. It is recommended to use MPFR in all new projects.

The documentation for these two libraries will be essential. I call those pages out here: GMP, MPFR.

Both of these libraries have been installed on Lemuria and on the cluster nodes. To install them on a Debian-like system, such as Ubuntu, use:

  $ sudo apt install libgmp-dev libmpfr-dev

On macOS using homebrew, the libraries can be installed using:

  $ brew install gmp mpfr

On Cygwin, use the setup program to search for libgmp-devel and libmpfr-devel.

To you help understand how to use GMP and MPFR, consider this sample program that computes e, the base of natural logarithms, using an infinite series in a manner similar to the series used by Ramanujan's formula. To compile this program using a command such as the follow.

  $ gcc -Wall -o compute-e compute-e.c -lgmp -lmpfr

Ramanujan's Formula

Srinivasa Ramanujan, an Indian mathematician from the early 20th century, proved the following result (taken from Wikipedia).

Ramanujan's Formula

We will call the variable k the term index variable. Notice that the first term when k = 0 is 1103. However, the latter terms (larger k) are smaller because of the k!4 in the denominator of each term.

This formula adds about eight decimal digits of precision with each term computed. Suppose we wanted to calculate one million decimal digits of π. That requires approximately 125,000 terms to be evaluated. The last term thus involves computing 125000!4, which is an exceedingly large number. That's where GMP comes in.

The numerator and denominator of each term can be represented by integers (using type mpz_t). However, those integers need to be converted to floating point values and divided using MPFR. The MPFR library requires that each floating point value (of type mpfr_t) be configured with the desired precision when it is initialized. One million decimal digits requires 1000000*log_2(10) = 3321928 bits of precision. This is important because MPFR requires that precision be stated in bits.

The Steps

Do not attempt to compute π to one million digits of precision right away. Move toward that goal in smaller steps.

  1. First, compile the sample program to ensure you have the necessary libraries installed. You can also use the sample program as a starting point for your π-computing program.

  2. Write a single threaded program that computes π to an adjustable number of digits. The desired precision can be declared as a constant in the program. Start small with, say, 100 digits of precision. You can likely find the correct result online to check your work.

  3. Increase the desired precision in steps until the program requires an inconvenient about of time to finish. Can you go as high as one million digits? Can you go higher? Use your "final" version of the serial program as your baseline.

  4. Using OpenMP, parallelize your program. Split the loop that computes the terms manually. Simply using omp parallel for won't work well (do you see why not?). The calculation of k!, and certain other factors, needs to be handled carefully.

  5. Measure the performance of the OpenMP version to see if the parallelization is helping.

  6. Finally, parallelize the program again using MPI to see if you can accelerate it further.

The Gauss-Legendre Algorithm

This section describes an alternative method for computing π using the Gauss-Legendre algorithm. It is provided for your information only. While I encourage you to experiment with this method, it is not a required part of the assignment.

Gauss-Legendre is an iterative method. In some ways it is easier to implement than Ramanujan's formula, but it is awkward to parallelize. However, it converges extremely quickly, approximately doubling the number of correct digits with each iteration! Thus, computing one million digits of π only requires about 20 iterations of the algorithm. The method is as shown below (taken from Wikipedia).

Gauss-Legendre Algorithm

It is possible to run the computations of an, bn, tn, and pn in parallel, but even that must be done with care to account for the thread safety properties of the MPFR library.

Submission

Submit ... along with ... . Pack all the files into a zip archive and submit that archive to Canvas.


Last Revised: 2024-03-14
© Copyright 2024 by Peter Chapin <peter.chapin@vermontstate.edu>