Initial version

master
klauspost 2015-06-19 16:31:24 +02:00
parent 8bebd4bfa9
commit 2b4171b9e6
7 changed files with 1396 additions and 0 deletions

191
galois.go Normal file
View File

@ -0,0 +1,191 @@
/**
* 8-bit Galois Field
* Copyright 2015, Klaus Post
* Copyright 2015, Backblaze, Inc. All rights reserved.
*/
package reedsolomon
const (
// The number of elements in the field.
fieldSize = 256
// The polynomial used to generate the logarithm table.
//
// There are a number of polynomials that work to generate
// a Galois field of 256 elements. The choice is arbitrary,
// and we just use the first one.
//
// The possibilities are: 29, 43, 45, 77, 95, 99, 101, 105,
//* 113, 135, 141, 169, 195, 207, 231, and 245.
generatingPolynomial = 29
)
var logTable = [fieldSize]byte{
0, 0, 1, 25, 2, 50, 26, 198,
3, 223, 51, 238, 27, 104, 199, 75,
4, 100, 224, 14, 52, 141, 239, 129,
28, 193, 105, 248, 200, 8, 76, 113,
5, 138, 101, 47, 225, 36, 15, 33,
53, 147, 142, 218, 240, 18, 130, 69,
29, 181, 194, 125, 106, 39, 249, 185,
201, 154, 9, 120, 77, 228, 114, 166,
6, 191, 139, 98, 102, 221, 48, 253,
226, 152, 37, 179, 16, 145, 34, 136,
54, 208, 148, 206, 143, 150, 219, 189,
241, 210, 19, 92, 131, 56, 70, 64,
30, 66, 182, 163, 195, 72, 126, 110,
107, 58, 40, 84, 250, 133, 186, 61,
202, 94, 155, 159, 10, 21, 121, 43,
78, 212, 229, 172, 115, 243, 167, 87,
7, 112, 192, 247, 140, 128, 99, 13,
103, 74, 222, 237, 49, 197, 254, 24,
227, 165, 153, 119, 38, 184, 180, 124,
17, 68, 146, 217, 35, 32, 137, 46,
55, 63, 209, 91, 149, 188, 207, 205,
144, 135, 151, 178, 220, 252, 190, 97,
242, 86, 211, 171, 20, 42, 93, 158,
132, 60, 57, 83, 71, 109, 65, 162,
31, 45, 67, 216, 183, 123, 164, 118,
196, 23, 73, 236, 127, 12, 111, 246,
108, 161, 59, 82, 41, 157, 85, 170,
251, 96, 134, 177, 187, 204, 62, 90,
203, 89, 95, 176, 156, 169, 160, 81,
11, 245, 22, 235, 122, 117, 44, 215,
79, 174, 213, 233, 230, 231, 173, 232,
116, 214, 244, 234, 168, 80, 88, 175,
}
/**
* Inverse of the logarithm table. Maps integer logarithms
* to members of the field. There is no entry for 255
* because the highest log is 254.
*
* This table was generated by `go run gentables.go`
*/
var expTable = []byte{0x1, 0x2, 0x4, 0x8, 0x10, 0x20, 0x40, 0x80, 0x1d, 0x3a, 0x74, 0xe8, 0xcd, 0x87, 0x13, 0x26, 0x4c, 0x98, 0x2d, 0x5a, 0xb4, 0x75, 0xea, 0xc9, 0x8f, 0x3, 0x6, 0xc, 0x18, 0x30, 0x60, 0xc0, 0x9d, 0x27, 0x4e, 0x9c, 0x25, 0x4a, 0x94, 0x35, 0x6a, 0xd4, 0xb5, 0x77, 0xee, 0xc1, 0x9f, 0x23, 0x46, 0x8c, 0x5, 0xa, 0x14, 0x28, 0x50, 0xa0, 0x5d, 0xba, 0x69, 0xd2, 0xb9, 0x6f, 0xde, 0xa1, 0x5f, 0xbe, 0x61, 0xc2, 0x99, 0x2f, 0x5e, 0xbc, 0x65, 0xca, 0x89, 0xf, 0x1e, 0x3c, 0x78, 0xf0, 0xfd, 0xe7, 0xd3, 0xbb, 0x6b, 0xd6, 0xb1, 0x7f, 0xfe, 0xe1, 0xdf, 0xa3, 0x5b, 0xb6, 0x71, 0xe2, 0xd9, 0xaf, 0x43, 0x86, 0x11, 0x22, 0x44, 0x88, 0xd, 0x1a, 0x34, 0x68, 0xd0, 0xbd, 0x67, 0xce, 0x81, 0x1f, 0x3e, 0x7c, 0xf8, 0xed, 0xc7, 0x93, 0x3b, 0x76, 0xec, 0xc5, 0x97, 0x33, 0x66, 0xcc, 0x85, 0x17, 0x2e, 0x5c, 0xb8, 0x6d, 0xda, 0xa9, 0x4f, 0x9e, 0x21, 0x42, 0x84, 0x15, 0x2a, 0x54, 0xa8, 0x4d, 0x9a, 0x29, 0x52, 0xa4, 0x55, 0xaa, 0x49, 0x92, 0x39, 0x72, 0xe4, 0xd5, 0xb7, 0x73, 0xe6, 0xd1, 0xbf, 0x63, 0xc6, 0x91, 0x3f, 0x7e, 0xfc, 0xe5, 0xd7, 0xb3, 0x7b, 0xf6, 0xf1, 0xff, 0xe3, 0xdb, 0xab, 0x4b, 0x96, 0x31, 0x62, 0xc4, 0x95, 0x37, 0x6e, 0xdc, 0xa5, 0x57, 0xae, 0x41, 0x82, 0x19, 0x32, 0x64, 0xc8, 0x8d, 0x7, 0xe, 0x1c, 0x38, 0x70, 0xe0, 0xdd, 0xa7, 0x53, 0xa6, 0x51, 0xa2, 0x59, 0xb2, 0x79, 0xf2, 0xf9, 0xef, 0xc3, 0x9b, 0x2b, 0x56, 0xac, 0x45, 0x8a, 0x9, 0x12, 0x24, 0x48, 0x90, 0x3d, 0x7a, 0xf4, 0xf5, 0xf7, 0xf3, 0xfb, 0xeb, 0xcb, 0x8b, 0xb, 0x16, 0x2c, 0x58, 0xb0, 0x7d, 0xfa, 0xe9, 0xcf, 0x83, 0x1b, 0x36, 0x6c, 0xd8, 0xad, 0x47, 0x8e, 0x1, 0x2, 0x4, 0x8, 0x10, 0x20, 0x40, 0x80, 0x1d, 0x3a, 0x74, 0xe8, 0xcd, 0x87, 0x13, 0x26, 0x4c, 0x98, 0x2d, 0x5a, 0xb4, 0x75, 0xea, 0xc9, 0x8f, 0x3, 0x6, 0xc, 0x18, 0x30, 0x60, 0xc0, 0x9d, 0x27, 0x4e, 0x9c, 0x25, 0x4a, 0x94, 0x35, 0x6a, 0xd4, 0xb5, 0x77, 0xee, 0xc1, 0x9f, 0x23, 0x46, 0x8c, 0x5, 0xa, 0x14, 0x28, 0x50, 0xa0, 0x5d, 0xba, 0x69, 0xd2, 0xb9, 0x6f, 0xde, 0xa1, 0x5f, 0xbe, 0x61, 0xc2, 0x99, 0x2f, 0x5e, 0xbc, 0x65, 0xca, 0x89, 0xf, 0x1e, 0x3c, 0x78, 0xf0, 0xfd, 0xe7, 0xd3, 0xbb, 0x6b, 0xd6, 0xb1, 0x7f, 0xfe, 0xe1, 0xdf, 0xa3, 0x5b, 0xb6, 0x71, 0xe2, 0xd9, 0xaf, 0x43, 0x86, 0x11, 0x22, 0x44, 0x88, 0xd, 0x1a, 0x34, 0x68, 0xd0, 0xbd, 0x67, 0xce, 0x81, 0x1f, 0x3e, 0x7c, 0xf8, 0xed, 0xc7, 0x93, 0x3b, 0x76, 0xec, 0xc5, 0x97, 0x33, 0x66, 0xcc, 0x85, 0x17, 0x2e, 0x5c, 0xb8, 0x6d, 0xda, 0xa9, 0x4f, 0x9e, 0x21, 0x42, 0x84, 0x15, 0x2a, 0x54, 0xa8, 0x4d, 0x9a, 0x29, 0x52, 0xa4, 0x55, 0xaa, 0x49, 0x92, 0x39, 0x72, 0xe4, 0xd5, 0xb7, 0x73, 0xe6, 0xd1, 0xbf, 0x63, 0xc6, 0x91, 0x3f, 0x7e, 0xfc, 0xe5, 0xd7, 0xb3, 0x7b, 0xf6, 0xf1, 0xff, 0xe3, 0xdb, 0xab, 0x4b, 0x96, 0x31, 0x62, 0xc4, 0x95, 0x37, 0x6e, 0xdc, 0xa5, 0x57, 0xae, 0x41, 0x82, 0x19, 0x32, 0x64, 0xc8, 0x8d, 0x7, 0xe, 0x1c, 0x38, 0x70, 0xe0, 0xdd, 0xa7, 0x53, 0xa6, 0x51, 0xa2, 0x59, 0xb2, 0x79, 0xf2, 0xf9, 0xef, 0xc3, 0x9b, 0x2b, 0x56, 0xac, 0x45, 0x8a, 0x9, 0x12, 0x24, 0x48, 0x90, 0x3d, 0x7a, 0xf4, 0xf5, 0xf7, 0xf3, 0xfb, 0xeb, 0xcb, 0x8b, 0xb, 0x16, 0x2c, 0x58, 0xb0, 0x7d, 0xfa, 0xe9, 0xcf, 0x83, 0x1b, 0x36, 0x6c, 0xd8, 0xad, 0x47, 0x8e}
/*
var expTable = [fieldSize * 2]int8{
1, 2, 4, 8, 16, 32, 64, -128,
29, 58, 116, -24, -51, -121, 19, 38,
76, -104, 45, 90, -76, 117, -22, -55,
-113, 3, 6, 12, 24, 48, 96, -64,
-99, 39, 78, -100, 37, 74, -108, 53,
106, -44, -75, 119, -18, -63, -97, 35,
70, -116, 5, 10, 20, 40, 80, -96,
93, -70, 105, -46, -71, 111, -34, -95,
95, -66, 97, -62, -103, 47, 94, -68,
101, -54, -119, 15, 30, 60, 120, -16,
-3, -25, -45, -69, 107, -42, -79, 127,
-2, -31, -33, -93, 91, -74, 113, -30,
-39, -81, 67, -122, 17, 34, 68, -120,
13, 26, 52, 104, -48, -67, 103, -50,
-127, 31, 62, 124, -8, -19, -57, -109,
59, 118, -20, -59, -105, 51, 102, -52,
-123, 23, 46, 92, -72, 109, -38, -87,
79, -98, 33, 66, -124, 21, 42, 84,
-88, 77, -102, 41, 82, -92, 85, -86,
73, -110, 57, 114, -28, -43, -73, 115,
-26, -47, -65, 99, -58, -111, 63, 126,
-4, -27, -41, -77, 123, -10, -15, -1,
-29, -37, -85, 75, -106, 49, 98, -60,
-107, 55, 110, -36, -91, 87, -82, 65,
-126, 25, 50, 100, -56, -115, 7, 14,
28, 56, 112, -32, -35, -89, 83, -90,
81, -94, 89, -78, 121, -14, -7, -17,
-61, -101, 43, 86, -84, 69, -118, 9,
18, 36, 72, -112, 61, 122, -12, -11,
-9, -13, -5, -21, -53, -117, 11, 22,
44, 88, -80, 125, -6, -23, -49, -125,
27, 54, 108, -40, -83, 71, -114,
// Repeat the table a second time, so multiply()
// does not have to check bounds.
1, 2, 4, 8, 16, 32, 64, -128,
29, 58, 116, -24, -51, -121, 19, 38,
76, -104, 45, 90, -76, 117, -22, -55,
-113, 3, 6, 12, 24, 48, 96, -64,
-99, 39, 78, -100, 37, 74, -108, 53,
106, -44, -75, 119, -18, -63, -97, 35,
70, -116, 5, 10, 20, 40, 80, -96,
93, -70, 105, -46, -71, 111, -34, -95,
95, -66, 97, -62, -103, 47, 94, -68,
101, -54, -119, 15, 30, 60, 120, -16,
-3, -25, -45, -69, 107, -42, -79, 127,
-2, -31, -33, -93, 91, -74, 113, -30,
-39, -81, 67, -122, 17, 34, 68, -120,
13, 26, 52, 104, -48, -67, 103, -50,
-127, 31, 62, 124, -8, -19, -57, -109,
59, 118, -20, -59, -105, 51, 102, -52,
-123, 23, 46, 92, -72, 109, -38, -87,
79, -98, 33, 66, -124, 21, 42, 84,
-88, 77, -102, 41, 82, -92, 85, -86,
73, -110, 57, 114, -28, -43, -73, 115,
-26, -47, -65, 99, -58, -111, 63, 126,
-4, -27, -41, -77, 123, -10, -15, -1,
-29, -37, -85, 75, -106, 49, 98, -60,
-107, 55, 110, -36, -91, 87, -82, 65,
-126, 25, 50, 100, -56, -115, 7, 14,
28, 56, 112, -32, -35, -89, 83, -90,
81, -94, 89, -78, 121, -14, -7, -17,
-61, -101, 43, 86, -84, 69, -118, 9,
18, 36, 72, -112, 61, 122, -12, -11,
-9, -13, -5, -21, -53, -117, 11, 22,
44, 88, -80, 125, -6, -23, -49, -125,
27, 54, 108, -40, -83, 71, -114,
}
*/
func galAdd(a, b byte) byte {
return a ^ b
}
func galSub(a, b byte) byte {
return a ^ b
}
// galMultiply multiplies to elements of the field.
func galMultiply(a, b byte) byte {
if a == 0 || b == 0 {
return 0
}
logA := int(logTable[a])
logB := int(logTable[b])
return expTable[logA+logB]
}
// galDivide is inverse of galMultiply.
func galDivide(a, b byte) byte {
if a == 0 {
return 0
}
if b == 0 {
panic("Argument 'divisor' is 0")
}
logA := int(logTable[a])
logB := int(logTable[b])
logResult := logA - logB
if logResult < 0 {
logResult += 255
}
return expTable[logResult]
}
// Computes a**n.
//
// The result will be the same as multiplying a times itself n times.
func galExp(a byte, n int) byte {
if n == 0 {
return 1
}
if a == 0 {
return 0
}
logA := logTable[a]
logResult := int(logA) * n
for logResult >= 255 {
logResult -= 255
}
return byte(expTable[logResult])
}

138
galois_test.go Normal file
View File

@ -0,0 +1,138 @@
/**
* Unit tests for Galois
*
* Copyright 2015, Klaus Post
* Copyright 2015, Backblaze, Inc.
*/
package reedsolomon
import (
"testing"
)
func TestAssociativity(t *testing.T) {
for i := 0; i < 256; i++ {
a := byte(i)
for j := 0; j < 256; j++ {
b := byte(j)
for k := 0; k < 256; k++ {
c := byte(k)
x := galAdd(a, galAdd(b, c))
y := galAdd(galAdd(a, b), c)
if x != y {
t.Fatal("add does not match:", x, "!=", y)
}
x = galMultiply(a, galMultiply(b, c))
y = galMultiply(galMultiply(a, b), c)
if x != y {
t.Fatal("multiply does not match:", x, "!=", y)
}
}
}
}
}
func TestIdentity(t *testing.T) {
for i := 0; i < 256; i++ {
a := byte(i)
b := galAdd(a, 0)
if a != b {
t.Fatal("Add zero should yield same result", a, "!=", b)
}
b = galMultiply(a, 1)
if a != b {
t.Fatal("Mul by one should yield same result", a, "!=", b)
}
}
}
func TestInverse(t *testing.T) {
for i := 0; i < 256; i++ {
a := byte(i)
b := galSub(0, a)
c := galAdd(a, b)
if c != 0 {
t.Fatal("inverse sub/add", c, "!=", 0)
}
if a != 0 {
b = galDivide(1, a)
c = galMultiply(a, b)
if c != 1 {
t.Fatal("inverse div/mul", c, "!=", 1)
}
}
}
}
func TestCommutativity(t *testing.T) {
for i := 0; i < 256; i++ {
a := byte(i)
for j := 0; j < 256; j++ {
b := byte(j)
x := galAdd(a, b)
y := galAdd(b, a)
if x != y {
t.Fatal(x, "!= ", y)
}
x = galMultiply(a, b)
y = galMultiply(b, a)
if x != y {
t.Fatal(x, "!= ", y)
}
}
}
}
func TestDistributivity(t *testing.T) {
for i := 0; i < 256; i++ {
a := byte(i)
for j := 0; j < 256; j++ {
b := byte(j)
for k := 0; k < 256; k++ {
c := byte(k)
x := galMultiply(a, galAdd(b, c))
y := galAdd(galMultiply(a, b), galMultiply(a, c))
if x != y {
t.Fatal(x, "!= ", y)
}
}
}
}
}
func TestExp(t *testing.T) {
for i := 0; i < 256; i++ {
a := byte(i)
power := byte(1)
for j := 0; j < 256; j++ {
x := galExp(a, j)
if x != power {
t.Fatal(x, "!=", power)
}
power = galMultiply(power, a)
}
}
}
func TestGalois(t *testing.T) {
// These values were copied output of the Python code.
if galMultiply(3, 4) != 12 {
t.Fatal("galMultiply(3, 4) != 12")
}
if galMultiply(7, 7) != 21 {
t.Fatal("galMultiply(7, 7) != 21")
}
if galMultiply(23, 45) != 41 {
t.Fatal("galMultiply(23, 45) != 41")
}
if galExp(2, 2) != 4 {
t.Fatal("galExp(2, 2) != 4")
}
if galExp(5, 20) != 235 {
t.Fatal("galExp(5, 20) != 235")
}
if galExp(13, 7) != 43 {
t.Fatal("galExp(13, 7) != 43")
}
}

