Frequency Spectrum of Pi

Pi is an irrational number. Truncated to 30 decimal digits, its representation is

3.14159 26535 89793 23846 26433 83279

Pi has been computed to more than trillion decimal digits. What I find really fascinating is that no pattern has been found among trillion+ digits. I spent a bit of time searching the Internet for the probability distribution of decimal fractions of Pi. I asked whether the decimal digits are distributed for example as White Noise. I didn't find anything immediately obvious, so I wrote a simple program to calculate the frequency spectrum of Pi.

If it does behave as White Noise, the spectrum should contain relatively equal amplitudes. I used the awesome fftw library for DFT. You can find the source code of my little silly program at the end of this post. I ran it on various number of decimal digits of Pi up to 1946 digits, and I always get what seems to be completely random spectrum, where minimum and maximum amplitudes differ up to around 100 times. It doesn't look like a white noise to me, or does it? It is more like a totally stochastic process, if there is such a thing. You can find the spectrum graph below. The amplitudes fluctuate between around 0.002 and 0.3+.

The average of all digits is 4.5. Considering the absence of any pattern, the 4.5 average means that the digits (from 0 to 9) are very evenly distributed. It should be easy to confirm by just counting the number of times a digit appears in the sequence.





The source code:

// pi.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <math.h>
#include <fftw/fftw3.h>

