src/rcont.c
int rcont2(double* matrix, int nr, const double* r, int nc, const double* c,
           double* ntot, bitgen_t* rstate) {
  int status = 0;
  if (*ntot == 0) // perform checks only once
    status = rcont_check(ntot, matrix, nr, r, nc, c);
  if (status != 0)
    return status;

  // jwork is folded into matrix using last row
  double* jwork = ptr(matrix, nr, nc, nr - 1, 0);
  for (int i = 0; i < nc; ++i) {
    jwork[i] = c[i];
  }

  double jc = *ntot;
  double ib = 0;
  // last row is not random due to constraint
  for (int l = 0; l < nr - 1; ++l) {
    double ia = r[l]; // first term
    if (ia == 0) {
      for (int i = 0; i < nc; ++i)
        *ptr(matrix, nr, nc, l, i) = 0;
      continue;
    }
    double ic = jc; // second term
    jc -= r[l];
    // last column is not random due to constraint
    for (int m = 0; m < nc - 1; ++m) {
      const double id = jwork[m]; // third term
      const double ie = ic; // eight term
      ic -= id;
      ib = ie - ia;
      // must be after ib calculation, which is used at the end
      if (c[m] == 0) {
        for (int i = 0; i < nr; ++i)
          *ptr(matrix, nr, nc, i, m) = 0;
        continue;
      }
      const double ii = ib - id; // forth term
      if (ie == 0) {
        for (int j = m; j < nc - 1; ++j)
          *ptr(matrix, nr, nc, l, j) = 0;
        ia = 0;
        break;
      }
      double z = random_standard_uniform(rstate);
      double nlm;
      l131: nlm = floor(ia * id / ie + 0.5);
      double x = exp(
          lgamma(ia + 1)
          + lgamma(ib + 1)
          + lgamma(ic + 1)
          + lgamma(id + 1)
          - lgamma(ie + 1)
          - lgamma(nlm + 1)
          - lgamma(id - nlm + 1)
          - lgamma(ia - nlm + 1)
          - lgamma(ii + nlm + 1)
      );
      if (x >= z) goto l160;
      double sumprb = x;
      double y = x;
      double nll = nlm;
      int lsp = 0;
      int lsm = 0;
      // increment entry at (l,m)
      double j;
      l140: j = (id - nlm) * (ia - nlm);
      if (j == 0) goto l156;
      nlm += 1;
      x *= j / (nlm * (ii + nlm));
      sumprb += x;
      if (sumprb >= z) goto l160;
      l150: if (lsm) goto l155;
      // decrement entry at (l,m)
      j = nll * (ii + nll);
      if (j == 0) goto l154;
      nll -= 1;
      y *= j / ((id - nll) * (ia - nll));
      sumprb += y;
      if (sumprb >= z) goto l159;
      if (!lsp) goto l140;
      goto l150;
      l154: lsm = 1;
      l155: if (!lsp) goto l140;
      z = random_standard_uniform(rstate) * sumprb;
      goto l131;
      l156: lsp = 1;
      goto l150;
      l159: nlm = nll;
      l160: *ptr(matrix, nr, nc, l, m) = nlm;
      ia -= nlm;
      jwork[m] -= nlm;
    }
    // compute entry in last column of matrix
    *ptr(matrix, nr, nc, l, nc-1) = ia;
  }
  // compute entries in last row of matrix
  // jwork is already last row of matrix, so nothing to be done up to nc - 2
  *ptr(matrix, nr, nc, nr-1, nc-1) = ib - *ptr(matrix, nr, nc, nr - 1, nc - 2);

  return 0;
}
/*
  Generate random two-way table with given marginal totals.

  Patefield's algorithm adapted from AS 159 Appl. Statist. (1981) vol. 30, no. 1.
  This algorithm has O(K log(N)) complexity in time for a table with K cells and N
  entries in total. It requires no extra space (apart from a few stack variables).

  The original FORTRAN code was hand-translated to C. Changes to the original:

  - The computation of a look-up table of log-factorials was replaced with
    calls to lgamma, which lifts the limitation that the code only works for tables
    with less then 5000 entries.
  - The data type of input and output arrays was changed to double to minimize
    type conversions. Users are responsible for passing only integral numbers.
  - The original implementation allocated a column vector JWORK, but this is not
    necessary. The vector can be folded into the last column of the output matrix.
  - The function uses Numpy's random number generator and distribution library.
  - The algorithm now handles zero entries in row or column vector. When a zero is
    encountered, the output matrix is filled with zeros along that row or column and
    the algorithm proceeds to the next entry.

  The argument ntot is used to detect whether the function has been run before and
  has to be zero initialised.
*/
Generate random two-way table with given marginal totals.
Patefield's algorithm adapted from AS 159 Appl. Statist.

