优化估计 Pi 函数 C++

Question

我写了一个程序，它使用蒙特卡洛方法近似 Pi。它工作正常，但我想知道我是否可以让它工作得更好更快，因为当从那个点插入类似 ~n = 100000000 和更大的东西时，需要一些时间来进行计算和打印结果。

我曾想过如何通过对 n 结果取中值来更好地近似它，但考虑到我的算法对大数字有多慢，我决定不这样做。

基本上，问题是：我怎样才能让这个功能运行得更快？

这是我目前得到的代码：

double estimate_pi(double n)
{
    int i;
    double x,y, distance;
    srand( (unsigned)time(0));
    double num_point_circle = 0;
    double num_point_total = 0;
    double final;

    for (i=0; i<n; i++)
    {
        x = (double)rand()/RAND_MAX;
        y = (double)rand()/RAND_MAX;
        distance = sqrt(x*x + y*y);
        if (distance <= 1)
        {
            num_point_circle+=1;
        }
        num_point_total+=1;
    }
    final = ((4 * num_point_circle) / num_point_total);
    return final;
}

Answer 1

一个明显的（小的）加速是去掉平方根运算。 sqrt(x*x + y*y)正好小于1，当x*x + y*y也小于1。所以你可以写：

double distance2 = x*x + y*y;
if (distance2 <= 1) {
    ...
}

这可能会有所收获，因为计算平方根是一项昂贵的操作，并且需要更多的时间（比加法慢 4-20 倍，具体取决于 CPU、体系结构、优化级别、. ..).

您还应该对计数的变量使用整数值，例如 num_point_circle、n 和 num_point_total。对它们使用 int、long 或 long long。

Monte Carlo算法也是并行的尴尬。因此，加速它的一个想法是运行具有多个线程的算法，每个线程处理 n 的一小部分，然后在最后总结圆内的点数。

此外，您可以尝试使用 SIMD 指令提高速度。

最后，还有更有效的计算 Pi 的方法。使用 Monte Carlo 方法，您可以进行数百万次迭代，并且只能获得几位数的精度。使用更好的算法可以计算数千、数百万或数十亿位数字。

例如你可以在网站上阅读 https://www.i4cy.com/pi/

Answer 2

首先，你应该修改代码，考虑迭代一定次数后估计值变化多少，达到满意精度时停止。

对于您的代码，迭代次数固定，编译器无法比您更好地优化。您可以改进的一件小事是在每次迭代中放弃计算平方根。你不需要它，因为如果 x <= 1.

sqrt(x) <= 1 为真

非常粗略地说，一旦您的 x,y 点均匀分布在四分之一圆中，您就可以从 Monte Carlo 方法中得到一个很好的估计。考虑到这一点，有一种更简单的方法可以使点均匀分布而不随机：使用均匀网格并计算圆内有多少点，圆外有多少点。

我希望它比 Monte Carlo 收敛得更好（当减小网格间距时）。但是，您也可以尝试使用它来加速您的 Monte Carlo 代码，方法是从从统一网格上的点获得的 num_point_circle 和 num_point_total 的计数开始，然后通过继续递增计数器随机分布点。

Answer 3

如何让这个功能运行得更快？不是更精确！

当Monte Carlo无论如何都是一个估计，最后的乘法是2,4,8的倍数等等你也可以做位运算。任何 if 语句都会使它变慢，因此试图摆脱它。但是，当我们增加 1 或什么都不增加 (0) 时，您可以摆脱该语句并将其简化为简单的数学运算，这应该会更快。当i在循环前初始化，我们在循环内递增，也可以是while循环

#include <time.h>
double estimate_alt_pi(double n){
    uint64_t num_point_total = 0;
    double x, y;
    srand( (unsigned)time(0));
    uint64_t num_point_circle = 0;
    while (num_point_total<n) {
        x = (double)rand()/RAND_MAX;
        y = (double)rand()/RAND_MAX;
        num_point_circle += (sqrt(x*x + y*y) <= 1);
        num_point_total++; // +1.0 we don't need 'i'
    }
    return ((num_point_circle << 2) / (double)num_point_total); // x<<2 == x * 4
}

