packedmatrix.h File Reference

Dense matrices over GF(2) represented as a bit field. More...

#include <math.h>
#include <assert.h>
#include <stdio.h>

Go to the source code of this file.

Data Structures
struct	mzd_t
	Dense matrices over GF(2). More...
Defines
#define	MZD_MUL_BLOCKSIZE   MIN(((int)sqrt((double)(4*CPU_L2_CACHE)))/2,2048)
	Matrix multiplication block-ing dimension.
#define	mzd_free_window   mzd_free
	Free a matrix window created with mzd_init_window.
#define	mzd_sub   mzd_add
	Same as mzd_add.
#define	_mzd_sub   _mzd_add
	Same as mzd_sub but without any checks on the input.
Functions
mzd_t *	mzd_init (const size_t r, const size_t c)
	Create a new matrix of dimension r x c.
void	mzd_free (mzd_t *A)
	Free a matrix created with mzd_init.
mzd_t *	mzd_init_window (const mzd_t *M, const size_t lowr, const size_t lowc, const size_t highr, const size_t highc)
	Create a window/view into the matrix M.
static void	mzd_row_swap (mzd_t *M, const size_t rowa, const size_t rowb)
	Swap the two rows rowa and rowb.
void	mzd_copy_row (mzd_t *B, size_t i, const mzd_t *A, size_t j)
	copy row j from A to row i from B.
void	mzd_col_swap (mzd_t *M, const size_t cola, const size_t colb)
	Swap the two columns cola and colb.
static void	mzd_col_swap_in_rows (mzd_t *M, const size_t cola, const size_t colb, const size_t start_row, const size_t stop_row)
	Swap the two columns cola and colb but only between start_row and stop_row.
static BIT	mzd_read_bit (const mzd_t *M, const size_t row, const size_t col)
	Read the bit at position M[row,col].
static void	mzd_write_bit (mzd_t *M, const size_t row, const size_t col, const BIT value)
	Write the bit value to position M[row,col].
void	mzd_print (const mzd_t *M)
	Print a matrix to stdout.
void	mzd_print_tight (const mzd_t *M)
	Print the matrix to stdout.
static void	mzd_row_add_offset (mzd_t *M, size_t dstrow, size_t srcrow, size_t coloffset)
	Add the rows sourcerow and destrow and stores the total in the row destrow, but only begins at the column coloffset.
void	mzd_row_add (mzd_t *M, const size_t sourcerow, const size_t destrow)
	Add the rows sourcerow and destrow and stores the total in the row destrow.
mzd_t *	mzd_transpose (mzd_t *DST, const mzd_t *A)
	Transpose a matrix.
mzd_t *	mzd_mul_naive (mzd_t *C, const mzd_t *A, const mzd_t *B)
	Naive cubic matrix multiplication.
mzd_t *	mzd_addmul_naive (mzd_t *C, const mzd_t *A, const mzd_t *B)
	Naive cubic matrix multiplication and addition.
mzd_t *	_mzd_mul_naive (mzd_t *C, const mzd_t *A, const mzd_t *B, const int clear)
	Naive cubic matrix multiplication with the pre-transposed B.
mzd_t *	_mzd_mul_va (mzd_t *C, const mzd_t *v, const mzd_t *A, const int clear)
	Matrix multiplication optimized for v*A where v is a vector.
void	mzd_randomize (mzd_t *M)
	Fill matrix M with uniformly distributed bits.
void	mzd_set_ui (mzd_t *M, const unsigned value)
	Set the matrix M to the value equivalent to the integer value provided.
int	mzd_gauss_delayed (mzd_t *M, const size_t startcol, const int full)
	Gaussian elimination.
int	mzd_echelonize_naive (mzd_t *M, const int full)
	Gaussian elimination.
BIT	mzd_equal (const mzd_t *A, const mzd_t *B)
	Return TRUE if A == B.
