# r2dtable()関数

• こちらで、固定した周辺度数のランダム分割表の発生関数r2dtable()を教えていただきました
• 発生のアルゴリズムはPatefield's algorithmというものだそうです。
• RではCで書かれたR_r2dtable()関数を呼び出しています
• veganというパッケージには、行列(分割表)のパーミュテーション処理(らしき)ものがあるようです。こちらを参照。
```SEXP
R_r2dtable(SEXP n, SEXP r, SEXP c)
{
int nr, nc, *row_sums, *col_sums, i, *jwork;
int n_of_samples, n_of_cases;
double *fact;
SEXP ans, tmp;

nr = length(r);
nc = length(c);

/* Note that the R code in r2dtable() also checks for missing and
negative values.
Should maybe do the same here ...
*/
if(!isInteger(n) || (length(n) == 0) ||
!isInteger(r) || (nr <= 1) ||
!isInteger(c) || (nc <= 1))
error(_("invalid arguments"));

n_of_samples = INTEGER(n)[0];
row_sums = INTEGER(r);
col_sums = INTEGER(c);

/* Compute total number of cases as the sum of the row sums.
Note that the R code in r2dtable() also checks whether this is
the same as the sum of the col sums.
Should maybe do the same here ...
*/
n_of_cases = 0;
jwork = row_sums;
for(i = 0; i < nr; i++)
n_of_cases += *jwork++;

/* Log-factorials from 0 to n_of_cases.
(I.e., lgamma(1), ..., lgamma(n_of_cases + 1).)
*/
fact = (double *) R_alloc(n_of_cases + 1, sizeof(double));
fact[0] = 0.;
for(i = 1; i <= n_of_cases; i++)
fact[i] = lgammafn((double) (i + 1));

jwork = (int *) R_alloc(nc, sizeof(int));

PROTECT(ans = allocVector(VECSXP, n_of_samples));

GetRNGstate();

for(i = 0; i < n_of_samples; i++) {
PROTECT(tmp = allocMatrix(INTSXP, nr, nc));
rcont2(&nr, &nc, row_sums, col_sums, &n_of_cases, fact,
jwork, INTEGER(tmp));
SET_VECTOR_ELT(ans, i, tmp);
UNPROTECT(1);
}

PutRNGstate();

UNPROTECT(1);

return(ans);
}
```