PI (=3.14...) Calculation Program based on FFT and AGM


Description

This is a package to calculate PI in huge digits. I programmed it for a benchmark of a FFT based multiple-precision routine.

Distribution

Files in the Package (pi_fftc6)

Makefile : make file - for 32bit CPU: digits <= 6,000,000,000
Makefile_64bit : make file - for 64bit CPU
Makefile_quad : make file - for 128bit FPU
config.h : machine depending macros
fftsgx.c : FFT Package in C - split-radix - use work areas
fftsg_hx.c : FFT Package in C - split-radix - use no work areas
pi_fftca.c : PI Calc. Program - standard version - use fft*gx.c"
pi_fftcs.c : PI Calc. Program - mem save version - use "fft*g_hx.c"
pi_fftcw.c : PI Calc. Program - mem swap version - use "fft*g_hx.c"
dgt_div.c : data converter - from pi.dat to Super_PI format
README.TXT : readme file
win32bin/ : DIR
win32bin/Makefile : make file - for intel C++ compiler
win32bin/pi_cs.exe : Win32 exec. - mem save version
win32bin/pi_cs_thread.exe : Win32 exec. - use FFT level threads
win32bin/pi_cw.exe : Win32 exec. - mem swap version
win32bin/pi_cw_thread.exe : Win32 exec. - use FFT level threads
win32bin/dgt_div.exe : data converter - pi.dat => pi.txt a printout format

To Compile (pi_fftc6)

Check macros in "config.h" and modify if necessary. FFT_ERROR_MARGIN is very impotant parameter. If FFT_ERROR_MARGIN is very small then efficiency will be bad. If FFT_ERROR_MARGIN >= 0.5 then it may calculate a wrong result.

Performance (pi_fftc5: old version)

Number of Floating Point Operations (pi_fftc5)

---- Relationship between Number of Digits and FFT Length ----
ndigit = nfft*log_10(R), R >= 10000 or 1000
R is a radix of multiple-precision format. R depends on DBL_ERROR_MARGIN and FFT+machine+compiler's tolerance.

Memory Use (pi_fftc5)

FFT Algorithm

AGM Algorithm

  ---- a formula based on the AGM (Arithmetic-Geometric Mean) ----
    c = sqrt(0.125);
    a = 1 + 3 * c;
    b = sqrt(a);
    e = b - 0.625;
    b = 2 * b;
    c = e - c;
    a = a + e;
    npow = 4;
    do {
        npow = 2 * npow;
        e = (a + b) / 2;
        b = sqrt(a * b);
        e = e - b;
        b = 2 * b;
        c = c - e;
        a = e + b;
    } while (e > SQRT_SQRT_EPSILON);
    e = e * e / 4;
    a = a + b;
    pi = (a * a - e - e / 2) / (a * c - e) / npow;
  ---- modification ----
    This is a modified version of Gauss-Legendre formula
    (by T.Ooura). It is faster than original version.

Reference

  1. E.Salamin, Computation of PI Using Arithmetic-Geometric Mean, Mathematics of Computation, Vol.30 1976.
  2. R.P.Brent, Fast Multiple-Precision Evaluation of Elementary Functions, J. ACM 23 1976.
  3. D.Takahasi, Y.Kanada, Calculation of PI to 51.5 Billion Decimal Digits on Distributed Memoriy Parallel Processors, Transactions of Information Processing Society of Japan, Vol.39 No.7 1998.
  4. T.Ooura, Improvement of the PI Calculation Algorithm and Implementation of Fast Multiple-Precision Computation, Information Processing Society of Japan SIG Notes, 98-HPC-74, 1998.

License


Links


Main Page / FFT Page