197 lines
5 KiB
C++
197 lines
5 KiB
C++
|
|
/*
|
|
* This file is part of Codecrypt.
|
|
*
|
|
* Codecrypt is free software: you can redistribute it and/or modify it
|
|
* under the terms of the GNU Lesser General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or (at
|
|
* your option) any later version.
|
|
*
|
|
* Codecrypt is distributed in the hope that it will be useful, but WITHOUT
|
|
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
|
|
* License for more details.
|
|
*
|
|
* You should have received a copy of the GNU Lesser General Public License
|
|
* along with Codecrypt. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#include "qd_utils.h"
|
|
|
|
#include <vector>
|
|
using namespace std;
|
|
|
|
/*
|
|
* we count on that all integers are sufficiently large.
|
|
* They should be, largest value occuring should be O(k*n) if initial vector is
|
|
* consisted only from {0,1}^n, and we don't usually have codes of this size.
|
|
*/
|
|
|
|
static void fwht (vector<int> x, vector<int>&r)
|
|
{
|
|
int bs, s;
|
|
s = x.size();
|
|
r.resize (s);
|
|
bs = s >> 1;
|
|
r.swap (x);
|
|
while (bs) {
|
|
x.swap (r);
|
|
for (uint i = 0; i < s; ++i) {
|
|
if ( (i / bs) & 1)
|
|
r[i] = x[i - bs] - x[i];
|
|
else
|
|
r[i] = x[i] + x[i + bs];
|
|
}
|
|
bs >>= 1;
|
|
}
|
|
}
|
|
|
|
//we expect correct parameter size and preallocated out.
|
|
void fwht_dyadic_multiply (const bvector& a, const bvector& b, bvector& out)
|
|
{
|
|
|
|
//lift everyting to Z.
|
|
vector<int> t, A, B;
|
|
uint i;
|
|
|
|
t.resize (a.size() );
|
|
A.resize (a.size() );
|
|
B.resize (a.size() );
|
|
|
|
for (i = 0; i < a.size(); ++i) t[i] = a[i];
|
|
fwht (t, A);
|
|
|
|
for (i = 0; i < b.size(); ++i) t[i] = b[i];
|
|
fwht (t, B);
|
|
|
|
//multiply diagonals to A
|
|
for (i = 0; i < A.size(); ++i) A[i] *= B[i];
|
|
fwht (A, t);
|
|
|
|
uint bitpos = a.size(); //no problem as a.size() == 1<<m == 2^m
|
|
for (i = 0; i < t.size(); ++i) out[i] = (t[i] & bitpos) ? 1 : 0;
|
|
}
|
|
|
|
bool qd_to_right_echelon_form (std::vector<std::vector<bvector> >&mat)
|
|
{
|
|
uint w = mat.size();
|
|
if (!w) return false;
|
|
uint h = mat[0].size();
|
|
if (!h) return false;
|
|
uint bs = mat[0][0].size();
|
|
|
|
uint i, j, k, l;
|
|
|
|
/*
|
|
* Inversion is done the quasi-dyadic way:
|
|
*
|
|
* - because for QD matrix m=delta(h) the product
|
|
* m*m = sum(h) * I, binary QD matrix m is either
|
|
* inversion of itself (m*m=I) or isn't invertible
|
|
* and m*m=0. sum(h), the "count of ones in QD
|
|
* signature mod 2", easily determines the result.
|
|
*
|
|
* - Using blockwise invertions/multiplications,
|
|
* gaussian elimination needed to invert the right
|
|
* square of H can be performed in O(m^2*block_count)
|
|
* matrix operations. Matrix operations are either
|
|
* addition (O(t) on QDs), multiplication(O(t log t)
|
|
* on QDs) or inversion (O(t), as shown above).
|
|
* Whole proces is therefore quite fast.
|
|
*
|
|
* Gaussian elimination on the QD signature should
|
|
* result in something like this: (for m=3, t=4)
|
|
*
|
|
* 1010 0101 1001 1000 0000 0000
|
|
* 0101 1100 1110 0000 1000 0000
|
|
* 0111 1110 0100 0000 0000 1000
|
|
*/
|
|
|
|
bvector tmp;
|
|
tmp.resize (bs);
|
|
for (i = 0; i < h; ++i) { //gauss step
|
|
//first, find a nonsingular matrix in the column
|
|
for (j = i; j < h; ++j)
|
|
if (mat[w - h + i][j]
|
|
.hamming_weight() % 2) break;
|
|
if (j >= h) //none found, die!
|
|
return false;
|
|
|
|
//bring it to correct position (swap it to i-th row)
|
|
if (j > i) for (k = 0; k < w; ++k)
|
|
mat[k][i].swap
|
|
(mat[k][j]);
|
|
|
|
//now normalize the row
|
|
for (j = i; j < h; ++j) {
|
|
l = mat [w - h + i]
|
|
[j].hamming_weight();
|
|
if (l == 0) continue; //zero is just okay :]
|
|
if (! (l % 2) ) //singular, make it regular by adding the i-th row
|
|
for (k = 0;
|
|
k < w;
|
|
++k)
|
|
mat[k][j].add
|
|
(mat[k][i]);
|
|
|
|
//now a matrix is regular, we can easily make it I.
|
|
//first, multiply the row
|
|
for (k = 0; k < w; ++k) {
|
|
//don't overwrite the matrix we're counting with
|
|
if (k == w - h + i) continue;
|
|
fwht_dyadic_multiply
|
|
(mat[w - h + i][j],
|
|
mat[k][j], tmp);
|
|
mat[k][j] = tmp;
|
|
}
|
|
//change the block on the diagonal
|
|
fwht_dyadic_multiply
|
|
(mat[w - h + i][j],
|
|
mat[w - h + i][j], tmp);
|
|
mat[w - h + i][j] = tmp;
|
|
|
|
//and zero the column below diagonal
|
|
if (j > i) for (k = 0; k < w; ++k)
|
|
mat[k][j].add
|
|
(mat[k][i]);
|
|
}
|
|
}
|
|
|
|
for (i = 0; i < h; ++i) { //jordan step
|
|
//normalize diagonal
|
|
for (k = 0; k < w - i; ++k) {
|
|
//we can safely rewrite the diagonal here (nothing's behind it)
|
|
fwht_dyadic_multiply
|
|
(mat[w - i - 1][h - i - 1],
|
|
mat[k][h - i - 1], tmp);
|
|
mat[k][h - i - 1] = tmp;
|
|
}
|
|
|
|
//now make zeroes above
|
|
for (j = i + 1; j < h; ++j) {
|
|
l = mat[w - i - 1]
|
|
[h - j - 1].hamming_weight();
|
|
if (l == 0) continue; //already zero
|
|
if (! (l % 2) ) { //nonsingular, fix it by adding diagonal
|
|
for (k = 0; k < w; ++k)
|
|
mat[k][h - j - 1].add
|
|
(mat[k][h - i - 1]);
|
|
}
|
|
for (k = 0; k < w - i; ++k) {
|
|
//overwrite is also safe here
|
|
fwht_dyadic_multiply
|
|
(mat[w - i - 1]
|
|
[h - j - 1],
|
|
mat[k][h - j - 1], tmp);
|
|
mat[k][h - j - 1] = tmp;
|
|
}
|
|
//I+I=0
|
|
for (k = 0; k < w; ++k)
|
|
mat[k][h - j - 1].add
|
|
(mat[k][h - i - 1]);
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|