我寫了一個程式，它使用蒙特卡洛方法逼近 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;
}
uj5u.com熱心網友回復：
一個明顯的（小）加速是擺脫平方根運算。
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為他們。
蒙特卡洛演算法也是令人尷尬的并行。所以加速它的一個想法是用多個執行緒運行演算法，每個執行緒處理 的一小部分n，然后最后總結圓內的點數。
此外，您可以嘗試使用 SIMD 指令提高速度。
最后，還有更有效的計算 Pi 的方法。使用蒙特卡羅方法，您可以進行數百萬次迭代，并且只能獲得幾位數的準確度。有了更好的演算法，就可以計算數千、數百萬或數十億的數字。
例如，您可以在網站https://www.i4cy.com/pi/上閱讀
uj5u.com熱心網友回復：
首先，您應該修改代碼以考慮在經過一定次數的迭代后估計會發生多少變化，并在達到令人滿意的準確度時停止。
對于您的代碼，在迭代次數固定的情況下，沒有太多需要優化的編譯器不能比您做得更好。您可以改進的一件小事是在每次迭代中都不再計算平方根。你不需要它，因為它sqrt(x) <= 1是真的 if x <= 1。
粗略地說，一旦您的x,y點均勻分布在四分之一圓中，您就可以從蒙特卡羅方法中得到很好的估計。考慮到這一點，有一種更簡單的方法可以使點均勻分布而沒有隨機性：使用均勻網格并計算圓內有多少點以及圓外有多少點。
我希望這比蒙特卡羅收斂得更好（當減小網格間距時）。但是，您也可以嘗試使用它來加速您的 Monte Carlo 代碼num_point_circle，方法num_point_total是從均勻網格上的點的計數開始，然后通過繼續隨機分布的點來增加計數器。
uj5u.com熱心網友回復：
隨著誤差隨著樣本數平方根的倒數而減少，要獲得一位數，您需要增加 100 倍的樣本。沒有任何微觀優化會顯著幫助。只是考慮到這種方法對于有效計算 π 沒有用處。
如果您堅持：https://en.wikipedia.org/wiki/Variance_reduction。
uj5u.com熱心網友回復：
如 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。
uj5u.com熱心網友回復：
我怎樣才能使這個功能更快地作業？不是更精確！
當 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 位數則更慢。所以整個蒙特卡羅演算法確實不值得深入挖掘，因為它只是關于速度同時保持某種可靠性。
從字面上看，沒有什么比使用固定數字作為 Pi 的 #define 或靜態更快的了，但這不是你的問題。
uj5u.com熱心網友回復：
從性能方面來看，您的代碼非常糟糕，因為：
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。順便說一句，num_point_total當你已經得到時，你為什么要計算n？還：
num_point_circle = (sqrt(x*x y*y) <= 1);
為什么是個壞主意sqrt？你知道1^2 = 1所以你可以簡單地做：
num_point_circle = (x*x y*y <= 1);
為什么不做連續計算？......所以你需要實作的是：
應用程式啟動時的狀態負載
在計時器或 OnIdle 事件中計算
所以在每次迭代中/甚至你會做 N 次 Pi 迭代（添加到一些全域總和和計數）
在應用程式退出時保存狀態
請注意 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
所以你應該添加狀態的保存/加載......
如前所述，有很多更好的方法可以像BBP一樣獲得 Pi
上面的代碼也使用了我的時間測量 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 值一次？迭代并停止一些關鍵的點擊，比如逃生。
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/368828.html
標籤：C 功能 循环 时间：2019-05-06 标签：c builder
上一篇：適當遞增