75
gentables.go Normal file
View File

@ -0,0 +1,75 @@
//+build ignore
package main
import (
"fmt"
)
var logTable = [fieldSize]int16{
-1, 0, 1, 25, 2, 50, 26, 198,
3, 223, 51, 238, 27, 104, 199, 75,
4, 100, 224, 14, 52, 141, 239, 129,
28, 193, 105, 248, 200, 8, 76, 113,
5, 138, 101, 47, 225, 36, 15, 33,
53, 147, 142, 218, 240, 18, 130, 69,
29, 181, 194, 125, 106, 39, 249, 185,
201, 154, 9, 120, 77, 228, 114, 166,
6, 191, 139, 98, 102, 221, 48, 253,
226, 152, 37, 179, 16, 145, 34, 136,
54, 208, 148, 206, 143, 150, 219, 189,
241, 210, 19, 92, 131, 56, 70, 64,
30, 66, 182, 163, 195, 72, 126, 110,
107, 58, 40, 84, 250, 133, 186, 61,
202, 94, 155, 159, 10, 21, 121, 43,
78, 212, 229, 172, 115, 243, 167, 87,
7, 112, 192, 247, 140, 128, 99, 13,
103, 74, 222, 237, 49, 197, 254, 24,
227, 165, 153, 119, 38, 184, 180, 124,
17, 68, 146, 217, 35, 32, 137, 46,
55, 63, 209, 91, 149, 188, 207, 205,
144, 135, 151, 178, 220, 252, 190, 97,
242, 86, 211, 171, 20, 42, 93, 158,
132, 60, 57, 83, 71, 109, 65, 162,
31, 45, 67, 216, 183, 123, 164, 118,
196, 23, 73, 236, 127, 12, 111, 246,
108, 161, 59, 82, 41, 157, 85, 170,
251, 96, 134, 177, 187, 204, 62, 90,
203, 89, 95, 176, 156, 169, 160, 81,
11, 245, 22, 235, 122, 117, 44, 215,
79, 174, 213, 233, 230, 231, 173, 232,
116, 214, 244, 234, 168, 80, 88, 175,
}
const (
// The number of elements in the field.
fieldSize = 256
// The polynomial used to generate the logarithm table.
//
// There are a number of polynomials that work to generate
// a Galois field of 256 elements. The choice is arbitrary,
// and we just use the first one.
//
// The possibilities are: 29, 43, 45, 77, 95, 99, 101, 105,
//* 113, 135, 141, 169, 195, 207, 231, and 245.
generatingPolynomial = 29
)
func main() {
t := generateExpTable()
fmt.Printf("var expTable = %#v\n", t)
}
/**
* Generates the inverse log table.
*/
func generateExpTable() []byte {
result := make([]byte, fieldSize*2-2)
for i := 1; i < fieldSize; i++ {
log := logTable[i]
result[log] = byte(i)
result[log+fieldSize-1] = byte(i)
}
return result
}

