1// (C) 1996-1998 Luigi Rizzo (luigi@iet.unipi.it)
2//     2009-2010 Jack Lloyd (lloyd@randombit.net)
3//     2011 Billy Brumley (billy.brumley@aalto.fi)
4//     2016-2017 Vivint, Inc. (jeff.wendling@vivint.com)
5//
6// Portions derived from code by Phil Karn (karn@ka9q.ampr.org),
7// Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu) and Hari
8// Thirumoorthy (harit@spectra.eng.hawaii.edu), Aug 1995
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15//    notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18//    notice, this list of conditions and the following disclaimer in the
19//    documentation and/or other materials provided with the
20//    distribution.
21//
22// THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
23// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
24// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
25// DISCLAIMED. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT,
26// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
27// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
28// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
30// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
31// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
32// POSSIBILITY OF SUCH DAMAGE.
33
34package infectious
35
36import (
37	"bytes"
38	"errors"
39)
40
41type pivotSearcher struct {
42	k    int
43	ipiv []bool
44}
45
46func newPivotSearcher(k int) *pivotSearcher {
47	return &pivotSearcher{
48		k:    k,
49		ipiv: make([]bool, k),
50	}
51}
52
53func (p *pivotSearcher) search(col int, matrix []byte) (int, int, error) {
54	if p.ipiv[col] == false && matrix[col*p.k+col] != 0 {
55		p.ipiv[col] = true
56		return col, col, nil
57	}
58
59	for row := 0; row < p.k; row++ {
60		if p.ipiv[row] {
61			continue
62		}
63
64		for i := 0; i < p.k; i++ {
65			if p.ipiv[i] == false && matrix[row*p.k+i] != 0 {
66				p.ipiv[i] = true
67				return row, i, nil
68			}
69		}
70	}
71
72	return 0, 0, errors.New("pivot not found")
73}
74
75func swap(a, b *byte) {
76	tmp := *a
77	*a = *b
78	*b = tmp
79}
80
81// TODO(jeff): matrix is a K*K array, row major.
82func invertMatrix(matrix []byte, k int) error {
83	pivot_searcher := newPivotSearcher(k)
84	indxc := make([]int, k)
85	indxr := make([]int, k)
86	id_row := make([]byte, k)
87
88	for col := 0; col < k; col++ {
89		icol, irow, err := pivot_searcher.search(col, matrix)
90		if err != nil {
91			return err
92		}
93
94		if irow != icol {
95			for i := 0; i < k; i++ {
96				swap(&matrix[irow*k+i], &matrix[icol*k+i])
97			}
98		}
99
100		indxr[col] = irow
101		indxc[col] = icol
102		pivot_row := matrix[icol*k:][:k]
103		c := pivot_row[icol]
104
105		if c == 0 {
106			return errors.New("singular matrix")
107		}
108
109		if c != 1 {
110			c = gf_inverse[c]
111			pivot_row[icol] = 1
112			mul_c := gf_mul_table[c][:]
113
114			for i := 0; i < k; i++ {
115				pivot_row[i] = mul_c[pivot_row[i]]
116			}
117		}
118
119		id_row[icol] = 1
120		if !bytes.Equal(pivot_row, id_row) {
121			p := matrix
122			for i := 0; i < k; i++ {
123				if i != icol {
124					c = p[icol]
125					p[icol] = 0
126					addmul(p[:k], pivot_row, c)
127				}
128				p = p[k:]
129			}
130		}
131
132		id_row[icol] = 0
133	}
134
135	for i := 0; i < k; i++ {
136		if indxr[i] != indxc[i] {
137			for row := 0; row < k; row++ {
138				swap(&matrix[row*k+indxr[i]], &matrix[row*k+indxc[i]])
139			}
140		}
141	}
142	return nil
143}
144
145func createInvertedVdm(vdm []byte, k int) {
146	if k == 1 {
147		vdm[0] = 1
148		return
149	}
150
151	b := make([]byte, k)
152	c := make([]byte, k)
153
154	c[k-1] = 0
155	for i := 1; i < k; i++ {
156		mul_p_i := gf_mul_table[gf_exp[i]][:]
157		for j := k - 1 - (i - 1); j < k-1; j++ {
158			c[j] ^= mul_p_i[c[j+1]]
159		}
160		c[k-1] ^= gf_exp[i]
161	}
162
163	for row := 0; row < k; row++ {
164		index := 0
165		if row != 0 {
166			index = int(gf_exp[row])
167		}
168		mul_p_row := gf_mul_table[index][:]
169
170		t := byte(1)
171		b[k-1] = 1
172		for i := k - 2; i >= 0; i-- {
173			b[i] = c[i+1] ^ mul_p_row[b[i+1]]
174			t = b[i] ^ mul_p_row[t]
175		}
176
177		mul_t_inv := gf_mul_table[gf_inverse[t]][:]
178		for col := 0; col < k; col++ {
179			vdm[col*k+row] = mul_t_inv[b[col]]
180		}
181	}
182}
183