Quellcode-Bibliothek vecgf2.c Sprache: C

/****************************************************************************
**
**  This file is part of GAP, a system for computational discrete algebra.
**
**  Copyright of GAP belongs to its developers, whose names are too numerous
**  to list here. Please refer to the COPYRIGHT file for details.
**
**  SPDX-License-Identifier: GPL-2.0-or-later
*/

#include "vecgf2.h"

#include "ariths.h"
#include "bits_intern.h"
#include "blister.h"
#include "bool.h"
#include "error.h"
#include "finfield.h"
#include "gvars.h"
#include "integer.h"
#include "lists.h"
#include "modules.h"
#include "opers.h"
#include "plist.h"
#include "precord.h"
#include "range.h"
#include "records.h"
#include "stats.h"
#include "vec8bit.h"

#include <gmp.h>

/****************************************************************************
**
*F * * * * * * * * * * * imported library variables * * * * * * * * * * * * *
*/

/****************************************************************************
**
*V  TYPE_LIST_GF2VEC  . . . . . . . . . . . . . . type of a GF2 vector object
*/
Obj TYPE_LIST_GF2VEC;

/****************************************************************************
**
*V  TYPE_LIST_GF2VEC_LOCKED. . . .  type of a mutable GF2 vector object
**                                          with locked representation
*/
Obj TYPE_LIST_GF2VEC_LOCKED;

/****************************************************************************
**
*V  TYPE_LIST_GF2VEC_IMM  . . . . . .  type of an immutable GF2 vector object
*/
Obj TYPE_LIST_GF2VEC_IMM;

/****************************************************************************
**
*V  TYPE_LIST_GF2VEC_IMM_LOCKED . . .  type of an immutable GF2 vector object
**                                          with locked representation
*/
Obj TYPE_LIST_GF2VEC_IMM_LOCKED;

/****************************************************************************
**
*V  TYPE_LIST_GF2MAT  . . . . . . . . . . . . . . type of a GF2 matrix object
*/
Obj TYPE_LIST_GF2MAT;

/****************************************************************************
**
*V  TYPE_LIST_GF2MAT_IMM  . . . . . .  type of an immutable GF2 matrix object
*/
Obj TYPE_LIST_GF2MAT_IMM;

/****************************************************************************
**
*V  IsGF2VectorRep  . . . . . . . . . . . . . . . . . . . . . . . . .  filter
*/
Obj IsGF2VectorRep;

/****************************************************************************
**
*V  GF2One  . . . . . . . . . . . . . . . . . . . . . . . . . . .  one of GF2
*/
static Obj GF2One;

/****************************************************************************
**
*V  GF2Zero . . . . . . . . . . . . . . . . . . . . . . . . . . . zero of GF2
*/
static Obj GF2Zero;

/****************************************************************************
**
*F * * * * * * * * * * * * arithmetic operations  * * * * * * * * * * * * * *
*/

static inline void AddGF2VecToGF2Vec(UInt * ptS, const UInt * ptV, UInt len)
{
    register UInt ct;
    ct = (len + BIPEB - 1) / BIPEB;
    while (ct--) {
        *ptS++ ^= *ptV++;
    }
}

/****************************************************************************
**
*F  AddCoeffsGF2VecGF2Vec( <sum>, <vec> ) . . . . . . . .  add <vec> to <sum>
**
**  `AddCoeffsGF2VecGF2Vec' adds  the   entries of <vec>  to <sum>.    If the
**  length  are not equal the  missing entries are  assumed  to be zero.  The
**  position of the rightmost Z(2) is returned.
*/

static UInt RightMostOneGF2Vec(Obj vec)
{
    UInt len;

    len = LEN_GF2VEC(vec);
    while (0 < len) {
        if (CONST_BLOCK_ELM_GF2VEC(vec, len) == 0)
            len = BIPEB * ((len - 1) / BIPEB);
        else if (BLOCK_ELM_GF2VEC(vec, len) & MASK_POS_GF2VEC(len))
            break;
        else
            len--;
    }
    return len;
}

static Obj AddCoeffsGF2VecGF2Vec(Obj sum, Obj vec)
{
    UInt *       ptS;
    const UInt * ptV;
    UInt         len;

    // get the length
    len = LEN_GF2VEC(vec);

    // grow <sum> is necessary
    if (LEN_GF2VEC(sum) < len) {
        ResizeWordSizedBag(sum, SIZE_PLEN_GF2VEC(len));
        SET_LEN_GF2VEC(sum, len);
    }

    // add <vec> to <sum>
    ptS = BLOCKS_GF2VEC(sum);
    ptV = CONST_BLOCKS_GF2VEC(vec);
    AddGF2VecToGF2Vec(ptS, ptV, len);
    return INTOBJ_INT(RightMostOneGF2Vec(sum));
}

static inline void
CopySection_GF2Vecs(Obj src, Obj dest, UInt smin, UInt dmin, UInt nelts)
{
    UInt         soff;
    UInt         doff;
    const UInt * sptr;
    UInt *       dptr;

    // switch to zero-based indices and find the first blocks and so on
    soff = (smin - 1) % BIPEB;
    doff = (dmin - 1) % BIPEB;
    sptr = CONST_BLOCKS_GF2VEC(src) + (smin - 1) / BIPEB;
    dptr = BLOCKS_GF2VEC(dest) + (dmin - 1) / BIPEB;

    CopyBits(sptr, soff, dptr, doff, nelts);
    return;
}

/****************************************************************************
**
*F  AddPartialGF2VecGF2Vec( <sum>, <vl>, <vr>, <n> )  . . . . . . partial sum
**
**  'AddPartialGF2VecGF2Vec' adds the entries  of <vl> and <vr> starting from
**  that block which is corresponding to the entry with number <n> and stores
**  the result in <sum>.
**
**  Note: The other entries are set to be zero. So use a higher value for <n>
**        only for vectors, which both have leading zero-entries.
**
**  You  can use  the parameter  <n> for  example for  an  gauss-algorithm on
**  gf2-matrices  can be  improved,  because when  using the gauss-algorithm,
**  you  know that  the leading entries of two vectors to be  added are equal
**  zero. If <n> = 1 all entries are added.
**
**  Note that the caller has to ensure, that <sum> is a gf2-vector with the
**  correct size.
*/
static Obj AddPartialGF2VecGF2Vec(Obj sum, Obj vl, Obj vr, UInt n)
{
    const UInt * ptL;       // bit field of <vl>
    const UInt * ptR;       // bit field of <vr>
    UInt *       ptS;       // bit field of <sum>
    UInt *       end;       // end marker
    UInt         len;       // length of the list
    UInt         offset;    // number of block to start adding
    UInt         x;

    // both operands lie in the same field
    len = LEN_GF2VEC(vl);
    if (len != LEN_GF2VEC(vr)) {
        ErrorMayQuit("Vector +: vectors must have the same length", 0, 0);
    }

    // calculate the offset for adding
    if (n == 1) {
        ptL = CONST_BLOCKS_GF2VEC(vl);
        ptR = CONST_BLOCKS_GF2VEC(vr);
        ptS = BLOCKS_GF2VEC(sum);
        end = ptS + ((len + BIPEB - 1) / BIPEB);
    }
    else {
        offset = (n - 1) / BIPEB;
        ptL = CONST_BLOCKS_GF2VEC(vl) + offset;
        ptR = CONST_BLOCKS_GF2VEC(vr) + offset;
        ptS = BLOCKS_GF2VEC(sum) + offset;
        end = ptS + ((len + BIPEB - 1) / BIPEB) - offset;
    }

    // loop over the entries and add
    if (vl == sum)
        while (ptS < end) {
            // maybe remove this condition
            if ((x = *ptR) != 0)
                *ptS = *ptL ^ x;
            ptL++;
            ptS++;
            ptR++;
        }
    else if (vr == sum)
        while (ptS < end) {
            // maybe remove this condition
            if ((x = *ptL) != 0)
                *ptS = *ptR ^ x;
            ptL++;
            ptS++;
            ptR++;
        }
    else
        while (ptS < end)
            *ptS++ = *ptL++ ^ *ptR++;

    return sum;
}

/****************************************************************************
**
*F  ProdGF2VecGF2Vec( <vl>, <vr> )  . . . . . . .  product of two GF2 vectors
**
**  'ProdVecGF2VecGF2' returns  the product of  the two GF2 vectors <vl> and
**  <vr>.   The product is  the folded sum of   the corresponding entries of
**  <vl> and <vr>.
**
**  'ProdVecGF2VecGF2' is an improved  version of the general multiplication,
**  which  does not  call 'PROD'  but uses bit  operations instead.   It will
**  always return either 'GF2One' or 'GF2Zero'.
*/
#ifdef SYS_IS_64_BIT

#define PARITY_BLOCK(m)                                                      \
    do {                                                                     \
        m = m ^ (m >> 32);                                                   \
        m = m ^ (m >> 16);                                                   \
        m = m ^ (m >> 8);                                                    \
        m = m ^ (m >> 4);                                                    \
        m = m ^ (m >> 2);                                                    \
        m = m ^ (m >> 1);                                                    \
    } while (0)

#else

#define PARITY_BLOCK(m)                                                      \
    do {                                                                     \
        m = m ^ (m >> 16);                                                   \
        m = m ^ (m >> 8);                                                    \
        m = m ^ (m >> 4);                                                    \
        m = m ^ (m >> 2);                                                    \
        m = m ^ (m >> 1);                                                    \
    } while (0)

