// To compile and run under Linux:
//   g++ -O3 -o bcd bcd.cpp -lm
//   ./bcd

// Adjustable input sizes and parameters:
#define n 30
#define k 11 // 0<=k<=n
#define w 8 // 0<=w<=n
#define k1 6 // 0<=k1
#define k2 5 // 0<=k2; k1+k2=k
#define p1 2 // 0<=p1<=k1
#define p2 1 // 0<=p2<=k2
#define q1 1 // 0<=q1
#define q2 2 // 0<=q2
#define l1 4 // q1<=l1
#define l2 5 // q2<=l2
// 0<=w-p1-p2-q1-q2<=n-k-l1-l2

#include <iostream>
using std::cout;
using std::ostream;
#include <stdlib.h>
#include <math.h>

long long cost = 0;
// See Sections 5 and 6 of paper for discussion of the cost model.

class bit {
  int x;
public:
  // Initialization, copying, comparison:
  bit() { x = 0; }
  bit(int y) { x = y; }
  bit& operator=(int y) { x = y; return *this; }
  bit& operator=(const bit &a) { x = a.x; return *this; }
  int iszero() { return x == 0; }
  int isnonzero() { return x != 0; }
  friend int operator==(const bit &a,const bit &b) {
    return a.x == b.x;
  }
  // Arithmetic:
  friend bit operator+(const bit &a,const bit &b) {
    cost += 1;
    return a.x ^ b.x;
  }
  friend bit operator*(const bit &a,const bit &b) {
    cost += 1;
    return a.x & b.x;
  }
  bit& operator+=(const bit &a) {
    cost += 1;
    this->x = this->x ^ a.x;
  }
  bit& operator*=(const bit &a) {
    cost += 1;
    this->x = this->x & a.x;
  }
  // Output:
  friend ostream& operator<<(ostream &os,const bit &b) {
    os << b.x;
    return os;
  }
} ;

#if l1+l2>30
error! l1+l2 too big
#endif

class lbits {
  int x; // 0 <= x < 2^(l1+l2)-1
public:
  // Initialization, comparison:
  lbits() {
    x = 0;
  }
  lbits(bit *a) {
    x = 0;
    for (int i = 0;i < l1+l2;++i)
      if (a[i].isnonzero())
        x |= (1 << i);
  }
  friend int operator==(const lbits &a,const lbits &b) {
    return a.x == b.x;
  }
  int index() {
    return x;
  }
  // Arithmetic:
  friend lbits operator+(const lbits &a,const lbits &b) {
    lbits result;
    cost += l1+l2;
    result.x = a.x ^ b.x;
    return result;
  }
  lbits plusbit(int pos) {
    lbits result;
    cost += 1;
    result.x = x ^ (1 << pos);
    return result;
  }
} ;

int check(int i)
{
  if (!i) {
    cout << "assertion failed\n";
    exit(111);
  }
}

int randommod(int m)
{
  return random() % m;
}

double binomial(int m,int p)
{
  double result = 1;
  for (int i = 0;i < p;++i) result = (result * (m - i)) / (i + 1);
  return result;
}

double L(int m,int p)
{
  double result = 0;
  for (int i = 0;i < p;++i) result += binomial(m,p-i);
  return result;
}

double Lprime(int m,int p)
{
  double result = 0;
  for (int i = 0;i < p;++i) result += binomial(m-i,p-i);
  return result;
}

int min(int a,int b)
{
  if (a < b) return a;
  return b;
}


