The Shortest Implementations of π

As I have established (see here and in [1]), this permutation actually has a very strong and very simple structure that is completely captured by the following system of equations:

where \(\alpha\) is a root of \(X^8+X^4+X^3+X^2+1\) so that it generates the multiplicative subgroup of the finite field with 256 elements. The functions \(s\) is a permutation of \(\{ 0, 1,...,14 \}\) defined in Table 1.

The function \(\kappa\) is maps \(\mathbb{F}_2^4\) to \(\mathbb{F}_2^8\) and is defined as \(\kappa(x) = \Lambda(x) \oplus 252\), where \(\Lambda\) is a linear function defined as follows:

If you are not too familiar with finite field arithmetic, a simple and readable C implementation is provided below.

Our best C implementation (see below) can be rearranged as follows using only newlines and spaces.

p(x){
  unsigned char
      *t="@`rFTDVbpPBvdtfR@\xacp?\xe2>4\xa6\xe9{z\xe3q5\xa7\xe8",
      a=2,l=0,b=17;
  while(x && (l++,a^x)) a=2*a^a/128*29;
  return l%b ? t[l%b]^t[b+l/b]^b : t[l/b]^188;
}

The newlines and spaces allow us to better see the logic.

• t stores the precomputation of both \(i \mapsto \kappa(16-i)\) and \(j \mapsto (\alpha^{17})^{s(j)}\) up to some constants.
• The variable \(b\) merely serves as a shorthand for the constant \(17\).
• The while loop computes the discrete log of x using a counter (l) that is incremented each time the variable a, which is initialized to \(\alpha\), is multiplied by \(\alpha\). Said multiplication is implemented as a Galois LFSR.
• The final return value uses a ternary operator to implement the logic of the if statement that is used to define the TKlog.

An even more readable implementation works as follows.

unsigned char p(unsigned char x){
    unsigned char
        s[]={1,221,146,79,147,153,11,68,214,215,78,220,152,10,69},
        k[]={0,32,50,6,20,4,22,34,48,16,2,54,36,52,38,18,0};    
    if(x) {
        unsigned char l=1, a=2;
        while(a!=x) {
            a=(a<<1)^(a>>7)*29;
            l++;
        }
        if (l%17) return 252^k[l%17]^s[l/17];
        else return 252^k[l/17];
    }
    else return 252;
}

In this implementation:

• s stores \((\alpha^{17})^{s(j)}\) for all \(0 \leq j < 15\),
• k stores the value of \(\kappa(16-i) \oplus \kappa(0)\) for all \(0 \leq i < 17\),
• (a<<1)^(a>>7)*29 is the multiplication by \(\alpha\) in the finite field, and
• 252 is \(\kappa(0)\).

Brian Shaft (832 bits)

The following program can be compiled into a 832 bits binary file using the Algol compiler for the Unisys MCP operating system.

Real    Procedure P(X);
   Value X;
   Integer X;
