您的代码在性能方面非常糟糕:
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 变量,要么将标题更改为 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 事件中的计算
所以在每次迭代中/甚至你都会对 Pi 进行 N 次迭代(添加到一些全局总和和计数)
-
在应用退出时保存状态
注意 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 值迭代并停止在某些键击(如转义)上。