Some finite field arithmetic; a simple Swift implementation

The title of this post started out as  “A Swift implementation of binary arithmetic in GF(2)^n” but I decided against it because, even though it’s more accurate, it defeats the purpose of the post which is to (hopefully) explain the principles in-code and in human-readable form…

For your perusal; here is some background, references, and more formal details;

http://xcore.github.io/doc_tips_and_tricks/crc.html

https://sites.google.com/site/ctxtree/crc/crc-binarydivision

https://en.wikipedia.org/wiki/Modular_arithmetic

http://jameskbeard.com/Temple/Data/Binary_Polynomial_Division.pdf

Now, to divide two binary numbers, modulo 2, you use the same technique as “long division” but I’ve not done long division by hand since a very long time and quite frankly can’t remember much of it so I am not going to use that excuse and instead present you directly with an algorithm, firstly in pseudo-code, then in Swift;

divideBinaryModulo2(N, D) where N>D
 let divisor = msb of D aligned with the msb of N by 
               shifting it up k bits
 let dividend = N
 let quotient = 0
 let remainder = 0
 for i = k-1..0 {
   if high bit of dividend is set
     set i’th bit of quotient to 1
     dividend = dividend XOR divisor

   shift dividend up 1 bit
 }
 remainder = dividend shifted back down k bits

The result is the quotient and remainder such that

quotient * D + remainder = N

Addition and subtraction is the same, modulo 2, and is implemented with the XOR operation so another way of writing this is

quotient * D XOR remainder = N

Where multiplication of A and B, modulo 2, is done by summing A multiplied by each term of B, i.e;

multiplyBinaryModulo2(A,B)
 result = 0 
 for each set bit in B {   
   let i = bit position
   result = result XOR (A shifted up by i)
}

Note; For the multiplication to be correct we also need to ensure that the result doesn’t go beyond 32 bits and the way to do this is to do the multiplication modulo a primitive polynomial of order 32 (i.e. one that uses 33 bits.)

Imagine we wanted to work within only 8- of our 32 -bits. Since addition is an XOR operation, assuming the terms are of order less than 8, it can never bring a result to overflow. Multiplication, however, could easily do that and the question then is; given the result of a polynomial multiply which overflows our 8 bits, what do we do…?

One approach might be to just AND it by 0x7f (0b1111111) to mask out the upper bits. This certainly works but it’s not going to give us the right results because, firstly, it’s a logical operation (not an arithmetic one that’s either a + a – or a *, which is required for this to remain inside the confines of a “ring” – which is a bit like a mathematical “group” a little more structure, like the concept of an “inverse”), secondly it’s not guaranteed to give a non-0 result and 0 isn’t a polynomial (it’s just a constant) and we don’t want it to show up in our results. Fixing the overflow also needs to give us a result which can be any of the possible polynomials we can fit in the 8 bits; i.e. it shouldn’t rule any out and thereby reduce our set of polynomials.

To the rescue; prime numbers. Or, to be more specific, prime polynomials. Something is prime if it’s not a product of other parts so a prime polynomial, or a primitive-, and irreducible  -polynomial can’t be written out as the product of other, lower order, polynomials. Another important property of prime polynomials (or numbers) is that they guarantee that a “multiplicative inverse” exist in a collection where we define multiplication as including a modulo on the result; i.e. a*b = a*b mod p, where p is the prime (number or polynomial).

To be more specific; since a prime number can’t be written as the product of two or more, smaller, numbers the product a*b can never be a multiple of that (or any other) prime number. This means that taking the modulo will never result in 0, since you only get 0 if the remainder is 0, i.e. if a*b is a multiple of p…The same logic applies to p being a prime (irreducible, primitive) polynomial.

So; to make our multiplications correct and mathematically consistent for polynomials, and to fit inside our example of 8 bits we need to follow every multiplication with a modulo (i.e. remainder) operation on the result by a prime polynomial. This prime polynomial can be picked relatively arbitrarily and there are tables for these (they can also be found using the same techniques used to find prime numbers. I’ve tried this with the Sieve of Eratosthenes and I’ll publish this addition to the code below as soon as I have time.) So, therefore; the code below is, strictly, not correct for multiplications.

Now for the Swift part;

I’ve created a small class GfPolynomial32 which contains the code for performing modular arithmetic on bit strings up to 32 bits in length. The reason for the “Polynomial” part in the name, btw, is because another way of representing these bit strings is as polynomials of a parameter x which have either a “1” or a “0” coefficient. I.e., given the bit string 10110 we can interpret this as;

 1*x^4+0*x^3+1*x^2+1*x+0

