C++のポアソン乱数

しばらくRでrpoisを使ってシミュレーションしていると、C++になったときどうやってポアソン乱数を発生させるんだっけと調べることになるので、備忘録。
平均lambdaのポアソン乱数。

#include
#include
#include
#include
　　　using namespace std;
//関数宣言
int rpois(double);

int main(){
　　　//出力ファイル
　　　ofstream out ("rpois.txt", ios::out);
　　　if(!out){
　　　　　　cerr << "rpois.txt could not open." << endl;
　　　　　　exit(1);
　　　}

　　　//パラメータ設定
　　　int i, ri, max = 10000;
　　　double lambda = 3.0;
　　　srand( (unsigned)time(NULL));

　　　for(i=0;i < max;i++){
　　　　　　ri = rpois(lambda);
　　　　　　out << ri << endl;
　　　}
　　　out.close();
　　　return 0;
}

//関数rpois
int rpois(double lambda){
　　　double r;
　　　int irpois = 0;
　　　r = (double)rand() / RAND_MAX;
　　　lambda = exp(lambda) * r;

　　　while(lambda > 1){
　　　　　　r = (double)rand() / RAND_MAX;
　　　　　　lambda *= r;
　　　　　　irpois++;
　　　}
　　　return irpois;
}

ちなみに正規乱数はBox-Muller法。
平均0、標準偏差1の正規乱数から、引数にしたがって平均と標準偏差を与える。

//関数rnorm
double rnorm(double mean, double sd){
　　　double r1, r2, irnorm1, irnorm2;
　　　r1 = (double)rand() / RAND_MAX;
　　　r2 = (double)rand() / RAND_MAX;

　　　irnorm1 = sqrt(-2.0 * log(r1)) * sin(2 * PI * r2);
//　　　irnorm2 = sqrt(-2.0 * log(r1)) * cos(2 * PI * r2);

　　　irnorm1 *= sd;
　　　irnorm1 += mean;
　　　return irnorm1;
}

ただし、最初のとこで、

#define PI 3.14159265

としてあります。
これは2つの一様乱数から2つの独立した正規乱数を発生できるのに、関数にすると1つしか戻り値として渡せないので、なんか損した気分。
指数乱数は引数がlambdaのとき、平均1/lambda、標準偏差1/lambdaになる。

//関数rexp
double rexp(double lambda){
　　　double r, irexp;
　　　r = (double)rand() / RAND_MAX;
　　　irexp = -log(r) / lambda;
　　　return irexp;
}