The original FORTRAN code was hand-translated to C. Changes to the original.

The computation of a look-up table of log-factorials was replaced with
calls to lgamma, which lifts the limitation that the code only works for tables
with less then 5000 entries.
The data type of input and output arrays was changed to double to minimize
type conversions. Users are responsible for passing only integral numbers.
The original implementation allocated a column vector JWORK, but this is not
necessary. The vector can be folded into the last column of the output matrix.
The function uses Numpy's random number generator and distribution library.
The algorithm now handles zero entries in row or column vector. When a zero is
encountered, the output matrix is filled with zeros along that row or column and
the algorithm proceeds to the next entry.

The argument ntot is used to detect whether the function has been run before and
has to be zero initialised.
Generate random two - way table with given marginal totals . Patefield ' s algorithm adapted from AS 159 Appl . Statist . The original FORTRAN code was hand - translated to C . Changes to the original . The computation of a look - up table of log - factorials was replaced with calls to lgamma which lifts the limitation that the code only works for tables with less then 5000 entries . The data type of input and output arrays was changed to double to minimize type conversions . Users are responsible for passing only integral numbers . The original implementation allocated a column vector JWORK but this is not necessary . The vector can be folded into the last column of the output matrix . The function uses Numpy ' s random number generator and distribution library . The algorithm now handles zero entries in row or column vector . When a zero is encountered the output matrix is filled with zeros along that row or column and the algorithm proceeds to the next entry . The argument ntot is used to detect whether the function has been run before and has to be zero initialised .
int rcont2(double* matrix, int nr, const double* r, int nc, const double* c,
           double* ntot, bitgen_t* rstate) {
  int status = 0;
  if (*ntot == 0) 
    status = rcont_check(ntot, matrix, nr, r, nc, c);
  if (status != 0)
    return status;
  double* jwork = ptr(matrix, nr, nc, nr - 1, 0);
  for (int i = 0; i < nc; ++i) {
    jwork[i] = c[i];
  }
  double jc = *ntot;
  double ib = 0;
  for (int l = 0; l < nr - 1; ++l) {
    double ia = r[l]; 
    if (ia == 0) {
      for (int i = 0; i < nc; ++i)
        *ptr(matrix, nr, nc, l, i) = 0;
      continue;
    }
    double ic = jc; 
    jc -= r[l];
    for (int m = 0; m < nc - 1; ++m) {
      const double id = jwork[m]; 
      const double ie = ic; 
      ic -= id;
      ib = ie - ia;
      if (c[m] == 0) {
        for (int i = 0; i < nr; ++i)
          *ptr(matrix, nr, nc, i, m) = 0;
        continue;
      }
      const double ii = ib - id; 
      if (ie == 0) {
        for (int j = m; j < nc - 1; ++j)
          *ptr(matrix, nr, nc, l, j) = 0;
        ia = 0;
        break;
      }
      double z = random_standard_uniform(rstate);
      double nlm;
      l131: nlm = floor(ia * id / ie + 0.5);
      double x = exp(
          lgamma(ia + 1)
          + lgamma(ib + 1)
          + lgamma(ic + 1)
          + lgamma(id + 1)
          - lgamma(ie + 1)
          - lgamma(nlm + 1)
          - lgamma(id - nlm + 1)
          - lgamma(ia - nlm + 1)
          - lgamma(ii + nlm + 1)
      );
      if (x >= z) goto l160;
      double sumprb = x;
      double y = x;
      double nll = nlm;
      int lsp = 0;
      int lsm = 0;
      double j;
      l140: j = (id - nlm) * (ia - nlm);
      if (j == 0) goto l156;
      nlm += 1;
      x *= j / (nlm * (ii + nlm));
      sumprb += x;
      if (sumprb >= z) goto l160;
      l150: if (lsm) goto l155;
      j = nll * (ii + nll);
      if (j == 0) goto l154;
      nll -= 1;
      y *= j / ((id - nll) * (ia - nll));
      sumprb += y;
      if (sumprb >= z) goto l159;
      if (!lsp) goto l140;
      goto l150;
      l154: lsm = 1;
      l155: if (!lsp) goto l140;
      z = random_standard_uniform(rstate) * sumprb;
      goto l131;
      l156: lsp = 1;
      goto l150;
      l159: nlm = nll;
      l160: *ptr(matrix, nr, nc, l, m) = nlm;
      ia -= nlm;
      jwork[m] -= nlm;
    }
    *ptr(matrix, nr, nc, l, nc-1) = ia;
  }
  *ptr(matrix, nr, nc, nr-1, nc-1) = ib - *ptr(matrix, nr, nc, nr - 1, nc - 2);
  return 0;
}
int rcont2 ( double * matrix , int nr , const double * r , int nc , const double * c , double * ntot , bitgen_t * rstate ) { int status = 0 ; if ( * ntot == 0 ) status = rcont_check ( ntot , matrix , nr , r , nc , c ) ; if ( status != 0 ) return status ; double * jwork = ptr ( matrix , nr , nc , nr - 1 , 0 ) ; for ( int i = 0 ; i < nc ; ++ i ) { jwork [ i ] = c [ i ] ; } double jc = * ntot ; double ib = 0 ; for ( int l = 0 ; l < nr - 1 ; ++ l ) { double ia = r [ l ] ; if ( ia == 0 ) { for ( int i = 0 ; i < nc ; ++ i ) * ptr ( matrix , nr , nc , l , i ) = 0 ; continue ; } double ic = jc ; jc -= r [ l ] ; for ( int m = 0 ; m < nc - 1 ; ++ m ) { const double id = jwork [ m ] ; const double ie = ic ; ic -= id ; ib = ie - ia ; if ( c [ m ] == 0 ) { for ( int i = 0 ; i < nr ; ++ i ) * ptr ( matrix , nr , nc , i , m ) = 0 ; continue ; } const double ii = ib - id ; if ( ie == 0 ) { for ( int j = m ; j < nc - 1 ; ++ j ) * ptr ( matrix , nr , nc , l , j ) = 0 ; ia = 0 ; break ; } double z = random_standard_uniform ( rstate ) ; double nlm ; l131 : nlm = floor ( ia * id / ie + 0.5 ) ; double x = exp ( lgamma ( ia + 1 ) + lgamma ( ib + 1 ) + lgamma ( ic + 1 ) + lgamma ( id + 1 ) - lgamma ( ie + 1 ) - lgamma ( nlm + 1 ) - lgamma ( id - nlm + 1 ) - lgamma ( ia - nlm + 1 ) - lgamma ( ii + nlm + 1 ) ) ; if ( x >= z ) goto l160 ; double sumprb = x ; double y = x ; double nll = nlm ; int lsp = 0 ; int lsm = 0 ; double j ; l140 : j = ( id - nlm ) * ( ia - nlm ) ; if ( j == 0 ) goto l156 ; nlm += 1 ; x *= j / ( nlm * ( ii + nlm ) ) ; sumprb += x ; if ( sumprb >= z ) goto l160 ; l150 : if ( lsm ) goto l155 ; j = nll * ( ii + nll ) ; if ( j == 0 ) goto l154 ; nll -= 1 ; y *= j / ( ( id - nll ) * ( ia - nll ) ) ; sumprb += y ; if ( sumprb >= z ) goto l159 ; if ( ! lsp ) goto l140 ; goto l150 ; l154 : lsm = 1 ; l155 : if ( ! lsp ) goto l140 ; z = random_standard_uniform ( rstate ) * sumprb ; goto l131 ; l156 : lsp = 1 ; goto l150 ; l159 : nlm = nll ; l160 : * ptr ( matrix , nr , nc , l , m ) = nlm ; ia -= nlm ; jwork [ m ] -= nlm ; } * ptr ( matrix , nr , nc , l , nc - 1 ) = ia ; } * ptr ( matrix , nr , nc , nr - 1 , nc - 1 ) = ib - * ptr ( matrix , nr , nc , nr - 1 , nc - 2 ) ; return 0 ; }
Generate random two-way table with given marginal totals.
Generate random two - way table with given marginal totals .
// perform checks only once // jwork is folded into matrix using last row // last row is not random due to constraint // first term // second term // last column is not random due to constraint // third term // eight term // must be after ib calculation, which is used at the end // forth term // increment entry at (l,m) // decrement entry at (l,m) // compute entry in last column of matrix // compute entries in last row of matrix // jwork is already last row of matrix, so nothing to be done up to nc - 2