4000 Stellen von Pi mit ATtiny2313

ATtiny2313 im DIL-Gehäuse

Als neulich im Foren nach einem

"rechenintensiven Beispielprogramm für den kleinen ATtiny2313 in C oder ASM"

gesucht wurde und als Vorschlag ein

"Pi auf 1000 Stellen"

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, denn nicht alle zu berechnenden Stellen sind gleichzeitig im RAM des µC speicherbar: bereits berechnete Stellen sollen ausgegeben werden, und an der Ausgabe-Schnittstelle sollen Ziffernweise die Nachkommastellen von π heraus purzeln.

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 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 dabei nicht in stupidem Auswerten bis zur gewünschten Genauigkeit, sondern es wird die n-te Nachkommastelle von π damit bestimmt ohne vorherige oder nachfolgende Ziffern der Darstellung zu bestimmen. Der einzige Wermutstropfen ist, daß die Darstellung zur Basis 16 erhalten wird und nicht wie gewohnt zur Basis 10 — aber immerhin liefert die Formel überhaupt erst einen Weg, um die Berechnung auf einer kleinen Hardware 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 etwa in der Wikipedia

Die ganze Aufgabe besteht also nur noch darin, den Algorithmus hinzuschreiben und zu compilieren — und zu hoffen, dass das erzeugte Programm in den kleinen 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. Das Programm kann daher auch als Stresstest für die Tools verstanden werden.

Als Sprache wählte ich C99 und als Compiler avr-gcc 4.7 oder neuer, so daß ein Executable erzeugen wird, das in die 2k Programmspeicher eines ATtiny2313 hineinpasst.

Der 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, dass das Programm auch für einen PC übersetzt und dort ausgeführt 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 brauchbare 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 Programmspeicher des ATtiny2313 zu klein.
• Der Algorithmus berechnet wie gesagt Nachkommastellen der Hexadezimaldarstellung von π. 4000 hex-Nachkommastellen entsprechen übrigens ca. 4800 Dezimalstellen: log 16 / log 10 ≈ 6/5.
• Die Ausgabe lässt sich etwa vergleichen mit "First 8366 Hex Digits of PI"
• Den Quelltext gibt es auch hier als bereits fertig kompiliertes Atmel Studio 6.1 SP2 Project (inkl. Hex und ELF): Medium:4000Pi.zip (60KB)

/* pi.c: Compute the first 4000 hexadecimal digits of pi
         on an AVR ATtiny2313 microcontroller. */
/*
    Language: C99
    Compiler: An AVR toolchain based on avr-gcc 4.7 or newer and AVR-Libc
              if you like to compile for AVR ATtiny2313, or a C99-compliant
              compiler if you like to run this program 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 and floating point
              implementation, 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 is 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 is avr-gcc 4.7+ together with AVR-Libc.

                 avr-gcc pi.c -o pi.elf $(CFLAGS) $(OPT) $(DEFS) $(LDFLAGS)
    
              with the following abbreviations for convenience:
    
              CFLAGS  = -std=c99 -mmcu=attiny2313
              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 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
              convenient with baud rates.  Depending on your 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 this time by 16.

    The time to get the first n digits of pi with this method is roughly

        T(n) = K * n^2

    so that, given the time for the computation of n digits, you can easily
    compute K and estimate how long the computation will take for other
    values of n.  K will depend on the clock speed you use, for example.
    
    Some software metrics for the AVR binary compiled as above:
    
    Program Size:  Less than 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);

#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 that allows for
// larger moduli compared to (a * b) % n.
//
// The 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.
//
// Moreover, ATtiny2313 has no hardware multiplyer or divider, anyway.
// We have to cope with 2048 bytes of program memory and are happy
// to avoid dragging in library routines if possible.

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.  Return a value in (-2,2) that
// is equivalent to s mod 1.

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 an 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 routine that writes 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 the 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__