在 mac 上用 Xcode 进行基准测试，看起来像。

extern uint64_t dispatch_benchmark(size_t count, void (^block)(void));

int main(int argc, const char * argv[]) {
    size_t count = 1000000;
    double input = 1222.52764523423;
    __block double resultA;
    __block double resultB;
    uint64_t t = dispatch_benchmark(count,^{
        resultA = estimate_pi(input);
    });
    uint64_t t2 = dispatch_benchmark(count,^{
        resultB = estimate_alt_pi(input);
    });
    fprintf(stderr,"estimate_____pi=%f time used=%llu\n",resultA, t);
    fprintf(stderr,"estimate_alt_pi=%f time used=%llu\n",resultB, t2);
    return 0;
}

~快 1.35 倍，或花费原来时间的~73%
但当给定的数字较低时，差异显着较小。而且整个算法仅适用于最大 7 位的输入，这是因为使用的数据类型。低于 2 位数甚至更慢。所以整个 Monte Carlo 算法确实不值得深入挖掘，因为它只是关于速度，同时保持某种可靠性。

从字面上看，没有什么比使用带有固定数字的 #define 或 static 作为 Pi 更快的了，但这不是你的问题。

Answer 4

由于误差随着样本数平方根的倒数而减小，要获得一位数，您需要一百倍的样本。没有微观优化会有很大帮助。只是认为这种方法对有效计算 π 没有用。

如果你坚持：https://en.wikipedia.org/wiki/Variance_reduction.

Answer 5

第一级改进是使用蛮力算法在规则网格上采样，如 463035818_is_not_a_number 所建议的那样。

下一级别的改进是为每个 'x'“绘制”一个半圆，计算必须低于 y = sqrt(radius*radius - x*x) 的点数。这将复杂度从 O(N^2) 降低到 O(N)。

在网格大小 == 半径为 10、100、1000、10000 等的情况下，每一步应该有大约一位数字。

rand() 函数的一个问题是，数字很快就会开始重复——使用规则网格和这个 O(N) 转换问题，我们甚至可以模拟大小为 2 的网格^32 在相当合理的时间内。

利用 Bresenham 的圆绘图算法的思想，我们甚至可以快速评估是否应选择候选 (x+1, y_previous) 或 (x+1, y_previous-1)，因为其中只有一个在第一个八分圆的圆圈内。其次，需要一些其他技巧来避免 sqrt。

Answer 6

从性能方面看，您的代码非常糟糕：

x = (double)rand()/RAND_MAX;
y = (double)rand()/RAND_MAX;

是在 int 和 double 之间转换，也使用整数除法...为什么不使用 Random() 呢？还有这个：

for (i=0; i<n; i++)

是个坏主意，因为 n 是 double 而 i 是 int 所以要么在开始时将 n 存储到 int 变量，要么改变 header 到 int n。顺便说一句，当你已经得到 n 时，你为什么还要计算 num_point_total ？还有：

num_point_circle += (sqrt(x*x + y*y) <= 1);

是个坏主意，为什么 sqrt？你知道 1^2 = 1 所以你可以简单地做：

num_point_circle += (x*x + y*y <= 1);

为什么不做连续计算呢？ ...所以你需要实现的是：

应用程序启动时的状态加载
在计时器或 OnIdle 事件中计算

所以在每个 iteration/even 中你将对 Pi 进行 N 次迭代（添加到一些全局总和和计数）
应用程序退出时保存状态

当心Monte Carlo Pi 计算收敛非常缓慢，一旦总和变得太大，你就会遇到浮点精度问题

这是我很久以前做的小例子，连续Monte Carlo...

表格 cpp:

//$$---- Form CPP ----
//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "Unit1.h"
#include "performance.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
int i=0,n=0,n0=0;
//---------------------------------------------------------------------------
void __fastcall TForm1::Idleloop(TObject *Sender, bool &Done)
    {
    int j;
    double x,y;
    for (j=0;j<10000;j++,n++)
        {
        x=Random(); x*=x;
        y=Random(); y*=y;
        if (x+y<=1.0) i++;
        }
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    tbeg();
    Randomize();
    Application->OnIdle = Idleloop;
    }
//-------------------------------------------------------------------------
void __fastcall TForm1::Timer1Timer(TObject *Sender)
    {
    double dt;
    AnsiString txt;
    txt ="ref = 3.1415926535897932384626433832795\r\n";
    if (n) txt+=AnsiString().sprintf("Pi  = %.20lf\r\n",4.0*double(i)/double(n));
    txt+=AnsiString().sprintf("i/n = %i / %i\r\n",i,n);
    dt=tend();
    if (dt>1e-100) txt+=AnsiString().sprintf("IPS = %8.0lf\r\n",double(n-n0)*1000.0/dt);
    tbeg(); n0=n;
    mm_log->Text=txt;
    }
//---------------------------------------------------------------------------

表格 h:

//$$---- Form HDR ----
//---------------------------------------------------------------------------

#ifndef Unit1H
#define Unit1H
//---------------------------------------------------------------------------
#include <Classes.hpp>
#include <Controls.hpp>
#include <StdCtrls.hpp>
#include <Forms.hpp>
#include <ExtCtrls.hpp>
//---------------------------------------------------------------------------
class TForm1 : public TForm
{
__published:    // IDE-managed Components
    TMemo *mm_log;
    TTimer *Timer1;
    void __fastcall Timer1Timer(TObject *Sender);
private:    // User declarations
public:     // User declarations
    __fastcall TForm1(TComponent* Owner);
void __fastcall TForm1::Idleloop(TObject *Sender, bool &Done);
};
//---------------------------------------------------------------------------
extern PACKAGE TForm1 *Form1;
//---------------------------------------------------------------------------
#endif

表格 dfm:

object Form1: TForm1
  Left = 0
  Top = 0
  Caption = 'Project Euler'
  ClientHeight = 362
  ClientWidth = 619
  Color = clBtnFace
  Font.Charset = OEM_CHARSET
  Font.Color = clWindowText
  Font.Height = 14
  Font.Name = 'System'
  Font.Pitch = fpFixed
  Font.Style = [fsBold]
  OldCreateOrder = False
  PixelsPerInch = 96
  TextHeight = 14
  object mm_log: TMemo
    Left = 0
    Top = 0
    Width = 619
    Height = 362
    Align = alClient
    ScrollBars = ssBoth
    TabOrder = 0
  end
  object Timer1: TTimer
    Interval = 100
    OnTimer = Timer1Timer
    Left = 12
    Top = 10
  end
end

所以你应该添加 save/load 状态 ...

如前所述，有很多更好的方法来获取 Pi，例如 BBP

上面的代码也使用了我的时间测量 heder，所以这里是：

//---------------------------------------------------------------------------
//--- Performance counter time measurement: 2.01 ----------------------------
//---------------------------------------------------------------------------
#ifndef _performance_h
#define _performance_h
//---------------------------------------------------------------------------
const int   performance_max=64;                 // push urovni
double      performance_Tms=-1.0,               // perioda citaca [ms]
            performance_tms=0.0,                // zmerany cas po tend [ms]
            performance_t0[performance_max];    // zmerane start casy [ms]
int         performance_ix=-1;                  // index aktualneho casu
//---------------------------------------------------------------------------
void tbeg(double *t0=NULL)  // mesure start time
    {
    double t;
    LARGE_INTEGER i;
    if (performance_Tms<=0.0)
        {
        for (int j=0;j<performance_max;j++) performance_t0[j]=0.0;
        QueryPerformanceFrequency(&i); performance_Tms=1000.0/double(i.QuadPart);
        }
    QueryPerformanceCounter(&i); t=double(i.QuadPart); t*=performance_Tms;
    if (t0) { t0[0]=t; return; }
    performance_ix++;
    if ((performance_ix>=0)&&(performance_ix<performance_max)) performance_t0[performance_ix]=t;
    }