main()
{
  check(k >= 0);
  check(k <= n);
  check(w >= 0);
  check(w <= n);
  check(k1 >= 0);
  check(k2 >= 0);
  check(k == k1 + k2);
  check(p1 >= 0);
  check(p1 <= k1);
  check(p2 >= 0);
  check(p2 <= k2);
  check(q1 >= 0);
  check(q2 >= 0);
  check(q1 <= l1);
  check(q2 <= l2);
  check(w-p1-p2-q1-q2 >= 0);
  check(w-p1-p2-q1-q2 <= n-k-l1-l2);
  check(p1 >= 1); // Extra restriction imposed by implementation.
  check(p2 >= 1); // Extra restriction imposed by implementation.

  long long cost1 = 0;
  long long cost2 = 0;
  long long cost3 = 0;
  long long cost4 = 0;
  long long cost5 = 0;
  long long iterations = 0;
  double gausscost = 0;
  double scost = 0;
  double tcost = 0;
  double matchcost = 0;
  double totaliter = 0;
  double totalitersquared = 0;
  long long challenges;
  long long loop;
  long long Smax = binomial(k1,p1) * binomial(l1,q1);
  long long Tmax = binomial(k2,p2) * binomial(l2,q2);
  int twol = 1 << (l1 + l2);
  int head[twol];
  for (int i = 0;i < twol;++i) head[i] = -1;

  for (challenges = 0;challenges < 10000;++challenges) {
    newchallenge:

    bit H[n-k][n]; // Random parity-check matrix.
    for (int i = 0;i < n-k;++i)
      for (int j = 0;j < n;++j)
        H[i][j] = randommod(2);
    // Paper explicitly assumes that H has full rank.
    // Instead of checking here we check as part of elimination later,
    // making the decoding algorithm robust against non-full-rank inputs.
    // If H does not have full rank then we go back to newchallenge;
    // statistics are collected only for valid (full-rank) inputs.
  
    bit target[n]; // Random weight-w target.
    for (int j = 0;j < n;++j) target[j] = 0;
    for (int i = 0;i < w;++i) {
      int j;
      do j = randommod(n); while(target[j].isnonzero());
      target[j] = 1;
    }
  
    bit s[n-k]; // Syndrome.
    for (int j = 0;j < n-k;++j) s[j] = 0;
    for (int i = 0;i < n-k;++i)
      for (int j = 0;j < n;++j)
        s[i] += H[i][j] * target[j];

    int printinput = 0;
    if (printinput) {
      for (int i = 0;i < n-k;++i) {
        for (int j = 0;j < n;++j) cout << H[i][j]; cout << "\n";
      }
      for (int i = 0;i < n-k;++i) cout << s[i]; cout << "\n";
    }
  
    int found = 0;
    for (loop = 0;!found;++loop) {
      cost1 = cost;
  
      // Step 1+4: Choose Z, compute U. Asymptotically negligible.
      // This implementation ignores most of the standard elimination speedups.
      // It does omit fruitless Gaussian-elimination steps, as in Stern;
      // see "finding an information set" in paper.
  
      int roworder[n-k];
      for (int j = 0;j < n-k;++j) {
        int i = randommod(j + 1);
        roworder[j] = roworder[i];
        roworder[i] = j;
      }
  
      bit U[n-k][n-k];
      for (int i = 0;i < n-k;++i)
        for (int j = 0;j < n-k;++j)
          U[i][j] = (i == j);
  
      bit UH[n-k][n];
      for (int i = 0;i < n-k;++i)
        for (int j = 0;j < n;++j)
          UH[i][j] = H[i][j];

      bit Us[n-k];
      for (int i = 0;i < n-k;++i) Us[i] = s[i];
  
      int Z[n];
      for (int j = 0;j < n;++j) {
        int i = randommod(j + 1);
        Z[j] = Z[i];
        Z[i] = j;
      }
      // This is not the final Z; elimination will remove columns.
  
      int rowcol[n-k];
      for (int r = 0;r < n-k;++r) {
        int i = roworder[r];
        int j;
        int jindex;
        if (loop == 0) {
          for (jindex = 0;jindex < n-r;++jindex)
            if (UH[i][Z[jindex]].isnonzero()) break;
          if (jindex == n-r) goto newchallenge; // H was not full rank
        }
        do {
          jindex = randommod(n-r);
          j = Z[jindex];
        } while (UH[i][j].iszero());
        Z[jindex] = Z[n-r - 1];
        rowcol[r] = j;
        for (int e = 0;e < n-k;++e)
          if (e != i && UH[e][j].isnonzero()) {
            // Row operations to eliminate UH[e][j]:
            UH[e][j] = 0;
            for (int f = 0;f < n-r-1;++f) {
              UH[e][Z[f]] += UH[i][Z[f]];
            }
            for (int f = 0;f <= r;++f) {
              U[e][roworder[f]] += U[i][roworder[f]];
            }
            Us[e] += Us[i];
          }
        int printelimination = 0;
        if (printelimination) {
          for (int i = 0;i < n-k;++i) {
            for (int j = 0;j < n-k;++j) cout << U[i][j];
            cout << " ";
            for (int j = 0;j < n;++j) cout << UH[i][j];
            cout << "\n";
          }
          cout << "\n";
        }
      }
  
      int checku = 0;
      if (checku) {
        for (int i = 0;i < n-k;++i)
          for (int j = 0;j < n;++j) {
            bit b = 0;
            for (int t = 0;t < n-k;++t)
              b += U[i][t] * H[t][j];
            check(b == UH[i][j]);
          }
        for (int i = 0;i < n-k;++i) {
          bit b = 0;
          for (int t = 0;t < n-k;++t)
            b += U[i][t] * s[t];
          check(b == Us[i]);
        }
        for (int i = 0;i < n-k;++i)
          for (int r = 0;r < n-k;++r)
            check(UH[i][rowcol[r]] == (i == roworder[r]));
      }

      // Z[0],...,Z[k-1] form the paper's set Z.
      // rowcol[0],...,rowcol[n-k-1] form the complement of Z.
      // roworder[0],...,roworder[n-k-1] are the corresponding rows.
      // Could permute Z[0],...,Z[k-1] into increasing order;
      // could permute rowcol and roworder, along with U and UH;
      // but the algorithm doesn't need this visual normalization.
  
      // Step 2: Partition Z into parts of sizes k1,k2.
      // We already adequately randomized the initial Z order,
      // so no need to further randomize here.
      int Fk1[k1];
      int Fk2[k2];
      for (int i = 0;i < k1;++i) Fk1[i] = Z[i];
      for (int i = 0;i < k2;++i) Fk2[i] = Z[k1 + i];
  
      // Step 3: Partition complement into parts of sizes l1,l2,n-k-l1-l2.
      // Again have adequate randomization already.
      int Fl1[l1];
      int Fl2[l2];
      int Fnkl1l2[n-k-l1-l2];
      for (int i = 0;i < l1;++i) Fl1[i] = rowcol[i];
      for (int i = 0;i < l2;++i) Fl2[i] = rowcol[l1 + i];
      for (int i = 0;i < n-k-l1-l2;++i) Fnkl1l2[i] = rowcol[l1 + l2 + i];

      // Step 4, continued: Write Z-indexed columns of UH as A1,A2.
      bit lprep[l1+l2];
      lbits A1[k];
      bit A2[k][n-k-l1-l2];
      for (int j = 0;j < k;++j) {
        for (int i = 0;i < l1+l2;++i) lprep[i] = UH[roworder[i]][Z[j]];
        A1[j] = lbits(lprep);
        for (int i = 0;i < n-k-l1-l2;++i) A2[j][i] = UH[roworder[l1+l2+i]][Z[j]];
      }

      // Step 5: Write Us as s1,s2.
      for (int i = 0;i < l1+l2;++i) lprep[i] = Us[roworder[i]];
      lbits s1(lprep);
      bit s2[n-k-l1-l2];
      for (int i = 0;i < n-k-l1-l2;++i) s2[i] = Us[roworder[l1+l2+i]];

      cost2 = cost;

      // Step 6: Compute set S of (A1x0+x1,x0,x1)
      // where x0 ranges over Fk1 with weight p1
      // and x1 ranges over Fl1 with weight q1.

      struct {
        lbits sum;
        int x0pos[p1];
        int x1pos[q1];
        int next;
      } S[Smax];
      int Slen = 0;

      int x0pos[p1]; // x0 has bits Fk1[x0pos[0]], Fk1[x0pos[1]], ...
      for (int j = 0;j < p1;++j) x0pos[j] = j;

      int x1pos[q1]; // x1 has bits Fl1[x1pos[0]], Fl1[x1pos[1]], ...
      for (int j = 0;j < q1;++j) x1pos[j] = j;

      lbits Spartialsums[p1+q1+1];
      Spartialsums[1] = A1[x0pos[0]];
      for (int j = 1;j < p1;++j)
        Spartialsums[1+j] = Spartialsums[j] + A1[x0pos[j]];
      for (int j = 0;j < q1;++j)
        Spartialsums[1+p1+j] = Spartialsums[p1+j].plusbit(x1pos[j]);

      for (;;) {
        S[Slen].sum = Spartialsums[p1+q1];
        for (int j = 0;j < p1;++j) S[Slen].x0pos[j] = x0pos[j];
        for (int j = 0;j < q1;++j) S[Slen].x1pos[j] = x1pos[j];
        S[Slen].next = head[S[Slen].sum.index()];
        head[S[Slen].sum.index()] = Slen;
        ++Slen;

        int j;
        for (j = p1+q1;j > 0;--j) {
          if (j > p1 && x1pos[j-1-p1] < l1-q1+j-1-p1) break;
          if (j <= p1 && x0pos[j-1] < k1-p1+j-1) break;
        }
        if (j == 0) break;
        for (;;) {
          if (j > p1) {
            ++x1pos[j-1-p1];
            Spartialsums[j] = Spartialsums[j-1].plusbit(x1pos[j-1-p1]);
          } else {
            ++x0pos[j-1];
            if (j == 1) Spartialsums[j] = A1[x0pos[0]];
	    else Spartialsums[j] = Spartialsums[j-1] + A1[x0pos[j-1]];
          }
          ++j;
          if (j > p1+q1) break;
          if (j == p1+1) x1pos[0] = -1;
          else if (j > p1) x1pos[j-1-p1] = x1pos[j-2-p1];
          else x0pos[j-1] = x0pos[j-2];
        }
      }
      check(Slen == Smax);

      cost3 = cost;

      // Step 7: Compute set T of (A1y0+y1+s,y0,y1)
      // where y0 ranges over Fk2 with weight p2
      // and y1 ranges over Fl2 with weight q2.

      struct {
        lbits sum;
        int y0pos[p2];
        int y1pos[q2];
      } T[Tmax];
      int Tlen = 0;

      int y0pos[p2]; // y0 has bits Fk2[y0pos[0]], Fk2[y0pos[1]], ...
      for (int j = 0;j < p2;++j) y0pos[j] = j;

      int y1pos[q2]; // y1 has bits Fl2[y1pos[0]], Fl2[y1pos[1]], ...
      for (int j = 0;j < q2;++j) y1pos[j] = j;

      lbits Tpartialsums[p2+q2+1];
      Tpartialsums[0] = s1;
      for (int j = 0;j < p2;++j)
        Tpartialsums[1+j] = Tpartialsums[j] + A1[k1 + y0pos[j]];
      for (int j = 0;j < q2;++j)
        Tpartialsums[1+p2+j] = Tpartialsums[p2+j].plusbit(l1 + y1pos[j]);

      for (;;) {
        T[Tlen].sum = Tpartialsums[p2+q2];
        for (int j = 0;j < p2;++j) T[Tlen].y0pos[j] = y0pos[j];
        for (int j = 0;j < q2;++j) T[Tlen].y1pos[j] = y1pos[j];
        ++Tlen;

        int j;
        for (j = p2+q2;j > 0;--j) {
          if (j > p2 && y1pos[j-1-p2] < l2-q2+j-1-p2) break;
          if (j <= p2 && y0pos[j-1] < k2-p2+j-1) break;
        }
        if (j == 0) break;
        for (;;) {
          if (j > p2) {
            ++y1pos[j-1-p2];
            Tpartialsums[j] = Tpartialsums[j-1].plusbit(l1 + y1pos[j-1-p2]);
          } else {
            ++y0pos[j-1];
            Tpartialsums[j] = Tpartialsums[j-1] + A1[k1 + y0pos[j-1]];
          }
          ++j;
          if (j > p2+q2) break;
          if (j == p2+1) y1pos[0] = -1;
          else if (j > p2) y1pos[j-1-p2] = y1pos[j-2-p2];
          else y0pos[j-1] = y0pos[j-2];
        }
      }
      check(Tlen == Tmax);

      cost4 = cost;

      // Step 8: For each (v,x0,x1) in S:
      //           For each (y0,y1) such that (v,y0,y1) in T:
      //             If wt(A2(x0+y0)+s2)=w-p1-p2-q1-q2:
      //               Output x0+y0+x1+y1+A2(x0+y0)+s2.

      for (int j = 0;j < Tlen;++j)
        for (int i = head[T[j].sum.index()];i >= 0;i = S[i].next) {
          int weight = 0;
          bit A2x0y0[n-k-l1-l2];
          int f;
          for (f = 0;f < n-k-l1-l2;++f) {
            A2x0y0[f] = s2[f];
            for (int u = 0;u < p1;++u) A2x0y0[f] += A2[S[i].x0pos[u]][f];
            for (int u = 0;u < p2;++u) A2x0y0[f] += A2[k1 + T[j].y0pos[u]][f];
            if (A2x0y0[f].isnonzero()) ++weight;
            if (weight > w-p1-p2-q1-q2) break;
          }
          if (f == n-k-l1-l2) {
            bit output[n];
            for (int u = 0;u < p1;++u) output[Fk1[S[i].x0pos[u]]] = 1;
            for (int u = 0;u < q1;++u) output[Fl1[S[i].x1pos[u]]] = 1;
            for (int u = 0;u < p2;++u) output[Fk2[T[j].y0pos[u]]] = 1;
            for (int u = 0;u < q2;++u) output[Fl2[T[j].y1pos[u]]] = 1;
            for (int u = 0;u < n-k-l1-l2;++u) output[Fnkl1l2[u]] = A2x0y0[u];
            int checkoutput = 0;
            if (checkoutput) {
              bit soutput[n-k];
              for (int u = 0;u < n-k;++u) soutput[u] = s[u];
              for (int v = 0;v < n-k;++v)
                for (int u = 0;u < n;++u)
                  soutput[v] += H[v][u] * output[u];
              for (int u = 0;u < n-k;++u) check(soutput[u].iszero());
            }
            int u;
            for (u = 0;u < n;++u) if (!(output[u] == target[u])) break;
            if (u == n) ++found;
          }
	}

      cost5 = cost;

      // Clean up and tally costs:
      for (int i = 0;i < Slen;++i)
        head[S[i].sum.index()] = -1;

      gausscost += cost2 - cost1;
      scost += cost3 - cost2;
      tcost += cost4 - cost3;
      matchcost += cost5 - cost4;
      iterations += 1;
    }

    totaliter += loop;
    totalitersquared += loop * (double) loop;
  }

  cout << "n = " << n << "\n";
  cout << "k = " << k << "\n";
  cout << "w = " << w << "\n";
  cout << "k1 = " << k1 << "\n";
  cout << "k2 = " << k2 << "\n";
  cout << "l1 = " << l1 << "\n";
  cout << "l2 = " << l2 << "\n";
  cout << "p1 = " << p1 << "\n";
  cout << "p2 = " << p2 << "\n";
  cout << "q1 = " << q1 << "\n";
  cout << "q2 = " << q2 << "\n";
  cout << "\n";

  cout << "challenges = " << challenges << "\n";
  cout << "iterations = " << iterations << "\n";
  cout << "\n";

  double iterpaper =
    binomial(n,w)
    / (binomial(n-k-l1-l2,w-p1-p2-q1-q2)
      *binomial(k1,p1)
      *binomial(k2,p2)
      *binomial(l1,q1)
      *binomial(l2,q2)
      );
  double iterpredicted = iterpaper;
  cout << "average iterations to target = " << totaliter / challenges
    << " +- " <<
    sqrt(totalitersquared / challenges - (totaliter / challenges) * (totaliter / challenges))
    << "\n";
  cout << "predicted = " << iterpredicted << "; ";
  cout << "actual/predicted = " << (totaliter / challenges) / iterpredicted << "\n";
  cout << "~68% estimate: 1 +- " <<
    sqrt(totalitersquared / challenges - (totaliter / challenges) * (totaliter / challenges))
    / sqrt(challenges) / iterpredicted << "\n";
  cout << "~99.7% estimate: 1 +- " <<
    3 * sqrt(totalitersquared / challenges - (totaliter / challenges) * (totaliter / challenges))
    / sqrt(challenges) / iterpredicted << "\n";
  cout << "paper = " << iterpaper << "; ";
  cout << "predicted/paper = " << iterpredicted / iterpaper << "\n";
  cout << "\n";

  double gausspredicted = (n-k-1)*(n-k)*(n+1)/2.0;
  double gausspaper = (n-k)*(n-k)*(n+k)/2.0;
  cout << "gauss cost/iteration = " << gausscost / iterations << "\n";
  cout << "predicted = " << gausspredicted << "; ";
  cout << "actual/predicted = " << (gausscost/iterations) / gausspredicted << "\n";
  cout << "paper = " << gausspaper << "; ";
  cout << "predicted/paper = " << gausspredicted / gausspaper << "\n";
  cout << "\n";

  double spredicted = 
    (l1+l2)*(Lprime(k1,p1)-(k1-p1+1))
    + min(1,q1)*binomial(k1,p1)*Lprime(l1,q1);
  double spaper = 
    (l1+l2)*(L(k1,p1)-k1)
    + min(1,q1)*binomial(k1,p1)*L(l1,q1);
  cout << "s cost/iteration = " << scost / iterations << "\n";
  cout << "predicted = " << spredicted << "; ";
  cout << "actual/predicted = " << (scost/iterations) / spredicted << "\n";
  cout << "paper = " << spaper << "; ";
  cout << "predicted/paper = " << spredicted / spaper << "\n";
  cout << "\n";

  double tpredicted = 
    (l1+l2)*(Lprime(k2,p2))
    + min(1,q2)*binomial(k2,p2)*Lprime(l2,q2);
  double tpaper = 
    (l1+l2)*(L(k2,p2))
    + min(1,q2)*binomial(k2,p2)*L(l2,q2);
  cout << "t cost/iteration = " << tcost / iterations << "\n";
  cout << "predicted = " << tpredicted << "; ";
  cout << "actual/predicted = " << (tcost/iterations) / tpredicted << "\n";
  cout << "paper = " << tpaper << "; ";
  cout << "predicted/paper = " << tpredicted / tpaper << "\n";
  cout << "\n";

  double matchpaper =
    2*(w-p1-p2-q1-q2+1)*(p1+p2)
    *binomial(k1,p1)*binomial(k2,p2)
    *binomial(l1,q1)*binomial(l2,q2)
    *1.0/twol;
  double matchpredicted = matchpaper;
  cout << "match cost/iteration = " << matchcost / iterations << "\n";
  cout << "predicted = " << matchpredicted << "; ";
  cout << "actual/predicted = " << (matchcost/iterations) / matchpredicted << "\n";
  cout << "paper = " << matchpaper << "; ";
  cout << "predicted/paper = " << matchpredicted / matchpaper << "\n";
  return 0;
}