cubefs

1
/**
2
 * Matrix Algebra over an 8-bit Galois Field
3
 *
4
 * Copyright 2015, Klaus Post
5
 * Copyright 2015, Backblaze, Inc.
6
 */
7

8
package reedsolomon
9

10
import (
11
	"errors"
12
	"fmt"
13
	"strconv"
14
	"strings"
15
)
16

17
// byte[row][col]
18
type matrix [][]byte
19

20
// newMatrix returns a matrix of zeros.
21
func newMatrix(rows, cols int) (matrix, error) {
22
	if rows <= 0 {
23
		return nil, errInvalidRowSize
24
	}
25
	if cols <= 0 {
26
		return nil, errInvalidColSize
27
	}
28

29
	m := matrix(make([][]byte, rows))
30
	for i := range m {
31
		m[i] = make([]byte, cols)
32
	}
33
	return m, nil
34
}
35

36
// NewMatrixData initializes a matrix with the given row-major data.
37
// Note that data is not copied from input.
38
func newMatrixData(data [][]byte) (matrix, error) {
39
	m := matrix(data)
40
	err := m.Check()
41
	if err != nil {
42
		return nil, err
43
	}
44
	return m, nil
45
}
46

47
// IdentityMatrix returns an identity matrix of the given size.
48
func identityMatrix(size int) (matrix, error) {
49
	m, err := newMatrix(size, size)
50
	if err != nil {
51
		return nil, err
52
	}
53
	for i := range m {
54
		m[i][i] = 1
55
	}
56
	return m, nil
57
}
58

59
// errInvalidRowSize will be returned if attempting to create a matrix with negative or zero row number.
60
var errInvalidRowSize = errors.New("invalid row size")
61

62
// errInvalidColSize will be returned if attempting to create a matrix with negative or zero column number.
63
var errInvalidColSize = errors.New("invalid column size")
64

65
// errColSizeMismatch is returned if the size of matrix columns mismatch.
66
var errColSizeMismatch = errors.New("column size is not the same for all rows")
67

68
func (m matrix) Check() error {
69
	rows := len(m)
70
	if rows <= 0 {
71
		return errInvalidRowSize
72
	}
73
	cols := len(m[0])
74
	if cols <= 0 {
75
		return errInvalidColSize
76
	}
77

78
	for _, col := range m {
79
		if len(col) != cols {
80
			return errColSizeMismatch
81
		}
82
	}
83
	return nil
84
}
85

86
// String returns a human-readable string of the matrix contents.
87
//
88
// Example: [[1, 2], [3, 4]]
89
func (m matrix) String() string {
90
	rowOut := make([]string, 0, len(m))
91
	for _, row := range m {
92
		colOut := make([]string, 0, len(row))
93
		for _, col := range row {
94
			colOut = append(colOut, strconv.Itoa(int(col)))
95
		}
96
		rowOut = append(rowOut, "["+strings.Join(colOut, ", ")+"]")
97
	}
98
	return "[" + strings.Join(rowOut, ", ") + "]"
99
}
100

101
// Multiply multiplies this matrix (the one on the left) by another
102
// matrix (the one on the right) and returns a new matrix with the result.
103
func (m matrix) Multiply(right matrix) (matrix, error) {
104
	if len(m[0]) != len(right) {
105
		return nil, fmt.Errorf("columns on left (%d) is different than rows on right (%d)", len(m[0]), len(right))
106
	}
107
	result, _ := newMatrix(len(m), len(right[0]))
108
	for r, row := range result {
109
		for c := range row {
110
			var value byte
111
			for i := range m[0] {
112
				value ^= galMultiply(m[r][i], right[i][c])
113
			}
114
			result[r][c] = value
115
		}
116
	}
117
	return result, nil
118
}
119

120
// Augment returns the concatenation of this matrix and the matrix on the right.
121
func (m matrix) Augment(right matrix) (matrix, error) {
122
	if len(m) != len(right) {
123
		return nil, errMatrixSize
124
	}
125

126
	result, _ := newMatrix(len(m), len(m[0])+len(right[0]))
127
	for r, row := range m {
128
		for c := range row {
129
			result[r][c] = m[r][c]
130
		}
131
		cols := len(m[0])
132
		for c := range right[0] {
133
			result[r][cols+c] = right[r][c]
134
		}
135
	}
136
	return result, nil
137
}
138

139
// errMatrixSize is returned if matrix dimensions are doesn't match.
140
var errMatrixSize = errors.New("matrix sizes do not match")
141

142
func (m matrix) SameSize(n matrix) error {
143
	if len(m) != len(n) {
144
		return errMatrixSize
145
	}
146
	for i := range m {
147
		if len(m[i]) != len(n[i]) {
148
			return errMatrixSize
149
		}
150
	}
151
	return nil
152
}
153