#endif

static Obj ProdGF2VecGF2Vec(Obj vl, Obj vr)
{
    const UInt * ptL;     // bit field of <vl>
    const UInt * ptR;     // bit field of <vr>
    UInt         lenL;    // length of the list
    UInt         lenR;    // length of the list
    UInt         len;     // minimum of the lengths
    UInt         nrb;     // number of whole blocks to use
    UInt         m;       // number of bits in a block
    UInt         n;       // number of bits in blist
    UInt         i;       // loop variable
    UInt         mask;    // bit selecting mask

    // both operands lie in the same field
    lenL = LEN_GF2VEC(vl);
    lenR = LEN_GF2VEC(vr);
    len = (lenL < lenR) ? lenL : lenR;

    if (len == 0) {
        ErrorMayQuit("Vector *: both vectors must have at least one entry",
                     (Int)0, (Int)0);
    }

    // loop over the entries and multiply
    ptL = CONST_BLOCKS_GF2VEC(vl);
    ptR = CONST_BLOCKS_GF2VEC(vr);
    nrb = len / BIPEB;
    n = 0;
    for (i = nrb; i > 0; i--) {
        m = (*ptL++) & (*ptR++);
        PARITY_BLOCK(m);
        n ^= m;
    }
    // now process the remaining bits

    mask = 1;
    for (i = 0; i < len % BIPEB; i++) {
        n ^= (mask & *ptL & *ptR) >> i;
        mask <<= 1;
    }

    return (n & 1) ? GF2One : GF2Zero;
}

/****************************************************************************
**
*F  ProdGF2VecGF2Mat( <vl>, <vr> )  . .  product of GF2 vector and GF2 matrix
**
**  'ProdGF2VecGF2Mat'  returns the product  of the  GF2 vector <vl>  and the
**  GF2 matrix  <vr>.   The product is   the sum of  the  rows of <vr>,  each
**  multiplied by the corresponding entry of <vl>.  Note that the  caller has
**  to ensure, that <vl> is a gf2-vector and <vr> is a gf2-matrix.
*/
static Obj ProdGF2VecGF2Mat(Obj vl, Obj vr)
{
    UInt         len;    // length of the list
    UInt         stop;
    UInt         col;     // length of the rows
    UInt         i;       // loop variables
    Obj          prod;    // product, result
    Obj          row1;    // top row of matrix
    UInt *       start;
    const UInt * ptL;
    UInt         mask;

    // both operands lie in the same field
    len = LEN_GF2VEC(vl);
    if (len > LEN_GF2MAT(vr))
        len = LEN_GF2MAT(vr);

    // make the result vector
    row1 = ELM_GF2MAT(vr, 1);
    col = LEN_GF2VEC(row1);
    NEW_GF2VEC(prod,
               (IS_MUTABLE_OBJ(vl) || IS_MUTABLE_OBJ(row1))
                   ? TYPE_LIST_GF2VEC
                   : TYPE_LIST_GF2VEC_IMM,
               col);

    // get the start and end block
    start = BLOCKS_GF2VEC(prod);
    ptL = CONST_BLOCKS_GF2VEC(vl);

    // loop over the vector
    for (i = 1; i <= len; ptL++) {

        // if the whole block is zero, get the next entry
        if (*ptL == 0) {
            i += BIPEB;
            continue;
        }

        // run through the block
        stop = i + BIPEB - 1;
        if (len < stop)
            stop = len;
        for (mask = 1; i <= stop; i++, mask <<= 1) {

            // if there is entry add the row to the result
            if ((*ptL & mask) != 0) {
                const UInt * ptRR = CONST_BLOCKS_GF2VEC(ELM_GF2MAT(vr, i));
                AddGF2VecToGF2Vec(start, ptRR, col);
            }
        }
    }

    return prod;
}

/****************************************************************************
**
*F  ProdGF2MatGF2Vec( <ml>, <vr> )  . .  product of GF2 matrix and GF2 vector
**
**  'ProdGF2MatGF2Vec'  returns the product  of the  GF2 matrix <ml>  and the
**  GF2 vector  <vr>.   The ith entry of the
**  product is the inner product of  the  ith row of <ml> with <vr>.
**  Note that the  caller has
**  to ensure, that <ml> is a GF2 matrix and <vr> is a GF2 vector.
*/
static Obj ProdGF2MatGF2Vec(Obj ml, Obj vr)
{
    UInt         len;     // length of the vector
    UInt         ln1;     // length of the rows of the mx
    UInt         ln2;     // length of the matrix
    const UInt * ptL;     // bit field of <ml>[j]
    const UInt * ptR;     // bit field of <vr>
    UInt         nrb;     // number of blocks in blist
    UInt         m;       // number of bits in a block
    UInt         n;       // number of bits in blist
    UInt         i;       // loop variable
    UInt         j;       // loop variable
    Obj          prod;    // result
    UInt         mask;    // a one bit mask

    // both operands lie in the same field
    len = LEN_GF2VEC(vr);
    ln2 = LEN_GF2MAT(ml);
    if (0 == ln2) {
        ErrorMayQuit("PROD: empty GF2 matrix * GF2 vector not allowed", 0, 0);
    }

    ln1 = LEN_GF2VEC(ELM_GF2MAT(ml, 1));
    if (len > ln1) {
        len = ln1;
    }

    // make the result vector
    NEW_GF2VEC(prod,
               (IS_MUTABLE_OBJ(ELM_GF2MAT(ml, 1)) || IS_MUTABLE_OBJ(vr))
                   ? TYPE_LIST_GF2VEC
                   : TYPE_LIST_GF2VEC_IMM,
               ln2);

    // loop over the entries and multiply
    nrb = len / BIPEB;
    for (j = 1; j <= ln2; j++) {
        ptL = CONST_BLOCKS_GF2VEC(ELM_GF2MAT(ml, j));
        ptR = CONST_BLOCKS_GF2VEC(vr);
        n = 0;
        for (i = 1; i <= nrb; i++) {
            m = (*ptL++) & (*ptR++);
            PARITY_BLOCK(m);
            n ^= m;
        }

        mask = 1;
        for (i = 0; i < len % BIPEB; i++) {
            n ^= (mask & *ptL & *ptR) >> i;
            mask <<= 1;
        }

        if (n & 1)
            BLOCK_ELM_GF2VEC(prod, j) |= MASK_POS_GF2VEC(j);
    }

    return prod;
}

/****************************************************************************
**
*F  ProdGF2MatGF2MatSimple( <ml>, <mr> ) . . . .  product of two GF2 matrices
*F  ProdGF2MatGF2MatAdvanced( <ml>, <mr>, <greaselevel>, <blocksize> )
**                                     . .  product of twp GF2 matrices
**
**  'ProdGF2MatGF2MatSimple' returns the product of the GF2 matrix <ml> and
**  the GF2 matrix <mr>. This simply calls ProdGF2VecGF2Mat once on each row.
**
**  ProdGF2MatGF2MatAdvanced uses the specified grease and blocking to
**  accelerate larger matrix multiplies. In this case, the matrix dimensions
**  must be compatible.
*/

static Obj ProdGF2MatGF2MatSimple(Obj ml, Obj mr)
{
    Obj  prod;
    UInt i;
    UInt len;
    Obj  row;
    Obj  rtype;
    len = LEN_GF2MAT(ml);
    prod = NewBag(T_POSOBJ, SIZE_PLEN_GF2MAT(len));
    SET_LEN_GF2MAT(prod, len);
    if (IS_MUTABLE_OBJ(ml) || IS_MUTABLE_OBJ(mr)) {
        SET_TYPE_POSOBJ(prod, TYPE_LIST_GF2MAT);
        if (IS_MUTABLE_OBJ(ELM_GF2MAT(ml, 1)) ||
            IS_MUTABLE_OBJ(ELM_GF2MAT(mr, 1)))
            rtype = TYPE_LIST_GF2VEC_LOCKED;
        else
            rtype = TYPE_LIST_GF2VEC_IMM_LOCKED;
    }
    else {
        SET_TYPE_POSOBJ(prod, TYPE_LIST_GF2MAT_IMM);
        rtype = TYPE_LIST_GF2VEC_IMM_LOCKED;
    }
    for (i = 1; i <= len; i++) {
        row = ProdGF2VecGF2Mat(ELM_GF2MAT(ml, i), mr);

        // Since I'm going to put this vector into a matrix, I must lock its
        // representation, so that it doesn't get rewritten over GF(2^k)
        SetTypeDatObj(row, rtype);
        SET_ELM_GF2MAT(prod, i, row);
        CHANGED_BAG(prod);
        TakeInterrupt();
    }
    return prod;
}

// Utility functions for the advanced matrix multiply code below

// extract nbits bits starting from position from in vector vptr
// return them as the nbits least significant bits in a UInt.
// Bits are always numbered least-significant first

static inline UInt getbits(const UInt * vptr, UInt from, UInt nbits)
{
    UInt wno = (from - 1) / BIPEB;
    UInt word1 = vptr[wno];
    UInt shift1 = (from - 1) % BIPEB;
    UInt lbit = shift1 + nbits;
    UInt word2;
    if (lbit <= BIPEB) {
        // range is all in one word
        word1 <<= BIPEB - lbit;
        word1 >>= BIPEB - nbits;
    }
    else {
        // range is split across two words
        word1 >>= shift1;
        lbit -= BIPEB;
        word2 = vptr[wno + 1];
        word2 <<= BIPEB - lbit;
        word2 >>= shift1 - lbit;
        word1 |= word2;
    }
    return word1;
}