282
matrix.go Normal file
View File

@ -0,0 +1,282 @@
/**
* Matrix Algebra over an 8-bit Galois Field
*
* Copyright 2015, Klaus Post
* Copyright 2015, Backblaze, Inc.
*/
package reedsolomon
import (
"errors"
"fmt"
"strconv"
"strings"
)
// byte[row][col]
type matrix [][]byte
// newMatrix returns a matrix of zeros.
func newMatrix(rows, cols int) (matrix, error) {
if rows <= 0 {
return nil, ErrInvalidRowSize
}
if cols <= 0 {
return nil, ErrInvalidColSize
}
m := matrix(make([][]byte, rows))
for i := range m {
m[i] = make([]byte, cols)
}
return m, nil
}
// NewMatrixData initializes a matrix with the given row-major data.
// Note that data is not copied from input.
func newMatrixData(data [][]byte) (matrix, error) {
m := matrix(data)
err := m.Check()
if err != nil {
return nil, err
}
return m, nil
}
// IdentityMatrix returns an identity matrix of the given size.
func identityMatrix(size int) (matrix, error) {
m, err := newMatrix(size, size)
if err != nil {
return nil, err
}
for i := range m {
m[i][i] = 1
}
return m, nil
}
var ErrInvalidRowSize = errors.New("Invalid row size")
var ErrInvalidColSize = errors.New("Invalid column size")
var ErrColSizeMismatch = errors.New("Column size is not the same for all rows")
func (m matrix) Check() error {
rows := len(m)
if rows <= 0 {
return ErrInvalidRowSize
}
cols := len(m[0])
if cols <= 0 {
return ErrInvalidColSize
}
for _, col := range m {
if len(col) != cols {
return ErrColSizeMismatch
}
}
return nil
}
// String returns a human-readable string of the matrix contents.
//
// Example: [[1, 2], [3, 4]]
func (m matrix) String() string {
var rowOut []string
for _, row := range m {
var colOut []string
for _, col := range row {
colOut = append(colOut, strconv.Itoa(int(col)))
}
rowOut = append(rowOut, "["+strings.Join(colOut, ", ")+"]")
}
return "[" + strings.Join(rowOut, ", ") + "]"
}
// Multiply multiplies this matrix (the one on the left) by another
// matrix (the one on the right) and returns a new matrix with the result.
func (m matrix) Multiply(right matrix) (matrix, error) {
if len(m[0]) != len(right) {
return nil, fmt.Errorf("Columns on left (%d) is different than rows on right (%d)", len(m[0]), len(right))
}
result, err := newMatrix(len(m), len(right[0]))
if err != nil {
return nil, err
}
for r, row := range result {
for c := range row {
var value byte
for i := range m[0] {
value ^= galMultiply(m[r][i], right[i][c])
}
result[r][c] = value
}
}
return result, nil
}
// Augment returns the concatenation of this matrix and the matrix on the right.
func (m matrix) Augment(right matrix) (matrix, error) {
if len(m) != len(right) {
return nil, ErrMatrixSize
}
result, _ := newMatrix(len(m), len(m[0])+len(right[0]))
for r, row := range m {
for c := range row {
result[r][c] = m[r][c]
}
cols := len(m[0])
for c := range right[0] {
result[r][cols+c] = right[r][c]
}
}
return result, nil
}
var ErrMatrixSize = errors.New("matrix sizes does not match")
func (m matrix) SameSize(n matrix) error {
if len(m) != len(n) {
return ErrMatrixSize
}
for i := range m {
if len(m[i]) != len(n[i]) {
return ErrMatrixSize
}
}
return nil
}
// Returns a part of this matrix. Data is copied.
func (m matrix) SubMatrix(rmin, cmin, rmax, cmax int) (matrix, error) {
result, err := newMatrix(rmax-rmin, cmax-cmin)
if err != nil {
return nil, err
}
// OPTME: If used heavily, use copy function to copy slice
for r := rmin; r < rmax; r++ {
for c := cmin; c < cmax; c++ {
result[r-rmin][c-cmin] = m[r][c]
}
}
return result, nil
}
// SwapRows Exchanges two rows in the matrix.
func (m matrix) SwapRows(r1, r2 int) error {
if r1 < 0 || len(m) <= r1 || r2 < 0 || len(m) <= r2 {
return ErrInvalidRowSize
}
m[r2], m[r1] = m[r1], m[r2]
return nil
}
// IsSquare will return true if the matrix is square
// and nil if the matrix is square
func (m matrix) IsSquare() bool {
if len(m) != len(m[0]) {
return false
}
return true
}
var ErrSingular = errors.New("matrix is singular")
var ErrNotSquare = errors.New("only square matrices can be inverted")
// Invert returns the inverse of this matrix.
// Returns ErrSingular when the matrix is singular and doesn't have an inverse.
// The matrix must be square, otherwise ErrNotSquare is returned.
func (m matrix) Invert() (matrix, error) {
if !m.IsSquare() {
return nil, ErrNotSquare
}
size := len(m)
work, err := identityMatrix(size)
if err != nil {
return nil, err
}
work, err = m.Augment(work)
if err != nil {
return nil, err
}
err = work.gaussianElimination()
if err != nil {
return nil, err
}
return work.SubMatrix(0, size, size, size*2)
}
func (m matrix) gaussianElimination() error {
rows := len(m)
columns := len(m[0])
// Clear out the part below the main diagonal and scale the main
// diagonal to be 1.
for r := 0; r < rows; r++ {
// If the element on the diagonal is 0, find a row below
// that has a non-zero and swap them.
if m[r][r] == 0 {
for rowBelow := r + 1; rowBelow < rows; rowBelow++ {
if m[rowBelow][r] != 0 {
m.SwapRows(r, rowBelow)
break
}
}
}
// If we couldn't find one, the matrix is singular.
if m[r][r] == 0 {
return ErrSingular
}
// Scale to 1.
if m[r][r] != 1 {
scale := galDivide(1, m[r][r])
for c := 0; c < columns; c++ {
m[r][c] = galMultiply(m[r][c], scale)
}
}
// Make everything below the 1 be a 0 by subtracting
// a multiple of it. (Subtraction and addition are
// both exclusive or in the Galois field.)
for rowBelow := r + 1; rowBelow < rows; rowBelow++ {
if m[rowBelow][r] != 0 {
scale := m[rowBelow][r]
for c := 0; c < columns; c++ {
m[rowBelow][c] ^= galMultiply(scale, m[r][c])
}
}
}
}
// Now clear the part above the main diagonal.
for d := 0; d < rows; d++ {
for rowAbove := 0; rowAbove < d; rowAbove++ {
if m[rowAbove][d] != 0 {
scale := m[rowAbove][d]
for c := 0; c < columns; c++ {
m[rowAbove][c] ^= galMultiply(scale, m[d][c])
}
}
}
}
return nil
}
// Create a Vandermonde matrix, which is guaranteed to have the
// property that any subset of rows that forms a square matrix
// is invertible.
func vandermonde(rows, cols int) (matrix, error) {
result, err := newMatrix(rows, cols)
if err != nil {
return nil, err
}
for r, row := range result {
for c := range row {
result[r][c] = galExp(byte(r), c)
}
}
return result, nil
}

