多线程之基于积分法与欧拉恒等式法的圆周率计算及OMP优化

文章目录

一、问题描述
二、积分法
- 算法推导
- 编程实现
- OMP优化
三、欧拉恒等式
- 算法推导
- 编程实现
- - 前期准备
  - - 加法
    - 减法
    - 乘法
    - 除法
  - 算法实现
- OMP优化
四、总结
- 积分法与欧拉恒等式法的对比
- OMP实现方式的对比

一、问题描述

分别采用积分法和欧拉恒等式计算π，对比两种方法
使用OMP实现以上两种方法，再进行对比

二、积分法

算法推导

首先我们知道arctan(x)arctan(x)arctan(x)的导数f′(x)=11+t2f'(x)=\frac{1}{1+t^2}f′(x)=1+t21，所以有：
arctan(x)=∫0x11+t2dt(1)arctan(x)=\int_0^x\frac{1}{1+t^2}dt\tag{1} arctan(x)=∫0x1+t21dt(1)
取x=1x=1x=1得到
π4=arctan(1)=∫0111+t2dt(2)\frac{π}{4}=arctan(1)=\int_0^1\frac{1}{1+t^2}dt\tag{2} 4π=arctan(1)=∫011+t21dt(2)
离散化后得到
π=∑n=0N4Δt1+(nΔt)2,Δt=1N,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>
using namespace std;int main(void)
{clock_t startTime = clock();const long long N = 1e9;double deltaT = 1 / (N * 1.0);double sum = 0;for (long long 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);
}

运行结果：

OMP优化

首先添加头文件*<omp.h>*，在for循环前添加

#pragma omp parallel for reduction(+:sum)

reduction(+:sum)即为变量sum指定一个操作符+，每个线程都会创建变量sum的私有拷贝，在OMP区域结束后将使用各个线程的私有拷贝的值通过指定的操作符进行迭代运算，并赋值给原来的变量sum。

运行结果为：

一段程序可由 3 部分组成:准备(setup)、计算（compute)和结束(finalization) 部分，当计算量足够大时，可以忽略准备和结束部分所占用的时间，所以以下指标的计算只针对计算部分：

相对加速比：
S(P)=4.8711.424≈3.461S(P)=\frac{4.871}{1.424}≈3.461 S(P)=1.4244.871≈3.461
并行化效率：
Einte=S(P)8≈0.428E_{inte}=\frac{S(P)}{8}≈0.428 Einte=8S(P)≈0.428

三、欧拉恒等式

接下来使用欧拉恒等式实现高精度π的计算

算法推导

