Rの拡張でライフゲームを作ってみる

わざわざC言語で書くほどのものではなく、単に行列の受渡しを確認しただけですが、Rの中で使えるライフゲームを作ってみたので公開します。
以下のソースコードを LifeGame.c と名づけて保存します。

#include <R.h>
#include <Rinternals.h>

SEXP LifeGame(SEXP m)
{
  SEXP  ans, dim;
  int  nor, noc, i, j, nol, x, y, r, c;
/* 第一引数が行列か確認 */
  if(!isMatrix(m)){
    error("A matrix is required for the first argument.");
  }
/* 行列の行数、列数を得る */
  dim = getAttrib(m, R_DimSymbol);
  nor = INTEGER(dim)[0]; /* 行数 */
  noc = INTEGER(dim)[1]; /* 列数 */
/* 戻り値の行列を初期化（注意：REALSXPではなくINTSXPを指定） */
  PROTECT(ans = allocMatrix(INTSXP, nor, noc));
/* 端は消去 */
  for(i=0; i<nor; i++){
    INTEGER(ans)[i] = 0;
    INTEGER(ans)[i + nor*(noc-1)] = 0;
  }
  for(j=0; j<noc; j++){
    INTEGER(ans)[nor*j] = 0;
    INTEGER(ans)[nor-1 + nor*j] = 0;
  }
/* 中央部分だけを処理 */
  for(i=1; i<nor-1; i++){
    for(j=1; j<noc-1; j++){
/* 周辺の8マスの状態を確認 */
      nol = 0;
      for(y=-1; y<=1; y++){
        for(x=-1; x<=1; x++){
          if(!y && !x)
            continue;
          if(0<INTEGER(m)[i+y + nor*(j+x)])
            nol++;
        }
      }
/* 周囲のマスの状態で生死を決定 */
      if(0<INTEGER(m)[i + nor*j])
        INTEGER(ans)[i + nor*j] = nol<=1 ? 0 : (nol>=4 ? 0 : 1);
/* 周囲のマスの状態で発生を決定 */
      else
        INTEGER(ans)[i + nor*j] = nol==3 ? 1 : 0;
    }
  }
  UNPROTECT(1);
  return(ans);
}

コンパイルしてDLL化します。

R CMD SHLIB LifeGame.c

Rを起動して、例えばグライダーを飛ばして遊びます。行列内の1と0で識別しているので、分かりづらいですね。

# ライブラリを呼び出す
dyn.load(paste("LifeGame", .Platform$dynlib.ext, sep = ""))
# Rの関数でラッピングして使う
lgame <- function(m){
  .Call("LifeGame", m)
}

# グライダー
mg <-matrix(as.integer(c(
1,1,1,
1,0,0,
0,1,0)), 3, 3, byrow=TRUE)
# 移動物体なので広めの行列を作る
m <- matrix(as.integer(rep(0,12*12)),12,12)
# 広い行列にグライダーを写す
for(r in 1:nrow(mg)){
  for(c in 1:ncol(mg)){
    m[r+nrow(m)-nrow(mg)-1,c+ncol(m)-ncol(mg)-1] <- mg[r,c]
  }
}
# 移動させてみる
for(i in 1:15){
  print(m)
  m <- lgame(m)
}