Begin
   Value Array S(1,221,146,79,147,153,11,68,214,215,78,220,152,10,69),
               K(0,32,50,6,20,4,22,34,48,16,2,54,36,52,38,18,0);
   Integer L,A;
   A:= 2;
   If X ^= 0 Then
   Begin
      For L:= 1 Step 1 While A ^= X Do
         A:=(Real(Boolean(A+A).[7:8] eqv ^Boolean(A.[7:1]*29))).[7:8];
      P:= Real((If L mod 17 = 0 Then
                   Boolean(K[L div 17])
                Else
                   Boolean(K[L mod 17]) eqv ^Boolean(S[L div 17]) eqv ^Boolean(252));
   End
   Else
      P:= 252;
End;

Xavier Bonnetain & Gaëtan Leurent (640 bits)

The machine code corresponding to this implementation is 80 bytes long and thus fits on 640 bits. It corresponds to a complete 80 bytes long program, with 32 bytes of data and 20 instructions (16 16-bit instructions and 4 32-bit instructions). More details are given in [2]. Xavier's webpage; Gaëtan's webpage.

        .syntax unified
        .text
        .align  4

        .global p
p:
        movs    r3, #0
        cbz     r0, .L2
        mov     r2, #0x1000000
        lsls    r0, #24
.L3:
        lsls    r2, r2, #1
        it      cs
        eorcs   r2, r2, #0x1d000000
        adds    r3, r3, #1
        cmp     r2, r0
        bne     .L3
.L2:
        movs    r0, #17
        udiv    r2, r3, r0
        mls     r3, r0, r2, r3
        adr     r1, .table_k
        cbz     r3, .LL
        ldrb    r3, [r1, r3]
        add     r1, r0
.LL:
        ldrb    r0, [r1, r2]
        eors    r0, r0, r3
        bx      lr
        .align  2

.table_k:
        .byte   -4
        .byte   -36
        .byte   -50
        .byte   -6
        .byte   -24
        .byte   -8
        .byte   -22
        .byte   -34
        .byte   -52
        .byte   -20
        .byte   -2
        .byte   -54
        .byte   -40
        .byte   -56
        .byte   -38
        .byte   -18
        .byte   -4
.table_s:
        .byte   1
        .byte   -35
        .byte   -110
        .byte   79
        .byte   -109
        .byte   -103
        .byte   11
        .byte   68
        .byte   -42
        .byte   -41
        .byte   78
        .byte   -36
        .byte   -104
        .byte   10
        .byte   69

The programs here may not be perfectly standard (see for instance the omission of int in function prototypes) but they should be successfully processed by your compiler.

Xavier Bonnetain (1106 bits)

This implementation is 158 ASCII characters long and thus fits on 1106 bits. It simplifies the previous implementation by removing the variable a containing the generator raised to power l. Instead, this programs looks for l such that \(x^{255-\ell} = 1\). Xavier's webpage.

p(x){unsigned char*k="@`rFTDVbpPBvdtfR@\xacp?\xe2>4\xa6\xe9{z\xe3q5\xa7\xe8",l=0,b=17;while(--l&&x^1)x=2*x^x/128*285;return l%b?k[l%b]^k[b+l/b]^b:k[l/b]^188;}

Xavier Bonnetain (1127 bits)

This implementation is 161 ASCII characters long and thus fits on 1127 bits. Xavier's webpage.

p(x){unsigned char*t="@`rFTDVbpPBvdtfR@\xacp?\xe2>4\xa6\xe9{z\xe3q5\xa7\xe8",a=2,l=0,b=17;while(x&&(l++,a^x))a=2*a^a/128*29;return l%b?t[l%b]^t[b+l/b]^b:t[l/b]^188;}

In this category, I put C programs that can be compiled by some compilers and successfully evaluate pi. The functions in this part are not portable or even standard C—Lovecraftian horrors really.

Xavier Bonnetain (1071 bits)

This implementation is 153 ASCII characters long and thus fits on 1071 bits. Xavier's webpage.

p(x){char*k="@`rFTDVbpPBvdtfR@\xacp?\xe2>4\xa6\xe9{z\xe3q5\xa7\xe8";int l=256,b=17;while(--l*x^l)x=2*x^x/128*285;return l%b?k[l%b]^k[b+l/b]^b:k[l/b]^188;}

Léo Perrin (1351 bits)

This implementation is 193 ASCII characters long, meaning that it fits on 1351 bits. It is based on Xavier's C implementation. Note that it wouldn't work with Python2.

def p(x):
 t=b'@`rFTDVbpPBvdtfR@\xacp?\xe2>4\xa6\xe9{z\xe3q5\xa7\xe8'
 l,d=255,17
 if x<2:return(252,238)[x]
 while 1^x:x,l=(x<<1)^(x>>7)*285,l-1
 return(188^t[l//d],d^t[l%d]^t[d+l//d])[l%d>0]