Calculate Pi using Monte Carlo












9












$begingroup$


#include <iostream>
#include <iomanip>


#ifdef USE_OLD_RAND

#include <stdlib.h>
inline double getRandDart() {return rand() * 1.0 / RAND_MAX;}

#else

#include <random>
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0,1);
inline double getRandDart() {return distribution(generator);}

#endif


// Monte Carlo Simulator to estimate the value of PI.
//
// If we have a circle with a radius of 1.
// Then the smallest square that encloses the circle as sides of length 2.
//
// Area of circle pi r^2 = pi
// Area of square 2.r.2.r = 4
//
// Ratio of overlapping area: pi/4
//
// If we throw darts randomly at a dart board (with an even distribution) and always hit the square.
// Then the ratio of darts falling into the circle should be pi/4 of the total number of darts thrown.
//
// pi/4 * countInSquare = countInCircle
//
// pi = 4 . countInCircle / countInSquare
//
// To simplify the maths.
// We will set the center point as 0,0 and only use the top right quadrant of the circle.
// We have 1/4 the size of the square and circle but the same maths still apply.
//
// A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
// A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)
//

int main()
{

long countInSquare = 0;
long countInCircle = 0;

for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
double x = getRandDart();
double y = getRandDart();

double d = (x * x) + (y * y);

countInSquare += 1;
countInCircle += (d >= 1.0) ? 0 : 1;

if (iteration % 10'000'000 == 0) {
std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";
}
}
std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / countInSquare) << "n";

}


Output:



> ./a.out
9990000000 3.14158
10000000000 3.14158


3.14158355









share|improve this question











$endgroup$












  • $begingroup$
    Your comment says that xx+yy==1 is not outside the circle (hence it is inside), but your code considers it outside. Which way is the correct interpretation?
    $endgroup$
    – 1201ProgramAlarm
    yesterday










  • $begingroup$
    @1201ProgramAlarm: Good catch. Stupid comments; why don't they compile so we can check the code matches the comments.
    $endgroup$
    – Martin York
    yesterday












  • $begingroup$
    There are languages where that's kind of true (but also means that you can get syntax errors with your comments)
    $endgroup$
    – Foon
    yesterday










  • $begingroup$
    @Foon Interesting! Do you have any examples?
    $endgroup$
    – Solomon Ucko
    yesterday










  • $begingroup$
    Spark is the first one that comes to mind (en.wikipedia.org/wiki/SPARK_(programming_language)) ... I was thinking Eiffel also did it, but actually, Eiffel's contracts are probably natively supported (Spark is basically a subset of Ada with additional things embedded in Ada compilers for the Spark checker to work with)
    $endgroup$
    – Foon
    yesterday
















9












$begingroup$


#include <iostream>
#include <iomanip>


#ifdef USE_OLD_RAND

#include <stdlib.h>
inline double getRandDart() {return rand() * 1.0 / RAND_MAX;}

#else

#include <random>
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0,1);
inline double getRandDart() {return distribution(generator);}

#endif


// Monte Carlo Simulator to estimate the value of PI.
//
// If we have a circle with a radius of 1.
// Then the smallest square that encloses the circle as sides of length 2.
//
// Area of circle pi r^2 = pi
// Area of square 2.r.2.r = 4
//
// Ratio of overlapping area: pi/4
//
// If we throw darts randomly at a dart board (with an even distribution) and always hit the square.
// Then the ratio of darts falling into the circle should be pi/4 of the total number of darts thrown.
//
// pi/4 * countInSquare = countInCircle
//
// pi = 4 . countInCircle / countInSquare
//
// To simplify the maths.
// We will set the center point as 0,0 and only use the top right quadrant of the circle.
// We have 1/4 the size of the square and circle but the same maths still apply.
//
// A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
// A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)
//

int main()
{

long countInSquare = 0;
long countInCircle = 0;

for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
double x = getRandDart();
double y = getRandDart();

double d = (x * x) + (y * y);

countInSquare += 1;
countInCircle += (d >= 1.0) ? 0 : 1;

if (iteration % 10'000'000 == 0) {
std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";
}
}
std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / countInSquare) << "n";

}


Output:



> ./a.out
9990000000 3.14158
10000000000 3.14158


3.14158355









share|improve this question











