lib/checksum/gf2.hpp

/*******************************************************************************
 * Copyright 2014 Trevor Robinson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 ******************************************************************************/
#include "int_types.h"
#include <algorithm> // std::swap

#ifdef _MSC_VER
#pragma warning(disable:4146) // unary minus operator applied to unsigned type, result still unsigned
#endif

namespace pulsar {

// Type trait for unsigned integers of at least N bytes
template<unsigned int N>
struct uint_bytes {
    enum {
        is_defined = 0
    };
};

template<class T>
struct defined_uint_bytes {
    enum {
        is_defined = 1
    };
    typedef T type;
};

template<> struct uint_bytes<1> : defined_uint_bytes<uint8_t> {
};
template<> struct uint_bytes<2> : defined_uint_bytes<uint16_t> {
};
template<> struct uint_bytes<3> : defined_uint_bytes<uint32_t> {
};
template<> struct uint_bytes<4> : defined_uint_bytes<uint32_t> {
};
template<> struct uint_bytes<5> : defined_uint_bytes<uint64_t> {
};
template<> struct uint_bytes<6> : defined_uint_bytes<uint64_t> {
};
template<> struct uint_bytes<7> : defined_uint_bytes<uint64_t> {
};
template<> struct uint_bytes<8> : defined_uint_bytes<uint64_t> {
};

// Type trait for unsigned integers of at least N bits
template<unsigned int N>
struct uint_bits : uint_bytes<(N + 7) / 8> {
    enum {
        bits = 8 * sizeof(typename uint_bits::type)
    };
};

// Bit vector of N bits; currently just exposes an unsigned integer
template<unsigned int N>
class bitvector {
    typedef typename uint_bits<N>::type type;
    type value;
 public:
    bitvector() {
    }
    bitvector(type value)
            : value(value) {
    }
    operator type&() {
        return value;
    }
    operator type() const {
        return value;
    }
};

// Bit matrix of M columns by N rows
template<unsigned int M, unsigned int N = M>
class bitmatrix {
    typedef bitvector<M> row;
    row value[N];
 public:
    bitmatrix() {
    }
    explicit bitmatrix(bool b) {
        if (b)
            identity();
        else
            null();
    }
    void null() {
        for (unsigned int i = 0; i < N; ++i)
            value[i] = 0;
    }
    void identity() {
        for (unsigned int i = 0; i < N; ++i)
            value[i] = i < M ? (const row) 1 << i : 0;
    }
    void lower_shift() {
        for (unsigned int i = 0; i < N; ++i)
            value[i] = i > 0 && i <= M ? (const row) 1 << (i - 1) : 0;
    }
    void upper_shift() {
        for (unsigned int i = 0; i < N; ++i)
            value[i] = i + 1 < M ? (const row) 1 << (i + 1) : 0;
    }
    operator bitvector<M>*() {
        return value;
    }
    operator const bitvector<M>*() const {
        return value;
    }
};

/*
 * Multiplies MxN matrix A by N-row vector B in GF(2).
 *
 * For M,N = 3:
 *
 *     | a b c |      | x |       | ax + by + cz |
 * A = | d e f |, B = | y |, AB = | dx + ey + fz |
 *     | g h i |      | z |       | gx + hy + iz |
 *
 * In GF(2), addition corresponds to XOR and multiplication to AND:
 *
 *      | (a & x) ^ (b & y) ^ (c & z) |
 * AB = | (d & x) ^ (e & y) ^ (f & z) |
 *      | (g & x) ^ (h & y) ^ (i & z) |
 *
 * Trading variable names for [row,column] indices:
 *
 * AB = (A[,0] & B[0]) ^ (A[,1] & B[1]) ^ (A[,2] & B[2]) ^ ...
 *
 * Assuming columns are represented as words and rows as bit offsets,
 * all rows of AB can be calculated in parallel:
 *
 * AB = (A[0] & -((B >> 0) & 1) ^ (A[1] & -((B >> 1) & 1) ^ ...
 */
template<unsigned int M, unsigned int N>
bitvector<M> mul(const bitmatrix<M, N>& a, const bitvector<N> b) {
    bitvector<M> result(0);
    for (unsigned int i = 0; i < N; ++i)
        result ^= a[i] & -((b >> i) & 1);
    return result;
}

/*
 * Multiplies MxN matrix A by NxP matrix B in GF(2).
 *
 * For M,N,P = 3:
 *
 *     | a b c |      | j k l |       | (aj + bm + cp) (ak + bn + cq) (al + bo + cr) |
 * A = | d e f |, B = | m n o |, AB = | (dj + em + fp) (dk + en + fq) (dl + eo + fr) |
 *     | g h i |      | p q r |       | (gj + hm + ip) (gk + hn + iq) (gl + ho + ir) |
 */
template<unsigned int M, unsigned int N, unsigned int P>
void mul(bitmatrix<M, P>& result, const bitmatrix<M, N>& a, const bitmatrix<N, P>& b) {
    for (unsigned int i = 0; i < P; i++)
        result[i] = mul(a, b[i]);
}

/*
 * Squares an NxN matrix in GF(2).
 */
template<unsigned int N>
void sqr(bitmatrix<N, N>& result, const bitmatrix<N, N>& a) {
    mul(result, a, a);
}

/*
 * Raises an NxN matrix to the power n in GF(2) by squaring.
 */
template<unsigned int N>
void pow(bitmatrix<N, N>& result, const bitmatrix<N, N>& a, uint64_t n) {
    result.identity();
    if (n > 0) {
        bitmatrix<N, N> square = a;
        bitmatrix<N, N> temp;
        bitmatrix<N, N> *ptemp = &temp, *psquare = &square, *presult = &result;
        for (;;) {
            if (n & 1) {
                mul(*ptemp, *presult, *psquare);
                std::swap(ptemp, presult);
            }
            if (!(n >>= 1))
                break;
            sqr(*ptemp, *psquare);
            std::swap(ptemp, psquare);
        }
        if (presult != &result)
            result = *presult;
    }
}

} // namespace pulsar