int	mzd_cmp (const mzd_t *A, const mzd_t *B)
	Return -1,0,1 if if A < B, A == B or A > B respectively.
mzd_t *	mzd_copy (mzd_t *DST, const mzd_t *A)
	Copy matrix A to DST.
mzd_t *	mzd_concat (mzd_t *C, const mzd_t *A, const mzd_t *B)
	Concatenate B to A and write the result to C.
mzd_t *	mzd_stack (mzd_t *C, const mzd_t *A, const mzd_t *B)
	Stack A on top of B and write the result to C.
mzd_t *	mzd_submatrix (mzd_t *S, const mzd_t *M, const size_t lowr, const size_t lowc, const size_t highr, const size_t highc)
	Copy a submatrix.
mzd_t *	mzd_invert_naive (mzd_t *INV, mzd_t *A, const mzd_t *I)
	Invert the matrix target using Gaussian elimination.
mzd_t *	mzd_add (mzd_t *C, const mzd_t *A, const mzd_t *B)
	Set C = A+B.
mzd_t *	_mzd_add (mzd_t *C, const mzd_t *A, const mzd_t *B)
	Same as mzd_add but without any checks on the input.
void	mzd_combine (mzd_t *DST, const size_t row3, const size_t startblock3, const mzd_t *SC1, const size_t row1, const size_t startblock1, const mzd_t *SC2, const size_t row2, const size_t startblock2)
	row3[col3:] = row1[col1:] + row2[col2:]
static word	mzd_read_bits (const mzd_t *M, const size_t x, const size_t y, const int n)
static void	mzd_xor_bits (const mzd_t *M, const size_t x, const size_t y, const int n, word values)
	XOR n bits from values to M starting a position (x,y).
static void	mzd_and_bits (const mzd_t *M, const size_t x, const size_t y, const int n, word values)
	AND n bits from values to M starting a position (x,y).
static void	mzd_clear_bits (const mzd_t *M, const size_t x, const size_t y, const int n)
	Clear n bits in M starting a position (x,y).
int	mzd_is_zero (mzd_t *A)
	Zero test for matrix.
void	mzd_row_clear_offset (mzd_t *M, const size_t row, const size_t coloffset)
	Clear the given row, but only begins at the column coloffset.
int	mzd_find_pivot (mzd_t *M, size_t start_row, size_t start_col, size_t *r, size_t *c)
	Find the next nonzero entry in M starting at start_row and start_col.
double	mzd_density (mzd_t *A, int res)
	Return the number of nonzero entries divided by nrows * ncols.
double	_mzd_density (mzd_t *A, int res, size_t r, size_t c)
	Return the number of nonzero entries divided by nrows * ncols considering only the submatrix starting at (r,c).
size_t	mzd_first_zero_row (mzd_t *A)
	Return the first row with all zero entries.

---

Detailed Description

Dense matrices over GF(2) represented as a bit field.

Author:

Gregory Bard <bard@fordham.edu>

Martin Albrecht <M.R.Albrecht@rhul.ac.uk>

---

Define Documentation

#define _mzd_sub   _mzd_add

Same as mzd_sub but without any checks on the input.

Parameters:

C	Preallocated difference matrix, may be NULL for automatic creation.
A	Matrix
B	Matrix

Ignores offset atrtribute of packedmatrix.

#define mzd_free_window   mzd_free

Free a matrix window created with mzd_init_window.

Parameters:

A

Matrix

#define MZD_MUL_BLOCKSIZE   MIN(((int)sqrt((double)(4*CPU_L2_CACHE)))/2,2048)

Matrix multiplication block-ing dimension.

Defines the number of rows of the matrix A that are processed as one block during the execution of a multiplication algorithm.

#define mzd_sub   mzd_add

Same as mzd_add.

Parameters:

C	Preallocated difference matrix, may be NULL for automatic creation.
A	Matrix
B	Matrix

Ignores offset atrtribute of packedmatrix.

---

Function Documentation