//---------------------------------------------------------------------------
void tpause(double *t0=NULL)    // stop counting time between tbeg()..tend() calls
    {
    double t;
    LARGE_INTEGER i;
    QueryPerformanceCounter(&i); t=double(i.QuadPart); t*=performance_Tms;
    if (t0) { t0[0]=t-t0[0]; return; }
    if ((performance_ix>=0)&&(performance_ix<performance_max)) performance_t0[performance_ix]=t-performance_t0[performance_ix];
    }
//---------------------------------------------------------------------------
void tresume(double *t0=NULL)   // resume counting time between tbeg()..tend() calls
    {
    double t;
    LARGE_INTEGER i;
    QueryPerformanceCounter(&i); t=double(i.QuadPart); t*=performance_Tms;
    if (t0) { t0[0]=t-t0[0]; return; }
    if ((performance_ix>=0)&&(performance_ix<performance_max)) performance_t0[performance_ix]=t-performance_t0[performance_ix];
    }
//---------------------------------------------------------------------------
double tend(double *t0=NULL)    // return duration [ms] between matching tbeg()..tend() calls
    {
    double t;
    LARGE_INTEGER i;
    QueryPerformanceCounter(&i); t=double(i.QuadPart); t*=performance_Tms;
    if (t0) { t-=t0[0]; performance_tms=t; return t; }
    if ((performance_ix>=0)&&(performance_ix<performance_max)) t-=performance_t0[performance_ix]; else t=0.0;
    performance_ix--;
    performance_tms=t;
    return t;
    }
//---------------------------------------------------------------------------
double tper(double *t0=NULL)    // return duration [ms] between tper() calls
    {
    double t,tt;
    LARGE_INTEGER i;
    if (performance_Tms<=0.0)
        {
        for (int j=0;j<performance_max;j++) performance_t0[j]=0.0;
        QueryPerformanceFrequency(&i); performance_Tms=1000.0/double(i.QuadPart);
        }
    QueryPerformanceCounter(&i); t=double(i.QuadPart); t*=performance_Tms;
    if (t0) { tt=t-t0[0]; t0[0]=t; performance_tms=tt; return tt; }
    performance_ix++;
    if ((performance_ix>=0)&&(performance_ix<performance_max))
        {
        tt=t-performance_t0[performance_ix];
        performance_t0[performance_ix]=t;
        }
    else { t=0.0; tt=0.0; };
    performance_ix--;
    performance_tms=tt;
    return tt;
    }
//---------------------------------------------------------------------------
AnsiString tstr()
    {
    AnsiString s;
    s=s.sprintf("%8.3lf",performance_tms); while (s.Length()<8) s=" "+s; s="["+s+" ms]";
    return s;
    }
//---------------------------------------------------------------------------
AnsiString tstr(int N)
    {
    AnsiString s;
    s=s.sprintf("%8.3lf",performance_tms/double(N)); while (s.Length()<8) s=" "+s; s="["+s+" ms]";
    return s;
    }
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

几秒钟后结果如下：

其中 IPS 是每秒迭代次数，i,n 是保存实际计算状态的全局变量，Pi 是实际近似值，ref 是用于比较的参考 Pi 值。当我在 OnIdle 中计算并在 OnTimer 中显示时，不需要任何锁，因为它都是单线程。为了更快的速度，您可以启动更多的计算线程，但是您需要添加多线程锁定机制并将全局状态设置为 volatile.

如果你正在做控制台应用程序（无表单），你仍然可以这样做，只需将你的代码转换为无限循环，重新计算并输出 Pi 值一次 ???迭代并在某些关键命中时停止，例如转义。

优化估计 Pi 函数 C++

Optimize estimating Pi function C++

c++

function

loops

c++builder