欧拉恒等式被认为是数学上最优美的公式之一，它将自然常数eee，圆周率πππ，虚数单位iii，自然数1,01,01,0这五个最基本的数字连接在一起：
eiπ+1=0(1)e^{iπ}+1=0\tag{1} eiπ+1=0(1)
我们可以由它计算出π，由该式容易得到
π=ln(−1)i=2ilni(2)π=\frac{ln(-1)}{i}=\frac{2}{i}lni\tag{2} π=iln(−1)=i2lni(2)
再由虚数的性质可以得到
lni=ln1+i1−i=ln(1+i)−ln(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)
首先对ln(1+x)ln(1+x)ln(1+x)做泰勒级数展开，有：
ln(1+x)=x−12x2+13x3−14x4+15x5−……(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=±ix=±ix=±i，得：
ln(1+i)=i+12−13i−14+15i+……(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)

ln(1−i)=−i+12+13i−14−15i+……(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)

由上式得到
π4=1−13+15−……(7)\frac{π}{4}=1-\frac{1}{3}+\frac{1}{5}-\dots\dots\tag{7} 4π=1−31+51−……(7)
这个公式是莱布尼茨级数，但是用其来求圆周率π的效率过低，迭代10万次才能精确到小数点后六位。所以继续变换lnilnilni：
KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ lni=ln\frac{(2…
最后可以得到如下计算圆周率的公式：
π4=(12+13)−13(123+133)+15(125+135)−……(9)\frac{π}{4}=(\frac{1}{2}+\frac{1}{3})-\frac{1}{3}(\frac{1}{2^3}+\frac{1}{3^3})+\frac{1}{5}(\frac{1}{2^5}+\frac{1}{3^5})-\dots\dots\tag{9} 4π=(21+31)−31(231+331)+51(251+351)−……(9)
用这个公式只需要取前4项就可以达到祖冲之的效果。当然还可以推出其他效率更高的公式，但是以下的程序都以这个公式为框架。

编程实现

前期准备

为了实现高精度π的计算，我们使用数组来存储π，所以对应地需要实现数组间的加减乘除操作

加法

/*** @brief 实现两个数组的加法* @param b            加数1，以及结果保存数组* @param a          加数2* @param n          结果位数* @return none*/
void calc_add(int b[], int a[], int n)
{int carry = 0;for (int i = n - 1; i >= 0; i--){b[i] = a[i] + b[i] + carry;carry = b[i] / 10;b[i] %= 10;}
}

减法

/*** @brief 实现两个数组的减法* @param b            被减数，以及结果保存数组* @param a          减数* @param n           结果位数* @return none*/
void calc_sub(int b[], int a[], int n)
{int carry = 0;for (int i = n - 1; i >= 0; i--){b[i] = b[i] - a[i] - carry;carry = b[i] < 0;b[i] = carry ? (10 + b[i]) : b[i];}
}

乘法

/*** @brief 实现两个数组的乘法* @param a            乘数1，以及结果保存数组* @param b          乘数2* @param n          结果位数* @return none*/
void calc_multi(int a[], int b[], int n)
{int* result = new int[2 * n]{ 0 };for (int i = 0; i < n; i++){int* c = new int[2 * n]{ 0 };for (int j = 0; j < n; j++){c[i + j] = a[i] * b[j];}calc_add(result, c, n + n / 10);delete[] c;}for (int i = 0; i < n; i++){a[i] = result[i];}delete[] result;
}

除法

除法操作与其他操作有所不同，只用一个数组来保存结果，并且限制分母范围为(−231,231−1)(-2^{31},2^{31}-1)(−231,231−1)

/*** @brief 用数组存储两个整型数据相除结果* @param result[]   结果保存数组* @param y           被除数* @param x          除数* @param n           结果位数* @return none*/
void calc_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=1N(−1)i−142i−1(122i−1+132i−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)
令
yk=(−1)k−142k−1(122k−1+132k−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=1Nyk(12)π=\sum_{k=1}^Ny_k\tag{12} π=k=1∑Nyk(12)
我们可以对式(11)进行拆分，令{ak=42k−1bk=122k−1+132k−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，则：
yk=(−1)k−1akbk(13)y_k=(-1)^{k-1}a_kb_k\tag{13} yk=(−1)k−1akbk(13)
因此可以分别计算ak,bka_k,b_kak,bk，然后将两项的结果相乘得到yky_kyk，再对yky_kyk求和得到π

观察ak,bka_k,b_kak,bk，发现bkb_kbk的分母呈指数级增长形势，当计算{yk}\begin{Bmatrix} y_k\end{Bmatrix}{yk}第16项时，bkb_kbk的分母将会超过上述除法操作的分母限制范围(−231,231−1)(-2^{31},2^{31}-1)(−231,231−1)，同时也为了便于之后的OMP优化，我们继续对bkb_kbk进行分解变换。
122k−1=12×12×⋯×12⏞2k−1=12n×12m×12m×⋯×12m⏞[2k−1m](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+[2k−1m]×m=2k−1n+[\frac{2k-1}{m}]\times m = 2k-1n+[m2k−1]×m=2k−1，且n<m<16n<m<16n<m<16。同理可以得到132k−1\frac{1}{3^{2k-1}}32k−11

const int N = 500;          // pi的精确位数
const int TIMES = 1000;        // 算法迭代次数，即yk的项数
int b[TIMES][N] = { {0} }; // 存放bk项，因为其占用空间较大，所以定义为全局变量int main(void)
{clock_t startTime = clock();int result[N] = { 0 }; // 存放最终结果int const TDN = 8;  // m 或者理解为并行的线程数int x[N] = { 0 };  // 存放1/2的m次方int y[N] = { 0 };  // 存放1/3的m次方int(*x1)[N] = new int[TDN][N]{ {0} }; // 存放1/2的2k-1次方int(*y1)[N] = new int[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优化

分别在计算bkb_kbk、akbka_kb_kakbk代码段前添加以下语句，展开for循环。

#pragma omp parallel for

运行程序后有如下结果

仿照积分法，同样计算欧拉恒等式法的OMP优化指标

相对加速比：
S(P)=10.5241.937≈5.433S(P)=\frac{10.524}{1.937}≈5.433 S(P)=1.93710.524≈5.433
并行化效率：
Eeular=S(P)8≈0.679E_{eular}=\frac{S(P)}{8}≈0.679 Eeular=8S(P)≈0.679
相比于未进行OMP优化时的10.5秒，优化后的执行时间显著减少，观察计算bkb_kbk、akbka_kb_kakbk的代码，将第一层for循环拆开，分配到不同的核上执行，由于循环体内的程序每一次循环时都和前一次循环无关，并且不同次数的循环修改的内存地址不同，因此将循环拆开后，仍能计算出正确的结果。这是典型的以空间换效率，在不同地址存放每一次循环的结果，这样就避免了数据竞争的情况。

四、总结

积分法与欧拉恒等式法的对比

首先对比上文中积分法和欧拉恒等式法的计算通式：

积分法

π=∑i=0N4N1+(iN)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=1N(−1)i−142i−1(122i−1+132i−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→+∞

积分法的收敛速度仅为1x2\frac{1}{x^2}x21，为二次收敛，而本文所列举的欧拉恒等式法为指数收敛，收敛速度极快。

OMP实现方式的对比

对比两者的实现方式，都使用了将最外层循环展开的方法，然而两者的并行化效率却不一样，Einte=0.428＜Eeular=0.679<1E_{inte}=0.428＜E_{eular}=0.679<1Einte=0.428＜Eeular=0.679<1。造成并行化效率小于1的原因有以下两方面：

每个并行体的计算量有差别
在分发并行的时候，系统需要消耗资源。

针对这两个因素，我们可以将程序移植到Linux系统（双ARM核）上来做对照实验。

1.积分法

无OMP优化
OMP优化

可以计算此时的并行化效率
Einte−linux=S(P)2≈1E_{inte-linux}=\frac{S(P)}{2}≈1 Einte−linux=2S(P)≈1

2.欧拉恒等式法

无OMP优化
OMP优化

同样可以计算：
Eeular−linux=S(P)2≈1E_{eular-linux}=\frac{S(P)}{2}≈1 Eeular−linux=2S(P)≈1

通过上述对比，我们可以发现，在Linux下两种方法的并行化效率大致相同，且都约等于1。可见两种方法的并形体计算量基本一致，区别在于操作系统分发并行任务时的开销不同。

接着分析两种方法的并行化效率不同的区别。

本次实验的处理器为四核八线程的英特尔酷睿i7-7700HQ，就是说除了四个核心所能处理的四线程外，它还拥有四个超线程。

超线程HT（Hyper-Threading）技术是在单个核心处理单元中集成两个逻辑处理单元，也就是一个实体内核（共享的运算单元）拥有两个逻辑内核（有各自独立的处理器状态）。而其余部分如ALU（整数运算单元）、FPU（浮点运算单元）、L2 Cache（二级缓存）则保持不变，这些部分是被分享的。

这样就可以看出差别，本次实验中积分法大部分运算为浮点运算，进行OMP加速时运行在同一个核心上的两个线程将会争夺一个FPU，所以其加速比不能达到或接近核心（超线程）数量。而欧拉恒等式法全为整型运算，运行在同一个核心上的两个线程能够分别使用ALU和FPU，不会出现竞争情况，所以其加速比大于实际核心数量。

完成程序在这里https://download.csdn.net/download/qq_42688495/12291963