// To avoid having a lot of arguments to the recursive getgreasedata function,
// we put the things that don't change in the recursive call into this
// structure

struct greaseinfo {
    UInt *        pgtags;
    UInt *        pgbuf;
    UInt          nblocks;
    UInt *        pgrules;
    const UInt ** prrows;
};

// Make if necessary the grease row for bits controlled by the data in g.
// Recursive so can't be inlined
static const UInt * getgreasedata(struct greaseinfo * g, UInt bits)
{
    UInt         x, y;
    const UInt * ps;
    UInt *       pd;
    const UInt * ps2;
    UInt         i;
    UInt *       pd1;
    switch (g->pgtags[bits]) {
    case 0:
        // Need to make the row
        x = g->pgrules[bits];
        y = bits ^ (1 << x);
        // make it by adding row x to grease vector indexed y
        ps = g->prrows[x];
        ps2 = getgreasedata(g, y);
        pd1 = g->pgbuf + (bits - 3) * g->nblocks;
        pd = pd1;
        // time critical inner loop
        for (i = g->nblocks; i > 0; i--)
            *pd++ = *ps++ ^ *ps2++;
        // record that we made it
        g->pgtags[bits] = 1;
        return pd1;

    case 1:
        // we've made this one already, so just return it
        return g->pgbuf + (bits - 3) * g->nblocks;

    case 2:
        // This one does not need making, bits actually
        // has just a single 1 bit in it
        return g->prrows[g->pgrules[bits]];
    }
    return (UInt *)0;    // can't actually get here; include the return to
                         // pacify compiler
}

static Obj
ProdGF2MatGF2MatAdvanced(Obj ml, Obj mr, UInt greasesize, UInt blocksize)
{
    Obj          prod;          // Product Matrix
    UInt         i, j, k, b;    // Loop counters
    UInt         gs;            // Actual level of grease for current block
    const UInt * rptr;          // Pointer to current row of ml
    UInt bits;    // current chunk of current row, for lookup in grease tables
    const UInt * v;          // pointer to computed grease vector
    UInt len, rlen, ilen;    // len = length of ml, ilen = row length of ml =
                             // length of mr, rlen = row length of mr
    Obj    row;    // current row of ml, or row of prod when it is being built
    Obj    rtype;              // type of rows of prod
    Obj    gbuf = (Obj)0;      // grease buffer
    Obj    gtags = (Obj)0;     // grease tags (whether that row is known yet
    Obj    grules = (Obj)0;    // rules for making new grease vectors
    UInt * pgrules;            // pointer to contents of grules
    UInt * pgtags = (UInt *)0;     // pointer to contents of gtags
    UInt * pgbuf = (UInt *)0;      // pointer to grease buffer
    UInt   nwords;                 // number of words in a row of mr
    UInt   glen;                   // 1 << greasesize
    UInt   bs;                     // actual size of current block
    UInt * pprow;                  // pointer into current row of prod
    Obj    lrowptrs;               // cache of direct pointers to rows of ml
    const UInt **     plrows;      // and a direct pointer to that cache
    Obj               rrowptrs;    // and for mr
    const UInt **     prrows;
    Obj               prowptrs;    // and for prod
    UInt **           pprows;
    struct greaseinfo g;

    len = LEN_GF2MAT(ml);
    row = ELM_GF2MAT(mr, 1);
    rlen = LEN_GF2VEC(row);
    ilen = LEN_GF2MAT(mr);
    nwords = NUMBER_BLOCKS_GF2VEC(row);

    // Make a zero product matrix
    prod = NewBag(T_POSOBJ, SIZE_PLEN_GF2MAT(len));
    SET_LEN_GF2MAT(prod, len);
    if (IS_MUTABLE_OBJ(ml) || IS_MUTABLE_OBJ(mr)) {
        SET_TYPE_POSOBJ(prod, TYPE_LIST_GF2MAT);
        if (IS_MUTABLE_OBJ(ELM_GF2MAT(ml, 1)) ||
            IS_MUTABLE_OBJ(ELM_GF2MAT(mr, 1)))
            rtype = TYPE_LIST_GF2VEC_LOCKED;
        else
            rtype = TYPE_LIST_GF2VEC_IMM_LOCKED;
    }
    else {
        SET_TYPE_POSOBJ(prod, TYPE_LIST_GF2MAT_IMM);
        rtype = TYPE_LIST_GF2VEC_IMM_LOCKED;
    }

    for (i = 1; i <= len; i++) {
        NEW_GF2VEC(row, rtype, rlen);
        SET_ELM_GF2MAT(prod, i, row);
        CHANGED_BAG(prod);
    }

    // Cap greasesize and blocksize by the actual length
    if (ilen < greasesize)
        greasesize = ilen;
    if (ilen < greasesize * blocksize)
        blocksize = (ilen + greasesize - 1) / greasesize;

    // calculate glen
    glen = 1 << greasesize;

    // Allocate memory

    lrowptrs = NewBag(T_DATOBJ, sizeof(UInt *) * len);
    rrowptrs = NewBag(T_DATOBJ, sizeof(UInt *) * ilen);
    prowptrs = NewBag(T_DATOBJ, sizeof(UInt *) * len);

    if (greasesize >= 2) {
        gbuf =
            NewBag(T_DATOBJ, sizeof(UInt) * nwords * (glen - 3) * blocksize);
        gtags = NewBag(T_DATOBJ, sizeof(UInt) * glen * blocksize);
        grules = NewBag(T_DATOBJ, sizeof(Int) * glen);

        // From here no garbage collections

        pgtags = (UInt *)ADDR_OBJ(gtags);
        pgrules = (UInt *)ADDR_OBJ(grules);
        pgbuf = (UInt *)ADDR_OBJ(gbuf);

        // Calculate the greasing rules
        for (j = 3; j < glen; j++)
            for (i = 0; i < greasesize; i++)
                if ((j & (1 << i)) != 0) {
                    pgrules[j] = i;
                    break;
                }
        for (j = 0; j < greasesize; j++)
            pgrules[1 << j] = j;

        // fill in some more bits of g
        g.pgrules = pgrules;
        g.nblocks = nwords;
    }

    // Take direct pointers to all the parts of all the matrices to avoid
    // multiple indirection overheads
    plrows = (const UInt **)ADDR_OBJ(lrowptrs);
    prrows = (const UInt **)ADDR_OBJ(rrowptrs);
    pprows = (UInt **)ADDR_OBJ(prowptrs);

    for (i = 0; i < len; i++) {
        plrows[i] = CONST_BLOCKS_GF2VEC(ELM_GF2MAT(ml, i + 1));
        pprows[i] = BLOCKS_GF2VEC(ELM_GF2MAT(prod, i + 1));
    }
    for (i = 0; i < ilen; i++)
        prrows[i] = CONST_BLOCKS_GF2VEC(ELM_GF2MAT(mr, i + 1));

    // OK, finally ready to start work
    // loop over blocks
    for (b = 1; b <= ilen; b += blocksize * greasesize) {
        // last block may be a small one
        bs = blocksize;
        if ((b + bs * greasesize) > ilen)
            bs = (ilen - b + greasesize) / greasesize;

        // If we're greasing, start afresh
        if (greasesize > 1) {
            for (k = 0; k < bs; k++) {
                for (j = 0; j < 1 << greasesize; j++)
                    pgtags[k * glen + j] = 0;
                // powers of 2 correspond to rows of mr
                for (j = 0; j < greasesize; j++)
                    pgtags[k * glen + (1 << j)] = 2;
            }
        }

        // For each block, we run through rows of ml & prod
        for (j = 1; j <= len; j++) {
            // get pointers
            rptr = plrows[j - 1];
            pprow = pprows[j - 1];

            // Now within the block, we have multiple grease-units, run
            // through them
            for (i = 0; i < bs; i++) {
                // start of current grease unit
                k = b + i * greasesize;

                // last unit of last block may be short
                gs = greasesize;
                if (k + gs > ilen)
                    gs = ilen - k + 1;

                // find the appropriate parts of grease tags
                // grease buffer and mr. Store in g

                if (gs > 1) {
                    g.pgtags = pgtags + glen * i;
                    g.pgbuf = pgbuf + (glen - 3) * nwords * i;
                    g.prrows = prrows + k - 1;
                }

                // get a chunk from a row of ml
                bits = getbits(rptr, k, gs);

                // 0 means nothing to do
                if (bits == 0)
                    continue;
                else if (bits == 1)    // handle this one specially to speed
                                       // up the greaselevel 1 case
                    v = prrows[k - 1];    // -1 is because k is 1-based index
                else
                    v = getgreasedata(
                        &g, bits);    // The main case
                                      // This function should be inlined
                AddGF2VecToGF2Vec(pprow, v, rlen);
            }
        }

        // Allow GAP to respond to Ctrl-C
        if (TakeInterrupt()) {
            // Might have been a garbage collection, reload everything
            if (greasesize >= 2) {
                pgtags = (UInt *)ADDR_OBJ(gtags);
                pgrules = (UInt *)ADDR_OBJ(grules);
                pgbuf = (UInt *)ADDR_OBJ(gbuf);
                // fill in some more bits of g
                g.pgrules = pgrules;
                g.nblocks = nwords;
            }
            plrows = (const UInt **)ADDR_OBJ(lrowptrs);
            prrows = (const UInt **)ADDR_OBJ(rrowptrs);
            pprows = (UInt **)ADDR_OBJ(prowptrs);
            for (i = 0; i < len; i++) {
                plrows[i] = CONST_BLOCKS_GF2VEC(ELM_GF2MAT(ml, i + 1));
                pprows[i] = BLOCKS_GF2VEC(ELM_GF2MAT(prod, i + 1));
            }
            for (i = 0; i < ilen; i++)
                prrows[i] = CONST_BLOCKS_GF2VEC(ELM_GF2MAT(mr, i + 1));
        }
    }
    return prod;
}