I’ve added some extension methods to create string representations of the bits either as conventional 1’s and 0’s or in the polynomial representation.

Here, without further ado, is the code (just remember that it has not been written for speed and performance, but rather for understandability)

//

//  Gf2Polynomial32.swift

//

//  Created by Jarl Ostensen on 07/01/02015.

//  Copyright (c) 2015 SonarJetLens. All rights reserved.

//

import Foundation

// Gf2 representation of a polynomial of maximum order 31

public class Gf2Polynomial32 {

    private var _value:UInt32;

    private let _valueMsbPos:UInt32 = 0;

    // return the position of the highest set bit in value

    private func msbPosOf(Value:UInt32) -> UInt32 {

        var pos:UInt32 = 0;

        var c = Value;

        while(c != 0) {

            c >>= 1;

            ++pos;

        }

        return pos – 1;

    }

    public init(val:UInt32) {

        _value = val;

                  _valueMsbPos = msbPosOf(_value);

    }

    // order of the polynomial (equals position of msb)

    public var order:UInt32 {

        get {

            return _valueMsbPos;

        }

    }

    // raw value

    public var value:UInt32 {

        get {

            return _value;

        }

    }

}

// ====================================== various operators on Gf2Polynomial32’s;

public func == (left:Gf2Polynomial32, right:Gf2Polynomial32) -> Bool {

    return left.value == right.value;

}

public func != (left:Gf2Polynomial32, right:Gf2Polynomial32) -> Bool {

    return left.value != right.value;

}

public func + (left:Gf2Polynomial32, right:Gf2Polynomial32) -> Gf2Polynomial32 {

    return Gf2Polynomial32(val:(left.value ^ right.value));

}

public func – (left:Gf2Polynomial32, right:Gf2Polynomial32) -> Gf2Polynomial32 {

    // NOTE: same as for +

    return Gf2Polynomial32(val:(left.value ^ right.value));

}

// multiplication

public func * (left:Gf2Polynomial32, right:Gf2Polynomial32) -> Gf2Polynomial32 {

    var b:UInt32 = right.value;

    let a:UInt32 = left.value;

    var result:UInt32 = 0;

    var shift:UInt32 = 0;

    while(b != 0) {

        if ( b&1 == 1 ) {

            result ^= (a << shift);

        }

        ++shift;

        b >>= 1;

    }

    return Gf2Polynomial32(val:result);

}

// division (returns a pair; quotient and remainder)

public func / (left:Gf2Polynomial32, right:Gf2Polynomial32) -> (Gf2Polynomial32,Gf2Polynomial32) {

    let shiftAlign = (left.order – right.order);

    let divisor = right.value << shiftAlign;

    let highBitMask:UInt32 = 1 << left.order;

    var q:UInt32 = 0;

    var dividend = left.value;

    var qShift = shiftAlign+1;

    do {

        if ( (dividend&highBitMask) != 0 ) {

            dividend ^= divisor;

            q ^= (1 << (qShift-1));

        }

        –qShift;

        dividend = (dividend << 1);

    } while( qShift != 0 );

    return (Gf2Polynomial32(val:q),Gf2Polynomial32(val:(dividend>>(shiftAlign+1))));

}

And some helpful extensions;

//

//  Gf2Polynomial32+Extension.swift

//

//  Created by Jarl Ostensen on 07/01/02015.

//  Copyright (c) 2015 SonarJetLens. All rights reserved.

//

import Foundation

extension Gf2Polynomial32 {

    // return a 1’s and 0’s representation

    func asBinaryString() -> String {

        if ( value == 0 ) {

            return “0”;

        }

        else {

            var result = “”;

            var c = value;

            while( c != 0 ) {

                if ( (c & 1) != 0 ) {

                    result = “1” + result;

                }

                else {

                    result = “0” + result;

                }

                c >>= 1;

            }

            return “0b” + result;

        }

    }

    // return a polynomial representation

    func asPolynomial() -> String {

        if ( value==0 ) {

            return “0”;

        }

        else {

            var result = “”;

            var c = value;

            var pos = 0;

            while( c != 0 ) {

                if ( (c & 1) != 0 ) {

                    if ( pos>0 ) {

                        result = “x^\(pos) + (result.utf16Count>0 ? ” + \(result) : “”);

                    }

                    else {

                        result = “1”;

                    }

                }

                c >>= 1;

                ++pos;

            }

            return result;

        }

    }

}

And finally; here is some example output from using this class;

N = 0b1100110111; x^9 + x^8 + x^5 + x^4 + x^2 + x^1 + 1

D = 0b10011; x^4 + x^1 + 1

