Calculating Pi in Parallel on GPU and CPU, using the Bailey Borwein Plouffe formula.

Remember the benchmarking program Super Pi? It is a prime example of a benchmark which fell out of use because it was single threaded in an increasingly threaded era. All modern CPUs released in the last few years have at least 2 cores. The latest enthusiast chips have 32 cores and 64 threads (AMD Threadripper 2990WX). While single thread performance does improve with each generation, to make any real use of the most powerful processors, parallel programming must be used.

GPUs on the other hand have many small compute processors which all put together are often more powerful than a CPU. But on a GPU each thread is slow on its own so they are totally useless for any task which is unable to scale to use all the available performance.

The video below is a demonstration of my implementation of the Bailey–Borwein–Plouffe formula¹.

This is able to calculate any specified hexadecimal place of Pi, within certain limitations related to computational precision. It is an excellent candidate for parallel implementation and the two programs I wrote, one for CPU, one for GPU, are able to scale to almost any number of threads. For calculating the 10,000,000th hexadecimal place, 40 million threads could be used in theory, although based on my experience overheads would mean that 40,000 or so might be better.

My GPU implementation uses 700 threads and the CPU implementation uses only 8 yet it is not far behind. This is because the GPU threads are much slower. GPUs perform at their best when working on shorter data types. The Double Precision (fp64) performance of my R9 290 is rated at 600 GFLOPS³ yet the Single Precision (fp32,float in C) performance is said to be 4.8 TFLOPS. A ratio of 8. However, despite being much newer the vast majority of current GPUs have ratios of 16(AMD) or 32(Nvidia) crippling double precison performance. Only workstation class cards have meaningful double precison performance.

For reasons better explained by D. Bailey² in his excellent paper on the formula, it is necessary to use fp64 to estimate hexadecimal places up to 10 million. With 128-bit precision places up to 100 trillion can be estimated. However, CPUs can only handle 80-bit (aka long double) natively in hardware and GPUs only fp64 in hardware. Libraries can be used to work with 128-bit precision but as these are not natively implemented in hardware they are as much as 100 time slower than 64-bit.

In time I will add more to this page about the formula itself and my implementation. Including download links once I have implemented the GPU algorithm in CUDA and compiled win10 versions.

References