$endgroup$












  • $begingroup$
    Your comment says that xx+yy==1 is not outside the circle (hence it is inside), but your code considers it outside. Which way is the correct interpretation?
    $endgroup$
    – 1201ProgramAlarm
    yesterday










  • $begingroup$
    @1201ProgramAlarm: Good catch. Stupid comments; why don't they compile so we can check the code matches the comments.
    $endgroup$
    – Martin York
    yesterday












  • $begingroup$
    There are languages where that's kind of true (but also means that you can get syntax errors with your comments)
    $endgroup$
    – Foon
    yesterday










  • $begingroup$
    @Foon Interesting! Do you have any examples?
    $endgroup$
    – Solomon Ucko
    yesterday










  • $begingroup$
    Spark is the first one that comes to mind (en.wikipedia.org/wiki/SPARK_(programming_language)) ... I was thinking Eiffel also did it, but actually, Eiffel's contracts are probably natively supported (Spark is basically a subset of Ada with additional things embedded in Ada compilers for the Spark checker to work with)
    $endgroup$
    – Foon
    yesterday














9












9








9


1



$begingroup$


#include <iostream>
#include <iomanip>


#ifdef USE_OLD_RAND

#include <stdlib.h>
inline double getRandDart() {return rand() * 1.0 / RAND_MAX;}

#else

#include <random>
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0,1);
inline double getRandDart() {return distribution(generator);}

#endif


// Monte Carlo Simulator to estimate the value of PI.
//
// If we have a circle with a radius of 1.
// Then the smallest square that encloses the circle as sides of length 2.
//
// Area of circle pi r^2 = pi
// Area of square 2.r.2.r = 4
//
// Ratio of overlapping area: pi/4
//
// If we throw darts randomly at a dart board (with an even distribution) and always hit the square.
// Then the ratio of darts falling into the circle should be pi/4 of the total number of darts thrown.
//
// pi/4 * countInSquare = countInCircle
//
// pi = 4 . countInCircle / countInSquare
//
// To simplify the maths.
// We will set the center point as 0,0 and only use the top right quadrant of the circle.
// We have 1/4 the size of the square and circle but the same maths still apply.
//
// A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
// A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)
//

int main()
{

long countInSquare = 0;
long countInCircle = 0;

for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
double x = getRandDart();
double y = getRandDart();

double d = (x * x) + (y * y);

countInSquare += 1;
countInCircle += (d >= 1.0) ? 0 : 1;

if (iteration % 10'000'000 == 0) {
std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";
}
}
std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / countInSquare) << "n";

}


Output:



> ./a.out
9990000000 3.14158
10000000000 3.14158


3.14158355









share|improve this question











$endgroup$




#include <iostream>
#include <iomanip>


#ifdef USE_OLD_RAND

#include <stdlib.h>
inline double getRandDart() {return rand() * 1.0 / RAND_MAX;}

#else

#include <random>
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(0,1);
inline double getRandDart() {return distribution(generator);}

#endif


// Monte Carlo Simulator to estimate the value of PI.
//
// If we have a circle with a radius of 1.
// Then the smallest square that encloses the circle as sides of length 2.
//
// Area of circle pi r^2 = pi
// Area of square 2.r.2.r = 4
//
// Ratio of overlapping area: pi/4
//
// If we throw darts randomly at a dart board (with an even distribution) and always hit the square.
// Then the ratio of darts falling into the circle should be pi/4 of the total number of darts thrown.
//
// pi/4 * countInSquare = countInCircle
//
// pi = 4 . countInCircle / countInSquare
//
// To simplify the maths.
// We will set the center point as 0,0 and only use the top right quadrant of the circle.
// We have 1/4 the size of the square and circle but the same maths still apply.
//
// A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
// A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)
//

int main()
{

long countInSquare = 0;
long countInCircle = 0;

for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
double x = getRandDart();
double y = getRandDart();

double d = (x * x) + (y * y);

countInSquare += 1;
countInCircle += (d >= 1.0) ? 0 : 1;

if (iteration % 10'000'000 == 0) {
std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";
}
}
std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / countInSquare) << "n";

}


Output:



> ./a.out
9990000000 3.14158
10000000000 3.14158


3.14158355






c++ numerical-methods






share|improve this question















share|improve this question













share|improve this question




share|improve this question








edited yesterday









200_success

130k17154419




130k17154419










asked yesterday









Martin YorkMartin York

73.7k488271




