Today (14th March 2022)1 is special not only because it’s Pi Day, but also because I have a new piece of hardware to play with. Specifically, I have a small development board with an AVR AT90USB162 microcontroller on it. The obvious path to take is to use it to calculate pi.

The finished product The finished product. Not particularly pretty but definitely interesting.

This is interesting for me because don’t know much C and I’m very new to this sort of microcontroller programming. Naturally, that means my code won’t always be optimal, tidy or follow best practices.

The Calculation

The method I chose to implement was the Bailey-Borwein-Plouffe formula. It’s a hexadecimal spigot algorithm for the digits of pi. This means that it can produce the n-th (base 16) digit of pi without needing to calculate any others. I decided that this would work well for my simple microcontoller because I could display each digit in binary with four LEDs attached to the limited output pins on the board.

The exact formula for pi is

BBP formula

but we care about a modified version of it which extracts a specific digit.

BBP digit formula

BBP digit formula subsidiary

(Equations found here)

The Implementation

The first of the two formulas can be translated fairly directly into C. We’ll define intFloor later as the floor operation so that we can just return the actual digit.

//Return the nth (0-based) hex decimal place of pi
int nthDigit(int n) {
	return (unsigned int)intFloor(16 * (4 * bbpSum(n, 1) - 2 * bbpSum(n, 4)
		- bbpSum(n, 5) - bbpSum(n, 6))) % 16;
}

The second one takes a little more thought and a little more code but, at its core, is just two loops to sum the required terms. The first one simply iterates over values of k and does some modular exponentation.

The other loop needs to keep track of a few extra values. I decided to create variables for the denominator and the reciprocal of the numerator (since this is always an integer in this summation). We also need a constant ‘epsilon’ value, which will serve as the minimum size of terms before we consider them insignificant (since the sum would be infinite otherwise).

//Calculate a summation from the BBP digit extraction formula (mod 1)
double bbpSum(int n, int j) {
	double s = 0;
	for(int k = 0; k <= n; k++) {
		s += (double)binModExp(16, n - k, 8 * k + j) / (8 * k + j);
		s -= intFloor(s);
	}
	long numRecip = 16;
	int denom = 8 * (n+1) + j;
	double eps = 1e-6;
	double frac;
	while((frac = 1.0 / (numRecip * denom)) >= eps) {
		s += frac;
		numRecip *= 16;
		denom += 8;
	}
	return s;
}

Now, we need to define the subsidiary functions that we used earlier. Firstly, intFloor should round a double down to the nearest int below it. This can be achieved by ‘truncating’ it (take the integer part), and subtracting one if the number was negative (since truncation rounds towards zero).

//Return largest int less than or equal to x
int intFloor(double x) {
	int f = (int)x;	//Truncate x
	if((double)f == x || x >= 0)
		return f;	//Truncation should have produced correct value
	else
		return f - 1;	//Adjust because truncation rounded the wrong way
}

The final function we need is modular exponentation. This can be implemented naively by repeatedly multiplying and modulo-ing but this is too slow given the number of times that it will be called. Instead binModExp implements binary exponentation, being careful to use long multiplication to avoid any overflows.

//Calculate a^b mod m by binary exponentiation
unsigned int binModExp(unsigned int a, unsigned int b, unsigned int m) {
	int res = 1;
	while(b > 0) {
		if(b % 2 == 1)
			res = ((unsigned long)res * a) % m;	//Use long to avoid overflow
		a = ((unsigned long)a * a) % m;
		b >>= 1;
	}
	return res;
}

Electronics and AVR Code

My crude breadboard setup is shown below. It’s nothing special but it’s good enough for this project. And before you ask, the PCB pins are connected to the breadboard by small pieces of bent wire. It’s the best solution I could come up with since I didn’t feel like soldering anything to them for this quick project. They don’t make very consistent contact and frequently need jiggling to make sure the LEDs light up.

The electronics The breadboard. Note that the least significant bit is on the left.

I have four LEDs attached to the pins 0-3 of port D to display each digit. To light these up appropriately, I defined updateDigit which (if n has changed), calculates the nth digit and writes it to PORTD. It also turns the in-built LED (pin 4) on during the calculation to indicate that something is happening.

void updateDigit() {
	if(oldN == n)
		return;

	oldN = n;

	PORTD &= ~0xf;

	if(oldN == 0) {	//Display 3 as the first digit
		PORTD |= 0x3;
		return;
	}

	//Calculate nth digit, turning in-built LED on during calculation
	PORTD |= (1 << PD4);
	int digit = nthDigit(oldN - 1);
	PORTD &= ~(1 << PD4);
	//Output digit (mod 16) onto PORTD pins 0-3
	PORTD |= digit & 0xF;
}