/****************************************************************************
**
*F  FuncPROD_GF2VEC_ANYMAT( <self>, <v>, <m>)
**
**  method to handle vector*plain list of GF2Vectors reasonably efficiently.
*/
static Obj FuncPROD_GF2VEC_ANYMAT(Obj self, Obj vec, Obj mat)
{
    Obj  res;
    UInt len;
    UInt len1;
    Obj  row1;
    UInt i;
    UInt block = 0;

    len = LEN_GF2VEC(vec);
    if (len > LEN_PLIST(mat))
        len = LEN_PLIST(mat);

    // Get the first row, to establish the size of the result
    row1 = ELM_PLIST(mat, 1);
    if (!IS_GF2VEC_REP(row1))
        return TRY_NEXT_METHOD;
    len1 = LEN_GF2VEC(row1);

    // create the result space
    NEW_GF2VEC(res,
               (IS_MUTABLE_OBJ(vec) || IS_MUTABLE_OBJ(row1))
                   ? TYPE_LIST_GF2VEC
                   : TYPE_LIST_GF2VEC_IMM,
               len1);

    // Finally, we start work
    for (i = 1; i <= len; i++) {
        if (i % BIPEB == 1)
            block = CONST_BLOCK_ELM_GF2VEC(vec, i);
        if (block & MASK_POS_GF2VEC(i)) {
            row1 = ELM_PLIST(mat, i);
            if (!IS_GF2VEC_REP(row1))
                return TRY_NEXT_METHOD;
            AddPartialGF2VecGF2Vec(res, res, row1, 1);
        }
    }
    return res;
}

/****************************************************************************
**
*F  InversePlistGF2VecsDesstructive( <list> )
**
**  This is intended to form the core of a method for InverseOp.
**  by this point it should be checked that list is a plain list of GF2
**  vectors of equal lengths.
*/
static Obj InversePlistGF2VecsDesstructive(Obj list)
{
    UInt         len;     // dimension
    Obj          inv;     // result
    Obj          row;     // row vector
    Obj          old;     // row from <mat>
    Obj          tmp;     // temporary
    UInt *       ptQ;     // data block of <row>
    const UInt * ptP;     // data block of source row
    const UInt * end;     // end marker
    const UInt * end2;    // end marker
    UInt         i;       // loop variable
    UInt         k;       // loop variable

    len = LEN_PLIST(list);

    // create the identity matrix
    tmp = NEW_PLIST(T_PLIST, len);
    for (i = len; 0 < i; i--) {
        NEW_GF2VEC(row, TYPE_LIST_GF2VEC, len);
        BLOCK_ELM_GF2VEC(row, i) |= MASK_POS_GF2VEC(i);
        SET_ELM_PLIST(tmp, i, row);
        CHANGED_BAG(tmp);
    }
    SET_LEN_PLIST(tmp, len);
    inv = tmp;

    // now start with ( id | mat ) towards ( inv | id )
    for (k = 1; k <= len; k++) {

        // find a nonzero entry in column <k>
        for (i = k; i <= len; i++) {
            row = ELM_PLIST(list, i);
            if (CONST_BLOCK_ELM_GF2VEC(row, k) & MASK_POS_GF2VEC(k))
                break;
        }
        if (i > len) {
            return Fail;
        }
        if (i != k) {
            row = ELM_PLIST(list, i);
            SET_ELM_PLIST(list, i, ELM_PLIST(list, k));
            SET_ELM_PLIST(list, k, row);
            row = ELM_PLIST(inv, i);
            SET_ELM_PLIST(inv, i, ELM_PLIST(inv, k));
            SET_ELM_PLIST(inv, k, row);
        }

        // clear entries
        old = ELM_PLIST(list, k);
        end = CONST_BLOCKS_GF2VEC(old) + ((len + BIPEB - 1) / BIPEB);
        for (i = 1; i <= len; i++) {
            if (i == k)
                continue;
            row = ELM_PLIST(list, i);
            if (CONST_BLOCK_ELM_GF2VEC(row, k) & MASK_POS_GF2VEC(k)) {

                // clear <mat>
                ptQ = &(BLOCK_ELM_GF2VEC(row, k));
                ptP = &(CONST_BLOCK_ELM_GF2VEC(old, k));
                while (ptP < end) {
                    *ptQ++ ^= *ptP++;
                }

                // modify <inv>
                row = ELM_PLIST(inv, i);
                ptQ = BLOCKS_GF2VEC(row);
                row = ELM_PLIST(inv, k);
                ptP = CONST_BLOCKS_GF2VEC(row);
                end2 = ptP + ((len + BIPEB - 1) / BIPEB);
                while (ptP < end2) {
                    *ptQ++ ^= *ptP++;
                }
            }
        }
        TakeInterrupt();
    }
    return inv;
}

/****************************************************************************
**
*F  InverseGF2Mat( <mat> )  . . . . . . . . . . . . . . inverse of GF2 matrix
**
**  This should be improved to work with mutable GF2 matrices
**
*/

static Obj InverseGF2Mat(Obj mat, UInt mut)
{
    UInt         len;    // dimension
    Obj          inv;    // result
    Obj          row;    // row vector
    Obj          tmp;    // temporary
    UInt         i;      // loop variable
    Obj          old;    // row from <mat>
    const UInt * ptQ;    // data block of <row>
    UInt *       ptP;    // data block of source row
    UInt *       end;    // end marker
    Obj          rtype;

    // make a structural copy of <mat> as list of GF2 vectors
    len = LEN_GF2MAT(mat);

    // special routes for very small matrices
    if (len == 0) {
        return CopyObj(mat, 1);
    }
    if (len == 1) {
        row = ELM_GF2MAT(mat, 1);
        if (CONST_BLOCKS_GF2VEC(row)[0] & 1) {
            return CopyObj(mat, 1);
        }
        else
            return Fail;
    }

    tmp = NEW_PLIST(T_PLIST, len);
    for (i = len; 0 < i; i--) {
        old = ELM_GF2MAT(mat, i);
        NEW_GF2VEC(row, TYPE_LIST_GF2VEC_IMM, len);
        ptQ = CONST_BLOCKS_GF2VEC(old);
        ptP = BLOCKS_GF2VEC(row);
        end = ptP + ((len + BIPEB - 1) / BIPEB);
        while (ptP < end)
            *ptP++ = *ptQ++;
        SET_ELM_PLIST(tmp, i, row);
        CHANGED_BAG(tmp);
    }
    SET_LEN_PLIST(tmp, len);
    inv = InversePlistGF2VecsDesstructive(tmp);
    if (inv == Fail)
        return inv;

    // convert list <inv> into a matrix
    ResizeBag(inv, SIZE_PLEN_GF2MAT(len));
    if (mut == 2 || (mut == 1 && IS_MUTABLE_OBJ(mat) &&
                     IS_MUTABLE_OBJ(ELM_GF2MAT(mat, 1))))
        rtype = TYPE_LIST_GF2VEC_LOCKED;
    else
        rtype = TYPE_LIST_GF2VEC_IMM_LOCKED;
    for (i = len; 0 < i; i--) {
        row = ELM_PLIST(inv, i);
        SET_TYPE_POSOBJ(row, rtype);
        SET_ELM_GF2MAT(inv, i, row);
    }
    SET_LEN_GF2MAT(inv, len);
    RetypeBag(inv, T_POSOBJ);
    SET_TYPE_POSOBJ(inv, (mut == 2 || (mut == 1 && IS_MUTABLE_OBJ(mat)))
                             ? TYPE_LIST_GF2MAT
                             : TYPE_LIST_GF2MAT_IMM);
    return inv;
}

/****************************************************************************
**
*F  ShallowCopyVecGF2( <vec> )
**
*/

Obj ShallowCopyVecGF2(Obj vec)
{
    Obj          copy;
    UInt         len;
    const UInt * ptrS;
    UInt *       ptrD;
    len = LEN_GF2VEC(vec);
    NEW_GF2VEC(copy, TYPE_LIST_GF2VEC, len);
    ptrS = CONST_BLOCKS_GF2VEC(vec);
    ptrD = BLOCKS_GF2VEC(copy);
    memcpy(ptrD, ptrS, NUMBER_BLOCKS_GF2VEC(vec) * sizeof(UInt));
    return copy;
}

/****************************************************************************
**
*F  SemiEchelonPlistGF2Vecs( <mat>, <transformations-needed> )
**
**  The matrix needs to have mutable rows, so it can't be a GF2 mat
**
**  This has changed. There should now be a method for mutable GF2mats as
**  well.
**
**  This function DOES NOT CHECK that the rows are all GF2 vectors
**
**  Does not copy the matrix, may destroy it, may include some
**  of the rows among the returned vectors
*/