mzd_t* _mzd_add	(	mzd_t *	C,
		const mzd_t *	A,
		const mzd_t *	B
	)

Same as mzd_add but without any checks on the input.

Parameters:

C	Preallocated sum matrix, may be NULL for automatic creation.
A	Matrix
B	Matrix

Ignores offset atrtribute of packedmatrix.

double _mzd_density	(	mzd_t *	A,
		int	res,
		size_t	r,
		size_t	c
	)

Return the number of nonzero entries divided by nrows * ncols considering only the submatrix starting at (r,c).

If res = 0 then 100 samples per row are made, if res > 0 the function takes res sized steps within each row (res = 1 uses every word).

Parameters:

A	Matrix
res	Resolution of sampling
r	Row to start counting
c	Column to start counting

mzd_t* _mzd_mul_naive	(	mzd_t *	C,
		const mzd_t *	A,
		const mzd_t *	B,
		const int	clear
	)

Naive cubic matrix multiplication with the pre-transposed B.

That is, compute C such that C == AB^t.

Parameters:

C	Preallocated product matrix.
A	Input matrix A.
B	Pre-transposed input matrix B.
clear	Whether to clear C before accumulating AB

mzd_t* _mzd_mul_va	(	mzd_t *	C,
		const mzd_t *	v,
		const mzd_t *	A,
		const int	clear
	)

Matrix multiplication optimized for v*A where v is a vector.

Parameters:

C	Preallocated product matrix.
v	Input matrix v.
A	Input matrix A.
clear	If set clear C first, otherwise add result to C.

mzd_t* mzd_add	(	mzd_t *	C,
		const mzd_t *	A,
		const mzd_t *	B
	)

Set C = A+B.

C is also returned. If C is NULL then a new matrix is created which must be freed by mzd_free.

Parameters:

C	Preallocated sum matrix, may be NULL for automatic creation.
A	Matrix
B	Matrix

Examples:

testsuite/test_multiplication.c.

mzd_t* mzd_addmul_naive	(	mzd_t *	C,
		const mzd_t *	A,
		const mzd_t *	B
	)

Naive cubic matrix multiplication and addition.

That is, compute C such that C == C + AB.

Parameters:

C	Preallocated product matrix.
A	Input matrix A.
B	Input matrix B.

Note:

Normally, if you will multiply several times by b, it is smarter to calculate bT yourself, and keep it, and then use the function called _mzd_mul_naive

static void mzd_and_bits	(	const mzd_t *	M,
		const size_t	x,
		const size_t	y,
		const int	n,
		word	values
	)		[inline, static]

AND n bits from values to M starting a position (x,y).

Parameters:

M	Source matrix.
x	Starting row.
y	Starting column.
n	Number of bits (<= RADIX);
values	Word with values;

static void mzd_clear_bits	(	const mzd_t *	M,
		const size_t	x,
		const size_t	y,
		const int	n
	)		[inline, static]

Clear n bits in M starting a position (x,y).

Parameters:

M	Source matrix.
x	Starting row.
y	Starting column.
n	Number of bits (<= RADIX);

int mzd_cmp	(	const mzd_t *	A,
		const mzd_t *	B
	)

Return -1,0,1 if if A < B, A == B or A > B respectively.

Parameters:

A	Matrix.
B	Matrix.

Note:

This comparison is not well defined mathematically and relatively arbitrary since elements of GF(2) don't have an ordering.

Ignores offset atrtribute of packedmatrix.

void mzd_col_swap	(	mzd_t *	M,
		const size_t	cola,
		const size_t	colb
	)

Swap the two columns cola and colb.

Parameters:

M	Matrix.
cola	Column index.
colb	Column index.

static void mzd_col_swap_in_rows	(	mzd_t *	M,
		const size_t	cola,
		const size_t	colb,
		const size_t	start_row,
		const size_t	stop_row
	)		[inline, static]

Swap the two columns cola and colb but only between start_row and stop_row.

Parameters:

M	Matrix.
cola	Column index.
colb	Column index.
start_row	Row index.
stop_row	Row index (exclusive).

void mzd_combine	(	mzd_t *	DST,
		const size_t	row3,
		const size_t	startblock3,
		const mzd_t *	SC1,
		const size_t	row1,
		const size_t	startblock1,
		const mzd_t *	SC2,
		const size_t	row2,
		const size_t	startblock2
	)

row3[col3:] = row1[col1:] + row2[col2:]

Adds row1 of SC1, starting with startblock1 to the end, to row2 of SC2, starting with startblock2 to the end. This gets stored in DST, in row3, starting with startblock3.

Parameters:

DST	destination matrix
row3	destination row for matrix dst
startblock3	starting block to work on in matrix dst
SC1	source matrix
row1	source row for matrix sc1
startblock1	starting block to work on in matrix sc1
SC2	source matrix
startblock2	starting block to work on in matrix sc2
row2	source row for matrix sc2

Ignores offset atrtribute of packedmatrix.

Todo:

this code is slow if offset!=0

mzd_t* mzd_concat	(	mzd_t *	C,
		const mzd_t *	A,
		const mzd_t *	B
	)

Concatenate B to A and write the result to C.

That is,

 [ A ], [ B ] -> [ A  B ] = C
 

The inputs are not modified but a new matrix is created.

Parameters:

C	Matrix, may be NULL for automatic creation
A	Matrix
B	Matrix

Note:

This is sometimes called augment.

Ignores offset atrtribute of packedmatrix.

mzd_t* mzd_copy	(	mzd_t *	DST,
		const mzd_t *	A
	)

Copy matrix A to DST.

Parameters:

DST	May be NULL for automatic creation.
A	Source matrix.

Ignores offset atrtribute of packedmatrix.

Examples:

testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.

void mzd_copy_row	(	mzd_t *	B,
		size_t	i,
		const mzd_t *	A,
		size_t	j
	)

copy row j from A to row i from B.

The offsets of A and B must match and the number of columns of A must be less than or equal to the number of columns of B.

Parameters:

B	Target matrix.
i	Target row index.
A	Source matrix.
j	Source row index.

double mzd_density	(	mzd_t *	A,
		int	res
	)

Return the number of nonzero entries divided by nrows * ncols.

If res = 0 then 100 samples per row are made, if res > 0 the function takes res sized steps within each row (res = 1 uses every word).

Parameters:

A	Matrix
res	Resolution of sampling

int mzd_echelonize_naive	(	mzd_t *	M,
		const int	full
	)

Gaussian elimination.

This will do Gaussian elimination on the matrix m. If full=FALSE, then it will do triangular style elimination, and if full=TRUE, it will do Gauss-Jordan style, or full elimination.

Parameters:

M	Matrix
full	Gauss-Jordan style or upper triangular form only.

Ignores offset atrtribute of packedmatrix.

See also:

mzd_echelonize_m4ri(), mzd_echelonize_pluq()

Examples:

testsuite/bench_elimination.c, and testsuite/test_elimination.c.

BIT mzd_equal	(	const mzd_t *	A,
		const mzd_t *	B
	)

Return TRUE if A == B.

Parameters:

A	Matrix
B	Matrix

Ignores offset atrtribute of packedmatrix.

Examples:

testsuite/test_elimination.c, and testsuite/test_multiplication.c.

int mzd_find_pivot	(	mzd_t *	M,
		size_t	start_row,
		size_t	start_col,
		size_t *	r,
		size_t *	c
	)

Find the next nonzero entry in M starting at start_row and start_col.

This function walks down rows in the inner loop and columns in the outer loop. If a nonzero entry is found this function returns 1 and zero otherwise.

If and only if a nonzero entry is found r and c are updated.

Parameters:

M	Matrix
start_row	Index of row where to start search
start_col	Index of column where to start search
r	Row index updated if pivot is found
c	Column index updated if pivot is found

size_t mzd_first_zero_row

