else if (mod2sparse_col(e1)<mod2sparse_col(e2))      { mod2sparse_insert(r,i,mod2sparse_col(e1));        e1 = mod2sparse_next_in_row(e1);      }      else      { mod2sparse_insert(r,i,mod2sparse_col(e2));        e2 = mod2sparse_next_in_row(e2);             }    }    while (!mod2sparse_at_end(e1))    { mod2sparse_insert(r,i,mod2sparse_col(e1));      e1 = mod2sparse_next_in_row(e1);    }    while (!mod2sparse_at_end(e2))    { mod2sparse_insert(r,i,mod2sparse_col(e2));      e2 = mod2sparse_next_in_row(e2);    }  }}/* MULTIPLY TWO SPARSE MOD2 MATRICES. */void mod2sparse_multiply ( mod2sparse *m1, 	/* Left operand of multiply */  mod2sparse *m2,	/* Right operand of multiply */  mod2sparse *r		/* Place to store result of multiply */){  mod2entry *e1, *e2;  int i, j, b;  if (mod2sparse_cols(m1)!=mod2sparse_rows(m2)    || mod2sparse_rows(m1)!=mod2sparse_rows(r)    || mod2sparse_cols(m2)!=mod2sparse_cols(r))  { fprintf (stderr,      "mod2sparse_multiply: Matrices have incompatible dimensions\n");    exit(1);  }  if (r==m1 || r==m2)  { fprintf(stderr,     "mod2sparse_multiply: Result matrix is the same as one of the operands\n");    exit(1);  }  mod2sparse_clear(r);  for (i = 0; i<mod2sparse_rows(m1); i++)  {     if (mod2sparse_at_end(mod2sparse_first_in_row(m1,i)))     { continue;    }    for (j = 0; j<mod2sparse_cols(m2); j++)    {       b = 0;      e1 = mod2sparse_first_in_row(m1,i);      e2 = mod2sparse_first_in_col(m2,j);      while (!mod2sparse_at_end(e1) && !mod2sparse_at_end(e2))      {         if (mod2sparse_col(e1)==mod2sparse_row(e2))        { b ^= 1;          e1 = mod2sparse_next_in_row(e1);          e2 = mod2sparse_next_in_col(e2);         }        else if (mod2sparse_col(e1)<mod2sparse_row(e2))        { e1 = mod2sparse_next_in_row(e1);        }        else        { e2 = mod2sparse_next_in_col(e2);               }      }      if (b)      { mod2sparse_insert(r,i,j);      }    }  }}/* MULTIPLY VECTOR BY SPARSE MATRIX. */void mod2sparse_mulvec( mod2sparse *m,	/* The sparse matrix, with M rows and N columns */  char *u,		/* The input vector, N long */  char *v		/* Place to store the result, M long */){  mod2entry *e;  int M, N;  int i, j;  M = mod2sparse_rows(m);  N = mod2sparse_cols(m);  for (i = 0; i<M; i++) v[i] = 0;  for (j = 0; j<N; j++)  { if (u[j])    { for (e = mod2sparse_first_in_col(m,j);           !mod2sparse_at_end(e);           e = mod2sparse_next_in_col(e))      { v[mod2sparse_row(e)] ^= 1;      }    }  }}/* COUNT ENTRIES IN A ROW. */int mod2sparse_count_row( mod2sparse *m,  int row){  mod2entry *e;  int count;  if (row<0 || row>=mod2sparse_rows(m))  { fprintf(stderr,"mod2sparse_count_row: row index out of bounds\n");    exit(1);  }  count = 0;  for (e = mod2sparse_first_in_row(m,row);       !mod2sparse_at_end(e);       e = mod2sparse_next_in_row(e))  { count += 1;  }  return count;}/* COUNT ENTRIES IN A COLUMN. */int mod2sparse_count_col( mod2sparse *m,  int col){  mod2entry *e;  int count;  if (col<0 || col>=mod2sparse_cols(m))  { fprintf(stderr,"mod2sparse_count_col: column index out of bounds\n");    exit(1);  }  count = 0;  for (e = mod2sparse_first_in_col(m,col);       !mod2sparse_at_end(e);       e = mod2sparse_next_in_col(e))  { count += 1;  }  return count;}/* ADD TO A ROW. */void mod2sparse_add_row( mod2sparse *m1,	/* Matrix containing row to add to */  int row1,		/* Index in this matrix of row to add to */  mod2sparse *m2,	/* Matrix containing row to add from */  int row2		/* Index in this matrix of row to add from */){  mod2entry *f1, *f2, *ft;  if (mod2sparse_cols(m1)<mod2sparse_cols(m2))  { fprintf (stderr,     "mod2sparse_add_row: row added to is shorter than row added from\n");    exit(1);  }  if (row1<0 || row1>=mod2sparse_rows(m1)    || row2<0 || row2>=mod2sparse_rows(m2))  { fprintf (stderr,"mod2sparse_add_row: row index out of range\n");    exit(1);  }  f1 = mod2sparse_first_in_row(m1,row1);  f2 = mod2sparse_first_in_row(m2,row2);  while (!mod2sparse_at_end(f1) && !mod2sparse_at_end(f2))  { if (mod2sparse_col(f1)>mod2sparse_col(f2))    { mod2sparse_insert(m1,row1,mod2sparse_col(f2));      f2 = mod2sparse_next_in_row(f2);    }    else    { ft = mod2sparse_next_in_row(f1);        if (mod2sparse_col(f1)==mod2sparse_col(f2))      { mod2sparse_delete(m1,f1);        f2 = mod2sparse_next_in_row(f2);      }      f1 = ft;    }  }  while (!mod2sparse_at_end(f2))  { mod2sparse_insert(m1,row1,mod2sparse_col(f2));    f2 = mod2sparse_next_in_row(f2);  }}/* ADD TO A COLUMN. */void mod2sparse_add_col( mod2sparse *m1,	/* Matrix containing column to add to */  int col1,		/* Index in this matrix of column to add to */  mod2sparse *m2,	/* Matrix containing column to add from */  int col2		/* Index in this matrix of column to add from */){  mod2entry *f1, *f2, *ft;  if (mod2sparse_rows(m1)<mod2sparse_rows(m2))  { fprintf (stderr,     "mod2sparse_add_col: Column added to is shorter than column added from\n");    exit(1);  }  if (col1<0 || col1>=mod2sparse_cols(m1)    || col2<0 || col2>=mod2sparse_cols(m2))  { fprintf (stderr,"mod2sparse_add_col: Column index out of range\n");    exit(1);  }  f1 = mod2sparse_first_in_col(m1,col1);  f2 = mod2sparse_first_in_col(m2,col2);  while (!mod2sparse_at_end(f1) && !mod2sparse_at_end(f2))  { if (mod2sparse_row(f1)>mod2sparse_row(f2))    { mod2sparse_insert(m1,mod2sparse_row(f2),col1);      f2 = mod2sparse_next_in_col(f2);    }    else    { ft = mod2sparse_next_in_col(f1);      if (mod2sparse_row(f1)==mod2sparse_row(f2))      { mod2sparse_delete(m1,f1);        f2 = mod2sparse_next_in_col(f2);      }      f1 = ft;    }  }  while (!mod2sparse_at_end(f2))  { mod2sparse_insert(m1,mod2sparse_row(f2),col1);    f2 = mod2sparse_next_in_col(f2);  }}/* FIND AN LU DECOMPOSITION OF A SPARSE MATRIX. */int mod2sparse_decomp( mod2sparse *A,	/* Input matrix, M by N */  int K,		/* Size of sub-matrix to find LU decomposition of */  mod2sparse *L,	/* Matrix in which L is stored, M by K */  mod2sparse *U,	/* Matrix in which U is stored, K by N */  int *rows,		/* Array where row indexes are stored, M long */  int *cols,		/* Array where column indexes are stored, N long */  mod2sparse_strategy strategy, /* Strategy to follow in picking rows/columns */  int abandon_number,	/* Number of columns to abandon at some point */  int abandon_when	/* When to abandon these columns */){    int *rinv, *cinv, *acnt, *rcnt;  mod2sparse *B;  int M, N;  mod2entry *e, *f, *fn, *e2;  int i, j, k, cc, cc2, cc3, cr2, pr;  int found, nnf;  M = mod2sparse_rows(A);  N = mod2sparse_cols(A);  if (mod2sparse_cols(L)!=K || mod2sparse_rows(L)!=M   || mod2sparse_cols(U)!=N || mod2sparse_rows(U)!=K)  { fprintf (stderr,      "mod2sparse_decomp: Matrices have incompatible dimensions\n");    exit(1);  }  if (abandon_number>N-K)  { fprintf(stderr,"Trying to abandon more columns than allowed\n");    exit(1);  }  rinv = chk_alloc (M, sizeof *rinv);  cinv = chk_alloc (N, sizeof *cinv);  if (abandon_number>0)  { acnt = chk_alloc (M+1, sizeof *acnt);  }  if (strategy==Mod2sparse_minprod)  { rcnt = chk_alloc (M, sizeof *rcnt);  }  mod2sparse_clear(L);  mod2sparse_clear(U);  /* Copy A to B.  B will be modified, then discarded. */  B = mod2sparse_allocate(M,N);  mod2sparse_copy(A,B);  /* Count 1s in rows of B, if using minprod strategy. */  if (strategy==Mod2sparse_minprod)  { for (i = 0; i<M; i++)     { rcnt[i] = mod2sparse_count_row(B,i);    }  }  /* Set up initial row and column choices. */  for (i = 0; i<M; i++) rows[i] = rinv[i] = i;  for (j = 0; j<N; j++) cols[j] = cinv[j] = j;   /* Find L and U one column at a time. */  nnf = 0;  for (i = 0; i<K; i++)  {     /* Choose the next row and column of B. */    switch (strategy)    {      case Mod2sparse_first:       {         found = 0;        for (k = i; k<N; k++)        { e = mod2sparse_first_in_col(B,cols[k]);          while (!mod2sparse_at_end(e))          { if (rinv[mod2sparse_row(e)]>=i)            { found = 1;              goto out_first;            }            e = mod2sparse_next_in_col(e);          }        }      out_first:        break;      }      case Mod2sparse_mincol:      {         found = 0;        for (j = i; j<N; j++)        { cc2 = mod2sparse_count_col(B,cols[j]);          if (!found || cc2<cc)          { e2 = mod2sparse_first_in_col(B,cols[j]);            while (!mod2sparse_at_end(e2))            { if (rinv[mod2sparse_row(e2)]>=i)              { found = 1;                cc = cc2;                e = e2;                k = j;                break;              }              e2 = mod2sparse_next_in_col(e2);            }          }        }        break;      }      case Mod2sparse_minprod:      {         found = 0;        for (j = i; j<N; j++)        { cc2 = mod2sparse_count_col(B,cols[j]);          e2 = mod2sparse_first_in_col(B,cols[j]);          while (!mod2sparse_at_end(e2))          { if (rinv[mod2sparse_row(e2)]>=i)            { cr2 = rcnt[mod2sparse_row(e2)];              if (!found || cc2==1 || (cc2-1)*(cr2-1)<pr)              { found = 1;                pr = cc2==1 ? 0 : (cc2-1)*(cr2-1);                e = e2;                k = j;              }            }            e2 = mod2sparse_next_in_col(e2);          }        }        break;      }      default:      { fprintf(stderr,"mod2sparse_decomp: Unknown stategy\n");        exit(1);      }    }    if (!found)     { nnf += 1;    }    /* Update 'rows' and 'cols'.  Looks at 'k' and 'e' found above. */    if (found)    {       if (cinv[mod2sparse_col(e)]!=k) abort();      cols[k] = cols[i];      cols[i] = mod2sparse_col(e);      cinv[cols[k]] = k;      cinv[cols[i]] = i;      k = rinv[mod2sparse_row(e)];      if (k<i) abort();      rows[k] = rows[i];      rows[i] = mod2sparse_row(e);      rinv[rows[k]] = k;      rinv[rows[i]] = i;    }    /* Update L, U, and B. */    f = mod2sparse_first_in_col(B,cols[i]);     while (!mod2sparse_at_end(f))    {       fn = mod2sparse_next_in_col(f);      k = mod2sparse_row(f);      if (rinv[k]>i)      { mod2sparse_add_row(B,k,B,mod2sparse_row(e));        if (strategy==Mod2sparse_minprod)         { rcnt[k] = mod2sparse_count_row(B,k);        }        mod2sparse_insert(L,k,i);      }      else if (rinv[k]<i)      { mod2sparse_insert(U,rinv[k],cols[i]);      }      else      { mod2sparse_insert(L,k,i);        mod2sparse_insert(U,i,cols[i]);      }      f = fn;    }    /* Get rid of all entries in the current column of B, just to save space. */    for (;;)    { f = mod2sparse_first_in_col(B,cols[i]);      if (mod2sparse_at_end(f)) break;      mod2sparse_delete(B,f);    }    /* Abandon columns of B with lots of entries if it's time for that. */    if (abandon_number>0 && i==abandon_when)    {       for (k = 0; k<M+1; k++)       { acnt[k] = 0;      }      for (j = 0; j<N; j++)       { k = mod2sparse_count_col(B,j);        acnt[k] += 1;      }      cc = abandon_number;      k = M;      while (acnt[k]<cc)      { cc -= acnt[k];        k -= 1;        if (k<0) abort();      }      cc2 = 0;      for (j = 0; j<N; j++)      { cc3 = mod2sparse_count_col(B,j);        if (cc3>k || cc3==k && cc>0)        { if (cc3==k) cc -= 1;          for (;;)          { f = mod2sparse_first_in_col(B,j);            if (mod2sparse_at_end(f)) break;            mod2sparse_delete(B,f);          }          cc2 += 1;        }      }      if (cc2!=abandon_number) abort();      if (strategy==Mod2sparse_minprod)      { for (j = 0; j<M; j++)         { rcnt[j] = mod2sparse_count_row(B,j);        }      }    }  }  /* Get rid of all entries in the rows of L past row K, after reordering. */  for (i = K; i<M; i++)  { for (;;)    { f = mod2sparse_first_in_row(L,rows[i]);      if (mod2sparse_at_end(f)) break;      mod2sparse_delete(L,f);    }  }  mod2sparse_free(B);  free(rinv);  free(cinv);  if (strategy==Mod2sparse_minprod) free(rcnt);  if (abandon_number>0) free(acnt);  return nnf;}/* SOLVE A LOWER-TRIANGULAR SYSTEM BY FORWARD SUBSTITUTION. */int mod2sparse_forward_sub( mod2sparse *L,	/* Matrix that is lower triangular after reordering */  int *rows,		/* Array of indexes (from 0) of rows for new order */  char *x,		/* Vector on right of equation, also reordered */  char *y		/* Place to store solution */){  int K, i, j, ii, b, d;  mod2entry *e;  K = mod2sparse_cols(L);  /* Make sure that L is lower-triangular, after row re-ordering. */  for (i = 0; i<K; i++)  { ii = rows ? rows[i] : i;    e = mod2sparse_last_in_row(L,ii);    if (!mod2sparse_at_end(e) && mod2sparse_col(e)>i)    { fprintf(stderr,        "mod2sparse_forward_sub: Matrix is not lower-triangular\n");      exit(1);    }  }  /* Solve system by forward substitution. */  for (i = 0; i<K; i++)  {     ii = rows ? rows[i] : i;    /* Look at bits in this row, forming inner product with partial        solution, and seeing if the diagonal is 1. */    d = 0;    b = 0;    for (e = mod2sparse_first_in_row(L,ii);          !mod2sparse_at_end(e);         e = mod2sparse_next_in_row(e))    {       j = mod2sparse_col(e);      if (j==i)      { d = 1;      }      else      { b ^= y[j];      }    }    /* Check for no solution if the diagonal isn't 1. */    if (!d && b!=x[ii])     { return 0;    }    /* Set bit of solution, zero if arbitrary. */    y[i] = b^x[ii];  }  return 1;}/* SOLVE AN UPPER-TRIANGULAR SYSTEM BY BACKWARD SUBSTITUTION. */int mod2sparse_backward_sub( mod2sparse *U,	/* Matrix that is upper triangular after reordering */  int *cols,		/* Array of indexes (from 0) of columns for new order */  char *y,		/* Vector on right of equation */  char *z		/* Place to store solution, also reordered */){  int K, i, j, ii, b, d;  mod2entry *e;  K = mod2sparse_rows(U);  /* Make sure that U is upper-triangular, after column re-ordering. */  for (i = 0; i<K; i++)  { ii = cols ? cols[i] : i;    e = mod2sparse_last_in_col(U,ii);    if (!mod2sparse_at_end(e) && mod2sparse_row(e)>i)    { fprintf(stderr,        "mod2sparse_backward_sub: Matrix is not upper-triangular\n");      exit(1);    }  }  /* Solve system by backward substitution. */  for (i = K-1; i>=0; i--)  {     ii = cols ? cols[i] : i;    /* Look at bits in this row, forming inner product with partial        solution, and seeing if the diagonal is 1. */    d = 0;    b = 0;    for (e = mod2sparse_first_in_row(U,i);          !mod2sparse_at_end(e);         e = mod2sparse_next_in_row(e))    {       j = mod2sparse_col(e);      if (j==ii)      { d = 1;      }      else      { b ^= z[j];      }    }    /* Check for no solution if the diagonal isn't 1. */    if (!d && b!=y[i])     { return 0;    }    /* Set bit of solution, zero if arbitrary. */    z[ii] = b^y[i];  }  return 1;}