【问题标题】:Optimize estimating Pi function C++优化估计 Pi 函数 C++
【发布时间】:2022-01-04 11:18:55
【问题描述】:

我编写了一个程序,它使用蒙特卡罗方法逼近 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;
}

【问题讨论】:

  • 请注意,有更快的收敛方法来估计 Pi。使用蒙特卡罗来估计 pi ​​是一个很好的练习,但仅此而已。无论如何,当然您可以将其用作练习以使其更快。您可以尝试使用较小的切片而不是四分之一圈
  • 非常小的改进.... a) num_point_total = n 所以不需要在循环中计算它。 b) 使用 ++ 而不是 +=1(尽管这可能已被编译器优化)。 c) 将 int 用于 num_point_circle - 如果您在 n= 1000 时看到这些变化带来的巨大变化,您会感到惊讶 - 您的机器上的“某个时间”是什么?
  • 您希望通过使该方法更快获得什么?无论如何,您将永远达到十位数的准确度,并且您已经知道该方法有效。那又怎样?
  • 尝试使用 OpenMP 进行并行计算。如果这样做,请考虑使用与rand() 不同的东西,它不是线程安全的。您也可以尝试用不同的东西替换rand(),因为它很慢。
  • @MichałWalenciak:正如我在回答中所说,让算法更快实际上是没有希望的。并行化,比如 64 个核心最多只能产生 8 倍的加速,甚至不足以获得一个新的精确数字。

标签: c++ function loops c++builder


【解决方案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_circlennum_point_total。为他们使用intlonglong long


蒙特卡洛算法也是令人尴尬的并行。因此,加速它的一个想法是使用多个线程运行算法,每个线程处理n 的一小部分,然后最后将圆内的点数相加。


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


最后,还有更有效的计算 Pi 的方法。 使用蒙特卡洛方法,您可以进行数百万次迭代,并且只能获得几位数的准确度。使用更好的算法,可以计算数千、数百万或数十亿位数。

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

【讨论】:

  • 仅使用简单的 Machin 公式 16 arctan(1/5)-4 arctan(1/239),您已经达到了大约 30 个项的双精度精度。
【解决方案2】:

首先,您应该修改代码以考虑在一定次数的迭代后估计会发生多少变化,并在达到令人满意的准确度时停止。

对于您的代码,在迭代次数固定的情况下,编译器无法比您做得更好的优化。您可以改进的一件小事是在每次迭代中放弃计算平方根。你不需要它,因为如果x &lt;= 1sqrt(x) &lt;= 1 为真。

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

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

【讨论】:

    【解决方案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 位数字甚至更慢。所以整个蒙特卡洛算法确实不值得深入挖掘,因为它只关心速度,同时保持某种可靠性。

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

    【讨论】:

      【解决方案4】:

      随着误差随着样本数平方根的倒数而减小,要获得一位数,您需要一百倍以上的样本。没有微优化将有很大帮助。只要考虑这种方法对有效计算 π 没有用处。

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

      【讨论】:

        【解决方案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。

        【讨论】:

          【解决方案6】:

          您的代码在性能方面非常糟糕:

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

          在 int 和 double 之间进行转换,也使用整数除法...为什么不使用 Random() 来代替?还有这个:

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

          是个坏主意,因为 ndoubleiint 所以要么在开始时将 n 存储到 int 变量,要么将标题更改为 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);
          

          为什么不做连续计算? ...所以你需要实现的是:

          1. 应用启动时加载状态

          2. 计时器或 OnIdle 事件中的计算

            所以在每次迭代中/甚至你都会对 Pi 进行 N 次迭代(添加到一些全局总和和计数)

          3. 在应用退出时保存状态

          注意 Monte Carlo Pi 计算收敛速度非常慢,一旦总和变得太大,您将遇到浮点精度问题

          这是我很久以前做的连续蒙特卡罗的小例子......

          表格 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
          

          所以你应该添加状态的保存/加载...

          如前所述,有很多更好的获取 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 值迭代并停止在某些键击(如转义)上。

          【讨论】:

            猜你喜欢
            • 2018-07-21
            • 1970-01-01
            • 1970-01-01
            • 2018-09-12
            • 1970-01-01
            • 2023-02-08
            • 2016-11-05
            • 2016-02-21
            • 1970-01-01
            相关资源
            最近更新 更多