(

mzd_t *

A

)

Return the first row with all zero entries.

If no such row can be found returns nrows.

Parameters:

A

Matrix

void mzd_free

(

mzd_t *

A

)

Free a matrix created with mzd_init.

Parameters:

A

Matrix

Examples:

testsuite/bench_elimination.c, testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.

int mzd_gauss_delayed	(	mzd_t *	M,
		const size_t	startcol,
		const int	full
	)

Gaussian elimination.

This will do Gaussian elimination on the matrix m but will start not at column 0 necc but at column startcol. If full=FALSE, then it will do triangular style elimination, and if full=TRUE, it will do Gauss-Jordan style, or full elimination.

Parameters:

M	Matrix
startcol	First column to consider for reduction.
full	Gauss-Jordan style or upper triangular form only.

Ignores offset atrtribute of packedmatrix.

mzd_t* mzd_init	(	const size_t	r,
		const size_t	c
	)

Create a new matrix of dimension r x c.

Use mzd_free to kill it.

Parameters:

r	Number of rows
c	Number of columns

Examples:

testsuite/bench_elimination.c, testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.

mzd_t* mzd_init_window	(	const mzd_t *	M,
		const size_t	lowr,
		const size_t	lowc,
		const size_t	highr,
		const size_t	highc
	)

Create a window/view into the matrix M.

A matrix window for M is a meta structure on the matrix M. It is setup to point into the matrix so M must not be freed while the matrix window is used.

This function puts the restriction on the provided parameters that all parameters must be within range for M which is not enforced currently .

Use mzd_free_window to free the window.

Parameters:

M	Matrix
lowr	Starting row (inclusive)
lowc	Starting column (inclusive)
highr	End row (exclusive)
highc	End column (exclusive)

mzd_t* mzd_invert_naive	(	mzd_t *	INV,
		mzd_t *	A,
		const mzd_t *	I
	)

Invert the matrix target using Gaussian elimination.

To avoid recomputing the identity matrix over and over again, I may be passed in as identity parameter.

Parameters:

INV	Preallocated space for inversion matrix, may be NULL for automatic creation.
A	Matrix to be reduced.
I	Identity matrix.

Ignores offset atrtribute of packedmatrix.

int mzd_is_zero

(

mzd_t *

A

)

Zero test for matrix.

Parameters:

A	Input matrix.

mzd_t* mzd_mul_naive	(	mzd_t *	C,
		const mzd_t *	A,
		const mzd_t *	B
	)

Naive cubic matrix multiplication.

That is, compute C such that C == AB.

Parameters:

C	Preallocated product matrix, may be NULL for automatic creation.
A	Input matrix A.
B	Input matrix B.

Note:

Normally, if you will multiply several times by b, it is smarter to calculate bT yourself, and keep it, and then use the function called _mzd_mul_naive

Examples:

testsuite/test_multiplication.c.

void mzd_print

(

const mzd_t *

M

)

Print a matrix to stdout.

The output will contain colons between every 4-th column.

Parameters:

M

Matrix

void mzd_print_tight

(

const mzd_t *

M

)

Print the matrix to stdout.

Parameters:

M

Matrix

void mzd_randomize

(

mzd_t *

M

)

Fill matrix M with uniformly distributed bits.

Parameters:

M

Matrix

Todo:

Allow the user to provide a RNG callback.

Ignores offset atrtribute of packedmatrix.

Examples:

testsuite/bench_elimination.c, testsuite/test_elimination.c, testsuite/test_multiplication.c, and testsuite/test_pluq.c.

static BIT mzd_read_bit	(	const mzd_t *	M,
		const size_t	row,
		const size_t	col
	)		[inline, static]

Read the bit at position M[row,col].

Parameters:

M	Matrix
row	Row index
col	Column index

Note:

No bounds checks whatsoever are performed.

Examples:

testsuite/test_pluq.c.