static UInt RNheads, RNvectors, RNcoeffs, RNrelns;

static Obj SemiEchelonListGF2Vecs(Obj mat, UInt TransformationsNeeded)
{
    UInt  nrows, ncols;
    UInt  i, j, h;
    Obj   heads, vectors, coeffs = 0, relns = 0;
    UInt  nvecs, nrels = 0;
    Obj   coeffrow = 0;
    Obj   row;
    UInt *rowp, *coeffrowp = 0;
    Obj   res;
    nrows = LEN_PLIST(mat);
    ncols = LEN_GF2VEC(ELM_PLIST(mat, 1));
    heads = NEW_PLIST(T_PLIST_CYC, ncols);
    SET_LEN_PLIST(heads, ncols);
    vectors = NEW_PLIST(T_PLIST_TAB_RECT, nrows);
    nvecs = 0;
    if (TransformationsNeeded) {
        coeffs = NEW_PLIST(T_PLIST_TAB_RECT, nrows);
        relns = NEW_PLIST(T_PLIST_TAB_RECT, nrows);
        nrels = 0;
    }
    for (i = 1; i <= ncols; i++)
        SET_ELM_PLIST(heads, i, INTOBJ_INT(0));
    for (i = 1; i <= nrows; i++) {
        row = ELM_PLIST(mat, i);
        if (TransformationsNeeded) {
            NEW_GF2VEC(coeffrow, TYPE_LIST_GF2VEC, nrows);
            BLOCK_ELM_GF2VEC(coeffrow, i) |= MASK_POS_GF2VEC(i);
        }

        // No garbage collection risk from here
        rowp = BLOCKS_GF2VEC(row);
        if (TransformationsNeeded)
            coeffrowp = BLOCKS_GF2VEC(coeffrow);
        for (j = 1; j <= ncols; j++) {
            h = INT_INTOBJ(ELM_PLIST(heads, j));
            if (h != 0) {
                if (rowp[(j - 1) / BIPEB] & MASK_POS_GF2VEC(j)) {
                    AddGF2VecToGF2Vec(
                        rowp, CONST_BLOCKS_GF2VEC(ELM_PLIST(vectors, h)),
                        ncols);
                    if (TransformationsNeeded)
                        AddGF2VecToGF2Vec(
                            coeffrowp,
                            CONST_BLOCKS_GF2VEC(ELM_PLIST(coeffs, h)), nrows);
                }
            }
        }
        j = 1;
        while (j <= ncols && !*rowp) {
            j += BIPEB;
            rowp++;
        }
        while (j <= ncols && !(*rowp & MASK_POS_GF2VEC(j)))
            j++;

        // garbage collection OK again after here
        if (j <= ncols) {
            SET_ELM_PLIST(vectors, ++nvecs, row);
            CHANGED_BAG(vectors);    // Could be an old bag by now. Max.
            SET_LEN_PLIST(vectors, nvecs);
            SET_ELM_PLIST(heads, j, INTOBJ_INT(nvecs));
            if (TransformationsNeeded) {
                SET_ELM_PLIST(coeffs, nvecs, coeffrow);
                CHANGED_BAG(coeffs);    // Could be an old bag by now. Max.
                SET_LEN_PLIST(coeffs, nvecs);
            }
        }
        else if (TransformationsNeeded) {
            SET_ELM_PLIST(relns, ++nrels, coeffrow);
            CHANGED_BAG(relns);    // Could be an old bag by now. Max.
            SET_LEN_PLIST(relns, nrels);
        }
        TakeInterrupt();
    }
    if (RNheads == 0) {
        RNheads = RNamName("heads");
        RNvectors = RNamName("vectors");
    }
    res = NEW_PREC(TransformationsNeeded ? 4 : 2);
    AssPRec(res, RNheads, heads);
    AssPRec(res, RNvectors, vectors);
    if (LEN_PLIST(vectors) == 0)
        RetypeBag(vectors, T_PLIST_EMPTY);
    if (TransformationsNeeded) {
        if (RNcoeffs == 0) {
            RNcoeffs = RNamName("coeffs");
            RNrelns = RNamName("relations");
        }
        AssPRec(res, RNcoeffs, coeffs);
        if (LEN_PLIST(coeffs) == 0)
            RetypeBag(coeffs, T_PLIST_EMPTY);
        AssPRec(res, RNrelns, relns);
        if (LEN_PLIST(relns) == 0)
            RetypeBag(relns, T_PLIST_EMPTY);
    }
    SortPRecRNam(res);
    return res;
}

/****************************************************************************
**
*F  UInt TriangulizeListGF2Vecs( <mat>, <clearup> ) -- returns the rank
**
**  Again should add a method to work with mutable GF2 matrices
**
*/

static UInt TriangulizeListGF2Vecs(Obj mat, UInt clearup)
{
    UInt         nrows;
    UInt         ncols;
    UInt         workcol;
    UInt         workrow;
    UInt         rank;
    Obj          row, row2;
    const UInt * rowp;
    UInt *       row2p;
    UInt         block;
    UInt         mask;
    UInt         j;
    nrows = LEN_PLIST(mat);
    ncols = LEN_GF2VEC(ELM_PLIST(mat, 1));
    rank = 0;

    // Nothing here can cause a garbage collection

    for (workcol = 1; workcol <= ncols; workcol++) {
        block = (workcol - 1) / BIPEB;
        mask = MASK_POS_GF2VEC(workcol);
        for (workrow = rank + 1;
             workrow <= nrows &&
             !(CONST_BLOCKS_GF2VEC(ELM_PLIST(mat, workrow))[block] & mask);
             workrow++)
            ;
        if (workrow <= nrows) {
            rank++;
            row = ELM_PLIST(mat, workrow);
            if (workrow != rank) {
                SET_ELM_PLIST(mat, workrow, ELM_PLIST(mat, rank));
                SET_ELM_PLIST(mat, rank, row);
            }
            rowp = CONST_BLOCKS_GF2VEC(row);
            if (clearup)
                for (j = 1; j < rank; j++) {
                    row2 = ELM_PLIST(mat, j);
                    row2p = BLOCKS_GF2VEC(row2);
                    if (row2p[block] & mask)
                        AddGF2VecToGF2Vec(row2p, rowp, ncols);
                }
            for (j = workrow + 1; j <= nrows; j++) {
                row2 = ELM_PLIST(mat, j);
                row2p = BLOCKS_GF2VEC(row2);
                if (row2p[block] & mask)
                    AddGF2VecToGF2Vec(row2p, rowp, ncols);
            }
        }
        TakeInterrupt();
    }
    return rank;
}

/****************************************************************************
**
*F * * * * * * * * * * * *  conversion functions  * * * * * * * * * * * * * *
*/

/****************************************************************************
**
*F  PlainGF2Vec( <list> ) . . .  . convert a GF2 vector into an ordinary list
**
**  'PlainGF2Vec' converts the GF2 vector <list> to a plain list.
*/
static Obj IsLockedRepresentationVector;

static void PlainGF2Vec(Obj list)
{
    Int  len;          // length of <list>
    UInt i;            // loop variable
    Obj  first = 0;    // first entry

    // check for representation lock
    if (True == DoFilter(IsLockedRepresentationVector, list))
        ErrorMayQuit("Cannot convert a locked GF2 vector into a plain list",
                     0, 0);

    // resize the list and retype it, in this order
    len = LEN_GF2VEC(list);

    RetypeBagSM(list, (len == 0) ? T_PLIST_EMPTY : T_PLIST_FFE);

    GROW_PLIST(list, (UInt)len);
    SET_LEN_PLIST(list, len);

    // keep the first entry because setting the second destroys the first
    if (len == 0)
        SET_ELM_PLIST(list, 1, 0);
    else
        first = ELM_GF2VEC(list, 1);

    // wipe out the first entry of the GF2 vector (which becomes the  second
    // entry of the plain list, in case the list has length 1.
    if (len == 1)
        SET_ELM_PLIST(list, 2, 0);

    // replace the bits by 'GF2One' or 'GF2Zero' as the case may be
    // this must of course be done from the end of the list backwards
    for (i = len; 1 < i; i--)
        SET_ELM_PLIST(list, i, ELM_GF2VEC(list, i));
    if (len != 0)
        SET_ELM_PLIST(list, 1, first);

    CHANGED_BAG(list);
}

/****************************************************************************
**
*F  PlainGF2Mat( <list> ) . . .  . convert a GF2 matrix into an ordinary list
**
**  'PlainGF2Mat' converts the GF2 matrix <list> to a plain list.
*/
static void PlainGF2Mat(Obj list)
{
    Int  len;    // length of <list>
    UInt i;      // loop variable

    // resize the list and retype it, in this order
    len = LEN_GF2MAT(list);
    RetypeBagSM(list, T_PLIST);
    SET_LEN_PLIST(list, len);

    // shift the entries to the left
    for (i = 1; i <= len; i++) {
        SET_ELM_PLIST(list, i, ELM_GF2MAT(list, i));
    }
    SHRINK_PLIST(list, len);
    CHANGED_BAG(list);
}