154
// SubMatrix returns a part of this matrix. Data is copied.
155
func (m matrix) SubMatrix(rmin, cmin, rmax, cmax int) (matrix, error) {
156
	result, err := newMatrix(rmax-rmin, cmax-cmin)
157
	if err != nil {
158
		return nil, err
159
	}
160
	// OPTME: If used heavily, use copy function to copy slice
161
	for r := rmin; r < rmax; r++ {
162
		for c := cmin; c < cmax; c++ {
163
			result[r-rmin][c-cmin] = m[r][c]
164
		}
165
	}
166
	return result, nil
167
}
168

169
// SwapRows Exchanges two rows in the matrix.
170
func (m matrix) SwapRows(r1, r2 int) error {
171
	if r1 < 0 || len(m) <= r1 || r2 < 0 || len(m) <= r2 {
172
		return errInvalidRowSize
173
	}
174
	m[r2], m[r1] = m[r1], m[r2]
175
	return nil
176
}
177

178
// IsSquare will return true if the matrix is square
179
// and nil if the matrix is square
180
func (m matrix) IsSquare() bool {
181
	return len(m) == len(m[0])
182
}
183

184
// errSingular is returned if the matrix is singular and cannot be inversed
185
var errSingular = errors.New("matrix is singular")
186

187
// errNotSquare is returned if attempting to inverse a non-square matrix.
188
var errNotSquare = errors.New("only square matrices can be inverted")
189

190
// Invert returns the inverse of this matrix.
191
// Returns ErrSingular when the matrix is singular and doesn't have an inverse.
192
// The matrix must be square, otherwise ErrNotSquare is returned.
193
func (m matrix) Invert() (matrix, error) {
194
	if !m.IsSquare() {
195
		return nil, errNotSquare
196
	}
197

198
	size := len(m)
199
	work, _ := identityMatrix(size)
200
	work, _ = m.Augment(work)
201

202
	err := work.gaussianElimination()
203
	if err != nil {
204
		return nil, err
205
	}
206

207
	return work.SubMatrix(0, size, size, size*2)
208
}
209

210
func (m matrix) gaussianElimination() error {
211
	rows := len(m)
212
	columns := len(m[0])
213
	// Clear out the part below the main diagonal and scale the main
214
	// diagonal to be 1.
215
	for r := 0; r < rows; r++ {
216
		// If the element on the diagonal is 0, find a row below
217
		// that has a non-zero and swap them.
218
		if m[r][r] == 0 {
219
			for rowBelow := r + 1; rowBelow < rows; rowBelow++ {
220
				if m[rowBelow][r] != 0 {
221
					err := m.SwapRows(r, rowBelow)
222
					if err != nil {
223
						return err
224
					}
225
					break
226
				}
227
			}
228
		}
229
		// If we couldn't find one, the matrix is singular.
230
		if m[r][r] == 0 {
231
			return errSingular
232
		}
233
		// Scale to 1.
234
		if m[r][r] != 1 {
235
			scale := galDivide(1, m[r][r])
236
			for c := 0; c < columns; c++ {
237
				m[r][c] = galMultiply(m[r][c], scale)
238
			}
239
		}
240
		// Make everything below the 1 be a 0 by subtracting
241
		// a multiple of it.  (Subtraction and addition are
242
		// both exclusive or in the Galois field.)
243
		for rowBelow := r + 1; rowBelow < rows; rowBelow++ {
244
			if m[rowBelow][r] != 0 {
245
				scale := m[rowBelow][r]
246
				for c := 0; c < columns; c++ {
247
					m[rowBelow][c] ^= galMultiply(scale, m[r][c])
248
				}
249
			}
250
		}
251
	}
252

253
	// Now clear the part above the main diagonal.
254
	for d := 0; d < rows; d++ {
255
		for rowAbove := 0; rowAbove < d; rowAbove++ {
256
			if m[rowAbove][d] != 0 {
257
				scale := m[rowAbove][d]
258
				for c := 0; c < columns; c++ {
259
					m[rowAbove][c] ^= galMultiply(scale, m[d][c])
260
				}
261

262
			}
263
		}
264
	}
265
	return nil
266
}
267

268
// Create a Vandermonde matrix, which is guaranteed to have the
269
// property that any subset of rows that forms a square matrix
270
// is invertible.
271
func vandermonde(rows, cols int) (matrix, error) {
272
	result, err := newMatrix(rows, cols)
273
	if err != nil {
274
		return nil, err
275
	}
276
	for r, row := range result {
277
		for c := range row {
278
			result[r][c] = galExp(byte(r), c)
279
		}
280
	}
281
	return result, nil
282
}
283