73.7k488271












  • $begingroup$
    Your comment says that xx+yy==1 is not outside the circle (hence it is inside), but your code considers it outside. Which way is the correct interpretation?
    $endgroup$
    – 1201ProgramAlarm
    yesterday










  • $begingroup$
    @1201ProgramAlarm: Good catch. Stupid comments; why don't they compile so we can check the code matches the comments.
    $endgroup$
    – Martin York
    yesterday












  • $begingroup$
    There are languages where that's kind of true (but also means that you can get syntax errors with your comments)
    $endgroup$
    – Foon
    yesterday










  • $begingroup$
    @Foon Interesting! Do you have any examples?
    $endgroup$
    – Solomon Ucko
    yesterday










  • $begingroup$
    Spark is the first one that comes to mind (en.wikipedia.org/wiki/SPARK_(programming_language)) ... I was thinking Eiffel also did it, but actually, Eiffel's contracts are probably natively supported (Spark is basically a subset of Ada with additional things embedded in Ada compilers for the Spark checker to work with)
    $endgroup$
    – Foon
    yesterday


















  • $begingroup$
    Your comment says that xx+yy==1 is not outside the circle (hence it is inside), but your code considers it outside. Which way is the correct interpretation?
    $endgroup$
    – 1201ProgramAlarm
    yesterday










  • $begingroup$
    @1201ProgramAlarm: Good catch. Stupid comments; why don't they compile so we can check the code matches the comments.
    $endgroup$
    – Martin York
    yesterday












  • $begingroup$
    There are languages where that's kind of true (but also means that you can get syntax errors with your comments)
    $endgroup$
    – Foon
    yesterday










  • $begingroup$
    @Foon Interesting! Do you have any examples?
    $endgroup$
    – Solomon Ucko
    yesterday










  • $begingroup$
    Spark is the first one that comes to mind (en.wikipedia.org/wiki/SPARK_(programming_language)) ... I was thinking Eiffel also did it, but actually, Eiffel's contracts are probably natively supported (Spark is basically a subset of Ada with additional things embedded in Ada compilers for the Spark checker to work with)
    $endgroup$
    – Foon
    yesterday
















$begingroup$
Your comment says that xx+yy==1 is not outside the circle (hence it is inside), but your code considers it outside. Which way is the correct interpretation?
$endgroup$
– 1201ProgramAlarm
yesterday




$begingroup$
Your comment says that xx+yy==1 is not outside the circle (hence it is inside), but your code considers it outside. Which way is the correct interpretation?
$endgroup$
– 1201ProgramAlarm
yesterday












$begingroup$
@1201ProgramAlarm: Good catch. Stupid comments; why don't they compile so we can check the code matches the comments.
$endgroup$
– Martin York
yesterday






$begingroup$
@1201ProgramAlarm: Good catch. Stupid comments; why don't they compile so we can check the code matches the comments.
$endgroup$
– Martin York
yesterday














$begingroup$
There are languages where that's kind of true (but also means that you can get syntax errors with your comments)
$endgroup$
– Foon
yesterday




$begingroup$
There are languages where that's kind of true (but also means that you can get syntax errors with your comments)
$endgroup$
– Foon
yesterday












$begingroup$
@Foon Interesting! Do you have any examples?
$endgroup$
– Solomon Ucko
yesterday




$begingroup$
@Foon Interesting! Do you have any examples?
$endgroup$
– Solomon Ucko
yesterday












$begingroup$
Spark is the first one that comes to mind (en.wikipedia.org/wiki/SPARK_(programming_language)) ... I was thinking Eiffel also did it, but actually, Eiffel's contracts are probably natively supported (Spark is basically a subset of Ada with additional things embedded in Ada compilers for the Spark checker to work with)
$endgroup$
– Foon
yesterday




$begingroup$
Spark is the first one that comes to mind (en.wikipedia.org/wiki/SPARK_(programming_language)) ... I was thinking Eiffel also did it, but actually, Eiffel's contracts are probably natively supported (Spark is basically a subset of Ada with additional things embedded in Ada compilers for the Spark checker to work with)
$endgroup$
– Foon
yesterday










2 Answers
2






active

oldest

votes


















11












$begingroup$

One could consider at least the following points:




  • Instead of including <stdlib.h>, I'd include <cstdlib>.


  • In getRandDart(), it might in this case be more readable to do static_cast<double>(rand()) / RAND_MAX; instead of multiplying by 1.0.


  • In the for loop, all of x, y and d can be const, so I'd make them const. This has the potential to protect the programmer from unintended mistakes, and can sometimes allow the compiler to optimize better.


  • When you increment by one (in countInSquare += 1;), it makes more sense to use the ++ operator, i.e., to just write ++countInSquare. This is more idiomatic and protects us from unintended mistakes: ++ conveys the meaning of increment (by one), whereas with += we might accidentally write += 2; and that would be perfectly valid (but not what we wanted).


  • Regardless of the above point, notice that during the for-loop, it holds that iteration == countInSquare. So strictly speaking, the variable countInSquare is unnecessary and could be replaced by just iteration when needed.


  • You could consider making the number of iterations and the second operand of the % operand constants to allow for easier modification and perhaps to slightly improve readability.


  • Instead of typing (4.0 * countInCircle / countInSquare) twice, we could make a function that takes the two variables as parameters. This would allow us to save some typing, and again to protect us from unintended mistakes.