/****************************************************************************
**
*F  ConvGF2Vec( <list> )  . . . . . . convert a list into a GF2 vector object
*/
static void ConvGF2Vec(Obj list)
{
    Int  len;      // logical length of the vector
    Int  i;        // loop variable
    UInt block;    // one block of the boolean list
    UInt bit;      // one bit of a block
    Obj  x;

    // already in the correct representation
    if (IS_GF2VEC_REP(list)) {
        return;
    }

    // Otherwise make it a plain list so that we will know where it keeps
    // its data -- could do much better in the case of GF(2^n) vectors that
    // actually lie over GF(2)

    if (IS_VEC8BIT_REP(list))
        PlainVec8Bit(list);
    else
        PLAIN_LIST(list);

    // change its representation
    len = LEN_PLIST(list);

    // We may have to resize the bag now because a length 1
    // plain list is shorter than a length 1 GF2VEC
    if (SIZE_PLEN_GF2VEC(len) > SIZE_OBJ(list))
        ResizeBag(list, SIZE_PLEN_GF2VEC(len));

    // now do the work
    block = 0;
    bit = 1;
    for (i = 1; i <= len; i++) {
        x = ELM_PLIST(list, i);
        if (x == GF2One)
            block |= bit;
        else if (x != GF2Zero) {
            // might be GF(2) elt written over bigger field
            if (EQ(x, GF2One))
                block |= bit;
            else if (!EQ(x, GF2Zero))
                ErrorMayQuit(
                    "COPY_GF2VEC: argument must be a list of GF2 elements",
                    0, 0);
        }

        bit = bit << 1;
        if (bit == 0 || i == len) {
            BLOCK_ELM_GF2VEC(list, i) = block;
            block = 0;
            bit = 1;
        }
    }

    // retype and resize bag
    ResizeWordSizedBag(list, SIZE_PLEN_GF2VEC(len));
    SET_LEN_GF2VEC(list, len);
    if (IS_PLIST_MUTABLE(list)) {
        SetTypeDatObj(list, TYPE_LIST_GF2VEC);
    }
    else {
        SetTypeDatObj(list, TYPE_LIST_GF2VEC_IMM);
    }
    RetypeBag(list, T_DATOBJ);
}

/****************************************************************************
**
*F  FuncCONV_GF2VEC( <self>, <list> ) . . . . . convert into a GF2 vector rep
*/
static Obj FuncCONV_GF2VEC(Obj self, Obj list)
{
    ConvGF2Vec(list);
    return 0;
}

/****************************************************************************
**
*F  NewGF2Vec( <list> )  . . . . . . convert a list into a GF2 vector object
**
**  This is a non-destructive counterpart of ConvGF2Vec
*/
static Obj NewGF2Vec(Obj list)
{
    Int  len;      // logical length of the vector
    Int  i;        // loop variable
    UInt block;    // one block of the boolean list
    UInt bit;      // one bit of a block
    Obj  x;
    Obj  res;    // resulting GF2 vector object

    // already in the correct representation
    if (IS_GF2VEC_REP(list)) {
        res = ShallowCopyVecGF2(list);
        if (!IS_MUTABLE_OBJ(list))
            SetTypeDatObj(res, TYPE_LIST_GF2VEC_IMM);
        return res;
    }

    if (!IS_LIST(list)) {
        ErrorMayQuit("COPY_GF2VEC: argument must be a list of GF2 elements",
                     0, 0);
    }
    if (!IS_PLIST(list)) {
        list = SHALLOW_COPY_OBJ(list);
        // TODO: if list is in 8bit rep, we could do better
        if (IS_VEC8BIT_REP(list))
            PlainVec8Bit(list);
        else
            PLAIN_LIST(list);
    }

    len = LEN_PLIST(list);
    NEW_GF2VEC(res, TYPE_LIST_GF2VEC, len);

    // now do the work
    block = 0;
    bit = 1;
    for (i = 1; i <= len; i++) {
        x = ELM_PLIST(list, i);
        if (x == GF2One)
            block |= bit;
        else if (x != GF2Zero) {
            // might be GF(2) elt written over bigger field
            if (EQ(x, GF2One))
                block |= bit;
            else if (!EQ(x, GF2Zero))
                ErrorMayQuit(
                    "COPY_GF2VEC: argument must be a list of GF2 elements",
                    0, 0);
        }

        bit = bit << 1;
        if (bit == 0 || i == len) {
            BLOCK_ELM_GF2VEC(res, i) = block;    // only changed list to res
            block = 0;
            bit = 1;
        }
    }

    // mutability should be inherited from the argument
    if (IS_PLIST_MUTABLE(list))
        SetTypeDatObj(res, TYPE_LIST_GF2VEC);
    else
        SetTypeDatObj(res, TYPE_LIST_GF2VEC_IMM);

    return res;
}

/****************************************************************************
**
*F  FuncCOPY_GF2VEC( <self>, <list> ) . . . . . convert into a GF2 vector rep
**
**  This is a non-destructive counterpart of FuncCONV_GF2VEC
*/
static Obj FuncCOPY_GF2VEC(Obj self, Obj list)
{
    list = NewGF2Vec(list);
    return list;
}

/****************************************************************************
**
*F FuncCONV_GF2MAT (<self>, <list> ) . . . convert into a GF2 matrix rep
**
** <list> should be a list of compressed GF2 vectors
**
*/
static Obj FuncCONV_GF2MAT(Obj self, Obj list)
{
    UInt len, i;
    Obj  tmp;
    UInt mut;
    len = LEN_LIST(list);
    if (len == 0)
        return (Obj)0;

    PLAIN_LIST(list);
    GROW_PLIST(list, len + 1);
    for (i = len; i > 0; i--) {
        tmp = ELM_PLIST(list, i);
        if (!IS_GF2VEC_REP(tmp)) {
            int j;
            for (j = i + 1; j <= len; j++) {
                tmp = ELM_PLIST(list, j + 1);
                SET_ELM_PLIST(list, j, tmp);
            }
            ErrorMayQuit("CONV_GF2MAT: argument must be a list of compressed "
                         "GF2 vectors",
                         0, 0);
        }
        SetTypeDatObj(tmp, IS_MUTABLE_OBJ(tmp) ? TYPE_LIST_GF2VEC_LOCKED
                                               : TYPE_LIST_GF2VEC_IMM_LOCKED);
        SET_ELM_PLIST(list, i + 1, tmp);
    }
    SET_ELM_PLIST(list, 1, INTOBJ_INT(len));
    mut = IS_PLIST_MUTABLE(list);
    RetypeBag(list, T_POSOBJ);
    SET_TYPE_POSOBJ(list, mut ? TYPE_LIST_GF2MAT : TYPE_LIST_GF2MAT_IMM);
    return (Obj)0;
}

/****************************************************************************
**
*F  FuncPLAIN_GF2VEC( <self>, <list> ) . . .  convert back into ordinary list
*/
static Obj FuncPLAIN_GF2VEC(Obj self, Obj list)
{
    if (!IS_GF2VEC_REP(list)) {
        RequireArgument(SELF_NAME, list, "must be a GF2 vector");
    }
    PlainGF2Vec(list);
    return 0;
}

/****************************************************************************
**
*F  revertbits -- utility function to reverse bit orders
*/

// A list of flip values for bytes (i.e. ..xyz -> zyx..)
static const UInt1 revertlist[] = {
    0,  128, 64, 192, 32, 160, 96,  224, 16, 144, 80, 208, 48, 176, 112, 240,
    8,  136, 72, 200, 40, 168, 104, 232, 24, 152, 88, 216, 56, 184, 120, 248,
    4,  132, 68, 196, 36, 164, 100, 228, 20, 148, 84, 212, 52, 180, 116, 244,
    12, 140, 76, 204, 44, 172, 108, 236, 28, 156, 92, 220, 60, 188, 124, 252,
    2,  130, 66, 194, 34, 162, 98,  226, 18, 146, 82, 210, 50, 178, 114, 242,
    10, 138, 74, 202, 42, 170, 106, 234, 26, 154, 90, 218, 58, 186, 122, 250,
    6,  134, 70, 198, 38, 166, 102, 230, 22, 150, 86, 214, 54, 182, 118, 246,
    14, 142, 78, 206, 46, 174, 110, 238, 30, 158, 94, 222, 62, 190, 126, 254,
    1,  129, 65, 193, 33, 161, 97,  225, 17, 145, 81, 209, 49, 177, 113, 241,
    9,  137, 73, 201, 41, 169, 105, 233, 25, 153, 89, 217, 57, 185, 121, 249,
    5,  133, 69, 197, 37, 165, 101, 229, 21, 149, 85, 213, 53, 181, 117, 245,
    13, 141, 77, 205, 45, 173, 109, 237, 29, 157, 93, 221, 61, 189, 125, 253,
    3,  131, 67, 195, 35, 163, 99,  227, 19, 147, 83, 211, 51, 179, 115, 243,
    11, 139, 75, 203, 43, 171, 107, 235, 27, 155, 91, 219, 59, 187, 123, 251,
    7,  135, 71, 199, 39, 167, 103, 231, 23, 151, 87, 215, 55, 183, 119, 247,
    15, 143, 79, 207, 47, 175, 111, 239, 31, 159, 95, 223, 63, 191, 127, 255
};

// Takes an UInt a on n bits and returns the Uint obtained by reverting the
// bits
static UInt revertbits(UInt a, Int n)
{
    UInt b, c;
    b = 0;
    while (n > 8) {
        c = a & 0xff;    // last byte
        a = a >> 8;
        b = b << 8;
        b += (UInt)revertlist[(UInt1)c];    // add flipped
        n -= 8;
    }
    // cope with the last n bits
    a &= 0xff;
    b = b << n;
    c = (UInt)revertlist[(UInt1)a];
    c = c >> (8 - n);
    b += c;
    return b;
}

