伪蒙特卡洛(Quasi-Monte Carlo, QMC)随机

#include <cmath>
#include <cstdio>
#include <vector>
#include <cassert>
#include <omp.h>
const int UP=100;
bool sieve[UP+100];
int primes[UP],top=0;
void init()
{
  for (int i=2;i<=UP;++i)
    if (!sieve[i])
    {
      primes[top++]=i;
      for (int j=i;j<=UP/i;++j)
        sieve[i*j]=true;
    }
}
std::vector<double> halton(long long i,const int &dim)
{
  assert(dim<=top);
  std::vector<double> prime_inv(dim,0),r(dim,0);
  std::vector<long long> t(dim,i);
  for (int j=0;j<dim;++j)
    prime_inv[j]=1.0/primes[j];
  auto f=[](const std::vector<long long> &t)->long long {
    long long ret=0;
    for (const auto &e:t)
      ret+=e;
    return ret;
  };
  for (;f(t)>0;)
    for (int j=0;j<dim;++j)
    {
      long long d=t[j]%primes[j];
      r[j]+=d*prime_inv[j];
      prime_inv[j]/=primes[j];
      t[j]/=primes[j];
    }
  return r;
}
double experiment(long long idx)
{
  std::vector<double> li=halton(idx,6);
  double area1=fabs((li.at(0)-li.at(2))*(li.at(1)-li.at(3)));
  double area2=fabs((li.at(0)-li.at(4))*(li.at(1)-li.at(5)));
  double area3=fabs((li.at(2)-li.at(4))*(li.at(3)-li.at(5)));
  double w=area1+area2+area3-std::max(std::max(area1,area2),area3)-std::min(std::min(area1,area2),area3);
  return w;
 }
const int BATCH=100000;
const int THREADS=40;
int main()
{
  init();
  double total=0;
  for (long long trial=0;;)
  {
    std::vector<double> li(THREADS,0);
    omp_set_dynamic(0);
    omp_set_num_threads(THREADS);
    #pragma omp parallel for
    for (long long thread=0;thread<THREADS;++thread)
    {
      for (long long i=0;i<BATCH;++i)
        li.at(thread)+=experiment(trial+thread*BATCH+i);
    }
    for (const auto &d:li)
      total+=d;
    trial+=THREADS*BATCH;
    printf("%lld: %.10f\n",trial,total/trial),fflush(stdout);
  }
  return 0;
}