Ad

the solution is incredibly slow :: try .. 999999999 that's far from long long limits, butt still it takes forever for our soulution to calculate result. monte carlo algorithm is more fun than just straight bruteforce
we can replace big chunk of multiplications with comparison with constants,cutting square with R-size diagonal.. and go crazy a bit further.

also x++ is a "C" idiom, and it works the same way as ++x for integers in C++, so why bother with "more effective" way from books ++x.. we already can traverse collections without ++ing iterators.. and pow is overkill for one multiplication.. you use it for allocating space in stack and nothing more.. we can use inline one for same purpose.

obviously tests were designed for certain number of dots, but we can greatly increase it thus getting better results, so my approach slaps alot more dots in uniformal grid and gets to pi alot closer on the same n and alot faster :) say random generator slapped dots in uniformal grid fashion all over the 1sq quadrant .. and we use this possibility to speed up our algorithm

Code
Diff
  • #include <random>
    
    std::uniform_real_distribution<long double> dis(0.,1.); // Default is double
    std::random_device rd;
    std::mt19937 gen(rd());
    
    inline long double ppow(const long double& x) {
      return x * x;
    }
    
    long double pi_estimate_yours(long long n, int seed=0) {
        gen.seed(seed);
        long long inside = 0;
        for (auto i=0; i!=n; i++) {
          if (ppow(dis(gen))+ppow(dis(gen)) < 1.){
            inside++;
          }
        }
        return (4.*inside)/n; 
    }
    
    long double divisor(long double num) {
        return sqrt((2.*sqrt(2.) * static_cast<double>(num)) / (2. - sqrt(2.)));
    }
    
    //
    // i 
    //
    long double pi_estimate(long double n, int seed = 0) {
      
        gen.seed(seed); // warnings
        long double div = divisor(n);
    
        long double R = 1.L;
        long double Rsq = R * R;
        long double A = sqrt(0.5L);
        long double Rbit = 1.L / div;
        std::cout << Rbit << std::endl;
        long double Rnum = div;
        long double Anum = A / Rbit;
    
        long double inside_accum = Anum * Anum + (Rnum - Anum) * Anum; // regular grid
        
    
    
        long double yisq = 0;
        long double subaccum = 0;
        for (long double yi = A; yi < R; yi += Rbit) {
            yisq = yi * yi;
            for (long double xi = A * (R - yi)/(R - A); xi < A; xi += Rbit) {
                if (yisq + xi * xi < Rsq) {
                    subaccum++;
                } else break;
            }
        }
        
        inside_accum += 2 * subaccum;
        
        long double pee = inside_accum * 4. / (Rnum * Rnum);
    
    
        return (pee);
    }
    • #include <random>
    • std::uniform_real_distribution<long double> dis(0.,1.); // Default is double
    • std::random_device rd;
    • std::mt19937 gen(rd());
    • long double pi_estimate(long long n, int seed=0) {
    • inline long double ppow(const long double& x) {
    • return x * x;
    • }
    • long double pi_estimate_yours(long long n, int seed=0) {
    • gen.seed(seed);
    • long long inside=0;
    • for (auto i=0; i!=n; ++i) if (std::pow(dis(gen),2)+std::pow(dis(gen),2) < 1.) ++inside;
    • long long inside = 0;
    • for (auto i=0; i!=n; i++) {
    • if (ppow(dis(gen))+ppow(dis(gen)) < 1.){
    • inside++;
    • }
    • }
    • return (4.*inside)/n;
    • }
    • long double divisor(long double num) {
    • return sqrt((2.*sqrt(2.) * static_cast<double>(num)) / (2. - sqrt(2.)));
    • }
    • //
    • // i
    • //
    • long double pi_estimate(long double n, int seed = 0) {
    • gen.seed(seed); // warnings
    • long double div = divisor(n);
    • long double R = 1.L;
    • long double Rsq = R * R;
    • long double A = sqrt(0.5L);
    • long double Rbit = 1.L / div;
    • std::cout << Rbit << std::endl;
    • long double Rnum = div;
    • long double Anum = A / Rbit;
    • long double inside_accum = Anum * Anum + (Rnum - Anum) * Anum; // regular grid
    • long double yisq = 0;
    • long double subaccum = 0;
    • for (long double yi = A; yi < R; yi += Rbit) {
    • yisq = yi * yi;
    • for (long double xi = A * (R - yi)/(R - A); xi < A; xi += Rbit) {
    • if (yisq + xi * xi < Rsq) {
    • subaccum++;
    • } else break;
    • }
    • }
    • inside_accum += 2 * subaccum;
    • long double pee = inside_accum * 4. / (Rnum * Rnum);
    • return (pee);
    • }