4000 Stellen von Pi mit ATtiny2313
Als neulich im Foren nach einem
- "rechenintensiven Beispielprogramm für den kleinen ATtiny2313 in C oder ASM"[1]
gesucht wurde und als Vorschlag ein
- "Pi auf 1000 Stellen"[2]
kam, wolle ich es genauer wissen. Auch wenn der Vorschlag eher als Scherz aufzufassen ist — sind C-Compiler inzwischen so gut, daß damit π auf tausende Stellen berechnet werden kann? Bei einem Assembler-Programm hätte ich da keine Bedenken; aber C?
Der AVR ATtiny2313 ist bekanntlich ein 8-Bit Mikrocontroller mit 2048 Bytes an Programmspeicher und 128 Bytes RAM. Die Aufgabe erfordert also eine spezielle Herangehensweise, wenn nicht alle zu berechnenden Stellen sind gleichzeitig im RAM des µC speicherber: bereits bereichnete Stellen sollen ausgegeben werden, und an der Ausgabe-Schnittstelle sollen Ziffernweise die Nachkomastellen von π herauspurzeln.
Die Formel
Das theoretisch-algorithmische Rüstzeug dafür ist lange bekannt und das Web ist voll davon: Eine 1995 von Simon Plouffe gefundene Reihendarstellung[3] für π:
- [math]\displaystyle{ \pi=\sum_{k=0}^\infty \frac1{16^k} \left(\frac4{8k+1}-\frac2{8k+4}-\frac{1}{8k+5}-\frac1{8k+6}\right) }[/math]
Die Anwendung der Darstellung besteht nun nicht in stupidem Auswerten bis die genwünschte Zahl an Nachkommstellen erreicht ist, sondern es wird die n-te Nachkommastelle von π damit bestimmt ohne vorherige oder nachfolgende Ziffern der Darstellung zu bestimmen. Der einzige Wehrmutstropfen ist, daß die Darstellung zue Basis 16 erhalten wird und nicht wie gewohnt zur Basis 10 — aber immerhin liefert die Formel überhaupt erst eine Weg, um die Berechnung auf einen kleinen µC mit 2k Programmspeicher ausführen zu können. Auch dieser Weg ist hinreichend oft im Internet erklärt — je nach Seite mehr oder weniger gut. Eine gute Erklärung findet sich in der Wikipedia[4]
Die ganze Aufgabe besteht also nur noch darin, den Algorithmus hinzuschreiben und zu compilieren — und zu hoffen, daß das erzeugte Programm in den kleines Programmspeicher eines ATtiny hineinpasst.
Gerade der letzte Punkt ist nicht selbstverständlich, denn die Anwendung der obigen Formel verlangt den Einsatz von Fließkomma-Arithmetik. Zwar ist auch eine Implementierung mittels Festkomma-Zahlen möglich, aber ich wollte mit den Bordmitteln der C-Sprache auskommen: Ohne Inline-Assembler Hacks, ohne umständliches Hin-und-Herwandeln und Skalieren, ohne Implementierung eigener Divisionroutinen, etc.
Als Sprache wählte ich C99 und als Compiler avr-gcc[5] 4.7. Tatsächlich ist das die einzige Version von avr-gcc, die im Zusammenspiel mit der AVR Libc[6] ein Executable erzeugen kann, das in die 2k Flash des ATtiny hineinpasst.
Quelltext: C99
Zur Implementierung ist hier nicht viel zu sagen; die Quellen sind ausführlich kommentiert. Nur noch einige Anmerkungen:
- Die Quellen sind so allgemein gehalten, daß das Programm auf für einen PC übersetzt und dort gestartet werden kann.
- Wird für AVR übersetzt, erfolgt die Ausgabe der berechneten Ziffern auf der UART-Schnittstelle. Für Compiler- und Linker-Optionen siehe die Quellkommentare.
- Bei einer Taktfrequenz von 20 MHz dauert die Berechnung der ersten 1000 hex-Ziffern etwa 3–4 Minuten. Verdoppelt man die Anzahl der Ziffern, ist die 4-fache Rechenzeit zu erwarten.
- Der Algorithmus liefert brauchbre Resultate für die ersten 4096 Nachkommastellen von π. Danach reichten die verwendeten 16-Bit int für ganze Zahlen nicht mehr aus, und für 32-Bit Zahlen ist der Speicher zu klein.
- Die Ausgabe lässt sich etwa vergleichen mit "First 8366 Hex Digits of PI"[7]
<c> /* pi.c: Compute the first 1000 hexadecimal digits of pi
on an AVR ATtiny2313 microcontroller. */
/*
Language: C99
Compiler: avr-gcc 4.7 if you like to compile for AVR ATtiny2313, or any
C99-compliant compiler if you like to run this pogram on a PC.
Hardware: AVR ATtiny2313 which is an 8-bit microcontroller with
128 Bytes of RAM and 2048 Bytes of program memory.
http://www.atmel.com/Images/doc2543.pdf (Manual, PDF, 4 MB)
http://www.atmel.com/Images/doc2543S.pdf (Summary, PDF, 500 kB)
Compile: You need a reasonably optimizing compiler, or otherwise
the program won't fit into the tiny silicon. Notice that
the code below is *vanilla* C with IEEE-754 floating point
arithmetic and without assembler or fixed point hacks.
The program in generic enough to run on a PC. If you use GCC,
just compile with the following command line and run pi.exe
gcc pi.c -o pi.exe -std=c99 -O2 -lm
The only compiler that produces code small enough for the
AVR ATtiny2313 hardware is avr-gcc 4.7
avr-gcc pi.c -o pi.elf $(CFLAGS) $(OPT) $(DEFS) $(LDFLAGS)
with the following abbreviations for convenience:
CFLAGS = -std=c99 -mmcu=attiny2313 -fno-dse
OPT = -Os -mcall-prologues -fno-split-wide-types
-fno-caller-saves -fno-tree-loop-optimize
LDFLAGS = -Wl,--relax -Wl,--gc-sections -lm
DEFS = -DF_CPU=22118400 -DBAUDRATE=115200 -DHOST_WINDOWS
For documentation of GCC see see http://gcc.gnu.org/onlinedocs
DEFS represents the UART setup. In my circuit I use a
22.1184 MHz quartz which is a bit of overclocking but is so
comfortable with baud rates. Depending on you configuration,
you will use other values for F_CPU and BAUDRATE, see below.
The first 1000 digits of pi take about 3 Minutes at 22 MHz.
For 4000 digits multiply the time by 16.
Some software metrics for the AVR binary compiled as above:
Program Size: 2000 of 2048 bytes
RAM Usage:
Static storage : 0 bytes
Static stack usage : ~ 50 bytes
Dynamic : 0 bytes
Total : 40% of 128 bytes
Coding Style: Stroustrup Indentation: 4 Space
- /
- include <stdlib.h>
- include <stdint.h>
- include <math.h>
// Factor out some platform/compiler dependencies
- if defined (__AVR__) && defined (__GNUC__)
// Running on AVR
- include <avr/io.h>
- include <avr/wdt.h>
- define cput(C) uart_putc(C)
static void uart_init (void); static void uart_putc (const char); int __attribute__((OS_main)) main (void);
- else // ! avr-gcc
// Running on PC
- include <stdio.h>
- define wdt_reset() (void) 0
- define uart_init() (void) 0
static void cput (char c) {
fputc (c, stdout); fflush (stdout);
}
- endif // avr-gcc
// return a * b mod n // // 0 <= a < n // 0 <= b < n // // This is a school-book implementation of multiplication. // // First, ATtiny2313 has no hardware multiplyer or divider, anyway. // We have to cope with 2048 bytes of program memory and // thus avoid dragging in library routines if possible. // // Second, this implementation widens the range of valid moduli from // \sqrt{1+UINT_MAX} to (1+UINT_MAX)/2. Or vice versa, it allows us // to use a smaller type for the modulus -- 16 bits in the AVR case.
static unsigned mod_mul (unsigned a, unsigned b, unsigned n) {
unsigned ab = 0;
while (1)
{
if (a % 2 == 1)
{
ab += b;
if (ab >= n)
ab -= n;
}
a = a / 2;
if (a == 0)
return ab;
b = b * 2;
if (b >= n)
b -= n;
}
}
// Return 16^k mod n // // Exponentiation with the same approach as above except // that we binary-expand the exponent instead of a factor.
static unsigned mod_pow16 (unsigned k, unsigned n) {
unsigned p = 1;
unsigned _16 = 16;
if (n == 1)
return 0;
while (_16 >= n)
_16 -= n;
while (1)
{
if (k % 2 == 1)
p = mod_mul (_16, p, n);
k = k / 2;
if (k == 0)
break;
_16 = mod_mul (_16, _16, n);
}
return p;
}
// Helper for the function below.
static float tame (float s) {
int8_t si = (int8_t) lrintf (s);
if (si <= -2 || si >= 2)
s -= si;
return s;
}
// The finite part mod 1 of sigma_j, i.e. partial sum where the exponent // of 16 is >= 0. By "mod 1" we always mean "up to some integer", // i.e. the result needs not to be normalized to [0,1).
static float sigma_a (unsigned n, uint8_t j) {
float s = 0.0f;
for (unsigned k = n-1; k+1 != 0; k--)
{
unsigned j_8k = j + 8*k;
s += mod_pow16 (n-k, j_8k) / (float) j_8k;
// Cut down the sum and don't let it grow too big.
// The bigger the number grows the less precision is
// left for the fractional part we are interested in.
s = tame (s);
wdt_reset();
}
return s;
}
- define MARGIN 10
static float sigma_b (unsigned n, uint8_t j) {
float s = 0;
float _16 = 1.0f;
for (unsigned k = n; k <= n + MARGIN; k++)
{
s += _16 / (8*k + j);
_16 /= 16;
wdt_reset();
}
return s;
}
// Compute an approximation of 16^n * sigma(j) mod 1 // where // // sigma_j = \sum_0^oo 1 / (16^k * (8k + j)) // // The sum is split into two parts: // sigma_a is the finite sum up to n. // sigma_b is the finite sum from n+1 to oo // and approximated by a sum from n+1 to n+MARGIN
float sigma (unsigned n, uint8_t j) {
return sigma_a (n, j) + sigma_b (n, j);
}
// Compute pi * 16^n up to some integer // using a Bailey-Borwein-Plouffe formula for pi: // // pi = 4*sigma_1 - 2*sigma_4 - sigma_5 - sigma_6 // // with the definition of sigma_j from above. All this // is explained very nicely in the French wikipedia at // http://fr.wikipedia.org/wiki/Formule_BBP // // For a proof define the power series // // sigma_j (x) = \sum x^{8k} / (8k + j) // // write the sum as integral and evaluate it at // x = sqrt(1/2), see http://www.pi314.net/eng/plouffe.php
float pi_n (unsigned n) {
float s = 0.0;
for (uint8_t i = 0; i < 4; i++)
{
// j[i] = 1, 4, 5, 6
uint8_t j = i ? i + 3 : 1;
// c[i] = 4, -2, -1, -1
int8_t c = -1;
if (i == 0) c = 4;
if (i == 1) c = -2;
s += c * sigma (n, j);
}
return s;
}
// We computed pi_n = 16^n * pi mod 1 // Get the first fractional hexadecimal digit by multiplying // with 16 and extracting digit 0 of the result. Voila!
static uint8_t pi_dig16 (unsigned n) {
return 15 & lrintf (floorf (16 * pi_n (n)));
}
// Map 0 <= n < 16 to its hexadecimal ASCII digit // '0', '1', ... 'F'
static uint8_t hexdigit (uint8_t n) {
n += '0'; return n > '9' ? n + 'A'-'0'-10 : n;
}
int main (void) {
uart_init();
cput ('\n');
cput ('3');
cput ('.');
// As explained above, 16-bit integers limitate us to moduli <= 2^15.
// The biggest modulus for n is 8n+6 so that for n >= 4096 we expect
// garbage from the implementation if 16-bit integers are used like
// with avr-gcc. In fact, we get garbage for n > 4100.
// It's not exacly 4095 because of the denominators in sigma_a that
// delay the garbage for some values of n.
// Easy going 4000 hex-digits of pi.
for (unsigned n = 0; n < 4000; n++)
{
// Print a line break after every 64 digits. That way the output
// can easily be compared with, e.g. the hexadecimal representation
// of pi from blowfish listed in "First 8366 Hex Digits of PI" from
// http://www.herongyang.com/Cryptography/
if (n % 64 == 0)
cput ('\n');
// Get the n+1-th hexadecimal digit of pi and
// output it as ASCII character.
cput (hexdigit (pi_dig16 (n)));
}
cput ('\n');
return 0;
}
/************************************************************************/ // We are running on ATtiny2313 bare metal, of course. // Provide some minimalist output routines that write to UART. /************************************************************************/
- ifdef __AVR__
// !!! You must have set the fuses appropriately to run the !!! // !!! controller at desired F_CPU. Defining F_CPU is needed !!! // !!! to get correct values to set up UART. !!! // !!! !!! // !!! DEFINING F_CPU WILL NOT CHANGE THE FREQUENCY! !!!
void uart_init (void) {
// Define F_CPU and BAUDRATE depending on your hardware setup. // For my setup, it's the followin defines in avr-gcc's command line: // -DF_CPU=22118400 -DBAUDRATE=115200
unsigned ubrr = -.6 + F_CPU / (8L * BAUDRATE);
UBRRH = ubrr >> 8; UBRRL = ubrr;
// Enable UART Transmitter, data mode 8N1, asynchronous
UCSRA = (1 << U2X) | (1 << TXC); UCSRB = (1 << TXEN); UCSRC = (1 << UCSZ1) | (1 << UCSZ0);
}
/////////////////////////////////////////////////////////////////
// Write one char to UART (non-buffered, blocking version)
void uart_putc (char c) {
while (!(UCSRA & (1 << UDRE)))
wdt_reset();
UDR = c;
- ifdef HOST_WINDOWS
if (c == '\n')
uart_putc ('\r');
- endif // HOST_WINDOWS
}
void exit (int x) {
(void) x;
while (1)
wdt_reset();
}
- endif // __AVR__
</c>
Weblinks
- ↑ Forum: AVR Atiny2313: "Suche rechenintensivs Beispielprogramm für ATtiny2313"
- ↑ Forum: AVR Atiny2313: "Pi auf 1000 Stellen"
- ↑ Für einen Beweis siehe The World of π: Proof: Formula BBP.
- ↑ Wikipédia: Formule BBP: Exploitation de la formule pour calculer les chiffres après la virgule de π
- ↑ GCC: The GNU Compiler Collection
- ↑ AVR Libc: Home Page
- ↑ www.herongyang.com: First 8366 Hex Digits of PI