SCIP in C++11 ― 3.1.2節

モンテカルロシミュレーションと問題3.5

まあこの問題でrandom-in-rangeを自前のListを使った実装にする意味は無い、

というかパフォーマンスがとても悪いし、

この辺のコードは、この本の後のほうで確か使わなかったと思うが、

一応set!の練習を兼ねて、Listを使った実装で。

問題3.6はメッセージパッシングの練習ではあるのだが、

このC++11のMersenne-twister乱数の実装ではあまり意味が無いのでパス。

----

const List randomUpdate(const List& dummy)

{

    std::random_device rd;

    std::mt19937 x(rd());

    std::uniform_int_distribution<> dis(1,1000);

    return(makeLeaf(dis(x)));

}

const function<List(void)> randSICP(void)

{

    const List x(makeLeaf(0));

    return([x](void){

         setInsert(x,randomUpdate(x));

         return(x);

     });

}

const List randomInRange(const List& low, const List& high)

{

    std::random_device rd;

    std::mt19937 gen(rd());

    std::uniform_real_distribution<> dis

     (value<double>(low), value<double>(high));

    return(makeLeaf(dis(gen)));

}

//---------abstraction barrier---------

const List gcd(const List& n1, const List& n2)

{

    if(n2==makeLeaf(0)){return(n1);}

    return(gcd(n2,n1%n2));

}

const double MonteCarlo

(const int trials, const function<bool(void)> experiment)

{

    function<double(int,int)> iter;

    iter=[&iter,trials,&experiment]

     (const int trialsRemaining,

      const int trialsPassed)

     {

         if(trialsRemaining==0){

             return(static_cast<double>(trialsPassed)

                    /static_cast<double>(trials));

         }

         if(experiment()){

             return(iter(trialsRemaining-1,trialsPassed+1));

         }

         return(iter(trialsRemaining-1,trialsPassed));

     };

    return(iter(trials,0));

}

const function<bool(void)> cesaroTest

=[](void){

    return(gcd(randSICP()(),randSICP()())==makeLeaf(1));

};

const double estimatePi(const int trials)

{

    return(sqrt(6.0

             /MonteCarlo(trials,cesaroTest)));

}

const function<bool(double,double,double,double)>

integralTest

(const function<double(double)>& f)

{

    return([f](const double x1,const double x2,

            const double y1,const double y2){

            const double x

                (value<double>

                 (randomInRange(makeLeaf(x1),makeLeaf(x2))));

            const double y

                (value<double>

                 (randomInRange(makeLeaf(y1),makeLeaf(y2))));

            return(f(x)>y);});

}

const double estimateIntegral

(const function<bool(double,double,double,double)>& P,

 const double x1,const double x2, const double y1,const double y2,

 const int trials)

{

    const function<bool(void)> predicateIntegration

     =[P,x1,x2,y1,y2](void){

     return(P(x1,x2,y1,y2));

    };

    return(MonteCarlo(trials,predicateIntegration));

}

int main(int argc, char** argv)

{

    cout<<"(estimate-pi 10000) = "<<estimatePi(10000)<<endl;

   

    cout<<endl<<"Excersize 3.5:"<<endl;

    const function<double(double)> circleFunc

     =[](const double x){return(std::sqrt(1.0-square(x)));};

    cout<<4.0*estimateIntegral

     (integralTest(circleFunc),

      0.0,1.0,0.0,1.0,10000)<<endl;

    return(0);

}

----

出力

----

(estimate-pi 10000) = 3.12729

Excersize 3.5:

3.1256