Matt Parker ( @standupmaths ) has the opportunity to calculate Pi on a lander/rover on the Moon with Astrobotic! He has promised NOT to use some of his trademarked ‘Terrible Python Code’ TM and since I’ve actually written spaceflight software I though of a great way to do this…with fixed point arithmetic!!
You see, you don’t actually want to use floating point on a spacecraft processor if at all possible. And the way this estimation of Pi will be calculated is using a ratio of the random numbers inside a unit circle vs. outside the unit circle but inside a unit square. Since the area of a circle is Pi*r2 and r==1, the area is Pi. And the square surrounding the circle is 2 units wide so the area is 2*2=4.
But really we can just do the +x/+y quadrant since it’s all symmetric. Then we take the ratio of inside the circle to the entire area, which if we do the math is Pi/4 for the circle and 1×1==1 for the square – very convenient! So our ratio of ‘inside’ to ‘all’ x4 gives us our estimation of Pi.
So 2 key problems:
- We need random numbers to do this estimation
- We need to not use floating point calculations if at all possible, and in fact need to minimize ALL calculations…
OK, first the random numbers. Matt’s idea is to use data from the rover itself so do these calculations. So how to get random numbers from spacecraft telemetry? Well there are several options.
First the spacecraft processor will – somewhere – have an analog to digital converter (ADC). The ADC is likely 12 bits, giving values from 0..4095 (212-1). These are used to measure many things on a spacecraft including voltages, currents, temperatures, magnetic fields, etc. and we know from experience there will be a little bit of noise in them. Voltages tend to be pretty stable, so not a great choice. Currents (amps) will, many times, show periodic behavior as instruments, motors, and other things like that tend to run in cycles, so I think currents are also not a great choice.
For my money, I’m betting temperatures will be the right answer, as they not only tend to move around a bit, but since you are measuring resistance across a thermistor (a resistor that changes resistance with temperature) it tends to be a bit noisy – by at least a couple of counts, maybe +/- 5 or 6 counts or if you look at the lower bits of the 12 bit measurement maybe the bottom 3 bits (0 to 7 counts total). Concatenate some of these measurements and maybe move the bits around and I think you can create a good random number generator up to 15 bits pretty easily. By using different SHIFT and OR operations you can probably create several decent random numbers from a couple of measurements.
And that brings us to:
Using fixed point calculations instead of floating point. In fact we can just do this calculation with integer math and no floating point at all!
If we treat our new 15 bit random numbers (above) as the fractional portion of a value from 0..1, then adding two of them together will produce a value of at max 16 bits with a value from 0..2. Multiplying two of these numbers together will generate a new value of 0..1 but now with up to 30 bits – better than 1 part in a Billion!
So how to calculate the distance a 2 dimensional point is from the origin? This will allow us to answer the question of whether or not it falls inside the unit circle. It’s the same calculation as the hypotenuse of a triangle: square root of x2 + y2. “But Wait!” you say…”There’s a square root calculations there! Doesn’t that need floating point?” You’re right – it does require floating point and quite a bit of it!
BUT: we don’t need it! The magic of the unit circle (radius 1) is that the square root of a number >= 1.0 is still greater than 1.0. And the square root of a number <1.0 is still less than 1.0 ! So while, yes, technically, the equation requires the square root, we don’t need to know the actual value! Just whether the result of the square root would be inside or outside the circle IF we took its square root!
So I bashed together some not-entirely-terrible C++ code to show this actually works:
// Calculation of PI using only fixed point arithmetic.
// Matt Parker's "Moon Pi" project
// this would really run as a subroutine to which the main
// spacecraft code would pass a value presumably the result
// of an analog-to-digital conversion on something like a temperature
// We could then build a 15 bit value using the 'noise' in the measurement
// likely the bottom 2-3 bits since they're generally noisy.
// Take the bottom 3 bits of 5 different measurements and SHIFT/OR them
// together to make a random 15 bit value (I'm using rand() here for simplicity)
// You can probably make multiple different randomized value by using different
// SHIFT/OR patterns to build multiple number...
// Emory Stagmer aka VAXHeadroom 29-June-2025
#include <iostream>
// I know C++ has a PI constant somewhere, but I was too bored trying to find it, so here's a
// 'pretty' accurate define for comparing the result (good enough for what we're doing here)
const double pi=3.14159265358979323846;
using namespace std;
// a million cycles gives a pretty good approximation - within 0.001 or so...
#define RND_PNTS 1000000
int main()
{
int i;
unsigned long int tmpx, tmpy, tmpsum;
int count_in=0, count_out=0, count_ones=0;
float ratio;
srand(time(NULL));
cout << "Pi Fixed Point Calculator! Using " << RND_PNTS << " points." << endl;
for ( i=0; i<RND_PNTS; i++ )
{
// calc a 16 bit number where the low order 15 bits are the fraction and 0x8000 == 1.0
tmpx = ((rand() + rand()) & 0x7fff) +1;
// squaring the 15 bit fractional number gives at most a 30 fractional bit number
// where if bits 30 (or 31) are 1 then this overflowed to >1.0
tmpx = tmpx * tmpx;
// same for Y
tmpy = ((rand() + rand()) & 0x7fff) +1;
tmpy = tmpy * tmpy;
// add X+Y
tmpsum = tmpx + tmpy;
// if this overflows, then the answer is > 1.0
// The ACTUAL equation is sqrt( x^2 + y^2 )...
// BUT: Since the radius is 1.0 and we only care about
// whether it's inside or outside the circle we don't need to
// run the sqrt calc since sqrt( any number >=1.0 ) >= 1.0 and
// sqrt( any number <1.0 ) < 1.0
// so we don't care about the actual value of the sqrt only whether
// it WOULD have been >= 1.0 had we calculated it.
// We might do something fancy here if the value is actually == 1.0 (0x40000000)
// since that's actually a 1 in a billion chance (or a little more...2^30)
// and actually represents a loss of precision i.e. we don't know if it would have
// fallen inside or outside the circle if we had more bits...
if(tmpsum & 0xC0000000)
{
count_out++;
}
else
{
count_in++;
}
if( (tmpsum & 0x3FFFFFFF) == 0)
{
count_ones++;
}
}
// The rest of this code would be run on the ground after the fact and DOES include a floating point calc
cout << "Out = " << count_out << endl;
cout << "In = " << count_in << endl;
cout << "Ones = " << count_ones << endl;
ratio = (double)count_in / (double)(count_out+count_in);
cout << "Ratio = " << ratio << endl;
cout << "4xRatio = " << ratio*4.0 << endl;
cout << "Error = " << pi - (ratio*4.0) << endl;
}
So that’s how I’d calculate Pi on the moon (or any spacecraft!).
YMMV…
VAX's Headroom 
