755 lines
		
	
	
		
			16 KiB
		
	
	
	
		
			C
		
	
	
		
			Executable File
		
	
	
	
	
			
		
		
	
	
			755 lines
		
	
	
		
			16 KiB
		
	
	
	
		
			C
		
	
	
		
			Executable File
		
	
	
	
	
| /* MOD2DENSE.C - Procedures for handling dense mod2 matrices. */
 | |
| 
 | |
| /* Copyright (c) 1995-2012 by Radford M. Neal.
 | |
|  *
 | |
|  * Permission is granted for anyone to copy, use, modify, and distribute
 | |
|  * these programs and accompanying documents for any purpose, provided
 | |
|  * this copyright notice is retained and prominently displayed, and note
 | |
|  * is made of any changes made to these programs.  These programs and
 | |
|  * documents are distributed without any warranty, express or implied.
 | |
|  * As the programs were written for research purposes only, they have not
 | |
|  * been tested to the degree that would be advisable in any important
 | |
|  * application.  All use of these programs is entirely at the user's own
 | |
|  * risk.
 | |
|  */
 | |
| 
 | |
| 
 | |
| /* NOTE:  See mod2dense.html for documentation on these procedures. */
 | |
| 
 | |
| 
 | |
| #include <stdlib.h>
 | |
| #include <stdio.h>
 | |
| #include <math.h>
 | |
| 
 | |
| #include "alloc.h"
 | |
| #include "intio.h"
 | |
| #include "mod2dense.h"
 | |
| 
 | |
| 
 | |
| /* ALLOCATE SPACE FOR A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| mod2dense *mod2dense_allocate 
 | |
| ( int n_rows, 		/* Number of rows in matrix */
 | |
|   int n_cols		/* Number of columns in matrix */
 | |
| )
 | |