99
matrix_test.go Normal file
View File

@ -0,0 +1,99 @@
/**
* Unit tests for Matrix
*
* Copyright 2015, Klaus Post
* Copyright 2015, Backblaze, Inc. All rights reserved.
*/
package reedsolomon
import (
"testing"
)
func TestMatrixIdentity(t *testing.T) {
m, err := identityMatrix(3)
if err != nil {
t.Fatal(err)
}
str := m.String()
expect := "[[1, 0, 0], [0, 1, 0], [0, 0, 1]]"
if str != expect {
t.Fatal(str, "!=", expect)
}
}
func TestMatricMultiply(t *testing.T) {
m1, err := newMatrixData(
[][]byte{
[]byte{1, 2},
[]byte{3, 4},
})
if err != nil {
t.Fatal(err)
}
m2, err := newMatrixData(
[][]byte{
[]byte{5, 6},
[]byte{7, 8},
})
if err != nil {
t.Fatal(err)
}
actual, err := m1.Multiply(m2)
if err != nil {
t.Fatal(err)
}
str := actual.String()
expect := "[[11, 22], [19, 42]]"
if str != expect {
t.Fatal(str, "!=", expect)
}
}
func TestMatrixInverse(t *testing.T) {
m, err := newMatrixData([][]byte{
[]byte{56, 23, 98},
[]byte{3, 100, 200},
[]byte{45, 201, 123},
})
if err != nil {
t.Fatal(err)
}
m, err = m.Invert()
if err != nil {
t.Fatal(err)
}
str := m.String()
expect := "[[175, 133, 33], [130, 13, 245], [112, 35, 126]]"
if str != expect {
t.Fatal(str, "!=", expect)
}
}
func TestMatrixInverse2(t *testing.T) {
m, err := newMatrixData([][]byte{
[]byte{1, 0, 0, 0, 0},
[]byte{0, 1, 0, 0, 0},
[]byte{0, 0, 0, 1, 0},
[]byte{0, 0, 0, 0, 1},
[]byte{7, 7, 6, 6, 1},
})
if err != nil {
t.Fatal(err)
}
m, err = m.Invert()
if err != nil {
t.Fatal(err)
}
str := m.String()
expect := "[[1, 0, 0, 0, 0]," +
" [0, 1, 0, 0, 0]," +
" [123, 123, 1, 122, 122]," +
" [0, 0, 1, 0, 0]," +
" [0, 0, 0, 1, 0]]"
if str != expect {
t.Fatal(str, "!=", expect)
}
}