/****************************************************************************
**
*F  Cmp_GF2Vecs( <vl>, <vr> )   compare GF2 vectors -- internal
**                                    returns -1, 0 or 1
*/
static Int Cmp_GF2VEC_GF2VEC(Obj vl, Obj vr)
{
    UInt         i;                  // loop variable
    const UInt * bl;                 // block of <vl>
    const UInt * br;                 // block of <vr>
    UInt         len, lenl, lenr;    // length of the list
    UInt         a, b, nb;

    // get and check the length
    lenl = LEN_GF2VEC(vl);
    lenr = LEN_GF2VEC(vr);
    nb = NUMBER_BLOCKS_GF2VEC(vl);
    a = NUMBER_BLOCKS_GF2VEC(vr);
    if (a < nb) {
        nb = a;
    }

    // check all blocks
    bl = CONST_BLOCKS_GF2VEC(vl);
    br = CONST_BLOCKS_GF2VEC(vr);
    for (i = nb; 1 < i; i--, bl++, br++) {
        // comparison is numeric of the reverted lists
        if (*bl != *br) {
            a = revertbits(*bl, BIPEB);
            b = revertbits(*br, BIPEB);
            if (a < b)
                return -1;
            else
                return 1;
        }
    }

    // The last block remains
    len = lenl;
    if (len > lenr) {
        len = lenr;
    }

    // are both vectors length 0?
    if (len == 0)
        return 0;

    // is there still a full block in common?
    len = len % BIPEB;
    if (len == 0) {
        a = revertbits(*bl, BIPEB);
        b = revertbits(*br, BIPEB);
    }
    else {
        a = revertbits(*bl, len);
        b = revertbits(*br, len);
    }

    if (a < b)
        return -1;
    if (a > b)
        return 1;

    // blocks still the same --left length must be smaller to be true
    if (lenr > lenl)
        return -1;
    if (lenl > lenr)
        return 1;

    return 0;
}

/****************************************************************************
**
*F  FuncEQ_GF2VEC_GF2VEC( <self>, <vl>, <vr> )   test equality of GF2 vectors
*/
static Obj FuncEQ_GF2VEC_GF2VEC(Obj self, Obj vl, Obj vr)
{
    // we can do this case MUCH faster if we just want equality
    if (LEN_GF2VEC(vl) != LEN_GF2VEC(vr))
        return False;
    return (Cmp_GF2VEC_GF2VEC(vl, vr) == 0) ? True : False;
}

/****************************************************************************
**
*F  FuncLEN_GF2VEC( <self>, <list> )  . . . . . . . .  length of a GF2 vector
*/
static Obj FuncLEN_GF2VEC(Obj self, Obj list)
{
    return INTOBJ_INT(LEN_GF2VEC(list));
}

/****************************************************************************
**
*F  FuncELM0_GF2VEC( <self>, <list>, <pos> )  . select an elm of a GF2 vector
**
**  'ELM0_GF2VEC'  returns the element at the  position  <pos> of the boolean
**  list <list>, or `Fail' if <list> has no assigned  object at <pos>.  It is
**  the  responsibility of  the caller to   ensure  that <pos> is  a positive
**  integer.
*/
static Obj FuncELM0_GF2VEC(Obj self, Obj list, Obj pos)
{
    UInt p = GetSmallInt(SELF_NAME, pos);
    if (LEN_GF2VEC(list) < p) {
        return Fail;
    }
    else {
        return ELM_GF2VEC(list, p);
    }
}

/****************************************************************************
**
*F  FuncELM_GF2VEC( <self>, <list>, <pos> ) . . select an elm of a GF2 vector
**
**  'ELM_GF2VEC' returns the element at the position <pos>  of the GF2 vector
**  <list>.   An  error  is signalled  if  <pos>  is  not bound.    It is the
**  responsibility of the caller to ensure that <pos> is a positive integer.
*/
static Obj FuncELM_GF2VEC(Obj self, Obj list, Obj pos)
{
    UInt p = GetSmallInt(SELF_NAME, pos);
    if (LEN_GF2VEC(list) < p) {
        ErrorMayQuit("List Element: [%d] must have an assigned value",
                     p, 0);
    }
    else {
        return ELM_GF2VEC(list, p);
    }
}

/****************************************************************************
**
*F  FuncELMS_GF2VEC( <self>, <list>, <poss> ) . . . sublist from a GF2 vector
**
**  'ELMS_GF2VEC' returns a new list containing  the elements at the position
**  given in    the   list  <poss> from  the   vector   <list>.   It  is  the
**  responsibility of the caller to ensure that <poss>  is dense and contains
**  only positive integers.  An error is signalled if an element of <poss> is
**  larger than the length of <list>.
*/
static Obj FuncELMS_GF2VEC(Obj self, Obj list, Obj poss)
{
    Obj elms;       // selected sublist, result
    Int lenList;    // length of <list>
    Int lenPoss;    // length of positions
    Int pos;        // position as integer
    Int inc;        // increment in a range
    Int i;          // loop variable
    Obj apos;

    // get the length of <list>
    lenList = LEN_GF2VEC(list);

    // general code for arbitrary lists, which are ranges
    if (!IS_RANGE(poss)) {

        // get the length of <positions>
        lenPoss = LEN_LIST(poss);

        // make the result vector
        NEW_GF2VEC(elms, TYPE_LIST_GF2VEC, lenPoss);

        // loop over the entries of <positions> and select
        for (i = 1; i <= lenPoss; i++) {

            // get next position

            apos = ELM0_LIST(poss, i);
            if (!apos || !IS_INTOBJ(apos))
                ErrorMayQuit("ELMS_GF2VEC: error at position %d in positions "
                             "list, entry must be bound to a small integer",
                             i, 0);
            pos = INT_INTOBJ(apos);
            if (lenList < pos) {
                ErrorMayQuit("List Elements: [%d] must have a value",
                             pos, 0);
            }

            // assign the element into <elms>
            if (ELM_GF2VEC(list, pos) == GF2One) {
                BLOCK_ELM_GF2VEC(elms, i) |= MASK_POS_GF2VEC(i);
            }
        }
    }

    // special code for ranges
    else {

        // get the length of <positions>, the first elements, and the inc.
        lenPoss = GET_LEN_RANGE(poss);
        pos = GET_LOW_RANGE(poss);
        inc = GET_INC_RANGE(poss);

        // check that no <position> is larger than <lenList>
        if (lenList < pos) {
            ErrorMayQuit("List Elements: [%d] must have a value", pos,
                         0);
        }
        if (lenList < pos + (lenPoss - 1) * inc) {
            ErrorMayQuit("List Elements: [%d] must have a value",
                         pos + (lenPoss - 1) * inc, 0);
        }

        // make the result vector
        NEW_GF2VEC(elms, TYPE_LIST_GF2VEC, lenPoss);

        // increment 1 ranges is a block copy
        if (inc == 1)
            CopySection_GF2Vecs(list, elms, pos, 1, lenPoss);

        // loop over the entries of <positions> and select
        else {
            for (i = 1; i <= lenPoss; i++, pos += inc) {
                if (ELM_GF2VEC(list, pos) == GF2One) {
                    BLOCK_ELM_GF2VEC(elms, i) |= MASK_POS_GF2VEC(i);
                }
            }
        }
    }

    return elms;
}

/****************************************************************************
**
*F  FuncASS_GF2VEC( <self>, <list>, <pos>, <elm> ) set an elm of a GF2 vector
**
**  'ASS_GF2VEC' assigns the element  <elm> at the position  <pos> to the GF2
**  vector <list>.
**
**  It is the responsibility of the caller  to ensure that <pos> is positive,
**  and that <elm> is not 0.
*/

static Obj FuncASS_GF2VEC(Obj self, Obj list, Obj pos, Obj elm)
{
    // check that <list> is mutable
    RequireMutable("List Assignment", list, "list");

    // get the position
    UInt p = GetSmallInt(SELF_NAME, pos);

    // if <elm> is Z(2) or 0*Z(2) and the position is OK, keep rep
    if (p <= LEN_GF2VEC(list) + 1) {
        if (LEN_GF2VEC(list) + 1 == p) {
            if (DoFilter(IsLockedRepresentationVector, list) == True)
                ErrorMayQuit("Assignment forbidden beyond the end of locked "
                             "GF2 vector",
                             0, 0);
            ResizeWordSizedBag(list, SIZE_PLEN_GF2VEC(p));
            SET_LEN_GF2VEC(list, p);
        }
        if (EQ(GF2One, elm)) {
            BLOCK_ELM_GF2VEC(list, p) |= MASK_POS_GF2VEC(p);
        }
        else if (EQ(GF2Zero, elm)) {
            BLOCK_ELM_GF2VEC(list, p) &= ~MASK_POS_GF2VEC(p);
        }
        else if (IS_FFE(elm) && CHAR_FF(FLD_FFE(elm)) == 2 &&
                 DEGR_FF(FLD_FFE(elm)) <= 8) {
            RewriteGF2Vec(list, SIZE_FF(FLD_FFE(elm)));
            ASS_VEC8BIT(list, pos, elm);
        }
        else {
            PlainGF2Vec(list);
            ASS_LIST(list, p, elm);
        }
    }
    else {
        PlainGF2Vec(list);
        ASS_LIST(list, p, elm);
    }
    return 0;
}