static word mzd_read_bits	(	const mzd_t *	M,
		const size_t	x,
		const size_t	y,
		const int	n
	)		[inline, static]

Get n bits starting a position (x,y) from the matrix M.

Parameters:

M	Source matrix.
x	Starting row.
y	Starting column.
n	Number of bits (<= RADIX);

void mzd_row_add	(	mzd_t *	M,
		const size_t	sourcerow,
		const size_t	destrow
	)

Add the rows sourcerow and destrow and stores the total in the row destrow.

Parameters:

M	Matrix
sourcerow	Index of source row
destrow	Index of target row

Note:

this can be done much faster with mzd_combine.

static void mzd_row_add_offset	(	mzd_t *	M,
		size_t	dstrow,
		size_t	srcrow,
		size_t	coloffset
	)		[inline, static]

Add the rows sourcerow and destrow and stores the total in the row destrow, but only begins at the column coloffset.

Parameters:

M	Matrix
dstrow	Index of target row
srcrow	Index of source row
coloffset	Column offset

void mzd_row_clear_offset	(	mzd_t *	M,
		const size_t	row,
		const size_t	coloffset
	)

Clear the given row, but only begins at the column coloffset.

Parameters:

M	Matrix
row	Index of row
coloffset	Column offset

static void mzd_row_swap	(	mzd_t *	M,
		const size_t	rowa,
		const size_t	rowb
	)		[inline, static]

Swap the two rows rowa and rowb.

Parameters:

M	Matrix
rowa	Row index.
rowb	Row index.

void mzd_set_ui	(	mzd_t *	M,
		const unsigned	value
	)

Set the matrix M to the value equivalent to the integer value provided.

Specifically, this function does nothing if value%2 == 0 and returns the identity matrix if value%2 == 1.

If the matrix is not square then the largest possible square submatrix is set to the identity matrix.

Parameters:

M	Matrix
value	Either 0 or 1

mzd_t* mzd_stack	(	mzd_t *	C,
		const mzd_t *	A,
		const mzd_t *	B
	)

Stack A on top of B and write the result to C.

That is,

 [ A ], [ B ] -> [ A ] = C
                 [ B ]
 

The inputs are not modified but a new matrix is created.

Parameters:

C	Matrix, may be NULL for automatic creation
A	Matrix
B	Matrix

Ignores offset atrtribute of packedmatrix.

mzd_t* mzd_submatrix	(	mzd_t *	S,
		const mzd_t *	M,
		const size_t	lowr,
		const size_t	lowc,
		const size_t	highr,
		const size_t	highc
	)

Copy a submatrix.

Note that the upper bounds are not included.

Parameters:

S	Preallocated space for submatrix, may be NULL for automatic creation.
M	Matrix
lowr	start rows
lowc	start column
highr	stop row (this row is not included)
highc	stop column (this column is not included)

mzd_t* mzd_transpose	(	mzd_t *	DST,
		const mzd_t *	A
	)

Transpose a matrix.

This function uses the fact that:

   [ A B ]T    [AT CT]
   [ C D ]  =  [BT DT] 
 

and thus rearranges the blocks recursively.

Parameters:

DST	Preallocated return matrix, may be NULL for automatic creation.
A	Matrix

static void mzd_write_bit	(	mzd_t *	M,
		const size_t	row,
		const size_t	col,
		const BIT	value
	)		[inline, static]

Write the bit value to position M[row,col].

Parameters:

M	Matrix
row	Row index
col	Column index
value	Either 0 or 1

Note:

No bounds checks whatsoever are performed.

Examples:

testsuite/bench_elimination.c, and testsuite/test_pluq.c.

static void mzd_xor_bits	(	const mzd_t *	M,
		const size_t	x,
		const size_t	y,
		const int	n,
		word	values
	)		[inline, static]

XOR n bits from values to M starting a position (x,y).

Parameters:

M	Source matrix.
x	Starting row.
y	Starting column.
n	Number of bits (<= RADIX);
values	Word with values;