391
reedsolomon.go Normal file
View File

@ -0,0 +1,391 @@
/**
* Reed-Solomon Coding over 8-bit values.
*
* Copyright 2015, Klaus Post
* Copyright 2015, Backblaze, Inc.
*/
package reedsolomon
import (
"errors"
"runtime"
"sync"
)
type ReedSolomon struct {
DataShardCount int
ParityShardCount int
TotalShardCount int
m matrix
parity [][]byte
}
type Shard []byte
type Shards []Shard
func New(dataShardCount, parityShardCount int) (*ReedSolomon, error) {
r := ReedSolomon{
DataShardCount: dataShardCount,
ParityShardCount: parityShardCount,
TotalShardCount: dataShardCount + parityShardCount,
}
// Start with a Vandermonde matrix. This matrix would work,
// in theory, but doesn't have the property that the data
// shards are unchanged after encoding.
vm, err := vandermonde(r.TotalShardCount, dataShardCount)
if err != nil {
return nil, err
}
// Multiple by the inverse of the top square of the matrix.
// This will make the top square be the identity matrix, but
// preserve the property that any square subset of rows is
// invertible.
top, _ := vm.SubMatrix(0, 0, dataShardCount, dataShardCount)
top, _ = top.Invert()
r.m, err = vm.Multiply(top)
r.parity = make([][]byte, parityShardCount)
for i := range r.parity {
r.parity[i] = r.m[dataShardCount+i]
}
return &r, err
}
var ErrTooFewShards = errors.New("too few shards given to encode")
// Encodes parity for a set of data shards.
// shards An array containing data shards followed by parity shards.
// Each shard is a byte array, and they must all be the same size.
func (r ReedSolomon) Encode(shards [][]byte) error {
if len(shards) != r.TotalShardCount {
return ErrTooFewShards
}
err := checkShards(shards, false)
if err != nil {
return err
}
// Get the slice of output buffers.
output := shards[r.DataShardCount:]
// Do the coding.
r.codeSomeShards(r.parity, shards[0:r.DataShardCount], output, r.ParityShardCount, len(shards[0]))
return nil
}
// Verify returns true if the parity shards contain the right data.
func (r ReedSolomon) Verify(shards [][]byte) (bool, error) {
if len(shards) != r.TotalShardCount {
return false, ErrTooFewShards
}
err := checkShards(shards, false)
if err != nil {
return false, err
}
// Slice of buffers being checked.
toCheck := shards[r.DataShardCount:]
// Do the checking.
return r.checkSomeShards(r.parity, shards[0:r.DataShardCount], toCheck, r.ParityShardCount, len(shards[0])), nil
}
// TODO: Cleanup Doc:
// Multiplies a subset of rows from a coding matrix by a full set of
// input shards to produce some output shards.
// 'matrixRows' is The rows from the matrix to use.
// 'inputs' An array of byte arrays, each of which is one input shard.
// The inputs array may have extra buffers after the ones
// that are used. They will be ignored. The number of inputs used is
// determined by the length of each matrix row.
// outputs Byte arrays where the computed shards are stored.
// The outputs array may also have extra, unused, elements
// at the end. The number of outputs computed, and the
// number of matrix rows used, is determined by
// outputCount.
// outputCount The number of outputs to compute.
func (r ReedSolomon) codeSomeShards(matrixRows, inputs, outputs [][]byte, outputCount, byteCount int) {
if runtime.GOMAXPROCS(0) > 1 {
r.codeSomeShardsP(matrixRows, inputs, outputs, outputCount, byteCount)
return
}
for iByte := 0; iByte < byteCount; iByte++ {
for iRow := 0; iRow < outputCount; iRow++ {
matrixRow := matrixRows[iRow]
var value byte
for c := 0; c < r.DataShardCount; c++ {
// note: manual inlining is slower
value ^= galMultiply(matrixRow[c], inputs[c][iByte])
}
outputs[iRow][iByte] = value
}
}
}
//
const splitSize = 256
// Perform the same as codeSomeShards, but split the workload into
// several goroutines.
func (r ReedSolomon) codeSomeShardsP(matrixRows, inputs, outputs [][]byte, outputCount, byteCount int) {
var wg sync.WaitGroup
left := byteCount
start := 0
for {
do := left
if do > splitSize {
do = splitSize
}
if do == 0 {
break
}
left -= do
wg.Add(1)
go func(start, stop int) {
for iByte := start; iByte < stop; iByte++ {
for iRow := 0; iRow < outputCount; iRow++ {
matrixRow := matrixRows[iRow]
var value byte
for c := 0; c < r.DataShardCount; c++ {
// note: manual inlining is slower
value ^= galMultiply(matrixRow[c], inputs[c][iByte])
}
outputs[iRow][iByte] = value
}
}
wg.Done()
}(start, start+do)
start += do
}
wg.Wait()
}
/* Not really faster
func (r ReedSolomon) codeSomeShards(matrixRows, inputs, outputs [][]byte, outputCount, byteCount int) {
matrixRows = matrixRows[0:outputCount]
for iByte := 0; iByte < byteCount; iByte++ {
for iRow, matrixRow := range matrixRows {
var value byte
matrixRow := matrixRow[0:r.DataShardCount]
for c, mr := range matrixRow {
// note: manual inlining is slower
value ^= galMultiply(mr, inputs[c][iByte])
}
outputs[iRow][iByte] = value
}
}
}
*/
func (r ReedSolomon) checkSomeShards(matrixRows, inputs, toCheck [][]byte, outputCount, byteCount int) bool {
if runtime.GOMAXPROCS(0) > 1 {
return r.checkSomeShardsP(matrixRows, inputs, toCheck, outputCount, byteCount)
}
for iByte := 0; iByte < byteCount; iByte++ {
for iRow := 0; iRow < outputCount; iRow++ {
matrixRow := matrixRows[iRow]
var value byte
for c := 0; c < r.DataShardCount; c++ {
// note: manual inlining is slower
value ^= galMultiply(matrixRow[c], inputs[c][iByte])
}
if toCheck[iRow][iByte] != value {
return false
}
}
}
return true
}
func (r ReedSolomon) checkSomeShardsP(matrixRows, inputs, toCheck [][]byte, outputCount, byteCount int) bool {
var wg sync.WaitGroup
left := byteCount
start := 0
failed := make(chan bool)
for {
do := left
if do > splitSize {
do = splitSize
}
if do == 0 {
break
}
left -= do
wg.Add(1)
go func(start, stop int) {
defer wg.Done()
for iByte := start; iByte < stop; iByte++ {
for iRow := 0; iRow < outputCount; iRow++ {
matrixRow := matrixRows[iRow]
var value byte
for c := 0; c < r.DataShardCount; c++ {
// note: manual inlining is slower
value ^= galMultiply(matrixRow[c], inputs[c][iByte])
}
if toCheck[iRow][iByte] != value {
close(failed)
return
}
}
}
}(start, start+do)
start += do
}
wg.Wait()
select {
case <-failed:
return false
default:
}
return true
}
var ErrShardNoData = errors.New("no shard data")
var ErrShardSize = errors.New("shard sizes does not match")
func checkShards(shards [][]byte, nilok bool) error {
if len(shards) == 0 {
return ErrShardNoData
}
size := shardSize(shards)
for _, shard := range shards {
if len(shard) != size {
if len(shard) != 0 || !nilok {
return ErrShardSize
}
}
}
return nil
}
func shardSize(shards [][]byte) int {
for _, shard := range shards {
if len(shard) != 0 {
return len(shard)
}
}
panic("No shards are bigger than 0")
}
// Reconstruct will recreate the missing shards, if possible.
//
// Given a list of shards, some of which contain data, fills in the
// ones that don't have data.
//
// The length of the array must be equal to TotalShardCount.
// You indicate that a shard is missing by setting it to nil.
//
// If there are too few shards to reconstruct the missing
// ones, ErrTooFewShards will be returned.
//
// The reconstructed shard set is complete, but integrity is not verified.
// Use the Verify function to check if data set is ok.
func (r ReedSolomon) Reconstruct(shards [][]byte) error {
if len(shards) != r.TotalShardCount {
return ErrTooFewShards
}
// Check arguments.
err := checkShards(shards, true)
if err != nil {
return err
}
shardSize := shardSize(shards)
// Quick check: are all of the shards present? If so, there's
// nothing to do.
numberPresent := 0
for i := 0; i < r.TotalShardCount; i++ {
if len(shards[i]) != 0 {
numberPresent += 1
}
}
if numberPresent == r.TotalShardCount {
// Cool. All of the shards data data. We don't
// need to do anything.
return nil
}
// More complete sanity check
if numberPresent < r.DataShardCount {
return ErrTooFewShards
}
// Pull out the rows of the matrix that correspond to the
// shards that we have and build a square matrix. This
// matrix could be used to generate the shards that we have
// from the original data.
//
// Also, pull out an array holding just the shards that
// correspond to the rows of the submatrix. These shards
// will be the input to the decoding process that re-creates
// the missing data shards.
subMatrix, err := newMatrix(r.DataShardCount, r.DataShardCount)
if err != nil {
return err
}
subShards := make([][]byte, r.DataShardCount)
subMatrixRow := 0
for matrixRow := 0; matrixRow < r.TotalShardCount && subMatrixRow < r.DataShardCount; matrixRow++ {
if len(shards[matrixRow]) != 0 {
for c := 0; c < r.DataShardCount; c++ {
subMatrix[subMatrixRow][c] = r.m[matrixRow][c]
}
subShards[subMatrixRow] = shards[matrixRow]
subMatrixRow++
}
}
// Invert the matrix, so we can go from the encoded shards
// back to the original data. Then pull out the row that
// generates the shard that we want to decode. Note that
// since this matrix maps back to the orginal data, it can
// be used to create a data shard, but not a parity shard.
dataDecodeMatrix, err := subMatrix.Invert()
if err != nil {
return err
}
// Re-create any data shards that were missing.
//
// The input to the coding is all of the shards we actually
// have, and the output is the missing data shards. The computation
// is done using the special decode matrix we just built.
outputs := make([][]byte, r.ParityShardCount)
matrixRows := make([][]byte, r.ParityShardCount)
outputCount := 0
for iShard := 0; iShard < r.DataShardCount; iShard++ {
if len(shards[iShard]) == 0 {
shards[iShard] = make([]byte, shardSize)
outputs[outputCount] = shards[iShard]
matrixRows[outputCount] = dataDecodeMatrix[iShard]
outputCount++
}
}
r.codeSomeShards(matrixRows, subShards, outputs[:outputCount], outputCount, shardSize)
// Now that we have all of the data shards intact, we can
// compute any of the parity that is missing.
//
// The input to the coding is ALL of the data shards, including
// any that we just calculated. The output is whichever of the
// data shards were missing.
outputCount = 0
for iShard := r.DataShardCount; iShard < r.TotalShardCount; iShard++ {
if len(shards[iShard]) == 0 {
shards[iShard] = make([]byte, shardSize)
outputs[outputCount] = shards[iShard]
matrixRows[outputCount] = r.parity[iShard-r.DataShardCount]
outputCount++
}
}
r.codeSomeShards(matrixRows, shards[:r.DataShardCount], outputs[:outputCount], outputCount, shardSize)
return nil
}