I then decided to be fancy (read: pretty basic) and define an interrupt handler to increment n when the in-built button (pin 7) is pressed. Since n is used both in an interrupt handler and the main loop, I defined it as volatile (the Internet said this was a good idea).

volatile unsigned int n = 0;
ISR(INT7_vect) {
	n++;
}

Finally, I wrote some somewhat boilerplate setup code in main,

DDRD |= 0x1F;	//Set PORTD pins 0-4 as output
EICRB |= (1 << ISC71) | (1 << ISC70);	//Do interrupt on button rising edge
EIMSK |= (1 << INT7);	//Enable interrupts for button
sei();	//Global enable interrupts

followed by an infinite loop which repeatedly updates the displayed digit and then tells the processor to go into a light sleep (idle sleep mode) so that it consumes less power until the next interrupt occurs.

while(1) {
	updateDigit();

	//Enter idle to wait for next interrupt
	set_sleep_mode(SLEEP_MODE_IDLE);
	sleep_enable();
	sleep_cpu();
	sleep_disable();
}

Verifying and Correcting the Digits

Testing with USB and LUFA

With all of this code written, I plugged compiled the code, plugged in the board and flashed the code.2 I reset the board and pushed the button a few times with tentative delight as the correct initial digits showed up on the LEDs. The code clearly works in a basic capacity.

The first few correct digits 3.243f. Not quite as well known as the famous 3.1415.

But how far can it go? At what point will it go wrong?

I found 8336 hex digits of pi online to compare mine against but, rather than tediously comparing these by hand (although I did this more than I care to admit), I decided to automate and human error-proof it. This lead me down a rabbit hole with LUFA, an AVR library for USB communication, that is unfortunately beyond the scope of this post. The code is on Github for anyone who’s interested but mostly consists of LUFA’s VirtualSerial example.

Having used this to output digits from the microcontroller to a file on my laptop, the result was that the digits are correct until the 709th, which is a 7 instead of a 6. The discrepencies only become more frequent after that.

avr-gcc and Portability

Strangely, the same code on my laptop generates these digits correctly (and, unsurprisingly, far quicker). It turns out that avr-gcc (used to compile the code for the microcontroller) defines some data types with a different number of bits. Specifically, int, long and double are all half the size. I managed to make these consistent by specifying int16_t and int32_t for int and long and float instead of double (these were the sizes being used by avr-gcc anyway).

Fixing Overflows and Rounding Errors

With the code now portable between my laptop and the microcontroller, I ran some experiments on my laptop with larger data types in various places and determined that the primary source of error was s (a float) in bbpSum. Presumably, when calculating later digits, the value gets too large for the fractional part (the bit we care about) to be stored precisely enough. To keep the value small (and therefore precise), I simply subtracted the integer part on every iteration of the first bbpSum loop using the intFloor function.

for(int16_t k = 0; k <= n; k++) {
	s += (float)binModExp(16, n - k, 8 * k + j) / (8 * k + j);
	s -= intFloor(s);
}

On my laptop, this generated 4133 correct digits. This is suspiciously close to 4096 (212). Further investigation revealed that this is the point at which a couple of expressions (8k+j and 8(n+1)+j) overflow their 16 bits. Casting and storing these values as uint32_t instead of int16_t fixes this and allows for an impressive 8198 correct digits with only one mistake at position 4329.

I investigated a bit more but the only easy improvement I could find was changing float back to double, which fixed digit 4329. Unfortunately, this isn’t portable to the microcontroller because avr-gcc defines double as an alias for float. Regardless, 8000 digits with one mistake is good enough for me given that you’d have to press the button thousands of times to find a mistake.

Testing the final code the microcontroller should produce identical results, given the measures we took to ensure portability. To somewhat verify this, I flashed the code and left it running overnight. It managed to calculate 2304 completely correct digits.3

Conclusion

This project touched on a variety of areas and taught me a surprising amount. I’m really happy with the result, even if I finished it a week late. The code is on Github, in case you have a compatible AVR board and want to try it yourself.

  1. Technicially, this is no longer today’s date since I’m a week late on this. The project took longer than I anticipated and I kept wanting to do more with it. ↩

  2. This is obviously a huge simplification. A lot of testing, debugging and recompiling happened to get to this point. ↩

  3. There is a chance that generating the digits by pressing the button would uncover some unrelated error involving e.g. interrupts but I don’t feel like testing that. ↩