/****************************************************************************
**
*F  FuncPLAIN_GF2MAT( <self>, <list> ) . . .  convert back into ordinary list
*/
static Obj FuncPLAIN_GF2MAT(Obj self, Obj list)
{
    PlainGF2Mat(list);
    return 0;
}

/****************************************************************************
**
*F  FuncASS_GF2MAT( <self>, <list>, <pos>, <elm> ) set an elm of a GF2 matrix
**
**  'ASS_GF2MAT' assigns the element  <elm> at the position  <pos> to the GF2
**  matrix <list>.
**
**  It is the responsibility of the caller  to ensure that <pos> is positive,
**  and that <elm> is not 0.
*/
static Obj FuncASS_GF2MAT(Obj self, Obj list, Obj pos, Obj elm)
{
    // check that <list> is mutable
    RequireMutable("List Assignment", list, "list");

    // get the position
    UInt p = GetSmallInt(SELF_NAME, pos);

    // if <elm> is a GF2 vector and the length is OK, keep the rep
    if (!IS_GF2VEC_REP(elm)) {
        PlainGF2Mat(list);
        ASS_LIST(list, p, elm);
    }
    else if (p == 1 && 1 >= LEN_GF2MAT(list)) {
        ResizeBag(list, SIZE_PLEN_GF2MAT(p));
        SetTypeDatObj(elm, IS_MUTABLE_OBJ(elm) ? TYPE_LIST_GF2VEC_LOCKED
                                               : TYPE_LIST_GF2VEC_IMM_LOCKED);
        SET_ELM_GF2MAT(list, p, elm);
        CHANGED_BAG(list);
    }
    else if (p > LEN_GF2MAT(list) + 1) {
        PlainGF2Mat(list);
        ASS_LIST(list, p, elm);
    }
    else if (LEN_GF2VEC(elm) == LEN_GF2VEC(ELM_GF2MAT(list, 1))) {
        if (LEN_GF2MAT(list) + 1 == p) {
            ResizeBag(list, SIZE_PLEN_GF2MAT(p));
            SET_LEN_GF2MAT(list, p);
        }
        SetTypeDatObj(elm, IS_MUTABLE_OBJ(elm) ? TYPE_LIST_GF2VEC_LOCKED
                                               : TYPE_LIST_GF2VEC_IMM_LOCKED);
        SET_ELM_GF2MAT(list, p, elm);
        CHANGED_BAG(list);
    }
    else {
        PlainGF2Mat(list);
        ASS_LIST(list, p, elm);
    }
    return 0;
}

/****************************************************************************
**
*F  FuncELM_GF2MAT( <self>, <mat>, <row> ) . . . select a row of a GF2 matrix
**
*/
static Obj FuncELM_GF2MAT(Obj self, Obj mat, Obj row)
{
    UInt r = GetSmallInt(SELF_NAME, row);
    if (LEN_GF2MAT(mat) < r) {
        ErrorMayQuit("row index %d exceeds %d, the number of rows", r,
                     LEN_GF2MAT(mat));
    }
    return ELM_GF2MAT(mat, r);
}

/****************************************************************************
**
*F  FuncSWAP_ROWS_GF2MAT( <self>, <mat>, <r1>, <r2> )
**
*/
static Obj FuncSWAP_ROWS_GF2MAT(Obj self, Obj mat, Obj row1, Obj row2)
{
    RequireMutable(SELF_NAME, mat, "mat");

    UInt r1 = GetSmallInt(SELF_NAME, row1);
    UInt r2 = GetSmallInt(SELF_NAME, row2);
    UInt m = LEN_GF2MAT(mat);
    if (m < r1) {
        ErrorMayQuit("row index %d exceeds %d, the number of rows", r1, m);
    }
    if (m < r2) {
        ErrorMayQuit("row index %d exceeds %d, the number of rows", r2, m);
    }

    Obj a = ELM_GF2MAT(mat, r1);
    Obj b = ELM_GF2MAT(mat, r2);

    SET_ELM_GF2MAT(mat, r1, b);
    SET_ELM_GF2MAT(mat, r2, a);

    return 0;
}

/****************************************************************************
**
*F  FuncSWAP_COLS_GF2MAT( <self>, <mat>, <c1>, <c2> )
**
*/
static Obj FuncSWAP_COLS_GF2MAT(Obj self, Obj mat, Obj col1, Obj col2)
{
    UInt c1 = GetSmallInt(SELF_NAME, col1);
    UInt c2 = GetSmallInt(SELF_NAME, col2);
    UInt m = LEN_GF2MAT(mat);
    if (m == 0)
        return 0;

    UInt n = LEN_GF2VEC(ELM_GF2MAT(mat, 1));
    if (n < c1) {
        ErrorMayQuit("column index %d exceeds %d, the number of columns", c1, n);
    }
    if (n < c2) {
        ErrorMayQuit("column index %d exceeds %d, the number of columns", c2, n);
    }

    for (UInt i = 1; i <= m; ++i) {
        Obj vec = ELM_GF2MAT(mat, i);
        if (!IS_MUTABLE_OBJ(vec)) {
            ErrorMayQuit("row %d is immutable", i, 0);
        }
        if (LEN_GF2VEC(vec) != n) {
            ErrorMayQuit("row length mismatch, %d versus %d", n, LEN_GF2VEC(vec));
        }

        Obj a = ELM_GF2VEC(vec, c1);
        Obj b = ELM_GF2VEC(vec, c2);
        if (a != b) {
            // one entry is 1, the other is 0; so to switch them,
            // we can simply flip each individually
            BLOCK_ELM_GF2VEC(vec, c1) ^= MASK_POS_GF2VEC(c1);
            BLOCK_ELM_GF2VEC(vec, c2) ^= MASK_POS_GF2VEC(c2);
        }
    }

    return 0;
}

/****************************************************************************
**
*F  FuncUNB_GF2VEC( <self>, <list>, <pos> ) . unbind position of a GF2 vector
**
**  'UNB_GF2VEC' unbind  the element at  the position  <pos> in  a GF2 vector
**  <list>.
**
**  It is the responsibility of the caller  to ensure that <pos> is positive.
*/
static Obj FuncUNB_GF2VEC(Obj self, Obj list, Obj pos)
{
    // check that <list> is mutable
    RequireMutable("List Unbind", list, "vector");

    if (DoFilter(IsLockedRepresentationVector, list) == True) {
        ErrorMayQuit("Unbind forbidden on locked GF2 vector", 0, 0);
    }

    // get the position
    UInt p = GetSmallInt(SELF_NAME, pos);

    // if we unbind the last position keep the representation
    if (LEN_GF2VEC(list) < p) {
        ;
    }
    else if (LEN_GF2VEC(list) == p) {
        ResizeWordSizedBag(list, SIZE_PLEN_GF2VEC(p - 1));
        SET_LEN_GF2VEC(list, p - 1);
    }
    else {
        PlainGF2Vec(list);
        UNB_LIST(list, p);
    }
    return 0;
}

/****************************************************************************
**
*F  FuncUNB_GF2MAT( <self>, <list>, <pos> ) . unbind position of a GF2 matrix
**
**  'UNB_GF2VEC' unbind  the element at  the position  <pos> in  a GF2 matrix
**  <list>.
**
**  It is the responsibility of the caller  to ensure that <pos> is positive.
*/
static Obj FuncUNB_GF2MAT(Obj self, Obj list, Obj pos)
{
    // check that <list> is mutable
    RequireMutable("List Unbind", list, "matrix");

    // get the position
    UInt p = GetSmallInt(SELF_NAME, pos);

    // if we unbind the last position keep the representation
    if (p > 1 && LEN_GF2MAT(list) < p) {
        ;
    }
    else if (LEN_GF2MAT(list) == p) {
        ResizeBag(list, SIZE_PLEN_GF2MAT(p - 1));
        SET_LEN_GF2MAT(list, p - 1);
    }
    else {
        PlainGF2Mat(list);
        UNB_LIST(list, p);
    }
    return 0;
}

/****************************************************************************
**
*F * * * * * * * * * * * * arithmetic operations  * * * * * * * * * * * * * *
*/

/****************************************************************************
**
*F  FuncZERO_GF2VEC( <self>, <mat> )  . . . . . . . . . . . . zero GF2 vector
**
**  return the zero vector over GF2 of the same length as <mat>.
*/
static Obj FuncZERO_GF2VEC(Obj self, Obj mat)
{
    Obj  zero;
    UInt len;

    // create a new GF2 vector
    len = LEN_GF2VEC(mat);
    NEW_GF2VEC(zero, TYPE_LIST_GF2VEC, len);
    return zero;
}

/****************************************************************************
**
*F  FuncZERO_GF2VEC_2( <self>, <len>) . . . . . . . . . zero GF2 vector
**
**  return the zero vector over GF2 of length <len>
*/
static Obj FuncZERO_GF2VEC_2(Obj self, Obj len)
{
    Obj zero;
    RequireNonnegativeSmallInt(SELF_NAME, len);
    NEW_GF2VEC(zero, TYPE_LIST_GF2VEC, INT_INTOBJ(len));
    return zero;
}

--> --------------------

--> maximum size reached

--> --------------------

quality93%

¤ Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.0.87Bemerkung: (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.