220
reedsolomon_test.go Normal file
View File

@ -0,0 +1,220 @@
/**
* Unit tests for ReedSolomon
*
* Copyright 2015, Klaus Post
* Copyright 2015, Backblaze, Inc. All rights reserved.
*/
package reedsolomon
import (
"math/rand"
"testing"
)
func TestEncoding(t *testing.T) {
perShard := 50000
r, err := New(10, 3)
if err != nil {
t.Fatal(err)
}
shards := make([][]byte, 13)
for s := range shards {
shards[s] = make([]byte, perShard)
}
rand.Seed(0)
for s := 0; s < r.DataShardCount; s++ {
fillRandom(shards[s])
}
err = r.Encode(shards)
if err != nil {
t.Fatal(err)
}
ok, err := r.Verify(shards)
if err != nil {
t.Fatal(err)
}
if !ok {
t.Fatal("Verification failed")
}
}
func TestReconstruct(t *testing.T) {
perShard := 50000
r, err := New(10, 3)
if err != nil {
t.Fatal(err)
}
shards := make([][]byte, 13)
for s := range shards {
shards[s] = make([]byte, perShard)
}
rand.Seed(0)
for s := 0; s < r.DataShardCount; s++ {
fillRandom(shards[s])
}
err = r.Encode(shards)
if err != nil {
t.Fatal(err)
}
shards[0] = nil
shards[7] = nil
shards[11] = nil
err = r.Reconstruct(shards)
if err != nil {
t.Fatal(err)
}
ok, err := r.Verify(shards)
if err != nil {
t.Fatal(err)
}
if !ok {
t.Fatal("Verification failed")
}
}
func TestOneEncode(t *testing.T) {
codec, err := New(5, 5)
if err != nil {
t.Fatal(err)
}
shards := make([][]byte, 10)
shards[0] = []byte{0, 1}
shards[1] = []byte{4, 5}
shards[2] = []byte{2, 3}
shards[3] = []byte{6, 7}
shards[4] = []byte{8, 9}
shards[5] = []byte{0, 0}
shards[6] = []byte{0, 0}
shards[7] = []byte{0, 0}
shards[8] = []byte{0, 0}
shards[9] = []byte{0, 0}
codec.Encode(shards)
if shards[5][0] != 12 || shards[5][1] != 13 {
t.Fatal("shard 5 mismatch")
}
if shards[6][0] != 10 || shards[6][1] != 11 {
t.Fatal("shard 6 mismatch")
}
if shards[7][0] != 14 || shards[7][1] != 15 {
t.Fatal("shard 7 mismatch")
}
if shards[8][0] != 90 || shards[8][1] != 91 {
t.Fatal("shard 8 mismatch")
}
if shards[9][0] != 94 || shards[9][1] != 95 {
t.Fatal("shard 9 mismatch")
}
ok, err := codec.Verify(shards)
if err != nil {
t.Fatal(err)
}
if !ok {
t.Fatal("did not verify")
}
shards[8][0] += 1
ok, err = codec.Verify(shards)
if err != nil {
t.Fatal(err)
}
if ok {
t.Fatal("verify did not fail as expected")
}
}
func fillRandom(b []byte) {
for i := range b {
b[i] = byte(rand.Int() & 0xff)
}
}
func benchmarkEncode(b *testing.B, dataShards, parityShards, shardSize int) {
r, err := New(dataShards, parityShards)
if err != nil {
b.Fatal(err)
}
shards := make([][]byte, r.TotalShardCount)
for s := range shards {
shards[s] = make([]byte, shardSize)
}
rand.Seed(0)
for s := 0; s < r.DataShardCount; s++ {
fillRandom(shards[s])
}
b.SetBytes(int64(shardSize * dataShards))
b.ResetTimer()
for i := 0; i < b.N; i++ {
err = r.Encode(shards)
if err != nil {
b.Fatal(err)
}
}
}
func BenchmarkEncode10_2_10000(b *testing.B) {
benchmarkEncode(b, 10, 2, 10000)
}
func BenchmarkEncode100_20_10000(b *testing.B) {
benchmarkEncode(b, 100, 20, 10000)
}
// Benchmark 10 data shards and 4 parity shards with 16MB each.
func BenchmarkEncode10_4_16M(b *testing.B) {
benchmarkEncode(b, 10, 4, 16*1024*1024)
}
func benchmarkVerify(b *testing.B, dataShards, parityShards, shardSize int) {
r, err := New(dataShards, parityShards)
if err != nil {
b.Fatal(err)
}
shards := make([][]byte, r.TotalShardCount)
for s := range shards {
shards[s] = make([]byte, shardSize)
}
rand.Seed(0)
for s := 0; s < r.DataShardCount; s++ {
fillRandom(shards[s])
}
err = r.Encode(shards)
if err != nil {
b.Fatal(err)
}
b.SetBytes(int64(shardSize * dataShards))
b.ResetTimer()
for i := 0; i < b.N; i++ {
_, err = r.Verify(shards)
if err != nil {
b.Fatal(err)
}
}
}
// Benchmark 10 data slices with 2 parity slices holding 10000 bytes each
func BenchmarkVerify10_2_10000(b *testing.B) {
benchmarkVerify(b, 10, 2, 10000)
}
// Benchmark 50 data slices with 5 parity slices holding 100000 bytes each
func BenchmarkVerify50_5_50000(b *testing.B) {
benchmarkVerify(b, 50, 5, 100000)
}
// Benchmark 10 data slices with 4 parity slices holding 16MB bytes each
func BenchmarkVerify10_4_16M(b *testing.B) {
benchmarkVerify(b, 10, 4, 16*1024*1024)
}