share|improve this answer











$endgroup$





















    3












    $begingroup$

    #include <random>
    std::default_random_engine generator;


    If you require portable and reproducible results from your random engine, consider specifying the exact random engine you want instead of relying on an implementation-defined alias.



    // GCC & Clang
    using minstd_rand0 = linear_congruential_engine<uint_fast32_t, 16807, 0, 2147483647>;
    using default_random_engine = minstd_rand0;

    // MSVC
    using default_random_engine = mt19937;




    std::uniform_real_distribution<double>  distribution(0,1);


    For distributions, the standard doesn't require the result to be reproducible and identical across all implementations. If you care for portability, you'll need to write/copy the specific distribution behavior you want.





    inline double getRandDart() {return distribution(generator);}

    // A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
    // A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)


    This is an example of where the comment doesn't match the function. The function getRandDart() doesn't return a dart. It returns a single magnitude instead of a Dart/Vec2d.





        for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {


    Did you intend for ten billion and one iterations (<=)?





        long     countInSquare   = 0;
    for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
    countInSquare += 1;
    if (iteration % 10'000'000 == 0) {
    std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";


    Do you need countInSquare? Both iteration and countInSquare are maintaining the same count.



    When you don't care for the size in which an integer can represent because the value is countable, simply use int. In this case, you went with long as you needed a larger integral type to hold a value that couldn't be represented by an int. Compiling in an environment where long is 32-bit (windows) would obviously be bad. In these cases, use a specific fixed width integer type from <cstdint>. auto also works in determining the correct type (64-bit long on gcc/clang, 64-bit long long on msvc).





            countInCircle   += (d >= 1.0) ? 0 : 1;


    Compilers are nice nowadays and optimized your ternary add.



            countInCircle   += (d < 1.0);




    You can tighten the inner loop by reorganizing the operations. Rather than checking every iteration to report, tile the iterations into groups that calculate between reports.



        constexpr auto totalIterations = 10'000'000'001;
    constexpr auto reportInterval = 10'000'000;
    constexpr auto reports = std::div(totalIterations, reportInterval);

    while (reports.quot--) {
    for (auto iteration = 0; iteration < reportInterval; ++iteration) {
    const auto x = getRandDart();
    const auto y = getRandDart();
    const auto d = x * x + y * y;
    countInCircle += (d < 1.0);
    }
    countInSquare += reportInterval;
    std::cout << countInSquare << " " << (4.0 * countInCircle / countInSquare) << "n";
    }

    while (reports.rem--) {
    const auto x = getRandDart();
    const auto y = getRandDart();
    const auto d = x * x + y * y;
    countInCircle += (d < 1.0);
    }

    std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / totalIterations) << "n";





    share|improve this answer









    $endgroup$













      Your Answer





      StackExchange.ifUsing("editor", function () {
      return StackExchange.using("mathjaxEditing", function () {
      StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
      StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["\$", "\$"]]);
      });
      });
      }, "mathjax-editing");

      StackExchange.ifUsing("editor", function () {
      StackExchange.using("externalEditor", function () {
      StackExchange.using("snippets", function () {
      StackExchange.snippets.init();
      });
      });
      }, "code-snippets");

      StackExchange.ready(function() {
      var channelOptions = {
      tags: "".split(" "),
      id: "196"
      };
      initTagRenderer("".split(" "), "".split(" "), channelOptions);

      StackExchange.using("externalEditor", function() {
      // Have to fire editor after snippets, if snippets enabled
      if (StackExchange.settings.snippets.snippetsEnabled) {
      StackExchange.using("snippets", function() {
      createEditor();
      });
      }
      else {
      createEditor();
      }
      });

      function createEditor() {
      StackExchange.prepareEditor({
      heartbeatType: 'answer',
      autoActivateHeartbeat: false,
      convertImagesToLinks: false,
      noModals: true,
      showLowRepImageUploadWarning: true,
      reputationToPostImages: null,
      bindNavPrevention: true,
      postfix: "",
      imageUploader: {
      brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
      contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
      allowUrls: true
      },
      onDemand: true,
      discardSelector: ".discard-answer"
      ,immediatelyShowMarkdownHelp:true
      });


      }
      });














      draft saved

      draft discarded


















      StackExchange.ready(
      function () {
      StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fcodereview.stackexchange.com%2fquestions%2f215784%2fcalculate-pi-using-monte-carlo%23new-answer', 'question_page');
      }
      );

      Post as a guest















      Required, but never shown

























      2 Answers
      2






      active

      oldest

      votes








      2 Answers
      2






      active

      oldest

      votes









      active

      oldest

      votes






      active

      oldest

      votes









      11












      $begingroup$

      One could consider at least the following points:




      • Instead of including <stdlib.h>, I'd include <cstdlib>.


      • In getRandDart(), it might in this case be more readable to do static_cast<double>(rand()) / RAND_MAX; instead of multiplying by 1.0.


      • In the for loop, all of x, y and d can be const, so I'd make them const. This has the potential to protect the programmer from unintended mistakes, and can sometimes allow the compiler to optimize better.


      • When you increment by one (in countInSquare += 1;), it makes more sense to use the ++ operator, i.e., to just write ++countInSquare. This is more idiomatic and protects us from unintended mistakes: ++ conveys the meaning of increment (by one), whereas with += we might accidentally write += 2; and that would be perfectly valid (but not what we wanted).


      • Regardless of the above point, notice that during the for-loop, it holds that iteration == countInSquare. So strictly speaking, the variable countInSquare is unnecessary and could be replaced by just iteration when needed.


      • You could consider making the number of iterations and the second operand of the % operand constants to allow for easier modification and perhaps to slightly improve readability.


      • Instead of typing (4.0 * countInCircle / countInSquare) twice, we could make a function that takes the two variables as parameters. This would allow us to save some typing, and again to protect us from unintended mistakes.







      share|improve this answer











      $endgroup$


















        11












        $begingroup$

        One could consider at least the following points:




        • Instead of including <stdlib.h>, I'd include <cstdlib>.


        • In getRandDart(), it might in this case be more readable to do static_cast<double>(rand()) / RAND_MAX; instead of multiplying by 1.0.


        • In the for loop, all of x, y and d can be const, so I'd make them const. This has the potential to protect the programmer from unintended mistakes, and can sometimes allow the compiler to optimize better.


        • When you increment by one (in countInSquare += 1;), it makes more sense to use the ++ operator, i.e., to just write ++countInSquare. This is more idiomatic and protects us from unintended mistakes: ++ conveys the meaning of increment (by one), whereas with += we might accidentally write += 2; and that would be perfectly valid (but not what we wanted).


        • Regardless of the above point, notice that during the for-loop, it holds that iteration == countInSquare. So strictly speaking, the variable countInSquare is unnecessary and could be replaced by just iteration when needed.


        • You could consider making the number of iterations and the second operand of the % operand constants to allow for easier modification and perhaps to slightly improve readability.


        • Instead of typing (4.0 * countInCircle / countInSquare) twice, we could make a function that takes the two variables as parameters. This would allow us to save some typing, and again to protect us from unintended mistakes.







        share|improve this answer











        $endgroup$
















          11












          11








          11





          $begingroup$

          One could consider at least the following points:




          • Instead of including <stdlib.h>, I'd include <cstdlib>.


          • In getRandDart(), it might in this case be more readable to do static_cast<double>(rand()) / RAND_MAX; instead of multiplying by 1.0.


          • In the for loop, all of x, y and d can be const, so I'd make them const. This has the potential to protect the programmer from unintended mistakes, and can sometimes allow the compiler to optimize better.


          • When you increment by one (in countInSquare += 1;), it makes more sense to use the ++ operator, i.e., to just write ++countInSquare. This is more idiomatic and protects us from unintended mistakes: ++ conveys the meaning of increment (by one), whereas with += we might accidentally write += 2; and that would be perfectly valid (but not what we wanted).


          • Regardless of the above point, notice that during the for-loop, it holds that iteration == countInSquare. So strictly speaking, the variable countInSquare is unnecessary and could be replaced by just iteration when needed.


          • You could consider making the number of iterations and the second operand of the % operand constants to allow for easier modification and perhaps to slightly improve readability.


          • Instead of typing (4.0 * countInCircle / countInSquare) twice, we could make a function that takes the two variables as parameters. This would allow us to save some typing, and again to protect us from unintended mistakes.







          share|improve this answer











          $endgroup$



          One could consider at least the following points:




          • Instead of including <stdlib.h>, I'd include <cstdlib>.


          • In getRandDart(), it might in this case be more readable to do static_cast<double>(rand()) / RAND_MAX; instead of multiplying by 1.0.


          • In the for loop, all of x, y and d can be const, so I'd make them const. This has the potential to protect the programmer from unintended mistakes, and can sometimes allow the compiler to optimize better.


          • When you increment by one (in countInSquare += 1;), it makes more sense to use the ++ operator, i.e., to just write ++countInSquare. This is more idiomatic and protects us from unintended mistakes: ++ conveys the meaning of increment (by one), whereas with += we might accidentally write += 2; and that would be perfectly valid (but not what we wanted).


          • Regardless of the above point, notice that during the for-loop, it holds that iteration == countInSquare. So strictly speaking, the variable countInSquare is unnecessary and could be replaced by just iteration when needed.


          • You could consider making the number of iterations and the second operand of the % operand constants to allow for easier modification and perhaps to slightly improve readability.


          • Instead of typing (4.0 * countInCircle / countInSquare) twice, we could make a function that takes the two variables as parameters. This would allow us to save some typing, and again to protect us from unintended mistakes.








          share|improve this answer














          share|improve this answer



          share|improve this answer








          edited yesterday

























          answered yesterday









          JuhoJuho

          1,336511




          1,336511

























              3












              $begingroup$

              #include <random>
              std::default_random_engine generator;


              If you require portable and reproducible results from your random engine, consider specifying the exact random engine you want instead of relying on an implementation-defined alias.



              // GCC & Clang
              using minstd_rand0 = linear_congruential_engine<uint_fast32_t, 16807, 0, 2147483647>;
              using default_random_engine = minstd_rand0;

              // MSVC
              using default_random_engine = mt19937;




              std::uniform_real_distribution<double>  distribution(0,1);


              For distributions, the standard doesn't require the result to be reproducible and identical across all implementations. If you care for portability, you'll need to write/copy the specific distribution behavior you want.





              inline double getRandDart() {return distribution(generator);}

              // A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
              // A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)


              This is an example of where the comment doesn't match the function. The function getRandDart() doesn't return a dart. It returns a single magnitude instead of a Dart/Vec2d.





                  for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {


              Did you intend for ten billion and one iterations (<=)?





                  long     countInSquare   = 0;
              for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
              countInSquare += 1;
              if (iteration % 10'000'000 == 0) {
              std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";


              Do you need countInSquare? Both iteration and countInSquare are maintaining the same count.



              When you don't care for the size in which an integer can represent because the value is countable, simply use int. In this case, you went with long as you needed a larger integral type to hold a value that couldn't be represented by an int. Compiling in an environment where long is 32-bit (windows) would obviously be bad. In these cases, use a specific fixed width integer type from <cstdint>. auto also works in determining the correct type (64-bit long on gcc/clang, 64-bit long long on msvc).





                      countInCircle   += (d >= 1.0) ? 0 : 1;


              Compilers are nice nowadays and optimized your ternary add.



                      countInCircle   += (d < 1.0);




              You can tighten the inner loop by reorganizing the operations. Rather than checking every iteration to report, tile the iterations into groups that calculate between reports.



                  constexpr auto totalIterations = 10'000'000'001;
              constexpr auto reportInterval = 10'000'000;
              constexpr auto reports = std::div(totalIterations, reportInterval);

              while (reports.quot--) {
              for (auto iteration = 0; iteration < reportInterval; ++iteration) {
              const auto x = getRandDart();
              const auto y = getRandDart();
              const auto d = x * x + y * y;
              countInCircle += (d < 1.0);
              }
              countInSquare += reportInterval;
              std::cout << countInSquare << " " << (4.0 * countInCircle / countInSquare) << "n";
              }

              while (reports.rem--) {
              const auto x = getRandDart();
              const auto y = getRandDart();
              const auto d = x * x + y * y;
              countInCircle += (d < 1.0);
              }

              std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / totalIterations) << "n";





              share|improve this answer









              $endgroup$


















                3












                $begingroup$

                #include <random>
                std::default_random_engine generator;


                If you require portable and reproducible results from your random engine, consider specifying the exact random engine you want instead of relying on an implementation-defined alias.



                // GCC & Clang
                using minstd_rand0 = linear_congruential_engine<uint_fast32_t, 16807, 0, 2147483647>;
                using default_random_engine = minstd_rand0;

                // MSVC
                using default_random_engine = mt19937;




                std::uniform_real_distribution<double>  distribution(0,1);


                For distributions, the standard doesn't require the result to be reproducible and identical across all implementations. If you care for portability, you'll need to write/copy the specific distribution behavior you want.





                inline double getRandDart() {return distribution(generator);}

                // A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
                // A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)


                This is an example of where the comment doesn't match the function. The function getRandDart() doesn't return a dart. It returns a single magnitude instead of a Dart/Vec2d.





                    for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {


                Did you intend for ten billion and one iterations (<=)?





                    long     countInSquare   = 0;
                for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
                countInSquare += 1;
                if (iteration % 10'000'000 == 0) {
                std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";


                Do you need countInSquare? Both iteration and countInSquare are maintaining the same count.



                When you don't care for the size in which an integer can represent because the value is countable, simply use int. In this case, you went with long as you needed a larger integral type to hold a value that couldn't be represented by an int. Compiling in an environment where long is 32-bit (windows) would obviously be bad. In these cases, use a specific fixed width integer type from <cstdint>. auto also works in determining the correct type (64-bit long on gcc/clang, 64-bit long long on msvc).





                        countInCircle   += (d >= 1.0) ? 0 : 1;


                Compilers are nice nowadays and optimized your ternary add.



                        countInCircle   += (d < 1.0);




                You can tighten the inner loop by reorganizing the operations. Rather than checking every iteration to report, tile the iterations into groups that calculate between reports.



                    constexpr auto totalIterations = 10'000'000'001;
                constexpr auto reportInterval = 10'000'000;
                constexpr auto reports = std::div(totalIterations, reportInterval);

                while (reports.quot--) {
                for (auto iteration = 0; iteration < reportInterval; ++iteration) {
                const auto x = getRandDart();
                const auto y = getRandDart();
                const auto d = x * x + y * y;
                countInCircle += (d < 1.0);
                }
                countInSquare += reportInterval;
                std::cout << countInSquare << " " << (4.0 * countInCircle / countInSquare) << "n";
                }

                while (reports.rem--) {
                const auto x = getRandDart();
                const auto y = getRandDart();
                const auto d = x * x + y * y;
                countInCircle += (d < 1.0);
                }

                std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / totalIterations) << "n";





                share|improve this answer









                $endgroup$
















                  3












                  3








                  3





                  $begingroup$

                  #include <random>
                  std::default_random_engine generator;


                  If you require portable and reproducible results from your random engine, consider specifying the exact random engine you want instead of relying on an implementation-defined alias.



                  // GCC & Clang
                  using minstd_rand0 = linear_congruential_engine<uint_fast32_t, 16807, 0, 2147483647>;
                  using default_random_engine = minstd_rand0;

                  // MSVC
                  using default_random_engine = mt19937;




                  std::uniform_real_distribution<double>  distribution(0,1);


                  For distributions, the standard doesn't require the result to be reproducible and identical across all implementations. If you care for portability, you'll need to write/copy the specific distribution behavior you want.





                  inline double getRandDart() {return distribution(generator);}

                  // A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
                  // A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)


                  This is an example of where the comment doesn't match the function. The function getRandDart() doesn't return a dart. It returns a single magnitude instead of a Dart/Vec2d.





                      for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {


                  Did you intend for ten billion and one iterations (<=)?





                      long     countInSquare   = 0;
                  for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
                  countInSquare += 1;
                  if (iteration % 10'000'000 == 0) {
                  std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";


                  Do you need countInSquare? Both iteration and countInSquare are maintaining the same count.



                  When you don't care for the size in which an integer can represent because the value is countable, simply use int. In this case, you went with long as you needed a larger integral type to hold a value that couldn't be represented by an int. Compiling in an environment where long is 32-bit (windows) would obviously be bad. In these cases, use a specific fixed width integer type from <cstdint>. auto also works in determining the correct type (64-bit long on gcc/clang, 64-bit long long on msvc).





                          countInCircle   += (d >= 1.0) ? 0 : 1;


                  Compilers are nice nowadays and optimized your ternary add.



                          countInCircle   += (d < 1.0);




                  You can tighten the inner loop by reorganizing the operations. Rather than checking every iteration to report, tile the iterations into groups that calculate between reports.



                      constexpr auto totalIterations = 10'000'000'001;
                  constexpr auto reportInterval = 10'000'000;
                  constexpr auto reports = std::div(totalIterations, reportInterval);

                  while (reports.quot--) {
                  for (auto iteration = 0; iteration < reportInterval; ++iteration) {
                  const auto x = getRandDart();
                  const auto y = getRandDart();
                  const auto d = x * x + y * y;
                  countInCircle += (d < 1.0);
                  }
                  countInSquare += reportInterval;
                  std::cout << countInSquare << " " << (4.0 * countInCircle / countInSquare) << "n";
                  }

                  while (reports.rem--) {
                  const auto x = getRandDart();
                  const auto y = getRandDart();
                  const auto d = x * x + y * y;
                  countInCircle += (d < 1.0);
                  }

                  std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / totalIterations) << "n";





                  share|improve this answer









                  $endgroup$



                  #include <random>
                  std::default_random_engine generator;


                  If you require portable and reproducible results from your random engine, consider specifying the exact random engine you want instead of relying on an implementation-defined alias.



                  // GCC & Clang
                  using minstd_rand0 = linear_congruential_engine<uint_fast32_t, 16807, 0, 2147483647>;
                  using default_random_engine = minstd_rand0;

                  // MSVC
                  using default_random_engine = mt19937;




                  std::uniform_real_distribution<double>  distribution(0,1);


                  For distributions, the standard doesn't require the result to be reproducible and identical across all implementations. If you care for portability, you'll need to write/copy the specific distribution behavior you want.





                  inline double getRandDart() {return distribution(generator);}

                  // A dart thrown has a random x/y value in the range 0->1 (top right quadrant).
                  // A dart is outside the circle if x^2 + y^2 > 1 (note 1^2 is 1)


                  This is an example of where the comment doesn't match the function. The function getRandDart() doesn't return a dart. It returns a single magnitude instead of a Dart/Vec2d.





                      for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {


                  Did you intend for ten billion and one iterations (<=)?





                      long     countInSquare   = 0;
                  for(long iteration = 0; iteration <= 10'000'000'000; ++iteration) {
                  countInSquare += 1;
                  if (iteration % 10'000'000 == 0) {
                  std::cout << iteration << " " << (4.0 * countInCircle / countInSquare) << "n";


                  Do you need countInSquare? Both iteration and countInSquare are maintaining the same count.



                  When you don't care for the size in which an integer can represent because the value is countable, simply use int. In this case, you went with long as you needed a larger integral type to hold a value that couldn't be represented by an int. Compiling in an environment where long is 32-bit (windows) would obviously be bad. In these cases, use a specific fixed width integer type from <cstdint>. auto also works in determining the correct type (64-bit long on gcc/clang, 64-bit long long on msvc).





                          countInCircle   += (d >= 1.0) ? 0 : 1;


                  Compilers are nice nowadays and optimized your ternary add.



                          countInCircle   += (d < 1.0);




                  You can tighten the inner loop by reorganizing the operations. Rather than checking every iteration to report, tile the iterations into groups that calculate between reports.



                      constexpr auto totalIterations = 10'000'000'001;
                  constexpr auto reportInterval = 10'000'000;
                  constexpr auto reports = std::div(totalIterations, reportInterval);

                  while (reports.quot--) {
                  for (auto iteration = 0; iteration < reportInterval; ++iteration) {
                  const auto x = getRandDart();
                  const auto y = getRandDart();
                  const auto d = x * x + y * y;
                  countInCircle += (d < 1.0);
                  }
                  countInSquare += reportInterval;
                  std::cout << countInSquare << " " << (4.0 * countInCircle / countInSquare) << "n";
                  }

                  while (reports.rem--) {
                  const auto x = getRandDart();
                  const auto y = getRandDart();
                  const auto d = x * x + y * y;
                  countInCircle += (d < 1.0);
                  }

                  std::cout << "nn" << std::setprecision(9) << (4.0 * countInCircle / totalIterations) << "n";






                  share|improve this answer












                  share|improve this answer



                  share|improve this answer










                  answered 18 hours ago









                  SnowhawkSnowhawk

                  5,51611229




                  5,51611229






























                      draft saved

                      draft discarded




















































                      Thanks for contributing an answer to Code Review Stack Exchange!


                      • Please be sure to answer the question. Provide details and share your research!

                      But avoid



                      • Asking for help, clarification, or responding to other answers.

                      • Making statements based on opinion; back them up with references or personal experience.


                      Use MathJax to format equations. MathJax reference.


                      To learn more, see our tips on writing great answers.




                      draft saved


                      draft discarded














                      StackExchange.ready(
                      function () {
                      StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fcodereview.stackexchange.com%2fquestions%2f215784%2fcalculate-pi-using-monte-carlo%23new-answer', 'question_page');
                      }
                      );

                      Post as a guest















                      Required, but never shown





















































                      Required, but never shown














                      Required, but never shown












                      Required, but never shown







                      Required, but never shown

































                      Required, but never shown














                      Required, but never shown












                      Required, but never shown







                      Required, but never shown







                      Popular posts from this blog

                      數位音樂下載

                      When can things happen in Etherscan, such as the picture below?

                      格利澤436b