Skip to content

Instantly share code, notes, and snippets.

@raylee

raylee/Reed-Solomon-encoder.md

Last active September 27, 2022 04:07
Show Gist options
  • Save raylee/b532d6a763d618fc2f29cb02799dc985 to your computer and use it in GitHub Desktop.
Save raylee/b532d6a763d618fc2f29cb02799dc985 to your computer and use it in GitHub Desktop.
Download ZIP
gf256 implements arithmetic over the Galois Field GF(256) and Reed-Solomon encoding. A mechanical transliteration of rsc.io/qr/gf256. ECCs about 115MiB/s/thread on a low-end virt or laptop.
Raw
gf256.c
// gf256.c
// gf256 implements arithmetic over the Galois Field GF(256) and Reed-
// Solomon encoding. This is a mostly mechanical transliteration of
// rsc.io/qr/gf256, lowering it from Go to C. As such, this remains:
// Copyright 2010 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the License-gf256 file.
// Any errors or infelicities mine. github.com/raylee
// For a primer on how this works, see https://research.swtch.com/field.
#include <execinfo.h>
#include <signal.h>
#include <stdarg.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define Unused(x) (void)(x)
typedef unsigned char byte;
// A Field represents an instance of GF(256) defined by a specific polynomial.
typedef struct Field {
byte log[256]; // log[0] is unused
byte exp[510];
} Field;
// An RSEncoder implements Reed-Solomon encoding
// over a given field using a given number of error correction bytes.
typedef struct RSEncoder {
Field *f;
int c;
byte *gen;
byte *lgen;
byte *p;
size_t p_len;
} RSEncoder;
Field NewField(int poly, int alpha);
RSEncoder NewRSEncoder(Field *f, int c);
void ECC(RSEncoder *rs, byte *data, size_t data_len, byte *check, size_t check_len);
static byte Add(Field *f, byte x, byte y);
static byte Exp(Field *f, int e);
static int Log(Field *f, byte x);
static byte Inv(Field *f, byte x);
static byte Mul(Field *f, byte x, byte y);
static int mul(int x, int y, int poly);
static unsigned nbit(int p);
static bool reducible(int p);
static void panic(char *fmt, ...);
static int polyDiv(int p, int q);
static unsigned zallocCount = 0;
// zalloc malloc()s len bytes and clears the buffer before returning it.
static void *zalloc(size_t len) {
void *m = malloc(len);
bzero(m, len);
zallocCount++;
return m;
}
// panic is equivalent to fprintf(stderr, fmt, ...) followed by raising SIGABRT.
static void panic(char *fmt, ...) {
// Print the supplied message.
va_list ap;
va_start(ap, fmt);
fprintf(stderr, "panic: ");
vfprintf(stderr, fmt, ap);
va_end(ap);
fprintf(stderr, "\n");
// Print a backtrace.
void *array[10];
size_t size = backtrace(array, 10);
char **strings = backtrace_symbols(array, size);
for (size_t i = 0; i < size; i++) {
fprintf(stderr, " bt: %s\n", strings[i]);
}
free(strings);
fflush(stderr);
raise(SIGABRT);
}
// NewField returns a new field corresponding to the polynomial poly
// and generator alpha. The Reed-Solomon encoding in QR codes uses
// polynomial 0x11d with generator 2.
//
// The choice of generator alpha only affects the Exp and Log operations.
Field NewField(int poly, int alpha) {
if (poly < 0x100 || poly >= 0x200 || reducible(poly)) {
panic("gf256: invalid polynomial: %d", poly);
}
Field f;
int x = 1;
for (int i = 0; i < 255; i++) {
if (x == 1 && i != 0) {
panic("gf256: invalid generator %d for polynomial %d", alpha, poly);
}
f.exp[i] = x;
f.exp[i+255] = x;
f.log[x] = i;
x = mul(x, alpha, poly);
}
f.log[0] = 255;
for (int i = 0; i < 255; i++) {
if (f.log[f.exp[i]] != i) {
panic("bad log");
}
if (f.log[f.exp[i+255]] != i) {
panic("bad log");
}
}
for (int i = 1; i < 256; i++) {
if (f.exp[f.log[i]] != i) {
panic("bad log");
}
}
return f;
}
// nbit returns the number of significant bits in p.
static unsigned nbit(int p) {
unsigned n = 0;
for (; p > 0; p >>= 1) {
n++;
}
return n;
}
// polyDiv divides the polynomial p by q and returns the remainder.
static int polyDiv(int p, int q) {
unsigned np = nbit(p), nq = nbit(q);
for (; np >= nq; np--) {
if ((p&(1<<(np-1))) != 0) {
p ^= q << (np - nq);
}
}
return p;
}
// mul returns the product x*y mod poly, a GF(256) multiplication.
static int mul(int x, int y, int poly) {
int z = 0;
while (x > 0) {
if ((x&1) != 0) {
z ^= y;
}
x >>= 1;
y <<= 1;
if ((y&0x100) != 0) {
y ^= poly;
}
}
return z;
}
// reducible reports whether p is reducible.
static bool reducible(int p) {
// Multiplying n-bit * n-bit produces (2n-1)-bit,
// so if p is reducible, one of its factors must be
// of np/2+1 bits or fewer.
int np = nbit(p);
for (int q = 2; q < 1<<(np/2+1); q++) {
if (polyDiv(p, q) == 0) {
return true;
}
}
return false;
}
// Add returns the sum of x and y in the field.
static byte Add(Field *f, byte x, byte y) {
Unused(f);
return x ^ y;
}
// Exp returns the base-α exponential of e in the field.
// If e < 0, Exp returns 0.
static byte Exp(Field *f, int e) {
if (e < 0) {
return 0;
}
return f->exp[e%255];
}
// Log returns the base-α logarithm of x in the field.
// If x == 0, Log returns -1.
static int Log(Field *f, byte x) {
if (x == 0) {
return -1;
}
return f->log[x];
}
// Inv returns the multiplicative inverse of x in the field.
// If x == 0, Inv returns 0.
static byte Inv(Field *f, byte x) {
if (x == 0) {
return 0;
}
return f->exp[255-f->log[x]];
}
// Mul returns the product of x and y in the field.
static byte Mul(Field *f, byte x, byte y) {
if (x == 0 || y == 0) {
return 0;
}
return f->exp[f->log[x]+f->log[y]];
}
static void generate(Field *f, int e, byte **gen, byte **lgen) {
// p = 1
byte *p = zalloc(e+1);
p[e] = 1;
for (int i = 0; i < e; i++) {
// p *= (x + Exp(i))
// p[j] = p[j]*Exp(i) + p[j+1].
byte c = Exp(f, i);
for (int j = 0; j < e; j++) {
p[j] = Mul(f, p[j], c) ^ p[j+1];
}
p[e] = Mul(f, p[e], c);
}
// lp = log p.
byte *lp = zalloc(e+1);
for (int i = 0; i < e+1; i++) {
byte c = p[i];
if (c == 0) {
lp[i] = 255;
} else {
lp[i] = Log(f, c);
}
}
*gen = p;
*lgen = lp;
return;
}
// NewRSEncoder returns a new Reed-Solomon encoder
// over the given field and number of error correction bytes.
RSEncoder NewRSEncoder(Field *f, int c) {
byte *gen, *lgen;
generate(f, c, &gen, &lgen);
return (RSEncoder){.f = f, .c = c, .gen = gen, .lgen = lgen, .p_len = 0};
}
// ECC writes to check the error correcting code bytes
// for data using the given Reed-Solomon parameters.
void ECC(RSEncoder *rs, byte *data, size_t data_len, byte *check, size_t check_len) {
if (check_len < (size_t)rs->c) {
panic("gf256: invalid check byte length");
}
if (rs->c == 0) {
return;
}
// The check bytes are the remainder after dividing
// data padded with c zeros by the generator polynomial.
// p = data padded with c zeros.
byte *p;
size_t n = data_len + rs->c;
size_t p_len;
if (rs->p_len >= n) {
p = rs->p;
p_len = rs->p_len;
} else {
p = zalloc(n);
p_len = n;
}
memcpy(p, data, data_len);
for (size_t i = data_len; i < p_len; i++) {
p[i] = 0;
}
// Divide p by gen, leaving the remainder in p[data_len:].
// p[0] is the most significant term in p, and
// gen[0] is the most significant term in the generator,
// which is always 1.
// To avoid repeated work, we store various values as
// lv, not v, where lv = log[v].
Field *f = rs->f;
byte *lgen = &rs->lgen[1];
for (size_t i = 0; i < data_len; i++) {
byte c = p[i];
if (c == 0) {
continue;
}
byte *q = &p[i+1];
byte *exp = &f->exp[f->log[c]];
for (int j = 0; j < rs->c; j++) {
byte lg = lgen[j];
if (lg != 255) { // lgen uses 255 for log 0
q[j] ^= exp[lg];
}
}
}
memcpy(check, &p[data_len], rs->c);
rs->p = p;
rs->p_len = p_len;
}
#if defined(TEST) || defined(BENCHMARK)
Field f;
void hexdump(FILE *fp, byte *data, size_t len) {
char *hex="0123456789abcdef";
for (size_t i=0; i<len; i++) {
fputc(hex[(data[i]&0xf0) >> 4], fp);
fputc(hex[data[i]&0x0f], fp);
if (i % 16 == 15) {
fputc('\n', fp);
continue;
}
if (i % 4 == 3) {
fputc(' ', fp);
}
fputc(' ', fp);
}
fputc('\n', fp);
}
void TestBasic() {
if (Exp(&f, 0) != 1 || Exp(&f, 1) != 2 || Exp(&f, 255) != 1) {
panic("bad Exp");
}
}
#define Len(x) (sizeof(x)/sizeof(x[0]))
void TestECC() {
byte data[] = {0x10, 0x20, 0x0c, 0x56, 0x61, 0x80, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11};
byte check[] = {0xa5, 0x24, 0xd4, 0xc1, 0xed, 0x36, 0xc7, 0x87, 0x2c, 0x55};
byte out[Len(check)];
RSEncoder rs = NewRSEncoder(&f, Len(check));
ECC(&rs, data, Len(data), out, Len(check));
if (memcmp(out, check, Len(check))!=0) {
fputs("have: ", stderr);
hexdump(stderr, out, Len(out));
fputs("want: ", stderr);
hexdump(stderr, check, Len(check));
panic("have (out) != want (check)");
}
}
/*
func TestLinear(t *testing.T) {
d1 := []byte{0x00, 0x00}
c1 := []byte{0x00, 0x00}
out := make([]byte, len(c1))
rs := NewRSEncoder(f, len(c1))
if rs.ECC(d1, out); !bytes.Equal(out, c1) {
t.Errorf("ECBytes(%x, %d) = %x, want 0", d1, len(c1), out)
}
d2 := []byte{0x00, 0x01}
c2 := make([]byte, 2)
rs.ECC(d2, c2)
d3 := []byte{0x00, 0x02}
c3 := make([]byte, 2)
rs.ECC(d3, c3)
cx := make([]byte, 2)
for i := range cx {
cx[i] = c2[i] ^ c3[i]
}
d4 := []byte{0x00, 0x03}
c4 := make([]byte, 2)
rs.ECC(d4, c4)
if !bytes.Equal(cx, c4) {
t.Errorf("ECBytes(%x, 2) = %x\nECBytes(%x, 2) = %x\nxor = %x\nECBytes(%x, 2) = %x",
d2, c2, d3, c3, cx, d4, c4)
}
}
func TestGaussJordan(t *testing.T) {
rs := NewRSEncoder(f, 2)
m := make([][]byte, 16)
for i := range m {
m[i] = make([]byte, 4)
m[i][i/8] = 1 << uint(i%8)
rs.ECC(m[i][:2], m[i][2:])
}
if false {
fmt.Printf("---\n")
for _, row := range m {
fmt.Printf("%x\n", row)
}
}
b := []uint{0, 1, 2, 3, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25, 26, 27}
for i := 0; i < 16; i++ {
bi := b[i]
if m[i][bi/8]&(1<<(7-bi%8)) == 0 {
for j := i + 1; ; j++ {
if j >= len(m) {
t.Errorf("lost track for %d", bi)
break
}
if m[j][bi/8]&(1<<(7-bi%8)) != 0 {
m[i], m[j] = m[j], m[i]
break
}
}
}
for j := i + 1; j < len(m); j++ {
if m[j][bi/8]&(1<<(7-bi%8)) != 0 {
for k := range m[j] {
m[j][k] ^= m[i][k]
}
}
}
}
if false {
fmt.Printf("---\n")
for _, row := range m {
fmt.Printf("%x\n", row)
}
}
for i := 15; i >= 0; i-- {
bi := b[i]
for j := i - 1; j >= 0; j-- {
if m[j][bi/8]&(1<<(7-bi%8)) != 0 {
for k := range m[j] {
m[j][k] ^= m[i][k]
}
}
}
}
if false {
fmt.Printf("---\n")
for _, row := range m {
fmt.Printf("%x", row)
out := make([]byte, 2)
if rs.ECC(row[:2], out); !bytes.Equal(out, row[2:]) {
fmt.Printf(" - want %x", out)
}
fmt.Printf("\n")
}
}
}
*/
// Now returns the current time.
static struct timespec Now() {
struct timespec ts;
timespec_get(&ts, TIME_UTC);
return ts;
}
// since returns how long it's been since then in ns.
static uint64_t since(struct timespec then) {
struct timespec now = Now();
return (now.tv_sec - then.tv_sec) * 1e9 + now.tv_nsec - then.tv_nsec;
}
void BenchmarkECC() {
byte data[] = {0x10, 0x20, 0x0c, 0x56, 0x61, 0x80, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11, 0x10, 0x20, 0x0c, 0x56, 0x61, 0x80, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11, 0xec, 0x11};
byte check[] = {0x29, 0x41, 0xb3, 0x93, 0x8, 0xe8, 0xa3, 0xe7, 0x63, 0x8f};
byte out[Len(check)];
int iterations = 1e6;
struct timespec start = Now();
RSEncoder rs = NewRSEncoder(&f, Len(check));
for (int i = 0; i < iterations; i++) {
ECC(&rs, data, Len(data), out, Len(out));
}
uint64_t elapsed = since(start);
double rate = Len(data)*iterations;
rate /= (double)elapsed / 1.0e9;
fprintf(stderr, "%.0f bytes / second\n", rate);
if (memcmp(out, check, Len(check))!=0) {
fputs("have: ", stderr);
hexdump(stderr, out, Len(out));
fputs("want: ", stderr);
hexdump(stderr, check, Len(check));
panic("have (out) != want (check)");
}
}
void TestGenerate() {
for (int i = 0; i < 256; i++) {
byte *l, *lg;
generate(&f, i, &l, &lg);
if (lg[0] != 0) {
panic("#%d: %x", i, lg);
}
}
}
void TestReducible() {
byte count[] = {1, 2, 3, 6, 9, 18, 30, 56, 99, 186}; // oeis.org/A1037
for (size_t i = 0; i < Len(count); i++) {
byte want = count[i];
int n = 0;
for (int p = 1 << (i+2); p < 1<<(i+3); p++) {
if (!reducible(p)) {
n++;
}
}
if (n != want) {
panic("#reducible(%d-bit) = %d, want %d", i+2, n, want);
}
}
}
bool generates(int alpha, int poly);
void TestExhaustive() {
for (int poly = 0x100; poly < 0x200; poly++) {
if (reducible(poly)) {
continue;
}
int alpha = 2;
while (!generates(alpha, poly)) {
alpha++;
}
Field f = NewField(poly, alpha);
for (int p = 0; p < 256; p++) {
for (int q = 0; q < 256; q++) {
byte fm = Mul(&f, p, q);
byte pm = mul(p, q, poly);
if (fm != pm) {
panic("NewField(%#x).Mul(%#x, %#x) = %#x, want %#x", poly, p, q, fm, pm);
}
}
}
}
}
bool generates(int alpha, int poly) {
int x = alpha;
for (int i = 0; i < 254; i++) {
if (x == 1) {
return false;
}
x = mul(x, alpha, poly);
}
return true;
}
int main(void) {
f = NewField(0x11d, 2); // x^8 + x^4 + x^3 + x^2 + x^0
TestBasic();
TestReducible();
TestExhaustive();
TestGenerate();
TestECC();
zallocCount = 0;
#if defined(BENCHMARK)
BenchmarkECC();
#endif
fprintf(stderr, "zalloc count %d\n", zallocCount);
}
#endif
// Testing (add -DBENCHMARK for speed tests):
// cc -DTEST -O3 -o gf_test -fsanitize=address -fno-omit-frame-pointer gf.c
// ./gf_test
Raw
~ LICENSE
Copyright (c) 2009 The Go Authors. All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the
distribution.
* Neither the name of Google Inc. nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment