首先我们知道
a
r
c
t
a
n
(
x
)
arctan(x)
arctan(x)的导数
f
′
(
x
)
=
1
1
+
t
2
f'(x)=\frac{1}{1+t^2}
f′(x)=1+t21,所以有:
a
r
c
t
a
n
(
x
)
=
∫
0
x
1
1
+
t
2
d
t
(1)
arctan(x)=\int_0^x\frac{1}{1+t^2}dt\tag{1}
arctan(x)=∫0x1+t21dt(1) 取
x
=
1
x=1
x=1得到
π
4
=
a
r
c
t
a
n
(
1
)
=
∫
0
1
1
1
+
t
2
d
t
(2)
\frac{π}{4}=arctan(1)=\int_0^1\frac{1}{1+t^2}dt\tag{2}
4π=arctan(1)=∫011+t21dt(2) 离散化后得到
π
=
∑
n
=
0
N
4
Δ
t
1
+
(
n
Δ
t
)
2
,
Δ
t
=
1
N
,
N
→
+
∞
(3)
π=\sum_{n=0}^N\frac{4\varDelta t}{1+(n\varDelta t)^2},\varDelta t=\frac{1}{N},N\to +\infty\tag{3}
π=n=0∑N1+(nΔt)24Δt,Δt=N1,N→+∞(3)
编程实现
这种方法的计算速度比较慢,计算结果的最小位数与
Δ
t
\varDelta t
Δt成线性关系,所以一般不用于计算高精度π。如果不追求高精度的话,用C++实现式(3)难度不大。
#include<cstdio>#include<ctime>#include<cmath>usingnamespace std;intmain(void){
clock_t startTime =clock();constlonglong N =1e9;double deltaT =1/(N *1.0);double sum =0;for(longlong i =0; i < N; i++){double x = i * deltaT;
sum +=4/(1+ x * x)* deltaT;}printf("%.12lf\n", sum);printf("The run time is: %.3fs\n",(clock()-(int)startTime)/1000.0);}
相对加速比:
S
(
P
)
=
4.871
1.424
≈
3.461
S(P)=\frac{4.871}{1.424}≈3.461
S(P)=1.4244.871≈3.461 并行化效率:
E
i
n
t
e
=
S
(
P
)
8
≈
0.428
E_{inte}=\frac{S(P)}{8}≈0.428
Einte=8S(P)≈0.428
三、欧拉恒等式
接下来使用欧拉恒等式实现高精度π的计算
算法推导
欧拉恒等式被认为是数学上最优美的公式之一,它将自然常数
e
e
e,圆周率
π
π
π,虚数单位
i
i
i,自然数
1
,
0
1,0
1,0这五个最基本的数字连接在一起:
e
i
π
+
1
=
0
(1)
e^{iπ}+1=0\tag{1}
eiπ+1=0(1) 我们可以由它计算出π,由该式容易得到
π
=
l
n
(
−
1
)
i
=
2
i
l
n
i
(2)
π=\frac{ln(-1)}{i}=\frac{2}{i}lni\tag{2}
π=iln(−1)=i2lni(2) 再由虚数的性质可以得到
l
n
i
=
l
n
1
+
i
1
−
i
=
l
n
(
1
+
i
)
−
l
n
(
1
−
i
)
(3)
lni=ln\frac{1+i}{1-i}=ln(1+i)-ln(1-i)\tag{3}
lni=ln1−i1+i=ln(1+i)−ln(1−i)(3) 首先对
l
n
(
1
+
x
)
ln(1+x)
ln(1+x)做泰勒级数展开,有:
l
n
(
1
+
x
)
=
x
−
1
2
x
2
+
1
3
x
3
−
1
4
x
4
+
1
5
x
5
−
…
…
(4)
ln(1+x)=x-\frac{1}{2}x^2+\frac{1}{3}x^3-\frac{1}{4}x^4+\frac{1}{5}x^5-\dots\dots\tag{4}
ln(1+x)=x−21x2+31x3−41x4+51x5−……(4) 取
x
=
±
i
x=±i
x=±i,得:
l
n
(
1
+
i
)
=
i
+
1
2
−
1
3
i
−
1
4
+
1
5
i
+
…
…
(5)
ln(1+i)=i+\frac{1}{2}-\frac{1}{3}i-\frac{1}{4}+\frac{1}{5}i+\dots\dots\tag{5}
ln(1+i)=i+21−31i−41+51i+……(5)
l
n
(
1
−
i
)
=
−
i
+
1
2
+
1
3
i
−
1
4
−
1
5
i
+
…
…
(6)
ln(1-i)=-i+\frac{1}{2}+\frac{1}{3}i-\frac{1}{4}-\frac{1}{5}i+\dots\dots\tag{6}
ln(1−i)=−i+21+31i−41−51i+……(6)
/**
* @brief 用数组存储两个整型数据相除结果
* @param result[] 结果保存数组
* @param y 被除数
* @param x 除数
* @param n 结果位数
* @return none
*/voidcalc_div(int result[],int y,int x,int n){for(int i =0; i < n; i++){
result[i]= y / x;
y = y % x;
y *=10;}}
算法实现
根据式(9)写出π的计算通式:
π
=
∑
i
=
1
N
(
−
1
)
i
−
1
4
2
i
−
1
(
1
2
2
i
−
1
+
1
3
2
i
−
1
)
,
N
→
+
∞
(10)
π=\sum_{i=1}^N(-1)^{i-1}\frac{4}{2i-1}(\frac{1}{2^{2i-1}}+\frac{1}{3^{2i-1}}),N\to +\infty\tag{10}
π=i=1∑N(−1)i−12i−14(22i−11+32i−11),N→+∞(10) 令
y
k
=
(
−
1
)
k
−
1
4
2
k
−
1
(
1
2
2
k
−
1
+
1
3
2
k
−
1
)
,
k
∈
Z
+
(11)
y_k=(-1)^{k-1}\frac{4}{2k-1}(\frac{1}{2^{2k-1}}+\frac{1}{3^{2k-1}}),k\in Z^+\tag{11}
yk=(−1)k−12k−14(22k−11+32k−11),k∈Z+(11) 则
π
=
∑
k
=
1
N
y
k
(12)
π=\sum_{k=1}^Ny_k\tag{12}
π=k=1∑Nyk(12) 我们可以对式(11)进行拆分,令
{
a
k
=
4
2
k
−
1
b
k
=
1
2
2
k
−
1
+
1
3
2
k
−
1
\begin{cases}a_k=\frac{4}{2k-1} \\b_k=\frac{1}{2^{2k-1}}+\frac{1}{3^{2k-1}}\end{cases}
{ak=2k−14bk=22k−11+32k−11,则:
y
k
=
(
−
1
)
k
−
1
a
k
b
k
(13)
y_k=(-1)^{k-1}a_kb_k\tag{13}
yk=(−1)k−1akbk(13) 因此可以分别计算
a
k
,
b
k
a_k,b_k
ak,bk,然后将两项的结果相乘得到
y
k
y_k
yk,再对
y
k
y_k
yk求和得到π
观察
a
k
,
b
k
a_k,b_k
ak,bk,发现
b
k
b_k
bk的分母呈指数级增长形势,当计算
{
y
k
}
\begin{Bmatrix} y_k\end{Bmatrix}
{yk}第16项时,
b
k
b_k
bk的分母将会超过上述除法操作的分母限制范围
(
−
2
31
,
2
31
−
1
)
(-2^{31},2^{31}-1)
(−231,231−1),同时也为了便于之后的OMP优化,我们继续对
b
k
b_k
bk进行分解变换。
1
2
2
k
−
1
=
1
2
×
1
2
×
⋯
×
1
2
⏞
2
k
−
1
=
1
2
n
×
1
2
m
×
1
2
m
×
⋯
×
1
2
m
⏞
[
2
k
−
1
m
]
(14)
\frac{1}{2^{2k-1}} =\overbrace{\frac{1}{2}\times \frac{1}{2}\times \dots \times \frac{1}{2}}^{2k-1} =\frac{1}{2^n}\times \overbrace{\frac{1}{2^m}\times\frac{1}{2^m}\times \dots \times \frac{1}{2^m}}^{[\frac{2k-1}{m}]}\tag{14}
22k−11=21×21×⋯×212k−1=2n1×2m1×2m1×⋯×2m1[m2k−1](14) 其中
n
+
[
2
k
−
1
m
]
×
m
=
2
k
−
1
n+[\frac{2k-1}{m}]\times m = 2k-1
n+[m2k−1]×m=2k−1,且
n
<
m
<
16
n<m<16
n<m<16。同理可以得到
1
3
2
k
−
1
\frac{1}{3^{2k-1}}
32k−11
constint N =500;// pi的精确位数constint TIMES =1000;// 算法迭代次数,即yk的项数int b[TIMES][N]={{0}};// 存放bk项,因为其占用空间较大,所以定义为全局变量intmain(void){
clock_t startTime =clock();int result[N]={0};// 存放最终结果intconst TDN =8;// m 或者理解为并行的线程数int x[N]={0};// 存放1/2的m次方int y[N]={0};// 存放1/3的m次方int(*x1)[N]=newint[TDN][N]{{0}};// 存放1/2的2k-1次方int(*y1)[N]=newint[TDN][N]{{0}};// 存放1/3的2k-1次方// 计算1/2的m次方calc_div(x,1,(int)pow(2,2* TDN), N);// 计算1/3的m次方calc_div(y,1,(int)pow(3,2* TDN), N);// 计算1/2、1/3的n次方for(int k =0; k < TDN; k++){calc_div(x1[k],1,(int)pow(2,2* k +1), N);calc_div(y1[k],1,(int)pow(3,2* k +1), N);}// 计算ak*bkfor(int k =0; k < TDN; k++){for(int i =0; i < TIMES / TDN; i++){// 计算bkcalc_add(b[i * TDN + k], x1[k], N);calc_add(b[i * TDN + k], y1[k], N);calc_multi(x1[k], x, N);calc_multi(y1[k], y, N);// 计算ak*bkint a[N]={0};// 存放akint t =2*(i * TDN + k)+1;calc_div(a,4, t, N);calc_multi(b[i * TDN + k], a, N);}}// 将最终结果相加/减for(int i =0; i < TIMES; i++){if(i %2==1)calc_sub(result, b[i], N);elsecalc_add(result, b[i], N);}printf("The run time is: %.3fs\n",(clock()-(int)startTime)/1000.0);// 打印piprintf("%d.", result[0]);for(int i =1; i < N; i++){printf("%d", result[i]);}printf("\n");delete[] x1;delete[] y1;}
运行程序有如下结果
OMP优化
分别在计算
b
k
b_k
bk、
a
k
b
k
a_kb_k
akbk代码段前添加以下语句,展开for循环。
#pragma omp parallel for
运行程序后有如下结果
仿照积分法,同样计算欧拉恒等式法的OMP优化指标
相对加速比:
S
(
P
)
=
10.524
1.937
≈
5.433
S(P)=\frac{10.524}{1.937}≈5.433
S(P)=1.93710.524≈5.433 并行化效率:
E
e
u
l
a
r
=
S
(
P
)
8
≈
0.679
E_{eular}=\frac{S(P)}{8}≈0.679
Eeular=8S(P)≈0.679 相比于未进行OMP优化时的10.5秒,优化后的执行时间显著减少,观察计算
b
k
b_k
bk、
a
k
b
k
a_kb_k
akbk的代码,将第一层for循环拆开,分配到不同的核上执行,由于循环体内的程序每一次循环时都和前一次循环无关,并且不同次数的循环修改的内存地址不同,因此将循环拆开后,仍能计算出正确的结果。这是典型的以空间换效率,在不同地址存放每一次循环的结果,这样就避免了数据竞争的情况。
四、总结
积分法与欧拉恒等式法的对比
首先对比上文中积分法和欧拉恒等式法的计算通式:
积分法
π
=
∑
i
=
0
N
4
N
1
+
(
i
N
)
2
,
N
→
+
∞
π=\sum_{i=0}^N\frac{\frac{4}{N}}{1+(\frac{i}{N})^2},N\to +\infty
π=i=0∑N1+(Ni)2N4,N→+∞
欧拉恒等式法
π
=
∑
i
=
1
N
(
−
1
)
i
−
1
4
2
i
−
1
(
1
2
2
i
−
1
+
1
3
2
i
−
1
)
,
N
→
+
∞
π=\sum_{i=1}^N(-1)^{i-1}\frac{4}{2i-1}(\frac{1}{2^{2i-1}}+\frac{1}{3^{2i-1}}),N\to +\infty
π=i=1∑N(−1)i−12i−14(22i−11+32i−11),N→+∞
积分法的收敛速度仅为
1
x
2
\frac{1}{x^2}
x21,为二次收敛,而本文所列举的欧拉恒等式法为指数收敛,收敛速度极快。
OMP实现方式的对比
对比两者的实现方式,都使用了将最外层循环展开的方法,然而两者的并行化效率却不一样,
E
i
n
t
e
=
0.428
<
E
e
u
l
a
r
=
0.679
<
1
E_{inte}=0.428<E_{eular}=0.679<1
Einte=0.428<Eeular=0.679<1。造成并行化效率小于1的原因有以下两方面:
每个并行体的计算量有差别
在分发并行的时候,系统需要消耗资源。
针对这两个因素,我们可以将程序移植到Linux系统(双ARM核)上来做对照实验。
1.积分法
无OMP优化
OMP优化
可以计算此时的并行化效率
E
i
n
t
e
−
l
i
n
u
x
=
S
(
P
)
2
≈
1
E_{inte-linux}=\frac{S(P)}{2}≈1
Einte−linux=2S(P)≈1
2.欧拉恒等式法
无OMP优化
OMP优化
同样可以计算:
E
e
u
l
a
r
−
l
i
n
u
x
=
S
(
P
)
2
≈
1
E_{eular-linux}=\frac{S(P)}{2}≈1
Eeular−linux=2S(P)≈1