| {
 | |
|   mod2dense *m;
 | |
|   int j;
 | |
| 
 | |
|   if (n_rows<=0 || n_cols<=0)
 | |
|   { fprintf(stderr,"mod2dense_allocate: Invalid number of rows or columns\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   m = chk_alloc (1, sizeof *m);
 | |
| 
 | |
|   m->n_rows = n_rows;
 | |
|   m->n_cols = n_cols;
 | |
| 
 | |
|   m->n_words = (n_rows+mod2_wordsize-1) >> mod2_wordsize_shift;
 | |
| 
 | |
|   m->col = chk_alloc (m->n_cols, sizeof *m->col);
 | |
| 
 | |
|   m->bits = chk_alloc(m->n_words*m->n_cols, sizeof *m->bits);
 | |
| 
 | |
|   for (j = 0; j<m->n_cols; j++)
 | |
|   { m->col[j] = m->bits + j*m->n_words;
 | |
|   }
 | |
| 
 | |
|   return m;
 | |
| }
 | |
| 
 | |
| 
 | |
| /* FREE SPACE OCCUPIED BY A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| void mod2dense_free
 | |
| ( mod2dense *m		/* Matrix to free */
 | |
| )
 | |
| { if (m)
 | |
|   { free(m->bits);
 | |
|     free(m->col);
 | |
|     free(m);
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* CLEAR A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| void mod2dense_clear
 | |
| ( mod2dense *r
 | |
| )
 | |
| {
 | |
|   int k, j;
 | |
| 
 | |
|   for (j = 0; j<mod2dense_cols(r); j++)
 | |
|   { for (k = 0; k<r->n_words; k++)
 | |
|     { r->col[j][k] = 0;
 | |
|     }
 | |
|   }
 | |
| 
 | |
| }
 | |
| 
 | |
| 
 | |
| /* COPY A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| void mod2dense_copy
 | |
| ( mod2dense *m,		/* Matrix to copy */
 | |
|   mod2dense *r		/* Place to store copy of matrix */
 | |
| )
 | |
| { 
 | |
|   int k, j;
 | |
| 
 | |
|   if (mod2dense_rows(m)>mod2dense_rows(r) 
 | |
|    || mod2dense_cols(m)>mod2dense_cols(r))
 | |
|   { fprintf(stderr,"mod2dense_copy: Destination matrix is too small\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   for (j = 0; j<mod2dense_cols(m); j++)
 | |
|   { for (k = 0; k<m->n_words; k++)
 | |
|     { r->col[j][k] = m->col[j][k];
 | |
|     }
 | |
|     for ( ; k<r->n_words; k++)
 | |
|     { r->col[j][k] = 0;
 | |
|     }
 | |
|   }
 | |
| 
 | |
|   for ( ; j<mod2dense_cols(r); j++)
 | |
|   { for (k = 0; k<r->n_words; k++)
 | |
|     { r->col[j][k] = 0;
 | |
|     }
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* COPY ROWS OF A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| void mod2dense_copyrows
 | |
| ( mod2dense *m,		/* Matrix to copy */
 | |
|   mod2dense *r,		/* Place to store copy of matrix */
 | |
|   int *rows		/* Indexes of rows to copy, from 0 */
 | |
| )
 | |
| { 
 | |
|   int i, j;
 | |
| 
 | |
|   if (mod2dense_cols(m)>mod2dense_cols(r))
 | |
|   { fprintf(stderr,
 | |
|       "mod2dense_copyrows: Destination matrix has fewer columns than source\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   mod2dense_clear(r);
 | |
| 
 | |
|   for (i = 0; i<mod2dense_rows(r); i++)
 | |
|   { if (rows[i]<0 || rows[i]>=mod2dense_rows(m))
 | |
|     { fprintf(stderr,"mod2dense_copyrows: Row index out of range\n");
 | |
|       exit(1);
 | |
|     }
 | |
|     for (j = 0; j<mod2dense_cols(m); j++)
 | |
|     { mod2dense_set(r,i,j,mod2dense_get(m,rows[i],j));
 | |
|     }
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* COPY COLUMNS OF A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| void mod2dense_copycols
 | |
| ( mod2dense *m,		/* Matrix to copy */
 | |
|   mod2dense *r,		/* Place to store copy of matrix */
 | |
|   int *cols		/* Indexes of columns to copy, from 0 */
 | |
| )
 | |
| { 
 | |
|   int k, j;
 | |
| 
 | |
|   if (mod2dense_rows(m)>mod2dense_rows(r))
 | |
|   { fprintf(stderr,
 | |
|       "mod2dense_copycols: Destination matrix has fewer rows than source\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   for (j = 0; j<mod2dense_cols(r); j++)
 | |
|   { if (cols[j]<0 || cols[j]>=mod2dense_cols(m))
 | |
|     { fprintf(stderr,"mod2dense_copycols: Column index out of range\n");
 | |
|       exit(1);
 | |
|     }
 | |
|     for (k = 0; k<m->n_words; k++)
 | |
|     { r->col[j][k] = m->col[cols[j]][k];
 | |
|     }
 | |
|     for ( ; k<r->n_words; k++)
 | |
|     { r->col[j][k] = 0;
 | |
|     }
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* PRINT A DENSE MOD2 MATRIX IN HUMAN-READABLE FORM. */
 | |
| 
 | |
| void mod2dense_print     
 | |
| ( FILE *f,
 | |
|   mod2dense *m
 | |
| )
 | |
| { 
 | |
|   int i, j;
 | |
| 
 | |
|   for (i = 0; i<mod2dense_rows(m); i++)
 | |
|   { for (j = 0; j<mod2dense_cols(m); j++)
 | |
|     { fprintf(f," %d",mod2dense_get(m,i,j));
 | |
|     }
 | |
|     fprintf(f,"\n");
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* WRITE A DENSE MOD2 MATRIX TO A FILE IN MACHINE-READABLE FORM.
 | |
| 
 | |
|    Data is written using intio_write, so that it will be readable on a machine
 | |
|    with a different byte-ordering.  At present, this assumes that the words 
 | |
|    used to pack bits into are no longer than 32 bits. */
 | |
| 
 | |
| int mod2dense_write     
 | |
| ( FILE *f, 
 | |
|   mod2dense *m
 | |
| )
 | |
| { 
 | |
|   int j, k;
 | |
| 
 | |
|   intio_write(f,m->n_rows);
 | |
|   if (ferror(f)) return 0;
 | |
| 
 | |
|   intio_write(f,m->n_cols);
 | |
|   if (ferror(f)) return 0;
 | |
| 
 | |
|   for (j = 0; j<mod2dense_cols(m); j++)
 | |
|   {
 | |
|     for (k = 0; k<m->n_words; k++)
 | |
|     { intio_write(f,m->col[j][k]);
 | |
|       if (ferror(f)) return 0;
 | |
|     }
 | |
|   }
 | |
| 
 | |
|   return 1;
 | |
| }
 | |
| 
 | |
| 
 | |
| /* READ A DENSE MOD2 MATRIX STORED IN MACHINE-READABLE FORM FROM A FILE. */
 | |
| 
 | |
| mod2dense *mod2dense_read  
 | |
| ( FILE *f
 | |
| )
 | |
| { 
 | |
|   int n_rows, n_cols;
 | |
|   mod2dense *m;
 | |
|   int j, k;
 | |
|   
 | |
|   n_rows = intio_read(f);
 | |
|   if (feof(f) || ferror(f) || n_rows<=0) return 0;
 | |
| 
 | |
|   n_cols = intio_read(f);
 | |
|   if (feof(f) || ferror(f) || n_cols<=0) return 0;
 | |
| 
 | |
|   m = mod2dense_allocate(n_rows,n_cols);
 | |
| 
 | |
|   for (j = 0; j<mod2dense_cols(m); j++)
 | |
|   {
 | |
|     for (k = 0; k<m->n_words; k++)
 | |
|     { m->col[j][k] = intio_read(f);
 | |
|       if (feof(f) || ferror(f)) 
 | |
|       { mod2dense_free(m);
 | |
|         return 0;
 | |
|       }
 | |
|     }
 | |
|   }
 | |
| 
 | |
|   return m;
 | |
| }
 | |
| 
 | |
| 
 | |
| /* GET AN ELEMENT FROM A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| int mod2dense_get  
 | |
| ( mod2dense *m, 	/* Matrix to get element from */
 | |
|   int row,		/* Row of element (starting with zero) */
 | |
|   int col		/* Column of element (starting with zero) */
 | |
| )
 | |
| {
 | |
|   if (row<0 || row>=mod2dense_rows(m) || col<0 || col>=mod2dense_cols(m))
 | |
|   { fprintf(stderr,"mod2dense_get: row or column index out of bounds\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   return mod2_getbit (m->col[col][row>>mod2_wordsize_shift], 
 | |
|                       row&mod2_wordsize_mask);
 | |
| }
 | |
| 
 | |
| 
 | |
| /* SET AN ELEMENT IN A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| void mod2dense_set 
 | |
| ( mod2dense *m, 	/* Matrix to modify element of */
 | |
|   int row,		/* Row of element (starting with zero) */
 | |
|   int col,		/* Column of element (starting with zero) */
 | |
|   int value		/* New value of element (0 or 1) */
 | |
| )
 | |
| { 
 | |
|   mod2word *w;
 | |
| 
 | |
|   if (row<0 || row>=mod2dense_rows(m) || col<0 || col>=mod2dense_cols(m))
 | |
|   { fprintf(stderr,"mod2dense_set: row or column index out of bounds\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   w = &m->col[col][row>>mod2_wordsize_shift];
 | |
| 
 | |
|   *w = value ? mod2_setbit1(*w,row&mod2_wordsize_mask) 
 | |
|              : mod2_setbit0(*w,row&mod2_wordsize_mask);
 | |
| }
 | |
| 
 | |
| 
 | |
| /* FLIP AN ELEMENT OF A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| int mod2dense_flip  
 | |
| ( mod2dense *m, 	/* Matrix to flip element in */
 | |
|   int row,		/* Row of element (starting with zero) */
 | |
|   int col		/* Column of element (starting with zero) */
 | |
| )
 | |
| {
 | |
|   mod2word *w;
 | |
|   int b;
 | |
| 
 | |
|   if (row<0 || row>=mod2dense_rows(m) || col<0 || col>=mod2dense_cols(m))
 | |
|   { fprintf(stderr,"mod2dense_flip: row or column index out of bounds\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   b = 1 ^ mod2_getbit (m->col[col][row>>mod2_wordsize_shift], 
 | |
|                        row&mod2_wordsize_mask);
 | |
| 
 | |
|   w = &m->col[col][row>>mod2_wordsize_shift];
 | |
| 
 | |
|   *w = b ? mod2_setbit1(*w,row&mod2_wordsize_mask) 
 | |
|          : mod2_setbit0(*w,row&mod2_wordsize_mask);
 | |
| 
 | |
|   return b;
 | |
| }
 | |
| 
 | |
| 
 | |
| /* COMPUTE THE TRANSPOSE OF A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| void mod2dense_transpose
 | |
| ( mod2dense *m,		/* Matrix to compute transpose of (left unchanged) */
 | |
|   mod2dense *r		/* Result of transpose operation */
 | |
| )
 | |
| {
 | |
|   mod2word w, v, *p;
 | |
|   int k1, j1, i2, j2;
 | |
| 
 | |
|   if (mod2dense_rows(m)!=mod2dense_cols(r) 
 | |
|    || mod2dense_cols(m)!=mod2dense_rows(r))
 | |
|   { fprintf(stderr,
 | |
|      "mod2dense_transpose: Matrices have incompatible dimensions\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   if (r==m)
 | |
|   { fprintf(stderr, 
 | |
|      "mod2dense_transpose: Result matrix is the same as the operand\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   mod2dense_clear(r);
 | |
| 
 | |
|   for (j1 = 0; j1<mod2dense_cols(m); j1++)
 | |
|   { 
 | |
|     i2 = j1 >> mod2_wordsize_shift;
 | |
|     v = 1 << (j1 & mod2_wordsize_mask);
 | |
| 
 | |
|     p = m->col[j1];
 | |
|     k1 = 0;
 | |
| 
 | |
|     for (j2 = 0; j2<mod2dense_cols(r); j2++)
 | |
|     { if (k1==0)
 | |
|       { w = *p++;
 | |
|         k1 = mod2_wordsize;
 | |
|       }
 | |
|       if (w&1)
 | |
|       { r->col[j2][i2] |= v;
 | |
|       }
 | |
|       w >>= 1;     
 | |
|       k1 -= 1;
 | |
|     }
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* ADD TWO DENSE MOD2 MATRICES. */
 | |
| 
 | |
| void mod2dense_add
 | |
| ( mod2dense *m1,	/* Left operand of add */
 | |
|   mod2dense *m2,	/* Right operand of add */
 | |
|   mod2dense *r		/* Place to store result of add */
 | |
| )
 | |
| {
 | |
|   int j, k;
 | |
| 
 | |
|   if (mod2dense_rows(m1)!=mod2dense_rows(r) 
 | |
|    || mod2dense_cols(m1)!=mod2dense_cols(r) 
 | |
|    || mod2dense_rows(m2)!=mod2dense_rows(r)
 | |
|    || mod2dense_cols(m2)!=mod2dense_cols(r))
 | |
|   { fprintf(stderr,"mod2dense_add: Matrices have different dimensions\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   for (j = 0; j<mod2dense_cols(r); j++)
 | |
|   { for (k = 0; k<r->n_words; k++)
 | |
|     { r->col[j][k] = m1->col[j][k] ^ m2->col[j][k];
 | |
|     }
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* MULTIPLY TWO DENSE MOD2 MATRICES. 
 | |
| 
 | |
|    The algorithm used runs faster if the second matrix (right operand of the
 | |
|    multiply) is sparse, but it is also appropriate for dense matrices.  This
 | |
|    procedure could be speeded up a bit by replacing the call of mod2dense_get
 | |
|    with in-line code that avoids division, but this doesn't seem worthwhile
 | |
|    at the moment. 
 | |
| */
 | |
| 
 | |
| void mod2dense_multiply 
 | |
| ( mod2dense *m1, 	/* Left operand of multiply */
 | |
|   mod2dense *m2,	/* Right operand of multiply */
 | |
|   mod2dense *r		/* Place to store result of multiply */
 | |
| )
 | |
| {
 | |
|   int i, j, k;
 | |
| 
 | |
|   if (mod2dense_cols(m1)!=mod2dense_rows(m2) 
 | |
|    || mod2dense_rows(m1)!=mod2dense_rows(r) 
 | |
|    || mod2dense_cols(m2)!=mod2dense_cols(r))
 | |
|   { fprintf(stderr,
 | |
|      "mod2dense_multiply: Matrices have incompatible dimensions\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   if (r==m1 || r==m2)
 | |
|   { fprintf(stderr,
 | |
|       "mod2dense_multiply: Result matrix is the same as one of the operands\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   mod2dense_clear(r);
 | |
| 
 | |
|   for (j = 0; j<mod2dense_cols(r); j++)
 | |
|   { for (i = 0; i<mod2dense_rows(m2); i++)
 | |
|     { if (mod2dense_get(m2,i,j))
 | |
|       { for (k = 0; k<r->n_words; k++)
 | |
|         { r->col[j][k] ^= m1->col[i][k];
 | |
|         }
 | |
|       }
 | |
|     }
 | |
|   }
 | |
| }
 | |
| 
 | |
| 
 | |
| /* SEE WHETHER TWO DENSE MOD2 MATRICES ARE EQUAL. */
 | |
| 
 | |
| int mod2dense_equal
 | |
| ( mod2dense *m1,
 | |
|   mod2dense *m2
 | |
| )
 | |
| {
 | |
|   int k, j, w;
 | |
|   mod2word m;
 | |
| 
 | |
|   if (mod2dense_rows(m1)!=mod2dense_rows(m2) 
 | |
|    || mod2dense_cols(m1)!=mod2dense_cols(m2))
 | |
|   { fprintf(stderr,"mod2dense_equal: Matrices have different dimensions\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   w = m1->n_words;
 | |
| 
 | |
|   /* Form a mask that has 1s in the lower bit positions corresponding to
 | |
|      bits that contain information in the last word of a matrix column. */
 | |
| 
 | |
|   m = (1 << (mod2_wordsize - (w*mod2_wordsize-m1->n_rows))) - 1;
 | |
|   
 | |
|   for (j = 0; j<mod2dense_cols(m1); j++)
 | |
|   {
 | |
|     for (k = 0; k<w-1; k++)
 | |
|     { if (m1->col[j][k] != m2->col[j][k]) return 0;
 | |
|     }
 | |
| 
 | |
|     if ((m1->col[j][k]&m) != (m2->col[j][k]&m)) return 0;
 | |
|   }
 | |
| 
 | |
|   return 1;
 | |
| }
 | |
| 
 | |
| 
 | |
| /* INVERT A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| int mod2dense_invert 
 | |
| ( mod2dense *m,		/* The matrix to find the inverse of (destroyed) */
 | |
|   mod2dense *r		/* Place to store the inverse */
 | |
| )
 | |
| {
 | |
|   mod2word *s, *t;
 | |
|   int i, j, k, n, w, k0, b0;
 | |
| 
 | |
|   if (mod2dense_rows(m)!=mod2dense_cols(m))
 | |
|   { fprintf(stderr,"mod2dense_invert: Matrix to invert is not square\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   if (r==m)
 | |
|   { fprintf(stderr, 
 | |
|       "mod2dense_invert: Result matrix is the same as the operand\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   n = mod2dense_rows(m);
 | |
|   w = m->n_words;
 | |
| 
 | |
|   if (mod2dense_rows(r)!=n || mod2dense_cols(r)!=n)
 | |
|   { fprintf(stderr,
 | |
|      "mod2dense_invert: Matrix to receive inverse has wrong dimensions\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   mod2dense_clear(r);
 | |
|   for (i = 0; i<n; i++) 
 | |
|   { mod2dense_set(r,i,i,1);
 | |
|   }
 | |
| 
 | |
|   for (i = 0; i<n; i++)
 | |
|   { 
 | |
|     k0 = i >> mod2_wordsize_shift;
 | |
|     b0 = i & mod2_wordsize_mask;
 | |
| 
 | |
|     for (j = i; j<n; j++) 
 | |
|     { if (mod2_getbit(m->col[j][k0],b0)) break;
 | |
|     }
 | |
| 
 | |
|     if (j==n) return 0;
 | |
| 
 | |
|     if (j!=i)
 | |
|     {
 | |
|       t = m->col[i];
 | |
|       m->col[i] = m->col[j];
 | |
|       m->col[j] = t;
 | |
| 
 | |
|       t = r->col[i];
 | |
|       r->col[i] = r->col[j];
 | |
|       r->col[j] = t;
 | |
|     }
 | |
| 
 | |
|     for (j = 0; j<n; j++)
 | |
|     { if (j!=i && mod2_getbit(m->col[j][k0],b0))
 | |
|       { s = m->col[j];
 | |
|         t = m->col[i];
 | |
|         for (k = k0; k<w; k++) s[k] ^= t[k];
 | |
|         s = r->col[j];
 | |
|         t = r->col[i];
 | |
|         for (k = 0; k<w; k++) s[k] ^= t[k];
 | |
|       }
 | |
|     }
 | |
|   }
 | |
| 
 | |
|   return 1;
 | |
| }
 | |
| 
 | |
| 
 | |
| /* INVERT A DENSE MOD2 MATRIX WITH ROWS & COLUMNS SELECTED FROM BIGGER MATRIX.*/
 | |
| 
 | |
| int mod2dense_invert_selected
 | |
| ( mod2dense *m,		/* Matrix from which to pick a submatrix to invert */
 | |
|   mod2dense *r,		/* Place to store the inverse */
 | |
|   int *rows,		/* Set to indexes of rows used and not used */
 | |
|   int *cols		/* Set to indexes of columns used and not used */
 | |
| )
 | |
| {
 | |
|   mod2word *s, *t;
 | |
|   int i, j, k, n, n2, w, k0, b0, c, R;
 | |
| 
 | |
|   if (r==m)
 | |
|   { fprintf(stderr, 
 | |
|       "mod2dense_invert_selected2: Result matrix is the same as the operand\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   n = mod2dense_rows(m);
 | |
|   w = m->n_words;
 | |
| 
 | |
|   n2 = mod2dense_cols(m);
 | |
| 
 | |
|   if (mod2dense_rows(r)!=n || mod2dense_cols(r)!=n2)
 | |
|   { fprintf(stderr,
 | |
| "mod2dense_invert_selected2: Matrix to receive inverse has wrong dimensions\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   mod2dense_clear(r);
 | |
| 
 | |
|   for (i = 0; i<n; i++)
 | |
|   { rows[i] = i;
 | |
|   }
 | |
| 
 | |
|   for (j = 0; j<n2; j++)
 | |
|   { cols[j] = j;
 | |
|   }
 | |
| 
 | |
|   R = 0;
 | |
|   i = 0;
 | |
| 
 | |
|   for (;;)
 | |
|   { 
 | |
|     while (i<n-R)
 | |
|     {
 | |
|       k0 = rows[i] >> mod2_wordsize_shift;
 | |
|       b0 = rows[i] & mod2_wordsize_mask;
 | |
| 
 | |
|       for (j = i; j<n2; j++) 
 | |
|       { if (mod2_getbit(m->col[cols[j]][k0],b0)) break;
 | |
|       }
 | |
| 
 | |
|       if (j<n2) break;
 | |
| 
 | |
|       R += 1;
 | |
|       c = rows[i];
 | |
|       rows[i] = rows[n-R];
 | |
|       rows[n-R] = c;
 | |
| 
 | |
|     }
 | |
| 
 | |
|     if (i==n-R) break;
 | |
| 
 | |
|     c = cols[j];
 | |
|     cols[j] = cols[i];
 | |
|     cols[i] = c;
 | |
| 
 | |
|     mod2dense_set(r,rows[i],c,1);
 | |
| 
 | |
|     for (j = 0; j<n2; j++)
 | |
|     { if (j!=c && mod2_getbit(m->col[j][k0],b0))
 | |
|       { s = m->col[j];
 | |
|         t = m->col[c];
 | |
|         for (k = 0; k<w; k++) s[k] ^= t[k];
 | |
|         s = r->col[j];
 | |
|         t = r->col[c];
 | |
|         for (k = 0; k<w; k++) s[k] ^= t[k];
 | |
|       }
 | |
|     }
 | |
| 
 | |
|     i += 1;
 | |
|   }
 | |
| 
 | |
|   for (j = n-R; j<n; j++)
 | |
|   { s = r->col[cols[j]];
 | |
|     for (k = 0; k<w; k++) s[k] = 0;
 | |
|   }
 | |
| 
 | |
|   return R;
 | |
| }
 | |
| 
 | |
| 
 | |
| /* FORCIBLY INVERT A DENSE MOD2 MATRIX. */
 | |
| 
 | |
| int mod2dense_forcibly_invert 
 | |
| ( mod2dense *m, 	/* The matrix to find the inverse of (destroyed) */
 | |
|   mod2dense *r,		/* Place to store the inverse */
 | |
|   int *a_row,		/* Place to store row indexes of altered elements */
 | |
|   int *a_col		/* Place to store column indexes of altered elements */
 | |
| )
 | |
| {
 | |
|   mod2word *s, *t;
 | |
|   int i, j, k, n, w, k0, b0;
 | |
|   int u, c;
 | |
| 
 | |
|   if (mod2dense_rows(m)!=mod2dense_cols(m))
 | |
|   { fprintf(stderr,
 | |
|       "mod2dense_forcibly_invert: Matrix to invert is not square\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   if (r==m)
 | |
|   { fprintf(stderr, 
 | |
|       "mod2dense_forcibly_invert: Result matrix is the same as the operand\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   n = mod2dense_rows(m);
 | |
|   w = m->n_words;
 | |
| 
 | |
|   if (mod2dense_rows(r)!=n || mod2dense_cols(r)!=n)
 | |
|   { fprintf(stderr,
 | |
|  "mod2dense_forcibly_invert: Matrix to receive inverse has wrong dimensions\n");
 | |
|     exit(1);
 | |
|   }
 | |
| 
 | |
|   mod2dense_clear(r);
 | |
|   for (i = 0; i<n; i++) 
 | |
|   { mod2dense_set(r,i,i,1);
 | |
|   }
 | |
| 
 | |
|   for (i = 0; i<n; i++)
 | |
|   { a_row[i] = -1;
 | |
|     a_col[i] = i;
 | |
|   }
 | |
| 
 | |
|   for (i = 0; i<n; i++)
 | |
|   { 
 | |
|     k0 = i >> mod2_wordsize_shift;
 | |
|     b0 = i & mod2_wordsize_mask;
 | |
| 
 | |
|     for (j = i; j<n; j++) 
 | |
|     { if (mod2_getbit(m->col[j][k0],b0)) break;
 | |
|     }
 | |
| 
 | |
|     if (j==n)
 | |
|     { j = i;
 | |
|       mod2dense_set(m,i,j,1);
 | |
|       a_row[i] = i;
 | |
|     }
 | |
| 
 | |
|     if (j!=i)
 | |
|     { 
 | |
|       t = m->col[i];
 | |
|       m->col[i] = m->col[j];
 | |
|       m->col[j] = t;
 | |
| 
 | |
|       t = r->col[i];
 | |
|       r->col[i] = r->col[j];
 | |
|       r->col[j] = t;
 | |
| 
 | |
|       u = a_col[i];
 | |
|       a_col[i] = a_col[j];
 | |
|       a_col[j] = u;
 | |
|     }
 | |
| 
 | |
|     for (j = 0; j<n; j++)
 | |
|     { if (j!=i && mod2_getbit(m->col[j][k0],b0))
 | |
|       { s = m->col[j];
 | |
|         t = m->col[i];
 | |
|         for (k = k0; k<w; k++) s[k] ^= t[k];
 | |
|         s = r->col[j];
 | |
|         t = r->col[i];
 | |
|         for (k = 0; k<w; k++) s[k] ^= t[k];
 | |
|       }
 | |
|     }
 | |
|   }
 | |
| 
 | |
|   c = 0;
 | |
|   for (i = 0; i<n; i++)
 | |
|   { if (a_row[i]!=-1)
 | |
|     { a_row[c] = a_row[i];
 | |
|       a_col[c] = a_col[i];
 | |
|       c += 1;
 | |
|     }
 | |
|   }
 | |
| 
 | |
|   return c;
 | |
| }
 | 
