Ticket #1505: removing_m4ri.patch

File removing_m4ri.patch, 70.9 kB (added by malb, 1 year ago)

  • /dev/null

    old new  
    1 /****************************************************************************** 
    2 * 
    3 *            M4RI: Method of the Four Russians Inversion 
    4 * 
    5 *       Copyright (C) 2007 Gregory Bard <gregory.bard@ieee.org>  
    6 * 
    7 *  Distributed under the terms of the GNU General Public License (GPL) 
    8 * 
    9 *    This code is distributed in the hope that it will be useful, 
    10 *    but WITHOUT ANY WARRANTY; without even the implied warranty of 
    11 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
    12 *    General Public License for more details. 
    13 * 
    14 *  The full text of the GPL is available at: 
    15 * 
    16 *                  http://www.gnu.org/licenses/ 
    17 ******************************************************************************/ 
    18  
    19 #include "brilliantrussian.h" 
    20 #include <stdlib.h> 
    21  
    22 int forceNonZero2PackedFlex(packedmatrix *m, int xstart, int xstop, int y) { 
    23   int i; 
    24  
    25   for (i=xstart; i<=xstop; i++) { 
    26     if (readPackedCell(m, i, y)==1) { 
    27       if (i!=xstart) rowSwapPacked(m, i, xstart); 
    28       return YES; 
    29     } 
    30   } 
    31    
    32   return NO; 
    33 } 
    34  
    35  
    36 int prepPackedFlex(packedmatrix *m, int ai, int k) { 
    37   int pc; /* pivot column */ 
    38   int tr; /* target row */ 
    39   int good; 
    40  
    41   int rank = 0; 
    42  
    43   for (pc=ai; pc<min(ai+k,m->cols); pc++) { 
    44     /* Step one, find a pivot row in this column.*/ 
    45     good=forceNonZero2PackedFlex(m, pc, min( ai+k*3-1, m->rows-1 ), pc); 
    46  
    47     if (good==NO) return rank; 
    48  
    49     for (tr=ai; tr<min(ai+k*3, m->rows); tr++) { 
    50       /* Step two, add this pivot row to other rows as needed. */ 
    51       if (tr==pc) continue; 
    52        
    53       if (readPackedCell(m, tr, pc)==0) continue; 
    54  
    55       rowAddPackedOffset(m, pc, tr, ai); 
    56     } 
    57     rank++; 
    58   } 
    59  
    60   return rank; 
    61 } 
    62  
    63 void combineFlex( packedmatrix * s1, int row1, int startblock1,  
    64                   packedmatrix * s2, int row2, int startblock2, 
    65                   packedmatrix * dest, int row3, int startblock3 ) { 
    66   int wide=s1->width - startblock1; 
    67   int i; 
    68  
    69   word *b1_ptr = s1->values + startblock1 + s1->rowswap[row1]; 
    70   word *b2_ptr = s2->values + startblock2 + s2->rowswap[row2]; 
    71   word *b3_ptr; 
    72  
    73   /* this is a quite likely case, and we treat is specially to ensure 
    74      cache register friendlyness. (Keep in mind that the x86 has only 
    75      four general purpose registers) */ 
    76   if( dest == s1 && row1 == row3 && startblock1 == startblock3) { 
    77  
    78     /* A fair amount of time is spent in iterating i, thus we lower 
    79        the burden a bit here. 
    80      */ 
    81     if(wide%2==0) { 
    82       for(i = wide>>1 ; i > 0 ; i--) { 
    83         *b1_ptr++ ^= *b2_ptr++; 
    84         *b1_ptr++ ^= *b2_ptr++; 
    85       } 
    86       return; 
    87  
    88     } else { 
    89       for(i = wide ; i > 0 ; i--) { 
    90         *b1_ptr++ ^= *b2_ptr++; 
    91       } 
    92       return; 
    93  
    94     } 
    95  
    96   } else { 
    97     b3_ptr = dest->values + startblock3 + dest->rowswap[row3]; 
    98  
    99     for(i = 0 ; i < wide ; i++) { 
    100       *b3_ptr++ = *b1_ptr++ ^ *b2_ptr++; 
    101     } 
    102     return; 
    103   } 
    104 } 
    105  
    106 void makeTablePackedFlex( packedmatrix *m, int ai, int k, 
    107                           packedmatrix *tablepacked, int *lookuppacked, int full) { 
    108   int homeblock= full ? 0 : ai/RADIX; 
    109   int i, rowneeded, id; 
    110   int twokay= TWOPOW(k); 
    111  
    112   lookuppacked[0]=0; 
    113  
    114   for (i=1; i<twokay; i++) { 
    115     rowneeded=codebook[k]->inc[i-1]+ai; 
    116  
    117     id=codebook[k]->ord[i]; 
    118  
    119     lookuppacked[id]=i; 
    120  
    121     combineFlex(          m, rowneeded, homeblock,  
    122                 tablepacked,       i-1, homeblock,  
    123                 tablepacked,         i, homeblock); 
    124   } 
    125 } 
    126  
    127  
    128 inline int getValueFlex(packedmatrix *m, int x, int y, int k) { 
    129   int truerow = m->rowswap[x]; 
    130   int block; 
    131   int spot; 
    132  
    133   word temp; 
    134  
    135   word *values = m->values; 
    136  
    137   /** 
    138    * there are two possible situations. Either all bits are in one 
    139    * word or they are spread across two words. 
    140    */ 
    141  
    142   if ( (y%RADIX + k -1 ) < RADIX ) { 
    143     /** 
    144      * everything happens in one word here 
    145      */ 
    146     temp =  values[ y / RADIX + truerow ]; // get the value 
    147     temp <<= y%RADIX; // clear upper bits 
    148     temp >>= RADIX - k; // clear lower bits and move to correct position. 
    149     return (int)temp; 
    150  
    151   } else {  
    152     /** 
    153      * two words are affected 
    154      */ 
    155     block = y / RADIX + truerow; // correct block 
    156     spot = (y + k ) % RADIX; // correct offset 
    157     // make room by shifting spot times to the right, and add stuff from the second word 
    158     temp = (values[block] << spot) | ( values[block + 1] >> (RADIX - spot) );  
    159     return ((int)temp & ((1<<k)-1)); // clear upper bits and return 
    160    } 
    161 } 
    162  
    163 void processRowPackedFlex(packedmatrix *m, int row, int homecol, int k, packedmatrix *tablepacked, int *lookuppacked) { 
    164   int blocknum=homecol/RADIX; 
    165  
    166   int value=getValueFlex(m, row, homecol, k); 
    167  
    168   int tablerow=lookuppacked[value]; 
    169  
    170   combineFlex(          m,      row, blocknum,  
    171               tablepacked, tablerow, blocknum,  
    172                         m,      row, blocknum); 
    173 } 
    174  
    175 void processPackedFlex(packedmatrix *m, int startrow, int stoprow, int startcol, int k, packedmatrix *tablepacked, int *lookuppacked) { 
    176   int i; 
    177   int blocknum=startcol/RADIX; 
    178   int value; 
    179   int tablerow; 
    180   word *b1_ptr,*b2_ptr; 
    181  
    182   // for optimization reasons we distinguish several cases here. 
    183  
    184   switch(m->width - startcol/RADIX) { 
    185  
    186   case 1: 
    187     // no loop needed as only one block is operated on. 
    188     for (i=startrow; i<=stoprow; i++) { 
    189       value=getValueFlex(m, i, startcol, k); 
    190       tablerow=lookuppacked[value]; 
    191       b1_ptr = m->values + blocknum + m->rowswap[i]; 
    192       b2_ptr = tablepacked->values + blocknum + tablepacked->rowswap[tablerow]; 
    193       *b1_ptr ^= *b2_ptr; 
    194     } 
    195     break; 
    196  
    197   case 2: 
    198     // two blocks, no loop 
    199     for (i=startrow; i<=stoprow; i++) { 
    200       value=getValueFlex(m, i, startcol, k); 
    201       tablerow=lookuppacked[value]; 
    202       b1_ptr = m->values + blocknum + m->rowswap[i]; 
    203       b2_ptr = tablepacked->values + blocknum + tablepacked->rowswap[tablerow]; 
    204       *b1_ptr++ ^= *b2_ptr++; 
    205       *b1_ptr ^= *b2_ptr; 
    206  
    207     } 
    208     break; 
    209  
    210   default: 
    211     // the real deal more than two blocks. 
    212     for (i=startrow; i<=stoprow; i++) { 
    213       processRowPackedFlex(m, i, startcol, k, tablepacked, lookuppacked); 
    214     } 
    215     break; 
    216   } 
    217  
    218 } 
    219  
    220 int doAByteColumnFlex(packedmatrix *m, int full, int k, int ai,  
    221                       packedmatrix *tablepacked, int *lookuppacked) { 
    222   int submatrixrank; 
    223    
    224   /* 
    225    * Stage 1: Denote the first column to be processed in a given 
    226    * iteration as a_i . Then, perform Gaussian elimination on the 
    227    * first 3k rows after and including the i-th row to produce an 
    228    * identity matrix in $a_{(i,i)} ... a_{(i+k-1),(i+k-1)}$ , and 
    229    * zeroes in $a_{(i+k),i} ... a_{(i+3k-1),(i+k-1)}$. 
    230    */ 
    231  
    232   submatrixrank=prepPackedFlex(m, ai, k); 
    233  
    234  
    235   if (submatrixrank!=k) return submatrixrank; 
    236  
    237   /* 
    238    * Stage 2: Construct a table consisting of the 2k binary strings of 
    239    * length k in a Gray Code.  Thus with only 2k vector additions, all 
    240    * possible linear combinations of these k rows have been 
    241    * precomputed. 
    242    */ 
    243  
    244   makeTablePackedFlex(m, ai, k, tablepacked, lookuppacked, 0); 
    245  
    246  
    247   /* 
    248    * Stage 3: One can rapidly process the remaining rows from i + 3k 
    249    * until row m (the last row) by using the table. For example, 
    250    * suppose the jth row has entries $a_{(j,i)} ... a_{(j,i+k-1)}$ in 
    251    * the columns being processed. Selecting the row of the table 
    252    * associated with this k-bit string, and adding it to row j will 
    253    * force the k columns to zero, and adjust the remaining columns 
    254    * from i + k to n in the appropriate way, as if Gaussian 
    255    * Elimination had been performed. 
    256   */ 
    257  
    258   processPackedFlex(m, ai+k*3, m->rows-1, ai, k, 
    259                     tablepacked, lookuppacked); 
    260  
    261   /* While the above form of the algorithm will reduce a system of 
    262    * boolean linear equations to unit upper triangular form, and thus 
    263    * permit a system to be solved with back substitution, the M4RI 
    264    * algorithm can also be used to invert a matrix, or put the system 
    265    * into reduced row echelon form (RREF). Simply run Stage 3 on rows 
    266    * 0 ... i - 1 as well as on rows i + 3k · · · m. This only affects 
    267    * the complexity slightly, changing the 2.5 coeffcient to 3 
    268    */ 
    269  
    270   if (full==YES) processPackedFlex(m, 0, ai-1, ai, k,  
    271                                    tablepacked, lookuppacked); 
    272  
    273   return submatrixrank; 
    274 } 
    275  
    276 int fourRussiansPackedFlex(packedmatrix *m, int full, int k,  
    277                        packedmatrix *tablepacked, int *lookuppacked) { 
    278   int i, submatrixrank; 
    279   int stop=min(m->rows, m->cols); 
    280   int lastokay=-1; 
    281  
    282   int rank = 0; 
    283  
    284   for (i=0; i<stop; i+=k) { 
    285     // not enough room for M4RI left. 
    286     if ( ((i+k*3-1)>=m->rows) || ((i+k-1)>=m->cols) ) { 
    287       return rank + gaussianPackedDelayed(m, lastokay+1, full); 
    288     } 
    289      
    290     submatrixrank=doAByteColumnFlex(m, full, k, i, tablepacked, lookuppacked); 
    291  
    292     if (submatrixrank!=k) { 
    293       // not full rank, use Gaussian elimination :-( 
    294       return rank + gaussianPackedDelayed(m, lastokay+1, full); 
    295     }  
    296  
    297     lastokay=i+k-1; 
    298  
    299     rank += submatrixrank; 
    300   } 
    301    
    302   return rank;  
    303 } 
    304  
    305 int simpleFourRussiansPackedFlex(packedmatrix *m, int full, int k) { 
    306   int size=m->cols; 
    307   int twokay=TWOPOW(k); 
    308  
    309   packedmatrix *mytable=createPackedMatrix(twokay, size); 
    310  
    311   int *mylookups=(int *)safeCalloc(twokay, sizeof(int)); 
    312  
    313   int rank=fourRussiansPackedFlex(m, full, k, mytable, mylookups); 
    314  
    315   free(mylookups); 
    316  
    317   destroyPackedMatrix(mytable); 
    318    
    319   return rank; 
    320 } 
    321  
    322 packedmatrix *invertPackedFlexRussian(packedmatrix *m,  
    323                                       packedmatrix *identity, int k) { 
    324   packedmatrix *big=concatPacked(m, identity); 
    325   int size=m->cols; 
    326   int twokay=TWOPOW(k); 
    327   int i; 
    328   packedmatrix *mytable=createPackedMatrix(twokay, size*2); 
    329  
    330   int *mylookups=(int *)safeCalloc(twokay, sizeof(int)); 
    331   int rank; 
    332   packedmatrix *answer; 
    333  
    334   rank=fourRussiansPackedFlex(big, YES, k, mytable, mylookups); 
    335  
    336   for(i=0; i < size; i++) { 
    337     if (!readPackedCell(big, i,i )) { 
    338       answer=NULL; 
    339       break; 
    340     } 
    341   } 
    342   if (i == size) 
    343     answer=copySubMatrixPacked(big, 0, size, size-1, size*2-1); 
    344  
    345   free(mylookups); 
    346   destroyPackedMatrix(mytable); 
    347   destroyPackedMatrix(big); 
    348    
    349   return answer; 
    350 } 
    351  
    352 packedmatrix *m4rmTransposePacked(packedmatrix *A, packedmatrix *B, int k) { 
    353   packedmatrix *AT, *BT, *CT, *C; 
    354    
    355   if(A->cols != B->rows) die("A cols need to match B rows"); 
    356    
    357   AT = transposePacked(A); 
    358   BT = transposePacked(B); 
    359  
    360   CT = m4rmPacked(BT,AT,k); 
    361  
    362   destroyPackedMatrix(AT); 
    363   destroyPackedMatrix(BT); 
    364  
    365   C = transposePacked(CT); 
    366   destroyPackedMatrix(CT); 
    367   return C; 
    368 } 
    369  
    370  
    371 packedmatrix *m4rmPacked(packedmatrix *A, packedmatrix *B, int k) { 
    372   int i,j; 
    373   int a,b,c; 
    374   unsigned int x; 
    375   packedmatrix *C; 
    376   int *lookuppacked; 
    377   packedmatrix *T; 
    378  
    379   if(A->cols != B->rows) die("A cols need to match B rows"); 
    380  
    381   a = A->rows; 
    382   b = A->cols; 
    383   c = B->cols; 
    384  
    385   T =createPackedMatrix(TWOPOW(k), c); 
    386   lookuppacked = (int *)safeCalloc(TWOPOW(k), sizeof(int)); 
    387  
    388   C = createPackedMatrix(a,c); 
    389  
    390   for(i=0 ; i < b/k ; i++) { 
    391  
    392     //Make a Gray Code table of all the 2^k linear combinations of the k rows of Bi . 
    393     //Call the xth row Tx . 
    394     makeTablePackedFlex( B, i*k, k, T, lookuppacked, 1); 
    395      
    396     for(j = 0; j<a ; j++) { 
    397  
    398       //Read the entries aj,(i-1)k+1 , aj,(i-1)k+2 , . . . , aj,(i-1)k+k . 
    399       //Let x be the k bit binary number formed by the concatenation of aj,(i-1)k+1 , . . . , aj,ik . 
    400       x = lookuppacked[getValueFlex(A, j, i*k, k)]; 
    401  
    402       //for h = 1, 2, . . . , c do 
    403       //    Calculate Cjh = Cjh + Txh. 
    404       combineFlex(C,j,0,  T,x,0,  C,j,0 ); 
    405     } 
    406   } 
    407  
    408   //handle rest 
    409   if (b%k) { 
    410     makeTablePackedFlex( B, b/k * k , b%k, T, lookuppacked, 1); 
    411      
    412     for(j = 0; j<a ; j++) { 
    413       x = lookuppacked[getValueFlex(A, j, i*k, b%k)]; 
    414       combineFlex(C,j,0, T,x,0,  C,j,0); 
    415     } 
    416   } 
    417   destroyPackedMatrix(T); 
    418   free(lookuppacked); 
    419   return C; 
    420 } 

  • /dev/null

    old new  
    1 #ifndef BRILLIANTRUSSIAN_H 
    2 #define BRILLIANTRUSSIAN_H 
    3  
    4 /****************************************************************************** 
    5 * 
    6 *            M4RI: Method of the Four Russians Inversion 
    7 * 
    8 *       Copyright (C) 2007 Gregory Bard <gregory.bard@ieee.org>  
    9 * 
    10 *  Distributed under the terms of the GNU General Public License (GPL) 
    11 * 
    12 *    This code is distributed in the hope that it will be useful, 
    13 *    but WITHOUT ANY WARRANTY; without even the implied warranty of 
    14 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
    15 *    General Public License for more details. 
    16 * 
    17 *  The full text of the GPL is available at: 
    18 * 
    19 *                  http://www.gnu.org/licenses/ 
    20 ******************************************************************************/ 
    21  
    22  
    23 #include <math.h> 
    24 #include <string.h> 
    25 #include "packedmatrix.h" 
    26  
    27 #define TWOPOW(i) (1<<(i)) 
    28  
    29 /**  
    30  * Finds a pivot row between xstart and xstop. The column where this 
    31  * pivot is search is y. Returns YES if such a pivot row was 
    32  * found. Also, the appropriate row is swapped to the top (== xstart). 
    33  * 
    34  * INPUT:  
    35  *     m -- matrix to operate on 
    36  *     xstart -- start row 
    37  *     xstop  -- stop row (including) 
    38  *     y -- column to read 
    39  * 
    40  * OUTPUT: 
    41  *     YES if a pivot row was found 
    42  */ 
    43  
    44 int forceNonZero2PackedFlex(packedmatrix *m, int xstart, int xstop, int y); 
    45  
    46  
    47 /***************************************************/ 
    48  
    49 /** 
    50  * Performs Gaussian elimination on a submatrix of 3k x k starting at 
    51  * point (homepoint, homepoint) of m. 
    52  * 
    53  * INPUT: 
    54  *     m -- matrix to operate on 
    55  *     homepoint -- row,col where to start 
    56  *     k -- the parameter k of M4RI 
    57  * 
    58  * OUTPUT: 
    59  *     returns rank of 3k x k submatrix. 
    60  */ 
    61  
    62 int prepPackedFlex(packedmatrix *m, int ai, int k); 
    63  
    64 /** Adds row1 of s1, starting with col1 to the end, to row2 of s2, 
    65  *  starting with col2 to the end. This gets stored in dest, in row3, 
    66  *  starting with col3  
    67  * 
    68  *  row3[col3:] = row1[col1:] + row2[col2:] 
    69  *  
    70  */ 
    71  
    72 void combineFlex( packedmatrix * s1, int row1, int startblock1,  
    73                   packedmatrix * s2, int row2, int startblock2, 
    74                   packedmatrix * dest, int row3, int startblock3 ); 
    75  
    76 /*  
    77  * Constructs all possible 2^k row combinations using the gray code 
    78  * table. 
    79  * 
    80  * INPUT: 
    81  *     m -- matrix to operate on 
    82  *     ai -- the starting position 
    83  *     k -- the k parameter of M4RI 
    84  *     tablepacked -- prealloced matrix of dimension 2^k x m->cols 
    85  *     lookuppacked -- prealloced table of length 2^k 
    86  *     full -- touch columns before ai? 
    87  * 
    88  */ 
    89  
    90 void makeTablePackedFlex( packedmatrix *m, int ai, int k, packedmatrix *tablepacked, int *lookuppacked, int full); 
    91  
    92 /** 
    93  * returns k bits/entries starting at m[x,y]. 
    94  * 
    95  * WARNING: Precondition is k < RADIX 
    96  */ 
    97  
    98 inline int getValueFlex(packedmatrix *m, int x, int y, int k); 
    99  
    100 /** 
    101  *  
    102  * Adds the correct row from tablepacked to the row 'row' in 'm' starting at homecol. 
    103  *  
    104  * INPUT: 
    105  *     m -- matrix to operate on 
    106  *     row -- the row which is operated on 
    107  *     homecol -- starting column for addition 
    108  *     k -- M4RI parameter, used for gray table lookup 
    109  *     tablepacked -- contains the correct row to be added 
    110  *     lookuptable -- contains row number to be addede 
    111  */ 
    112  
    113 void processRowPackedFlex(packedmatrix *m, int row, int homecol, int k, packedmatrix *tablepacked, int *lookuppacked); 
    114  
    115 /*  
    116  *  Iterates proccessPackedRowFlex from startrow to stoprow. 
    117  */ 
    118  
    119 void processPackedFlex(packedmatrix *m, int startrow, int stoprow, int startcol, int k, packedmatrix *tablepacked, int *lookuppacked); 
    120  
    121 /** 
    122  * This is the actual heart of the M4RI algorithm.  
    123  * 
    124  */ 
    125  
    126 int doAByteColumnFlex(packedmatrix *m, int full, int k, int ai, packedmatrix *tablepacked, int *lookuppacked); 
    127  
    128 /*  
    129  * 
    130  * 
    131  */ 
    132  
    133 int fourRussiansPackedFlex(packedmatrix *m, int full, int k, packedmatrix *tablepacked, int *lookuppacked); 
    134  
    135 /** 
    136  * Top level function to call for M4RI 
    137  */ 
    138  
    139 int simpleFourRussiansPackedFlex(packedmatrix *m, int full, int k); 
    140  
    141  
    142 /** 
    143  * Inverts the matrix m using the M4RI algorithm. To avoid recomputing 
    144  * the identity matrix over and over again, I may be passed in as 
    145  * identity parameter. 
    146  */ 
    147  
    148 packedmatrix *invertPackedFlexRussian(packedmatrix *m, packedmatrix *identity, int k); 
    149  
    150  
    151 packedmatrix *m4rmPacked(packedmatrix *A, packedmatrix *B, int k); 
    152  
    153 packedmatrix *m4rmTransposePacked(packedmatrix *A, packedmatrix *B, int k); 
    154  
    155 #endif //BRILLIANTRUSSIAN_H 

  • /dev/null

    old new  
    1 /** 
    2  * gray code generation used by the M4RI algorithm. 
    3  *  
    4  * AUTHOR: malb 
    5  */ 
    6  
    7 /****************************************************************************** 
    8 * 
    9 *            M4RI: Method of the Four Russians Inversion 
    10 * 
    11 *       Copyright (C) 2007 Gregory Bard <gregory.bard@ieee.org>  
    12 *       Copyright (C) 2007 Martin Albrecht <malb@informatik.uni-bremen.de>  
    13 * 
    14 *  Distributed under the terms of the GNU General Public License (GPL) 
    15 * 
    16 *    This code is distributed in the hope that it will be useful, 
    17 *    but WITHOUT ANY WARRANTY; without even the implied warranty of 
    18 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
    19 *    General Public License for more details. 
    20 * 
    21 *  The full text of the GPL is available at: 
    22 * 
    23 *                  http://www.gnu.org/licenses/ 
    24 ******************************************************************************/ 
    25  
    26 #include "grayflex.h" 
    27 #include <stdlib.h> 
    28 #include <stdio.h> 
    29   
    30 code **codebook; 
    31 const int MAXKAY = 16; 
    32  
    33 void printBitString(int number, int length){ 
    34   int i; 
    35   for(i=length-1 ; i>=0; i--) { 
    36     ((1<<i) & number) ? printf("1") : printf("0"); 
    37   } 
    38   printf("\n"); 
    39 } 
    40  
    41 int swap_bits(int v,int length) { 
    42  unsigned int r = v; // r will be reversed bits of v; first get LSB of v 
    43  int s = length - 1; // extra shift needed at end 
    44   
    45  for (v >>= 1; v; v >>= 1) 
    46     {    
    47       r <<= 1; 
    48       r |= v & 1; 
    49       s--; 
    50     } 
    51  r <<= s; 
    52  return r; 
    53 } 
    54  
    55 int grayCode(int number, int length) { 
    56   int lastbit = 0; 
    57   int res = 0; 
    58   int i,bit; 
    59   for(i=length-1; i>=0;  i--) { 
    60     bit = number & (1<<i); 
    61     res |= ((lastbit>>1) ^ bit);  
    62     lastbit = bit; 
    63   }; 
    64   return swap_bits(res,length) & ((1<<length)-1); 
    65   //return res; 
    66 } 
    67  
    68 void buildCodeFlex(int *ord, int *inc, int length) { 
    69   int i,j; 
    70  
    71   // this one is easy. 
    72   for(i=0 ; i < TWOPOW(length) ; i++) { 
    73     ord[i] = grayCode(i,length); 
    74   } 
    75  
    76   for(i = length ; i>0 ; i--) { 
    77     for(j=1 ; j < TWOPOW(i) + 1 ; j++) { 
    78       inc[j *TWOPOW(length-i) -1 ] = length - i; 
    79     } 
    80   } 
    81 } 
    82  
    83 void buildAllCodes() { 
    84   int k; 
    85   codebook=calloc(MAXKAY+1, sizeof(code *)); 
    86  
    87   for(k=1 ; k<MAXKAY+1; k++) { 
    88     codebook[k] = (code *)calloc(sizeof(code),1); 
    89     codebook[k]->ord =(int *)calloc(TWOPOW(k),sizeof(int)); 
    90     codebook[k]->inc =(int *)calloc(TWOPOW(k),sizeof(int)); 
    91     buildCodeFlex(codebook[k]->ord, codebook[k]->inc, k); 
    92   } 
    93 } 
    94  
    95 void destroyAllCodes() { 
    96   int i; 
    97   for(i=1; i<MAXKAY+1; i++) { 
    98     free(codebook[i]->inc); 
    99     free(codebook[i]->ord); 
    100     free(codebook[i]); 
    101   } 
    102   free(codebook); 
    103 } 

  • /dev/null

    old new  
    1 #ifndef GRAYFLEX_H 
    2 #define GRAYFLEX_H 
    3  
    4 /** 
    5  * gray code generation used by the M4RI algorithm. 
    6  *  
    7  * AUTHOR: malb 
    8  */ 
    9  
    10 /****************************************************************************** 
    11 * 
    12 *            M4RI: Method of the Four Russians Inversion 
    13 * 
    14 *       Copyright (C) 2007 Gregory Bard <gregory.bard@ieee.org>  
    15 *       Copyright (C) 2007 Martin Albrecht <malb@informatik.uni-bremen.de>  
    16 * 
    17 *  Distributed under the terms of the GNU General Public License (GPL) 
    18 * 
    19 *    This code is distributed in the hope that it will be useful, 
    20 *    but WITHOUT ANY WARRANTY; without even the implied warranty of 
    21 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
    22 *    General Public License for more details. 
    23 * 
    24 *  The full text of the GPL is available at: 
    25 * 
    26 *                  http://www.gnu.org/licenses/ 
    27 ******************************************************************************/ 
    28   
    29 struct codestruct { 
    30   int *ord; 
    31   int *inc; 
    32 }; 
    33  
    34 typedef struct codestruct /*renamed as*/ code; 
    35  
    36 extern code **codebook; 
    37  
    38 #define TWOPOW(i) (1<<(i)) 
    39  
    40 void printBitString(int number, int length); 
    41  
    42 /** 
    43  * swaps length bits in v naively. 
    44  * 
    45  * WARNING: Uppper bits of return value may contain garbage after 
    46  * operation. 
    47  */ 
    48 int swap_bits(int v,int length); 
    49  
    50 /** 
    51  * Returns the 'number'-th gray code entry for a gray code of length 
    52  * $2^{length}$. 
    53  *  
    54  * INPUT: 
    55  *     number -- index in the gray code table 
    56  *     length -- length of the gray code 
    57  * 
    58  * OUTPUT: 
    59  *      number-th gray code entry 
    60  * 
    61  * AUTHOR: malb 
    62  * 
    63  * THANKS: Soroosh Yazdani explained the repeated sum idea to me. 
    64  *  
    65  */ 
    66  
    67 int grayCode(int number, int length); 
    68  
    69 /** 
    70  * Fills in 'ord' and 'inc' with gray code data for a gray code of 
    71  * length $2^{length}$. 
    72  * 
    73  * INPUT: 
    74  *    ord -- will hold gray code data, must be preallocated with correct size 
    75  *    inc -- will hold some increment data, must be preallocated with correct size 
    76  * 
    77  * AUTHOR: malb 
    78  * 
    79  * THANKS: Robert Miller had the idea for a non-recursive implementation. 
    80  * 
    81  */ 
    82  
    83 void buildCodeFlex(int *ord, int *inc, int length); 
    84  
    85 void buildAllCodes(); 
    86  
    87 void destroyAllCodes(); 
    88  
    89 #endif //GRAYFLEX_H 

  • /dev/null

    old new  
    1 /** 
    2  * Include this file if you which to use this M4RI implementation and configure as you wish. 
    3  *  
    4  * AUTHOR: malb 
    5  */ 
    6  
    7 /****************************************************************************** 
    8 * 
    9 *   M4RI: Method of the Four Russians Inversion 
    10 * 
    11 *       Copyright (C) 2006, 2007 Gregory Bard <bardg@math.umd.edu> 
    12 * 
    13 *  Distributed under the terms of the GNU General Public License (GPL) 
    14 * 
    15 *    This code is distributed in the hope that it will be useful, 
    16 *    but WITHOUT ANY WARRANTY; without even the implied warranty of 
    17 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
    18 *    General Public License for more details. 
    19 * 
    20 *  The full text of the GPL is available at: 
    21 * 
    22 *                  http://www.gnu.org/licenses/ 
    23 ******************************************************************************/ 
    24  
    25 #if ! defined M4RI_H 
    26 #define M4RI_H 
    27  
    28 #include <stdio.h> 
    29 #include <stdlib.h> 
    30 #include <math.h> 
    31 #include "watch.h" 
    32 #include "packedmatrix.h" 
    33 #include "brilliantrussian.h" 
    34  
    35 #endif //M4RI_H 

  • /dev/null

    old new  
    1 /****************************************************************************** 
    2 * 
    3 *            M4RI: Method of the Four Russians Inversion 
    4 * 
    5 *       Copyright (C) 2007 Gregory Bard <gregory.bard@ieee.org>  
    6 * 
    7 *  Distributed under the terms of the GNU General Public License (GPL) 
    8 * 
    9 *    This code is distributed in the hope that it will be useful, 
    10 *    but WITHOUT ANY WARRANTY; without even the implied warranty of 
    11 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
    12 *    General Public License for more details. 
    13 * 
    14 *  The full text of the GPL is available at: 
    15 * 
    16 *                  http://www.gnu.org/licenses/ 
    17 ******************************************************************************/ 
    18  
    19 #include <stdio.h> 
    20 #include <stdlib.h> 
    21 #include "matrix.h" 
    22  
    23 /***************************************************/ 
    24 /***************************************************/ 
    25 void die(char *errormessage) { 
    26   /*This function prints the error message and exits 
    27     the program.*/ 
    28  
    29   fprintf(stderr, "\a%s\n", errormessage); 
    30   exit(1); 
    31 }  
    32  
    33 /***************************************************/ 
    34 /***************************************************/ 
    35 int min ( int a, int b) { 
    36   if (a<=b) return a; 
    37   else return b; 
    38 } 
    39  
    40 /***************************************************/ 
    41 /***************************************************/ 
    42  
    43 /* MEMLEAK, use free */ 
    44 void *safeCalloc( int count, int size ) { 
    45   /* this function calls calloc with the given inputs,  
    46      but dies with an error message if a NULL is returned */ 
    47  
    48   void *newthing=calloc( count, size ); 
    49   if (newthing==NULL) { 
    50     die("calloc returned NULL"); 
    51     return NULL; /* unreachable. */ 
    52   } 
    53   else return newthing; 
    54 } 
    55  
    56 /***************************************************/ 
    57 /***************************************************/ 
    58  
    59 /* MEMLEAK, use free */ 
    60 void *safeMalloc( int count, int size ) { 
    61   /* this function calls malloc with the inputs, which are 
    62      to be provided in calloc notation. If the result is 
    63      NULL, the program dies with an error message.*/ 
    64  
    65   void *newthing=malloc( count*size ); 
    66   if (newthing==NULL) { 
    67     die("malloc returned NULL"); 
    68     return NULL; /* unreachable */ 
    69   } 
    70   else return newthing; 
    71 } 
    72  
    73 /***************************************************/ 
    74 /***************************************************/ 
    75  
    76 /* MEMLEAK (see destroyMatrix) */ 
    77 matrix *createMatrix( int rows, int cols) { 
    78   matrix *newmatrix=(matrix *)safeMalloc(1,sizeof(matrix)); 
    79    
    80   bit *newcells=(bit *)safeCalloc(rows*cols,sizeof(bit)); 
    81   int *newrowswap=(int *)safeMalloc(rows, sizeof(int)); 
    82   int *newcolswap=(int *)safeMalloc(cols, sizeof(int)); 
    83   int i; 
    84  
    85   newmatrix->cells=newcells; 
    86   newmatrix->rowswap=newrowswap; 
    87   newmatrix->colswap=newcolswap; 
    88  
    89   newmatrix->rows=rows; 
    90   newmatrix->cols=cols; 
    91  
    92   for (i=0; i<rows; i++) { newmatrix->rowswap[i]=i; } 
    93   for (i=0; i<cols; i++) { newmatrix->colswap[i]=i; } 
    94  
    95   return newmatrix; 
    96 } 
    97  
    98 /***************************************************/ 
    99 /***************************************************/ 
    100 void destroyMatrix(matrix *condemned) { 
    101   if (condemned==NULL) { 
    102     printf("\aTried to free a NULL data structure.\n"); 
    103     return; 
    104   } 
    105   free(condemned->cells); 
    106   free(condemned->rowswap); 
    107   free(condemned->colswap); 
    108   free(condemned); 
    109   return; 
    110 } 
    111  
    112 /***************************************************/ 
    113 /***************************************************/ 
    114 void swapRow( matrix *m, int a, int b) { 
    115   int temp=m->rowswap[a]; 
    116   m->rowswap[a]=m->rowswap[b]; 
    117   m->rowswap[b]=temp; 
    118 } 
    119  
    120 /***************************************************/ 
    121 /***************************************************/ 
    122 void swapCol( matrix *m, int a, int b) { 
    123   int temp=m->colswap[a]; 
    124   m->colswap[a]=m->colswap[b]; 
    125   m->colswap[b]=temp; 
    126 } 
    127  
    128 /***************************************************/ 
    129 /***************************************************/ 
    130  
    131 void addTrueCell(matrix *thematrix, int srow, int scol, int drow, int dcol) { 
    132   int sindex=srow * (thematrix->cols) + scol; 
    133   int dindex=drow * (thematrix->cols) + dcol; 
    134  
    135   thematrix->cells[dindex]= 
    136     (thematrix->cells[dindex]+thematrix->cells[sindex]) % 2; 
    137 } 
    138  
    139 /***************************************************/ 
    140 /***************************************************/ 
    141  
    142 void writeTrueCell(matrix *thematrix, int row, int col, bit value) { 
    143   int index=row * (thematrix->cols) + col; 
    144  
    145   thematrix->cells[index]=value; 
    146 } 
    147  
    148 /***************************************************/ 
    149 /***************************************************/ 
    150  
    151 void addCell(matrix *thematrix, int srow, int scol, int drow, int dcol) { 
    152   int newsrow=thematrix->rowswap[srow]; 
    153   int newscol=thematrix->colswap[scol]; 
    154   int newdrow=thematrix->rowswap[drow]; 
    155   int newdcol=thematrix->colswap[dcol]; 
    156  
    157   addTrueCell(thematrix, newsrow, newscol, newdrow, newdcol); 
    158 } 
    159  
    160 /***************************************************/ 
    161 /***************************************************/ 
    162  
    163 void writeCell(matrix *thematrix, int row, int col, bit value) { 
    164   int newrow=thematrix->rowswap[row]; 
    165   int newcol=thematrix->colswap[col]; 
    166  
    167   writeTrueCell(thematrix, newrow, newcol, value); 
    168 } 
    169 /***************************************************/ 
    170 /***************************************************/ 
    171  
    172 bit readTrueCell(matrix *thematrix, int row, int col) { 
    173   int index=row * (thematrix->cols) + col; 
    174  
    175   return thematrix->cells[index]; 
    176 } 
    177 /***************************************************/ 
    178 /***************************************************/ 
    179  
    180 bit readCell(matrix *thematrix, int row, int col) { 
    181   int newrow=thematrix->rowswap[row]; 
    182   int newcol=thematrix->colswap[col]; 
    183  
    184   return readTrueCell(thematrix, newrow, newcol); 
    185 } 
    186  
    187 /***************************************************/ 
    188 /***************************************************/ 
    189  
    190 void zeroMatrix(matrix *thematrix) { 
    191   int i, j; 
    192  
    193   for (i=0; i<thematrix->rows; i++) { 
    194      for (j=0; j<thematrix->cols; j++) { 
    195        writeCell(thematrix, i, j, 0); 
    196     } 
    197   } 
    198 } 
    199  
    200 /***************************************************/ 
    201 /***************************************************/ 
    202  
    203 void printMatrix(matrix *thematrix) { 
    204   int i, j; 
    205  
    206   for (i=0; i<thematrix->rows; i++) { 
    207     printf("[ "); 
    208     for (j=0; j<thematrix->cols; j++) { 
    209       printf("%u ", readCell(thematrix, i,j) ); 
    210     } 
    211     printf("]\n"); 
    212   } 
    213  
    214   printf("\n"); 
    215 } 
    216  
    217 /***************************************************/ 
    218 /***************************************************/ 
    219  
    220 bit coinFlip() { 
    221  
    222  
    223   if (rand() < RAND_MAX/2) { 
    224     return 0; 
    225   }  else { 
    226     return 1; 
    227   } 
    228 } 
    229  
    230 /***************************************************/ 
    231 /***************************************************/ 
    232  
    233 void fillRandomMatrix(matrix *thematrix) { 
    234   int i, j; 
    235  
    236   for (i=0; i<thematrix->rows; i++) { 
    237      for (j=0; j<thematrix->cols; j++) { 
    238        writeCell(thematrix, i, j, coinFlip()); 
    239     } 
    240   } 
    241 } 
    242  
    243 /***************************************************/ 
    244 /***************************************************/ 
    245  
    246 void setIdentityMatrix(matrix *thematrix) { 
    247   int i; 
    248   int limit=min(thematrix->rows, thematrix->cols); 
    249  
    250   zeroMatrix(thematrix); 
    251  
    252   for (i=0; i<limit; i++) { 
    253     writeCell(thematrix, i, i, 1); 
    254   } 
    255 } 
    256  
    257 /***************************************************/ 
    258 /***************************************************/ 
    259  
    260 /* MEMLEAK, use destroyMatrix */ 
    261 matrix *concatMatrix(matrix *left, matrix *right) { 
    262   int i,j,newwidth; 
    263   matrix *newmatrix; 
    264  
    265   if (left->rows != right->rows) { 
    266     printf("Attempted to concat matrix of height %d with height %d.",  
    267            left->rows, right->rows); 
    268     die("Terminating"); 
    269   } 
    270  
    271   newwidth= left->cols + right->cols; 
    272  
    273   newmatrix=createMatrix(left->rows, newwidth); 
    274  
    275   for (i=0; i<left->rows; i++) { 
    276     for (j=0; j<newwidth; j++) { 
    277       if (j<left->cols) { 
    278         writeCell(newmatrix, i, j, readCell(left, i, j)); 
    279       } else { 
    280         writeCell(newmatrix, i, j, readCell(right, i, j-left->cols)); 
    281       } 
    282     } 
    283   } 
    284  
    285   return newmatrix;  
    286 } 
    287  
    288 /***************************************************/ 
    289 /***************************************************/ 
    290  
    291 int forceNonZero(matrix *m, int x, int y) { 
    292   int i; 
    293   int j=y; 
    294  
    295   //for (j=y; j<m->cols; j++) { 
    296     for (i=x; i<m->rows; i++) { 
    297       if (readCell(m, i, j)==1) { 
    298         if (i!=x) swapRow(m, i, x); 
    299         if (j!=y) swapCol(m, j, y); 
    300         return YES; 
    301       } 
    302     } 
    303     //} 
    304  
    305   return NO; 
    306 } 
    307  
    308 /***************************************************/ 
    309 /***************************************************/ 
    310  
    311 void rowAdd(matrix *m, int rowa, int rowb) { 
    312   int i; 
    313  
    314   for (i=0; i<m->cols; i++) { 
    315     /*    writeCell(m, rowb, i, 
    316           (readCell(m, rowa, i) + readCell(m, rowb, i)) % 2 );*/ 
    317     addCell(m, rowa, i, rowb, i); 
    318   } 
    319 } 
    320  
    321 /***************************************************/ 
    322 /***************************************************/ 
    323  
    324 int delayedGaussian(matrix *m, int firstcol, int full) { 
    325   int pc; /* pivot column */ 
    326   int tr; /* target row */ 
    327   int stop=min(m->rows, m->cols); 
    328   int good; 
    329  
    330   for (pc=firstcol; pc<stop; pc++) { 
    331     /* Step one, find a pivot row in this column.*/ 
    332     good=forceNonZero(m, pc, pc); 
    333     if (good==NO) return NO; 
    334  
    335     for (tr=( full ? 0 : pc+1 ); tr<m->rows; tr++) { 
    336       /* Step two, add this pivot row to other rows as needed. */ 
    337       if (tr==pc) continue; 
    338        
    339       if (readCell(m, tr, pc)==0) continue; 
    340  
    341       rowAdd(m, pc, tr); 
    342     } 
    343   } 
    344  
    345   return YES; 
    346 } 
    347  
    348 int gaussian(matrix *m, int full) { return delayedGaussian(m,0, full); } 
    349  
    350 /***************************************************/ 
    351 /***************************************************/ 
    352  
    353 /* MEMLEAK, use destroyMatrix() */ 
    354 matrix *copySubMatrix(matrix *m, int rlow, int clow, int rhigh, int chigh) { 
    355   int i,j, newrows, newcols; 
    356   matrix *newmatrix; 
    357  
    358   if (rlow>rhigh) die("Bad args to subMatrix.\n"); 
    359   if (clow>chigh) die("Bad args to subMatrix.\n"); 
    360   if (clow>m->cols) die("Bad args to subMatrix.\n"); 
    361   if (chigh>m->cols) die("Bad args to subMatrix.\n"); 
    362   if (rlow>m->rows) die("Bad args to subMatrix.\n"); 
    363   if (rhigh>m->rows) die("Bad args to subMatrix.\n"); 
    364  
    365   newrows=rhigh-rlow+1; 
    366   newcols=chigh-clow+1; 
    367  
    368   newmatrix=createMatrix(newrows, newcols); 
    369  
    370   for (i=0; i<newrows; i++) { 
    371     for (j=0; j<newcols; j++) { 
    372         writeCell(newmatrix, i, j, readCell(m, i+rlow, j+clow)); 
    373     } 
    374   } 
    375  
    376   return newmatrix;  
    377 }   
    378  
    379 /***************************************************************/ 
    380  
    381 /* MEMLEAK, use destroyMatrix */ 
    382 matrix *cloneMatrix( matrix *a ) { 
    383   return copySubMatrix( a, 0, 0, (a->rows)-1, (a->cols)-1); 
    384 } 
    385  
    386 /***************************************************************/ 
    387 /***************************************************************/ 
    388  
    389 /* MEMLEAK, use destroyMatrix */ 
    390 matrix *matrixTimesMatrix(matrix *a, matrix *b) { 
    391   int i,j,k, total; 
    392   matrix *newmatrix; 
    393  
    394   if (a->cols!=b->rows) die("Incompatible Dimensions in Matrix Multiplic."); 
    395  
    396   newmatrix=createMatrix(a->rows, b->cols); 
    397  
    398   for (i=0; i<a->rows; i++) { 
    399     for (j=0; j<b->cols; j++) { 
    400       total=0; 
    401       for (k=0; k<a->cols; k++) { 
    402         total+= readCell(a, i, k) * readCell(b, k, j); 
    403       } 
    404 </