q = 0b110110; x^5 + x^4 + x^2 + x^1

r = 0b1101; x^3 + x^2 + 1

q*D     = x^9 + x^8 + x^5 + x^4 + x^3 + x^1

q*D + r = x^9 + x^8 + x^5 + x^4 + x^2 + x^1 + 1

Rule 30 – a 1D cellular automaton implemented using Processing

Processing is a simple but very powerful Java based scripting environment that allows one to quickly implement graphics (2D and 3D) prototypes. There are some pretty impressive visualization examples out there for you to feast your eyes on and mine is a simple contribution but hopefully interesting.

Cellular automaton are the basis for Stephen Wolfram’s theories in his book “A New Kind Of Science”. On the surface they are deceptively simple but, just like fractals, they can generate highly complex detail and seem to be able to replicate natural phenomena.

I’ve implemented a simple 1D cellular automaton in the script included here based on information on Wolfram’s Mathworld site.

The image this code generates is a pretty pyramid shape of generations of the bits generated by a rule (more on that below) called “Rule 30” which, according to Wolfram, generates a random sequence of bits. It’s actually used as the basis for Mathematica’s random number generator. It also looks quite pretty:
Rule 30

The algorithm is very simple, and I think (hope) it’s pretty well commented and should be easy to understand but let me just outline it for reference:

  • Start with a sequence of bits. It can be initialized to whatever you want but I’ve set it to all 0 except 1 bit in the “center” of the bit string (width of screen in pixels). This bit acts as a seed
  • Read each bit and select an output bit depending on it’s nearest neighbours (left,right) using a simple FSM rule. I’ve encoded the famous “Rule 30” but it could be any other rule…try it!
  • Render the line…and go to the next line

It really is that simple. The code does a little bit of fancy double buffering of the line rendered and the output line and the rest is basic use of Processing’s incredibly simple and powerful functionality.

Enjoy…


/*
Elementary Cellular Automaton (1-D) implementation
We render "generations" of the automaton, one per screen line, starting with at the top with a "seed".
Each generation (line) is rendered by testing each pixel's neighbours and selecting the appropriate rule.

*/

// the rule determining if a pixels is “on” or “off” depending on the state of its neighbours
int[] rule;
// the
int[] grid;
int WIDTH = 800;

void setup() {
size(WIDTH,WIDTH/2); //< half height to terminate nicely since we don’t deal with pixels at the borders
background(0);
noLoop(); //< only need to render one frame

rule = new int[2*2*2]; //< 8 rules
grid = new int[width*2]; //< double buffer; holds this and last generation

// initialize first generation with 1 pixel in the center
for ( int n=0; n < width; ++n ) {
grid[n] = (n == (width/2) ? 1 : 0);
}

// “rule 30”
// This rule is actually used to generate random numbers in Mathematica
rule[0] = 0;
rule[1] = 0;
rule[2] = 0;
rule[3] = 1;
rule[4] = 1;
rule[5] = 1;
rule[6] = 1;
rule[7] = 0;
}

void draw() {

loadPixels();

// double buffer selector
int sel = 1;
for( int generation = 1; generation < height; ++generation ) {

// double buffer selector: alternates between 1 and 0 (every other line)
sel ^= 1;
// read from last generation line
int read_offs = sel*width;
// write to this generation line
int write_offs = (sel^1)*width;

// offset into screen buffer of current line
int pixel_offs = generation*width;

// generate the image, line by line (generation by generation)
for ( int n = 1; n < (width-1); ++n ) {

// read out this and 2 neighbouring values
int n1 = grid[read_offs+(n-1)];
int n2 = grid[read_offs+n];
int n3 = grid[read_offs+(n+1)];
int value;

//simple state machine to select the appropriate rule
if ( n1==1 ) {
if ( n2==1 ) {
if ( n3==1 ) {
// 1,1,1
value = rule[0];
}
else {
// 1,1,0
value = rule[1];
}
}
else {
if ( n3==1 ) {
// 1,0,1
value = rule[2];
}
else {
// 1,0,0
value = rule[3];
}
}
}
else {
if ( n2==1 ) {
if ( n3==1 ) {
// 0,1,1
value = rule[4];
}
else {
// 0,1,0
value = rule[5];
}
}
else {
if ( n3==1 ) {
// 0,0,1
value = rule[6];
}
else {
// 0,0,0
value = rule[7];
}
}
}

// update the current “write buffer” for the next round
grid[write_offs+n] = value;
// draw and map so we get a fade out towards the bottom
pixels[pixel_offs+n] = color(value*map(generation, 1,height, 255,20));
}
}

updatePixels();

}