forループ

Rcpp で R コードを最適化する

Romain François

Consulting Datactive, ThinkR

C++ forループの4つの要素

初期化
継続条件
インクリメント
本体

forループ - 初期化

ループの開始時に実行される処理：

for( init ;  ;  ){

}

forループ - 継続条件

ループを続けるかどうかを制御する論理条件：

for(  ; condition ;  ){

}

forループ - インクリメント

各反復の終わりに実行される処理：

for(  ;  ; increment ){

}

forループ - 本体

各反復で実行される処理。ループの本体。

for(  ;  ;  ){
    body
}

典型的なforループ

for (int i=0; i<n; i++ ){
    // some code using i
}

典型的なforループ

for (int i=0; ; ){

}

典型的なforループ

for (int i=0; i<n; ){

}

典型的なforループ

for (int i=0; i<n; i++){

}

例：最初のn個の整数の合計

// [[Rcpp::export]]
int nfirst( int n ){
    if( n < 0 ) {
        stop( "n must be positive, I see n=%d", n ) ;
    }

    int result = 0 ;
    for( int i=0; i<n; i++){
        result = result + (i+1) ;
    }

    return result ;
}

forループからの脱出

// [[Rcpp::export]]
int nfirst( int n ){
    if( n < 0 ) {
        stop( "n must be positive, I see n=%d", n ) ;
    }    

    int result = 0 ;
    for( int i=0; i<n; i++){
        if( i == 13 ){
            Rprintf( "I cannot handle that, I am superstitious" ) ;
            break ;
        }
        result = result + (i+1) ;
    }

    return result ;
}

平方根を求めるニュートン反復法

$\sqrt{S}$ を求めることは、$f(x) = x^2 - S$ の根を求めることと同じです。

反復式は次のようになります：

$$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} = x_n - \frac{x^2_n - S}{2 x_n} = \frac{1}{2} \left( x_n + \frac{S}{x_n} \right) $$

アルゴリズム：

初期値 $x_0$ を設定する
上記の式を指定回数繰り返して $x$ を更新する

C++でのニュートン法

$$ x_{n+1} = \frac{1}{2} \left( x_n + \frac{S}{x_n} \right) $$

擬似コードに変換すると：

int n = ... // 反復回数
double res = ... // 初期化

for( int i=0; i<n; i++){
    // res の値を更新する
    // すなわち x_{n} から x_{n+1} を計算する
    res = ( res + S / res ) / 2.0 ; 
}

return res ;

では、練習しましょう！

Rcpp で R コードを最適化する