//1946 digits of pi
const int pi_fract[]= {

1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3,3,8,3,2,7,9,5,0,2,8,8,4,1,9,7,1,6,9,3,9,9,3,7,5,1,0,5,8,2,0,9,7,4,9,4,4,5,9,2,3,0,7,8,1,6,4,0,6,2,8,6,2,0,8
,9,9,8,6,2,8,0,3,4,8,2,5,3,4,2,1,1,7,0,6,7,9,8,2,1,4,8,0,8,6,5,1,3,2,8,2,3,0,6,6,4,7,0,9,3,8,4,4,6,0,9,5,5,0,5,8,2,2,3,1,7,2,5,3,5,9,4,0,8,1,2,8,4,8,1,1,1,7
,4,5,0,2,8,4,1,0,2,7,0,1,9,3,8,5,2,1,1,0,5,5,5,9,6,4,4,6,2,2,9,4,8,9,5,4,9,3,0,3,8,1,9,6,4,4,2,8,8,1,0,9,7,5,6,6,5,9,3,3,4,4,6,1,2,8,4,7,5,6,4,8,2,3,3
,7,8,6,7,8,3,1,6,5,2,7,1,2,0,1,9,0,9,1,4,5,6,4,8,5,6,6,9,2,3,4,6,0,3,4,8,6,1,0,4,5,4,3,2,6,6,4,8,2,1,3,3,9,3,6,0,7,2,6,0,2,4,9,1,4,1,2,7,3,7,2,4,5,8,7,0,0,6
,6,0,6,3,1,5,5,8,8,1,7,4,8,8,1,5,2,0,9,2,0,9,6,2,8,2,9,2,5,4,0,9,1,7,1,5,3,6,4,3,6,7,8,9,2,5,9,0,3,6,0,0,1,1,3,3,0,5,3,0,5,4,8,8,2,0,4,6,6,5,2,1,3,8,4,1,4,6
,9,5,1,9,4,1,5,1,1,6,0,9,4,3,3,0,5,7,2,7,0,3,6,5,7,5,9,5,9,1,9,5,3,0,9,2,1,8,6,1,1,7,3,8,1,9,3,2,6,1,1,7,9,3,1,0,5,1,1,8,5,4,8,0,7,4,4,6,2,3,7,9,9,6,2,7,4,9
,5,6,7,3,5,1,8,8,5,7,5,2,7,2,4,8,9,1,2,2,7,9,3,8,1,8,3,0,1,1,9,4,9,1,2,9,8,3,3,6,7,3,3,6,2,4,4,0,6,5,6,6,4,3,0,8,6,0,2,1,3,9,4,9,4,6,3,9,5,2,2,4,7,3,7,1,9,0
,7,0,2,1,7,9,8,6,0,9,4,3,7,0,2,7,7,0,5,3,9,2,1,7,1,7,6,2,9,3,1,7,6,7,5,2,3,8,4,6,7,4,8,1,8,4,6,7,6,6,9,4,0,5,1,3,2,0,0,0,5,6,8,1,2,7,1,4,5,2,6,3,5,6,0,8,2,7
,7,8,5,7,7,1,3,4,2,7,5,7,7,8,9,6,0,9,1,7,3,6,3,7,1,7,8,7,2,1,4,6,8,4,4,0,9,0,1,2,2,4,9,5,3,4,3,0,1,4,6,5,4,9,5,8,5,3,7,1,0,5,0,7,9,2,2,7,9,6,8,9,2,5,8,9,2,3
,5,4,2,0,1,9,9,5,6,1,1,2,1,2,9,0,2,1,9,6,0,8,6,4,0,3,4,4,1,8,1,5,9,8,1,3,6,2,9,7,7,4,7,7,1,3,0,9,9,6,0,5,1,8,7,0,7,2,1,1,3,4,9,9,9,9,9,9,8,3,7,2,9,7,8,0,4,9
,9,5,1,0,5,9,7,3,1,7,3,2,8,1,6,0,9,6,3,1,8,5,9,5,0,2,4,4,5,9,4,5,5,3,4,6,9,0,8,3,0,2,6,4,2,5,2,2,3,0,8,2,5,3,3,4,4,6,8,5,0,3,5,2,6,1,9,3,1,1,8,8,1,7,1,0,1,0
,0,0,3,1,3,7,8,3,8,7,5,2,8,8,6,5,8,7,5,3,3,2,0,8,3,8,1,4,2,0,6,1,7,1,7,7,6,6,9,1,4,7,3,0,3,5,9,8,2,5,3,4,9,0,4,2,8,7,5,5,4,6,8,7,3,1,1,5,9,5,6,2,8,6,3,8,8,2
,3,5,3,7,8,7,5,9,3,7,5,1,9,5,7,7,8,1,8,5,7,7,8,0,5,3,2,1,7,1,2,2,6,8,0,6,6,1,3,0,0,1,9,2,7,8,7,6,6,1,1,1,9,5,9,0,9,2,1,6,4,2,0,1,9,8,9,3,8,0,9,5,2,5,7,2,0,1
,0,6,5,4,8,5,8,6,3,2,7,8,8,6,5,9,3,6,1,5,3,3,8,1,8,2,7,9,6,8,2,3,0,3,0,1,9,5,2,0,3,5,3,0,1,8,5,2,9,6,8,9,9,5,7,7,3,6,2,2,5,9,9,4,1,3,8,9,1,2,4,9,7,2,1,7,7,5
,2,8,3,4,7,9,1,3,1,5,1,5,5,7,4,8,5,7,2,4,2,4,5,4,1,5,0,6,9,5,9,5,0,8,2,9,5,3,3,1,1,6,8,6,1,7,2,7,8,5,5,8,8,9,0,7,5,0,9,8,3,8,1,7,5,4,6,3,7,4,6,4,9,3,9,3,1,9
,2,5,5,0,6,0,4,0,0,9,2,7,7,0,1,6,7,1,1,3,9,0,0,9,8,4,8,8,2,4,0,1,2,8,5,8,3,6,1,6,0,3,5,6,3,7,0,7,6,6,0,1,0,4,7,1,0,1,8,1,9,4,2,9,5,5,5,9,6,1,9,8,9,4,6,7,6,7
,8,3,7,4,4,9,4,4,8,2,5,5,3,7,9,7,7,4,7,2,6,8,4,7,1,0,4,0,4,7,5,3,4,6,4,6,2,0,8,0,4,6,6,8,4,2,5,9,0,6,9,4,9,1,2,9,3,3,1,3,6,7,7,0,2,8,9,8,9,1,5,2,1,0,4,7,5,2
,1,6,2,0,5,6,9,6,6,0,2,4,0,5,8,0,3,8,1,5,0,1,9,3,5,1,1,2,5,3,3,8,2,4,3,0,0,3,5,5,8,7,6,4,0,2,4,7,4,9,6,4,7,3,2,6,3,9,1,4,1,9,9,2,7,2,6,0,4,2,6,9,9,2,2,7,9,6
,7,8,2,3,5,4,7,8,1,6,3,6,0,0,9,3,4,1,7,2,1,6,4,1,2,1,9,9,2,4,5,8,6,3,1,5,0,3,0,2,8,6,1,8,2,9,7,4,5,5,5,7,0,6,7,4,9,8,3,8,5,0,5,4,9,4,5,8,8,5,8,6,9,2,6,9,9,5
,6,9,0,9,2,7,2,1,0,7,9,7,5,0,9,3,0,2,9,5,5,3,2,1,1,6,5,3,4,4,9,8,7,2,0,2,7,5,5,9,6,0,2,3,6,4,8,0,6,6,5,4,9,9,1,1,9,8,8,1,8,3,4,7,9,7,7,5,3,5,6,6,3,6,9,8,0,7
,4,2,6,5,4,2,5,2,7,8,6,2,5,5,1,8,1,8,4,1,7,5,7,4,6,7,2,8,9,0,9,7,7,7,7,2,7,9,3,8,0,0,0,8,1,6,4,7,0,6,0,0,1,6,1,4,5,2,4,9,1,9,2,1,7,3,2,1,7,2,1,4,7,7,2,3,5,0
,1,4,1,4,4,1,9,7,3,5,6,8,5,4,8,1,6,1,3,6,1,1,5,7,3,5,2,5,5,2,1,3,3,4,7,5,7,4,1,8,4,9,4,6,8,4,3,8,5,2,3,3,2,3,9,0,7,3,9,4,1,4,3,3,3,4,5,4,7,7,6,2,4,1,6,8,6,2
,5,1,8,9,8,3,5,6,9,4,8,5,5,6,2,0,9,9,2,1,9,2,2,2,1,8,4,2,7,2,5,5,0,2,5,4,2,5,6,8,8,7,6,7,1,7,9,0,4,9,4,6,0,1,6,5,3,4,6,6,8,0,4,9,8,8,6,2,7,2,3,2,7,9,1,7,8,6
,0,8,5,7,8,4,3,8,3,8,2,7,9,6,7,9,7,6,6,8,1,4,5,4,1,0,0,9,5,3,8,8,3,7,8,6,3,6,0,9,5,0,6,8,0,0,6,4,2,2,5,1,2,5,2,0,5,1,1,7,3,9,2,9,8,4,8,9,6,0,8,4,1,2,8,4,8,8
,6,2,6,9,4,5,6,0,4,2,4,1,9,6,5,2,8,5,0,2,2,2,1,0,6,6,1,1,8,6,3,0,6,7,4,4,2,7,8,6,2,2,0,3,9,1,9,4,9,4,5,0,4,7,1,2,3,7,1,3,7,8,6,9,6,0,9,5,6,3,6,4,3,7,1,9,1

};

int main(int argc, char* argv[])
{
int IN_SIZE = 0;
double* in;
fftw_complex *out;

IN_SIZE = sizeof(pi_fract)/sizeof(pi_fract[0]);
std::cout << "Size: " << IN_SIZE << std::endl;

// fill allocated array
in = (double*) fftw_malloc(sizeof(double) * IN_SIZE);
for( int i = 0; i < IN_SIZE; ++i )
{
in[i] = pi_fract[i];
}

// calculate average
double av = 0;
for( int i = 0; i < IN_SIZE; ++i )
{
av += in[i];
}
av /= (double)IN_SIZE;
std::cout << "Avg: " << av << std::endl;

// subtract average
for( int i = 0; i < IN_SIZE; ++i )
{
in[i] -= av;
}

// DFT transform
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * IN_SIZE);
fftw_plan p = fftw_plan_dft_r2c_1d(IN_SIZE, in, out, FFTW_ESTIMATE);
fftw_execute(p);

// print results
for( int i = 0; i < IN_SIZE/2+1; ++i )
{
double a = 2.0*sqrt(out[i][0]*out[i][0] + out[i][1]*out[i][1])/IN_SIZE;
std::cout << a << std::endl;
}

// cleanup
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);

return 0;
}