cython m1ri base case

#  M1RI Base Case.
#
#  This is a standalone demo which implements 64x64 matrix multiplication
#  in a bitsliced representation of GF(3).
#
#  Copyright (c) 2008-2012 Tomas Boothby
#  Distributed under the terms of the GNU General Public License (GPL3)
#
#  NOT FIT FOR ANY USE

cdef struct vec3:
    long s
    long u

cdef inline void v3_add(vec3 *r,vec3 *a,vec3 *b):
    cdef long t

    r.u = b.s ^ a.u
    r.s = b.u ^ a.s
    r.s = r.s & r.u
    r.u = r.u ^ a.s
    t = b.u ^ a.u
    r.u = t | r.u

cdef inline void v3_sub(vec3 *r, vec3 *a, vec3 *b):
    cdef long t

    r.s = b.u ^ a.u
    r.u = b.s ^ a.s
    r.u = r.u | r.s
    r.s = r.s ^ b.s
    t = b.u ^ a.s
    r.s = t & r.s


cdef inline void v3_mul(vec3 *r, vec3 *a, vec3 *b):
    r.u = b.u & a.u
    r.s = b.s ^ a.s
    r.s = r.s & r.u


cdef inline void v3_iadd(vec3 *r, vec3 *a):
    cdef long t

    t = a.u ^ r.s
    r.s = a.s ^ r.u
    r.u = a.u ^ r.u
    r.s = r.s & t
    t = t ^ a.s
    r.u = t | r.u


cdef inline void v3_isub(vec3 *r, vec3 *a):
    cdef long t

    r.u = a.u ^ r.u
    t = r.u | r.s
    t = t ^ a.s
    r.s = a.u ^ r.s
    r.s = r.s & t
    r.u = t | r.u

cdef inline void combine4(vec3 *table, vec3 *input):
    cdef vec3 t, a, b, c, d
    t.s = t.u = 0
    a = input[0]
    b = input[1]
    c = input[2]
    d = input[3]

    table[0] = t
    table[1] = a
    table[2] = b
    table[4] = c
    table[8] = d

    v3_add(&t,&c,&d)
    table[12] = t

    v3_add(&t,&b,&c)
    table[6] = t
    v3_iadd(&t,&d)
    table[14] = t
    v3_isub(&t,&c)
    table[10] = t

    v3_add(&t,&a,&b)
    table[3] = t
    v3_iadd(&t,&d)
    table[11] = t
    v3_iadd(&t,&c)
    table[15] = t
    v3_isub(&t,&d)
    table[7] = t
    v3_isub(&t,&b)
    table[5] = t
    v3_iadd(&t,&d)
    table[13] = t
    v3_isub(&t,&c)
    table[9] = t

cdef inline void combine5(vec3 *table, vec3 *input):
    cdef vec3 e, *t4
    cdef int i

    combine4(table, input)
    e = input[4]
    t4 = table+16
    table[16] = e

    for i from 0 < i < 16:
        v3_add(t4 + i, table+i, &e)

cdef inline void combine6(vec3 *table, vec3 *input):
    cdef vec3 e, *t4
    cdef vec3 f, *t5
    cdef int i

    combine4(table, input)
    e = input[4]
    t4 = table+16
    table[16] = e

    f = input[5]
    t5 = table+32
    table[32] = f

    for i from 0 < i < 16:
        v3_add(t4 + i, table+i, &e)

    for i from 0 < i < 32:
        v3_add(t5 + i, table+i, &f)

cdef void mul_64(vec3 *R, vec3 *A, vec3 *B):
    cdef int i
    cdef vec3 t1, t2, r1, r2, a
    cdef long v1, v2

    cdef vec3 tables6[4][64]
    cdef vec3 tables5[8][32]
    for i from 0 <= i < 4:
        combine6(tables6[i], B + 0 + 6*i)
    for i from 0 <= i < 8:
        combine5(tables5[i], B + 24 + 5*i)
    for i from 0 <= i < 64:
        a = A[i]
        v2 = a.s
        v1 = a.u ^ v2
        r1 = tables6[0][v1&63];                    v1 >>= 6;
        r2 = tables6[0][v2&63];                    v2 >>= 6;
        t1 = tables6[1][v1&63]; v3_iadd(&r1, &t1); v1 >>= 6;
        t2 = tables6[1][v2&63]; v3_iadd(&r2, &t2); v2 >>= 6;
        t1 = tables6[2][v1&63]; v3_iadd(&r1, &t1); v1 >>= 6;
        t2 = tables6[2][v2&63]; v3_iadd(&r2, &t2); v2 >>= 6;
        t1 = tables6[3][v1&63]; v3_iadd(&r1, &t1); v1 >>= 6;
        t2 = tables6[3][v2&63]; v3_iadd(&r2, &t2); v2 >>= 6;
        t1 = tables5[0][v1&31]; v3_iadd(&r1, &t1); v1 >>= 5;
        t2 = tables5[0][v2&31]; v3_iadd(&r2, &t2); v2 >>= 5;
        t1 = tables5[1][v1&31]; v3_iadd(&r1, &t1); v1 >>= 5;
        t2 = tables5[1][v2&31]; v3_iadd(&r2, &t2); v2 >>= 5;
        t1 = tables5[2][v1&31]; v3_iadd(&r1, &t1); v1 >>= 5;
        t2 = tables5[2][v2&31]; v3_iadd(&r2, &t2); v2 >>= 5;
        t1 = tables5[3][v1&31]; v3_iadd(&r1, &t1); v1 >>= 5;
        t2 = tables5[3][v2&31]; v3_iadd(&r2, &t2); v2 >>= 5;
        t1 = tables5[4][v1&31]; v3_iadd(&r1, &t1); v1 >>= 5;
        t2 = tables5[4][v2&31]; v3_iadd(&r2, &t2); v2 >>= 5;
        t1 = tables5[5][v1&31]; v3_iadd(&r1, &t1); v1 >>= 5;
        t2 = tables5[5][v2&31]; v3_iadd(&r2, &t2); v2 >>= 5;
        t1 = tables5[6][v1&31]; v3_iadd(&r1, &t1); v1 >>= 5;
        t2 = tables5[6][v2&31]; v3_iadd(&r2, &t2); v2 >>= 5;
        t1 = tables5[7][v1&31]; v3_iadd(&r1, &t1);
        t2 = tables5[7][v2&31]; v3_iadd(&r2, &t2);

        v3_isub(&r1, &r2)
        R[i] = r1