目录

一、PI的计算方法

二、代码说明

1、类定义

1）class fftpack：

2）class longfloat：

3) class fft：

4）class refiner：

5）class pulser：

6）Class P2B3：

7）Class xthread：

8）Class xtask、Class xqueue：

9）Class mulwise

三、优化设计及版本历史

初始版本

2021年4月7日：

1、longfloat内部序列使用二进制代替10进制

2、   优化绝对值小于1的longfloat的有效长度

3、   使用P2B3模型，优化非迭代公式（适用于method < 30的所有PI计算公式）的计算方式

4、使用Karatsuba机制，提升长乘法的效率和精度

5、使用gcc/g++嵌入式汇编指令优化除法

6、其他

未来优化方向

1、使用SSE/SSE2、AVX指令集优化

2、使用并行计算优化

3、进制转换优化

4、增加乘数FFT缓存

2021年4月14日：

1、在调试新功能中发现数据或者内存异常，从而发现了之前版本隐藏较深的BUG，这些BUG在之前的版本运行中未暴露出来，实属意外，目前已经修正：

2、进制转换算法的实现和优

3、修改了pulser的实现机制

2021年4月16日：

2021年4月18日

1、发现并修正了程序中存在的个别缺陷和隐患

2、发现Borwein的9次方迭代求PI公式（Method=30）可以通过公式变换，减少1次大数乘法运算：

3、对所有使用refiner进行迭代精度控制的函数进行微调，从之前在每轮迭代结束时检查精度是否达标，前置到迭代过程中的某一步，发现精度达标即退出循环，这样可跳过后面无需执行的步骤。

2021年4月22日

一、PI的计算方法

程序中实现了4类14种PI的计算方法，每种方法都实测过，确保计算结果的准确性，所有方法至少试验过1~10万位的求解，输出有效数字完全相同，并与y-cruncher等标准程序输出结果进行了对照。

  1. Method 1-7为类Machin公式，也就是反正切函数的Taylor级数；程序中实现了传统累加法和P2B3模型法（P2B3模型的含义在后面章节有专门介绍）。
  2. Method 10-16为类BBP级数，包括标准BBP级数、Fabrice Bellard级数，程序实现了传统累加法和P2B3/BBP模型法。
  3. Method 20-21为类阶乘级数，业界求PI最常使用的Ramanujan公式和Chudnovsky公式，程序实现了传统递进累加法和P2B3模型法。使用了P2B3模型的这两个公式在所有方法中表现出最快的计算性能。
  4. Method 30-36为理论收敛速度最快的Borwein迭代法，包括9次收敛迭代、2次收敛迭代、4次收敛迭代、3次迭代、5次迭代、16次迭代。尽管理论收敛速度差异很大，但由于迭代步骤的复杂程度，实测中运行最快的是2次迭代法（Method=31），4次迭代性能相似（Method=32）、9次迭代和16次迭代实际运行速度不佳（Method=30、36）。

以下为各种求PI方法的实测性能，注意观察在相同位数情况下各方法的耗时对比，以及不同算法在位数增大10倍时的耗时变化情况，以对各种算法的复杂度有清晰的认识。部分算法因其性能不足以应对太高位数的计算，因此未进行实测。例如对于采用累加求和方式计算的级数，因其算法复杂度为O(n2)，也就是说N=100万所消耗的时间为N=10万时的100倍，因此当N为100万时，只实际测试了Ramanujan、Chudnovsky两种最快的方法（其10万位耗时在2秒左右，预测其计算100万位耗时约为几分钟）；表中的计算时间单位为秒，该时间不包括最终结果转化为10进制及打印输出的时间。同时，该耗时为在未进行并发优化的单线程模式下的运行时间。

Method

类别

方法

选项

10万

100万

1000万

1

Taylor级数

类Machin公式

P2B3

1.6

26

-

累加

11.4

-

-

2

P2B3

2.0

31

-

累加

9.0

-

-

3

P2B3

1.8

25

-

累加

7.7

-

-

4

P2B3

2.1

31

-

累加

12.4

-

-

5

P2B3

1.7

26

-

累加

8.9

-

-

6

P2B3

2.1

27

-

累加

8.6

-

-

7

P2B3

2.0

28

-

累加

8.8

-

-

10

BBP分数级数

BBP(16)

P2B3

4.3

63

-

累加

54.7

-

-

11

Fabrice Bellard(1024)

P2B3

3.1

55

-

累加

38.9

-

-

12

BBP(256)

P2B3

2.6

46

-

累加

55.7

-

-

13

BBP(256)

P2B3

2.5

47

-

累加

58.6

-

-

14

BBP(4096)

P2B3

2.6

44

-

累加

57.1

-

-

15

BBP(4096)

P2B3

3.9

67

-

累加

76.3

-

-

16

BBP(65536)

P2B3

3.3

50

-

累加

56.8

-

-

20

阶乘级数

Ramanujan

P2B3

0.7

11

214

累加

2.3

274

-

21

Chudnovsky

P2B3

0.5

6.5

179

累加

1.8

188

-

30

牛顿迭代

Borwein(9)

-

3.1

61

31

Gauss -Legendre(2)

-

2.3

33

32

Borwein(4)

-

2.4

31

33

Borwein(2)

-

3.7

65

34

Borwein(3)

-

2.3

41

35

Borwein(5)

-

4.1

43

36

Borwein(16)

-

2.5

44

二、代码说明

1、类定义

1）class fftpack：

用于将longfloat序列转换为适合FFT计算的序列；配合fft类使用将该类从fft中独立出来，是之前曾经想同时引入FFT和NTT算法，二者均需要使用相似的unpack/repack过程，将longfloat的序列转换为适合FFT/NTT计算的序列（unpack），并在计算完成后对结果进行逆变换（repack）。

2）class longfloat：

高精度浮点数类，使用科学计数三元法表示：<sign>A0.A1A2A3...An * U<exp>；sign为符号，A0/A1/A2为有效位，U为位基数（231），exp为指数该类中所有实例化对象拥有相同的最大精度（prec），即小数点后二进制位个数该类重载了主要的加、减、乘、除运算符，以及流输出运算符，并实现了乘方（sqr）、开方（sqrt）、阶乘（factorial）等函数该类的部分函数，包含了多种实现机制，根据需要灵活选择，程序中默认开启的是实测性能最高、且不损失精度的机制，其他机制的实现代码保留作为参考和结果校验。例如：

  1. 大数倒数（reciprocal， 1 / longfloat）：实现了试商法和牛顿迭代法。其中试商法即常规的用1作为被除数去除以该大数；程序根据该大数的位数，优选牛顿迭代法。
  2. 大数除法（longfloat / longfloat）：实现了试商法和乘倒数法（+dm/-dm）；其中试商法的每位试商，实现了最大值法(div_try_bit_max)、二分法（div_try_bit_bin）、倍增法（div_try_bit_dbl）。正常情况下优选乘倒数法，试商法仅用于除数长度很短的情形；试商所用机制为最大值法
  3. 大数开方（sqrt）：实现了试商法（division_sqrt，仅用于开平方，+qd/-qd）、正牛顿迭代法（iterative_sqrt）、倒数牛顿迭代法（inverse_sqrt）。正常情况下优选倒数牛顿迭代法
  4. 大数乘法：实现了基本算法（long_mul）和FFT算法（fft_mul），未强行指定（+-mul2fft）时根据位长自动选择所用机制，并使用Karatsuba算法在位数超过阈值时进行自动拆分-合并

3) class fft：

用于实现FFT变换、FFT乘法；考虑到应用中所有大数具有类似位长精度，该类尽量保存先前计算出的单位根和共轭根，并利用单位根的循环性质，在后续的所有大数乘法尽可能重复利用。

另外，标准的FFT乘法需要对2个乘数分别做1次FFT正变换，对结果做1次FFT逆变换。由于FFT变换是O(N*log2(N))复杂度的，是FFT乘法的性能关键，程序中使用了耦合-解耦机制，将两个乘数耦合在1组复数中，只做1次FFT正变换，然后解耦合得到各自正变换的值，这样可有效提升性能

之前的十进制版本中，还尝试使用过多模数FFT，但其复杂度、准确性和性能都不理想，故放弃

4）class refiner：

用于牛顿迭代法中的精度渐增控制。牛顿迭代法是逐步求精的，每步迭代出的值越来越精确，因此在起初的迭代步骤中，并不需要保留太多的位数，可以只使用较小的精度；而随着迭代次数增加，则逐步增加存储的数字位数以达到所需精度。这样可大大缩短前面若干次迭代的运算速度在用牛顿迭代法求倒数（reciprocal）、方根（iterative_sqrt）、方根倒数（inverse_sqrt）时，均使用了该精度控制机制，实测证明有效。

但是但是，对于3个Borwein迭代求PI的公式（Method = 30 / 31 /32），虽然使用了refiner，但实际并未采用渐进精度机制，而是一开始就使用全精度，否则不能保证迭代的收敛性，目前尚未有经过严格推敲的合理解释，推测原因是可能初始值已经包含了决定最终精度的重要信息，如果初始时截断精度，则会破坏迭代的收敛速度甚至导致无法收敛。未来可以进一步理清这个原因，给出更为科学严谨的结论。

5）class pulser：

用于在冗长的计算过程中，以适当的时间间隔输出中间结果信息，以给予调试者正向反馈。pulser的实现先后演进了3个版本：自动步速控制、独立定时器线程、定时器信号触发（最终版本）自动步速控制是在程序中每轮循环调用check()时，测算两次时间间隔之间的循环次数（即步速），根据最近的步速确定下一次需要触发输出的时机。

独立定时器线程是由单独的一个线程按秒进行计时，调用者通过检查两次之间计时器的间隔决定是否需要触发输出。定时器信号触发是利用Unix/Linux系统的信号机制（signal），定时触发回调函数进行计时，调用者通过检查两次计时器间隔决定是否需要触发输出。所有这些机制的目的都是避免调用者在主循环中频繁调用系统函数获取当前系统时间。

6）Class P2B3：

对适用的级数公式实施P2B3/BBP计算模型的基类。关于P2B3模型的介绍在后面叙述。P2B3基类实现了二分法计算P(a,b)、Q（a,b）、R（a,b）的算法，由每个子类通过calc纯虚函数提供单项P(b)、Q(b)、R(b)的计算方法即可。

Machin级数、BBP级数、Ramanujan级数、Chudnovsky级数均从该类继承出子类，完成自身的主体部分计算。

7）Class xthread：

并发多线程控制类。每个线程对象可单独执行并发任务，并提供线程间的同步互斥控制接口。主线程对象作为并发入口，负责根任务的拆分，并和其他线程一起执行所拆分处的子任务，并在所有子任务完成后进行结果的合并。任务的管理通过下面的任务队列类进行。

8）Class xtask、Class xqueue：

并发任务队列管理类。并发任务通过加入共享的任务队列，由各线程对象提取并执行。队列按照FIFO机制进行管理。

9）Class mulwise

简单乘除法合并优化类。当需要将一个longfloat型的大数连续乘以或者除以若干个整数时，通过将这些乘数/除数合并成大的整数，减少longfloat参与的乘法/除法次数。在Ramanujan级数、Chudnovsky级数的累加求和实现机制中被应用以提高计算性能。

三、优化设计及版本历史

初始版本

程序的初始版本设计如下：

  1. longfloat使用十进制的科学计数法表示一个高精度浮点数，每个单元（4字节整型）存储N位十进制位（N=9），计算过程按照十进制进位体系，实现了加、减、乘、除基本运算，以及乘方、开方、倒数等函数。
  2. 实现了以十进制为基底的FFT乘法，优化大数乘法。
  3. 实现了7种Machin公式、2种BBP级数、2种阶乘级数（Ramanujan、Chudnovsky）和3种Borwein迭代公式求PI的方法，并进行了验证。

2021年4月7日：

1、longfloat内部序列使用二进制代替10进制

1）先前序列中每项代表N位十进制数。N（ucdbits）默认为9，即基数（ucgrade）为109，为不超过32位字长的最大10的幂

2）新的序列中每项代表N位二进制数。N（ucbits）默认为31，即基数（ucgrade）为231，为不超过32位字长的最大2的幂

3）由于基数231 > 109，使用新序列在空间存储效率上略有提升，提升比例约为4% ）。原理上也可以使用N=32位，可以进一步提升空间效率（3%），但因为次数基数232超过32位字长，在编程中需要多处注意更换变量类型，否则更容易出现数值异常，故未采用

5）使用二进制表示方式后，内部计算过程中均按十六进制方式显示longfloat常量/变量；在显示最终结果时，调用convert10将其转换为十进制进行显示

6）增加了convert10()、convert16()接口，用于切换进制

7）FFT乘法也同时变更为2进制基数，FFT项从先前的4个十进制位（104=10000），变化为14个二进制位（214=16384），在FFT存储效率上提升5%，叠加之前的序列项效率提升，FFT乘法空间效率提升9%左右

8）使用二进制表示方式后，对于乘数、除数为2的整数次幂的乘除运算（longfloat vs int），可以使用二进制移位操作，大大提升计算效率，特别适用于2个BBP的求PI公式（Method = 10、11，每项中均会除以16或1024的N次幂）

2、   优化绝对值小于1的longfloat的有效长度

1）在Taylor级数（Machin公式）、BBP公式、Ramanujan公式、Chudnovsky公式的传统计算过程中，后面累加的项的数值指数级减小，对最终结果的影响主要是最后几位

2）这种情况下，无需保留太多有效数字，可以认为小数点后的前导0（即exp的绝对值）也是有效位数，计算后面项的除法（longfloat / int）中需要计算的位数越来越少，计算速度越来越快

3）优化方法为在longfloat的加减乘除函数中，根据指数（exp） + 最大精度（dbmax）选取操作数和结果的有效位数（eb）

4）采用该机制，可以成倍提高上述公式的计算效率。计算位数越多时效果越明显，往往能提高几倍的性能

5）可以通过命令行参数+clz/-clz打开或者关闭该优化选项，以进行对比和校验。原理上说，使用-clz能更精确，但实测中发现打开选项并未影响实际最终计算精度

3、   使用P2B3模型，优化非迭代公式（适用于method < 30的所有PI计算公式）的计算方式

1）Taylor级数（Machin公式）、BBP公式、Ramanujan公式、Chudnovsky公式的传统计算方法均是逐项累加法，各累加项的计算过程主要使用短乘法（longfloat * int）、短除法（longfloat / int），这种方式编程简单，但需要执行大量的除法指令

2）P2B3模型，将公式变换为先计算几个大数P、Q、R（P2B3模型）或P、Q（BBP模型），然后用P、Q、R进行几次乘、除运算得到最终结果

例如计算：1/2 + 1/3 + 1/4，传统是执行3次短除法（longfloat / int）和2次长加法（longfloat + longfloat），变换后为计算：(P = 3*4+2*4+2*3)/(Q = 2*3*4)，执行5次长乘法（longfloat * longfloat）、2次长加法（longfloat + longfloat）、1次长除法（longfloat / longfloat）

其中长乘法采用FFT进行优化；长除法a/b采用a * (1/b)进行优化（参见wise_div），而倒数1/b计算又可采用牛顿迭代法转换为长乘法（参见iterative_reciprocal函数），因此将传统计算中主要耗时的除法指令，转换为了以经过FFT优化的乘法运算为主

3）P、Q、R的计算过程，采用二分法模型进行递归，递归算法如下。可将计算复杂度大大降低

P(a,b) = P(a,m)Q(m,b) + P(m,b)R(a,m)

Q(a,b) = Q(a,m)Q(m,b)

R(a,b) = R(a,m)R(m,b)

其中m为a和b的中间值

4）程序中，上述所有PI的计算方法都同时实现了传统累加模型（accum）和P2B3模型（P2B3），以进行对比和校验；可以通过命令行开关+pb/-pb进行设置

5）根据测试结果，同一公式适用P2B3后，其运行速度比累加模型提高数十倍。例如对于非迭代求PI公式中的王者Chudnovsky公式，使用传统累加方式计算10万位的速度约为11.6秒，采用P2B3后速度仅为0.4秒（这里不包括最后转换为10进制表示形式），速度提高30倍

6）程序中使用抽象基类P2B3，实现二分递归机制，后面的每个实际公式均从该基类派生，实现各自的eval函数，计算P(b-1,b)、Q(b-1,b)、R(b-1,b)，以便未来扩展支持其他适用的公式

4、使用Karatsuba机制，提升长乘法的效率和精度

1）之前版本的长乘法中，根据乘法操作数的长度是否大于某个阈值，决定采用基本乘法（long_mul）还是FFT乘法（fft_mul）机制。该阈值根据实测经验定为1024

2）调测中遇到两个问题：当长度太大时（求解几百万位），使用基本乘法性能固然很低（复杂度为O(N2)），但使用FFT会带来精度损失。

3）原因是FFT中主要使用double浮点数运算，必须精确控制好不产生精度丢失，之前只采用了限制FFT项的基数合理（13-14bit），但项数对FFT精度也有影响，实测发现，在FFT基数为214时，乘法操作数的长度超过3万多时会造成计算偏差

4）为平衡性能和精度，采用分段计算方法，当操作数长度超过FFT阈值时，使用二分法将操作数拆分为2个较小的段进行计算，再对结果进行合并，基本原理为：

（a1 + a2*E） * (b1 + b2*E ) = a1b1 + (a1b2 + a2b1)*E + a2b2*E2

其中E为拆分处对应的指数（E = (231)-K，K为拆分点）。这样将原来的一个超阈值的长乘法，转换为阈值内的4次长乘法+3次长加法。计算复杂度由 O(2N*log2(2N))变为4*O(N*log2N) + 3*O'(N)，虽然保证了计算准确性，但计算性能不增反减

  1. 进一步公式变换如下：

（a1+a2*E）*(b1+b2*E ) =a1b1*(1-E)+(a1+a2)(b1+b2)*E - a2b2*(E-E2)

其中只包含3次降阶的长乘法和6次长加减法，性能可提高约25%，在长度较大时，还可以对其中的3次2/N长度的乘法进一步拆分为9次4/N长度的乘法，每拆分1次，性能提高25%；最终实际性能甚至可优于拆分前的性能

6）程序的wise_mul中包含了长乘法的算法控制代码

5、使用gcc/g++嵌入式汇编指令优化除法

1）使用asm汇编指令执行long long / int的计算，一次得到整数商和余数，余数用于序列中下一项的计算。可以将耗时的除法运算量减小一半，提升此类计算公式30%-40%的性能

6、其他

1）mulwise机制已经在先前版本引入，目的是合并乘数/除数以减少实际需要执行的短乘法/短除法（longfloat vs int）的次数，在Ramanujan公式、Chudnovsky传统累加计算中有明显的提升效果，约提高30%左右性能

2）中间版本尝试过引入NTT（数论变换）替代FFT，以用整数运算代替double浮点运算，实现性能和精度提升。但实测NTT的效果很不理想，一方面涉及到大量的64位整数除法和求余运算，另外由

于计算机整数的表示范围，很容易发生运算溢出，其不损失精度时的基数和项数阈值都远小于FFT。尽管有人尝试通过多模数NTT来解决此问题，但其引入的编程和计算复杂度都大大升高，最终放弃

3）在未实现P2B3模型优化之前，各种求PI的算法中，当位数小于10万时，最快的是Chudnovsky公式（Method=21），计算10万位10秒左右；当位数为几十万、上百万时，只能采用Borwein迭代法（Method=30/31/32），其中100万位时用Gause-Legendre迭代公式用时60秒左右，Chudnovsky公式需要约3分钟

Borwein迭代的优势是当求解位数很大的时候，因其非线性收敛特性，所需的迭代次数很少；缺点是每步迭代的计算比较复杂，都需要平方、开方等，因此在10万位以下的性能优势不明显

Chudnovsky公式及之前的所有公式都是线性收敛的，优点是后续项使用前导项递推，计算相对简单；缺点是当要求解的位数增加10倍时，其所需时间增加100倍

但使用P2B3模型优化后，Chudnovsky公式重新成为目前最快的计算方法，计算100万位大致耗时16秒左右

未来优化方向

1、使用SSE/SSE2、AVX指令集优化

1）理论上，使用SSE、AVX指令后，可以一次指令计算128/256位的浮点数运算，即2个或者4个double浮点数的计算，也可以用于32位/64位整数的定点运算

2）实际上，为了能有效利用其提升效率，程序的内存数据组织方式、对齐方式、存取方式均需要进行较大的调整，否则实际提升效果不会明显，甚至还会带来额外的性能损失。这是因为先需要将native的int、double类型操作数加载到SSE/AVX寄存器，计算完成后再存回native操作数；

3）容易想到的可用点为FFT的复数乘法、加法运算，实现实部、虚部的同时乘加；但编写的简单测试程序发现，性能没有明显提升，故放弃

4）longfloat的加法（减法）理论上也可以用SSE指令加速，但因为longfloat加法涉及到进位处理、指数对齐，并且长度不一定是mm128的整数，因此还需要一系列的额外运算。在各个求PI公式的计算中，占据核心运算量的都是FFT长乘法（Machin公式、BBP公式、Ramanujan公式、Chudnovsky公式均已经转换为P2B3模型后，也以长乘法为主），而FFT的核心计算则为复数（即双double浮点数）运算，长加、长减所占用的比例很少，不值得进行SSE优化

3）测试的虚拟主机上不支持AVX指令，只能使用128位的SSE。这样一个复数乘法运算，仍需执行2次mm128的SSE乘法指令，整体性能提升有限

2、使用并行计算优化

1）并行计算可以在多CPU、多核处理器环境下，充分调动CPU的利用率

2）并行计算的两个优化方向：

A、针对具体求PI公式优化：将原公式拆解为2个或更多个独立的互不依赖的分支公式，每个分支公式用1个线程分别求解，最后将各个分支公式计算结果合并

优点：有针对性，在较高粒度启动并行任务，可以一开始就创建好多个分支，并让其连贯运行到底，CPU的并行利用效率更高

缺点：每个公式需要单独设计和实现并行优化策略；另外，有的公式（如Chudnovsky）的传统实现机制为根据上一项推导出下一项，改为并行后，必然会引入一些无效或者重复计算量（例如，虽然可以用2个分支分别计算奇偶子序列，但每个分支的前后项推导过程变的更复杂），不仅程序会更复杂，实际效果会有折扣，一般都会有1 + 1 < 2，如1 + 1 = 1.5的实际性能

B、针对基本运算进行优化：对longfloat的基本操作（加、减、乘、除等）进行并行优化，将长序列的运算分解为N个短序列的同性质运算，最后将运算结果合并

优点：在底层细粒度的优化，可以作用于所有上层算法

缺点：优化层次太低，会导致频繁创建和结束细粒度并发任务，带来额外的性能损失；基本操作并行优化并不容易，如加法涉及到进位，分支结果合并时可能导致还需要多执行一遍O(N)的合并操作；除法则更难，因为前一项的余数是决定下一项结果的关键，只能串行执行

3）由于大数乘法占据绝大部分的运行时间，因此未来仍重点考虑对大数乘进行并行优化，初步思路有2个：

A、在wise_mul中已经使用Karatsuba机制，拆分为3个子段的乘法，再进行合并，这3个子段乘法可以考虑并行计算

B、在fft_mul中，FFT变换和逆变换因为具有串行逻辑不能分拆合并，但FFT变换后的复数向量点积运算是可以拆为多个独立计算任务的

4）实现并行优化时，需要考虑一套较好的细粒度子任务控制机制，尽量减少大量派生、消亡线程带来的额外开销。这套机制已经设计好，但在之前针对大数乘法进行测试时，实测效果不明显，故未将其代码纳入。

3、进制转换优化

目前longfloat实现了convert10和convert16函数，完成十进制和二进制表示法的双向转换。转换函数已经使用指数平移手段，使整个转换过程都只用乘法实现（包括长乘法和短乘法）。但目前在针对超大数列转换的运行速度还是不理想，最终结果转换为十进制所消耗的时间甚至远远大于前面整个求解过程的时间。例如用Chudnovsky公式计算PI的1000万位时，计算用时7分钟，但转换用时多达45分钟！！因为原理上每转换出1位，需要执行1次大数乘法，其复杂度为O(N2)，未来可考察有否进一步优化机制

2021.4.14补记：4.14版本已经利用分段法实现了进制转换的优化，在计算较高的位数时，可大大缩短耗时，时间复杂度从O(n2)降低为O(n)，1000万位长度的浮点数转换为10进制缩短为3分钟。

4、增加乘数FFT缓存

由于在最终的各种优化后的计算公式以及进制转换中，大数乘法是其中耗时最大的部分，因此对大数乘的优化也是最关键的部分，上面提到的SSE、并行计算等都是围绕大数乘法引入的CPU优化机制：指令集和多线程。

由于大数乘法的关键是FFT，而FFT的主要性能集中在FFT变换和逆变换，所有的longfloat形式的操作数，需要首先进行FFT变换为复数形式，计算完成后再逆变换为longfloat形式，如何减少不必要的变换就是关键。可以考虑机制包括：

（1）尽量缓存最近的几个使用频度较高的longfloat的FFT变换结果，未来检测到是同一操作数时，直接从缓存中得到其FFT变换结果，节省1次变换操作,这种缓存在乘方运算、牛顿迭代公式、进制变换过程中都能发挥出实际的效用，因为会多次乘以同一个大数。但是这里需要一个复杂的控制机制，考虑到：

  1. 在缓存的FFT结果和longfloat浮点数之间建立有效的映射关系，而又尽可能对上层计算过程透明，避免其深度参与底层计算控制。因为某个longfloat浮点数存储的结果、位长会在计算过程中发生变化，此时其缓存中的结果不再有效，如果每次变化都去通知缓存管理机构会造成代码逻辑过于复杂和晦涩，但靠缓存管理机构自己去判别longfloat是否发生变化也比较困难。这里的思路可以参考前面的“addr归位机制”的实现
  2. 由于FFT乘法需要对2个乘数进行FFT变换，即执行2次FFT变换，缓存的目的要减少变换的次数为1或者0；但实际两个乘数都在缓存中的概率很低，因此优化效果也就是将FFT次数减为1。而当前版本中，已经实现了通过“耦合-解耦”机制只执行1次FFT即得到两个乘数的FFT结果。即使加上缓存，除非两个乘数都在缓存中，否则还是至少需要1次FFT，对性能的提升效果极为有限（当然可以节省的是解耦合过程中的O(n)的复数加减运算，但其相对于变换本身的复杂度应该是微不足道的）

(2) 对于需要连续进行大数乘的过程（如P2B3中的P、Q、R的计算，以及乘方运算），对中间的乘积结果，不实施FFT逆变换、甚至不还原为longfloat形式，直接将点积结果以复数形式表示，下一次计算时直接使用。

FFT逆变换的时间复杂度为O(n*log2n)，但因为涉及除以N,以及进位处理，其耗时长于正向FFT。而FFT复数还原为longfloat的时间复杂度为O(n)。可见，如果能减少FFT逆变换次数，应该可以大大提高此类计算的速度，因为下次计算时的FFT正向变换也省略了。而省略还原过程也能在一定程度上提升效率，并节省中间过程的存储空间。但具体实现时仍需考虑如下几点：

  1. 程序中对于常规循环迭代的计算过程，已经充分考虑存储空间复用，减少释放和分配开销，故不会有空间效率；对于递归迭代的计算过程如P2B3，因中间结果均为临时动态的，因此减少中间存储空间的开销是有意义的，当然由于递归深度一般都是O(log2n）,这点提升可能从整体上也不明显
  2. 尽管FFT点积结果不实施逆变换等待重用听起来是个很振奋人心的性能提升思路，但这里有一个大的陷阱，就是FFT复数里double类型的精度限制。double类型在计算过程中会损失有效精度是毋庸置疑的，之所以能用它来进行大数乘法，是精确测算过，加上了FFT位长限制和Karatsuba分段计算，确保只进行1次FFT变换和点积的情况下，结果仍在double的有效精度范围内。但若直接对此结果再进行点积运算，则几乎肯定会带来精度损失，从而导致计算出的结果不准确。所以FFT点积结果必须逆变换后，重新进行FFT才能被再次使用，以保证精度。这个问题是暂时无法克服的。
  3. 既然FFT逆变换是不可避免的，缓存机制能解决的也就是节省转换为longfloat形式的时间开销，这点对整体性能提升应该很有限。如果要实现，则需要考虑如何对上层计算过程比较透明，不导致代码的过于复杂。基本思路是改造longfloat的内部结构，使之可以用整数、复数两种形式之一（或同时）存在，只在必要的时候才进行不同形式之间的转换。

2021.4.16补记：在最新版本中实现了对最近的乘数的FFT结果的缓存（对点积结果的FFT形式没有缓存，原因在上面已经分析了，没有实际意义）。实现机理与上面的思路基本相同：

  1. 在class fft中开辟缓存fftcache，用longfloat的地址addr和长度eb映射到FFT缓存项
  2. 执行FFT乘法时先查找缓存，没有找到再进行FFT变换，并将计算结果缓存起来
  3. 使用LRU机制，保留最近使用的若干个结果，超出时自动淘汰最老的缓存项
  4. 当longfloat::update()被调用时，即意味着即将对其中的数据进行修改，此时通知缓存管理器淘汰其FFT缓存即可，无需在其他地方逐一修改。

经测试，对于不同的算法，缓存确实能达到10%-50%的命中率，但是不幸的是，这对最终性能却几乎没有提升，因此这个优化机制是一个花瓶性质的东西。

经过跟踪和逻辑分析，原因其实和前面分析过的一样，因为运算中几乎不可能出现两个乘数都是完全和先前一样的情况（其中必有1个乘数是最近的运算产生的新数），因此即使每次均有1个乘数被命中，仍然需要对另一个乘数执行1次FFT变换。这样每次大数乘法执行1次FFT是必不可少的。而代码早已实现的交织耦合-解耦（fft::couple/decouple）机制，本来就能做到对2个乘数只做1次FFT，因此缓存起到的作用几乎可以忽略不计了。当然最新版本仍然保留了这个实现的代码，毕竟这也是一种经典的有意义的探索和尝试，如同前面说的递除法进行进制转换一样，虽然不用，但是作为整个算法体系的一部分存在；代码里默认将fft::caching开关置为false，不开启FFT缓存，以节省缓存查找、保存所额外的时间和空间开销。

2021年4月14日：

1、在调试新功能中发现数据或者内存异常，从而发现了之前版本隐藏较深的BUG，这些BUG在之前的版本运行中未暴露出来，实属意外，目前已经修正：

(1) 未及时处理longfloat的addr指针归位：addr指针为longfloat的第一个有效位地址，其在运算过程中可能向前或向后偏离初始位若干位置，以便快速处理去前导0、增加进位的情况（避免大块内存搬迁，直接移动首有效位地址即可）；之前在check函数中，会对addr积累的偏离进行修正，以使addr的偏移控制在很小的范围内；程序其他位置均未处理归位问题，但也未擅自改变addr指针位置，加上运算过程中因为不断调用check检查计算结果，因此未出现意外情况。

本次增加的进制转换在实现时，利用了直接修改addr的方式获得一个longfloat的小数部分（a.addr + a.exp + 1即为a的小数部分第一位），从而无需内存拷贝即可创建出1个新的longfloat数，这样做的好处是在运算中加快速度，但引入了addr与初始位置偏离大大，无法及时归位的隐患，此时若该浮点数为只读不写时不会产生致命问题，但若对该浮点数进行修改时则会发生异常。

修正方式有两种思路：一是提供显性归位函数，在各函数中对于要存放计算结果的longfloat，显式调用该函数进行归位，这种方式需要大量修改代码；第二种思路是修改update函数，当准备对一个longfloat进行修改、并且其先前存储数据不需要保留时，进行自动归位，这样修改量最小且有效。程序最终采用第二种思路，经测试不再有先前的异常

(2) longfloat::swap()中处理交换的2个对象时，如果对象本身引用其他对象（即采用浅拷贝赋值，ref指向另一个对象，与对方共享同一片存储区域，这种机制可以最大减少内存拷贝。例如：c = b = a * 1，我们只需要将b和c均指向a的同一片内存即可，不需要为b、c分配新空间并进行拷贝，直至后续对a、b、c其中之一进行修改时，才会执行真正的拷贝，并解除先前的引用关系），之前的处理是直接通过深拷贝解除原对象的内存共享，但这样会引起一次大的内存拷贝，降低效率。

改进机制为当存在对象引用时，直接调用unlink将原始对象的引用解除，并调relink将该引用转移到交换的对象上去，这样不会涉及内存拷贝。但是在该机制刚编码运行时，出现了异常，原因是当交换的2个对象存在彼此引用、或者引用同一个第三者时（即a.ref == b或b.ref == a或a.ref == b.ref），应做特殊处理，无需额外执行解除、转移操作。

(3)  longfloat::unsigned_sub()函数，执行c = a - b时，如果c与a或b为同一对象，则c原来的数据区域应该保留不能舍弃。函数中之前仅判断了c与 a相同的情况，对于c与b相同未做处理，之前运行中可能缺少这样的调用，BUG没有暴露出来；而Machin公式的P2B3实现中有这样的运算逻辑，计算结果发生异常。经修正后结果完全正常。

2、进制转换算法的实现和优

（1）先前版本已经实现了用小数乘法实现的2->10进制转换函数，其实现机理为：

H0.H1H2H3...Hn  * Ex ==> 0.B0B1B2...Bm  * Ry

这里E为2进制单元（一个32位字长）的基底，程序中默认为231；R为10进制单元的基底，程序中默认为109

即先通过乘以R的y次方将原来的浮点数H变换为小于1的浮点数B，之后将该浮点数B每次乘以R，得到的整数部分即为R进制形式的下一位：

0.B0B1B2...Bm *  R ==> C0.C1C2...Cn

C0即为转换结果的最左1位，然后去掉该整数位后，得到0.C1C2C3...Cn，再重复执行上述乘以R、取整操作，即可得到其余的位数。

递乘转换法也可称为小数部分转换法，原来用于对一个二进制数的小数部分进行转换，但经过上面的* Ry变换，实际上也可以针对带整数部分的任何整数/浮点数进行转换

（2）与递乘转换法对应的，也可以用递除的方式进行转换，其机理为：

H0.H1H2H3...Hn  * Ex ==> B0B1B2...Bm * Ry

这里的B为一个不带小数部分的大整数，然后不断用B去整除R，每次得到的余数（0 ~ R-1）即为从右到左的下一位数字：

B0B1B2...Bm / R = C0C1C2...Cm（商） 。。。D（余数）

D即为转换结果的“个位”，之后用商C0C1C2...Cm去整除R，得到下一位“十位”，以此类推，直至商为0为止。最后将结果逆序输出即可

（3）程序中实现了递乘、递除两种机制，以便进行算法对比和校验，分别为：convert10_MC、convert10_DC

（4）算法性能分析：在递乘转换法中，首先要执行一次* Ry的操作，这里包含几次大数乘法（含R的乘方），采用非FFT乘法时，为O(n2)，采用FFT后，为O(n*log2(n))；其后每计算1位，都要执行一次大数*R的乘法运算，每次的计算量为O(n)，则总的计算量为：

O(n2) + O(n*log2(n))，其中基本都是整数乘和浮点乘（FFT）运算。可见当n较大时，主要的性能开销是在后面的递乘上，当n=100万时，大概需要执行5000亿次整数乘法和1亿次浮点数乘法，后者的开销可忽略不计

递除转换法的性能与之类似，不同之处在于每计算1位，需要执行n次整数除法指令，因此除法计算量为O（n2）。虽然与递乘法都是相同量级的，但除法指令本身慢得多，所用时间约为递乘法的9~10倍左右，而且最后结果还需要一次逆序变换。因此程序默认使用递乘法代替递除法，可通过命令行开关进行控制（+/-cvm，使用递乘法/递除法）

（5）性能优化：O（n2）数量级的乘除运算，在n很大的时候耗时急剧增加。例如实测当n=100万时，使用递乘法完成一次转换约需要20分钟以上，使用递除法完成转换需要180分钟，这甚至已经大大超出了经过优化后整个求PI过程所花的时间。

对进制转换的优化思路是采用分而治之（Divide&Conquer）的思想，将大数切割成若干段，每段分别转换，最后将各段的转换结果进行合并。假设可将长度为n的原大数切割为4段，每段分别转换的结果进行简单拼接就能得到最终结果，那么切割前的计算量为：n2，切割后的计算量为4 * (n/4)2 = n2/4，可降低为原来的1/4，当然，这里需要加上切割、合并本身需要的额外计算量。因为转换前的大数是二进制基底，不能简单按照其长度进行分割，需要确保分隔后的各段之间维持10机制的幂关系，因此计算方式如下（以递乘法为例，递除法与之类似）：

第一步：仍是将大数转换为小于1的浮点数：

H0.H1H2H3...Hn  * Ex ==> 0.B0B1B2...Bm  * Ry

第二步：将0.B0B1B2...Bm * Rm ==> C0C1...Ck  . Ck+1Ck+2...Cn

结果包含整数部分C0...Ck和小数部分0.Ck+1...Cn，其中整数部分使用递乘法转换出的m位数字即为最终结果的最左m位数字，小数部分使用递乘法转换出的后面若干位数字。由于整数和小数部分的长度比n小很多，其计算量也大为减少；切割过程中使用的* Rm 包含几次大数乘法，可使用FFT使之仅有O（n*log2(n)）。这样总的计算量降低为：

O（n * log22(n)）

当n很大时，其较之先前的计算时间大大减少。

当然，由于FFT在n较大时才能发挥出其相对于简单大数相乘（long_mul）的性能优势，加之分段切割、组装会带来额外的代价，因此并非在所有情况下都把n分到最小的长度为止。应选择一个N值，所切割出的每段所包含的十进制位数不小于N，对这样的N位段直接使用递乘法转换，不再进行切割。

根据实际测试效果，在使用递乘法时，若采用的最小分段长度为10万~13万之间时，使用分割法比不分割有明显性能提升，小于这个长度时，分割计算所花的时间反而大于直接对原数进行递乘转换。以n=1000万为例，使用分割+递乘的转换时间仅略大于3分钟，大大小于先前的20分钟。

递除法也可以使用类似的分割机制优化，时间可以缩短为先前的

1/6 ~ 1/7左右，n=100万时，可从180分钟缩短为25~30分钟

（6）程序最终实现并保留了4种转换机制，可通过命令行开关进行控制，默认为分段递乘：

命令行

函数

机制

性能系数

-cvm -cvd

Convert10_DC_nodivide

直接递除

1

+cvm -cvd

Convert10_MC_nodivide

直接递乘

9

-cvm +cvd

Convert10_DC_divide

分段递除

7

+cvm +cvd

Convert10_MC_divide

分段递乘

60

3、修改了pulser的实现机制

pulser为心跳信号触发对象，供一个执行时间较长的过程内周期性向运行者提供进度指示。该耗时过程内部应能拆解为细分的、循环的操作，例如级数求和、牛顿迭代、进制转换等。一般在过程内每执行1次循环(称为1步)调用1次pulser::check()，如果返回为true，则程序可提供中间进度信息的显示，打印中间结果。

之前的pulser实现为采用步速控制机理，自动测算每次心跳之间的步数，决定下1次触发心跳的时机，并根据计算结果进行步速修正，确保两次心跳输出的时间间隔基本（非绝对）均匀。

修改后的机制，通过启动1个专门的计时器线程，简化每个pulser对象的判断逻辑，无需测算步速，只要上次计时器的值与本次之差大于某个间隔，就触发心跳输出。

计时器线程在创建第一个pulser对象时启动，不断计时，更新计时器当前值。当所有pulser对象消失后，其延时10秒，若此间再无新的pulser对象创建，则线程自动退出。整个程序结束时，建议调用pulser::terminate()通知计时器线程尽快退出。

2021年4月16日：

1、试验并实现了对FFT结果的缓存与查找，以期望减少执行FFT的开销，实现方法和效果在上面有详细补记说明，总之一句话，理想很丰满，现实很骨感。

2、发现了Borwein的4次迭代求PI公式（Method=32）可以通过公式变换，减少1次大数乘法运算：

原来的迭代公式为：

y₂ = (1 - (1 - y₁4)1/4) / (1 + (1 - y₁4)1/4)

a₂ = a₁*(1+y₂)4 - 22n+1*y₂(1+y₂+y₂²)

其中对于新的a值即a₂的计算，之前要做1个（1+y₂）的4次方、1个乘以a₁，1个y₂的平方，以及1个乘以y₂，共需执行5次大数乘法（4次方至少需要执行2次大数乘法，乘以22n+1用二进制移位法快速完成），现在可以等价变换如下：

a₂ = a₁*((1+y₂)²)² - 22n+1*y₂((1+y₂)² - y₂)

可以看出，减号左边的（1+y₂）²计算结果可以暂存下来（记为x₂）被右边括号内重复使用，这样最后只需要做4次大数乘法即可：x₂ = （1+y₂）²，x₂²，a1*x₂²，y₂*(x₂ - y₂)

经过该优化后，性能提高6%左右，计算100万位耗时从33秒降为31秒。

3、发现并修正了加法/减法函数unsigned_add()/unsigned_sub()中一个隐藏的BUG，当函数调用为如下形式：a = a + b、a = a - b，即计算结果放回其中的一个操作数时，程序之前做了优化处理：当b的指数小于a（例如1 + 0.000006），那么从右到左加（减）至b的最高有效位（即6）时，且未产生进位（借位），则结束循环，a的其余位无需计算维持原值。但这里少了个必要的操作，就是把该位到a的最右有效位之间的位置填补0，避免这段区域存放之前的垃圾数据，否则会导致计算结果出现随机偏差。

该问题在运行Borwein-9次方迭代（Method=30）的过程中发现，当n大于200万位时，牛顿迭代法求3次方根倒数时出现不收敛的现象，最后计算出的数值大大偏离正确值。

2021年4月18日

1、发现并修正了程序中存在的个别缺陷和隐患

（1）在执行大数乘法时，即使不开启FFT，也应该使用Karatsuba机制进行分段乘，这样每分割1次，可将乘法运算量降为先前的3/4，直至分到最小长度1024时才执行基本乘法long_mul。当n为100万时，乘法计算量从1万1千亿次降为620亿次，性能约提高20倍。该问题的发现和修改来源于下面第（2）个问题发现

（2）在用Borwein的9次方迭代公式（Method = 30）计算PI的200万位以上时，发现计算结果出现误差，此时其他公式均结果正常。起初怀疑是这个公式自身的实现逻辑问题，追踪发现是其迭代到第4次时倒数公式不收敛、迭代第8次立方根倒数公式不收敛，从而出现计算偏差。由于这两个公式均是采用牛顿迭代法实现的，其中主要运算为大数乘法，因此初步判断又是FFT误差所致，但要证明这点，必须证明在不开启FFT的情况下，公式能正常收敛。因此时位数已经很大，不开启FFT耗时很大，而之前程序中未开启FFT时直接使用long_mul进行大数乘法，执行1次乘法就需要20分钟，几乎无法做验证。这才有了上面（1）中对基本乘法的优化。改进后虽然仍比开启FFT慢得多，但执行时间亦可忍受，也证明了不开启FFT确实能正确收敛。

对FFT误差的唯一修正方式，就是减小乘数长度或者FFT基数，最终采用将乘数最大长度减半的方法，在wise_mul中从判别a、b的最大长度不超过阈值，改为a、b的长度之和不超过阈值。修改后执行结果即完全正常了。

2、发现Borwein的9次方迭代求PI公式（Method=30）可以通过公式变换，减少1次大数乘法运算：

原公式中v的迭代公式如下：

v₂ = t₂² + t₂*u₂ + u₂²

其中涉及到3次大数乘法和2次加减法，可以等价变换为：

v₂ = （t₂ + u₂）² - t₂*u₂

只需要执行2次大数乘法和2次加减法

由于Borwein 9次方迭代公式远比其他公式复杂，该优化对整体性能提升幅度很小，在n=100万时，可将整体计算时间从63秒提升到61秒。

3、对所有使用refiner进行迭代精度控制的函数进行微调，从之前在每轮迭代结束时检查精度是否达标，前置到迭代过程中的某一步，发现精度达标即退出循环，这样可跳过后面无需执行的步骤。

进行该优化，是因为根据运行输出中间结果，发现许多情况下最后1次迭代结果与前1次完全或者几乎一样，这表明其实在上一轮迭代时实际已经达到了所需精度，只是按照上次的判断逻辑认为还需要下一轮迭代。

前置的最佳位置需要根据每个公式的特点设定，举例说明：

Gauss-Legendre迭代公式：

a₂ = (a₁+b₁)/2

b₂ = √(a₁b₁)

t₂ = t₁ - p₁(a₁-a₂)²

p₂ = 2p₁

当a₂和a₁的误差小于精度要求时，即可满足要求而退出。之前是在最后的步骤：p₂ = 2p₁执行完后才比较a₂和a₁，观察这个公式可以发现，其实在第一步a₂ = (a₁+b₁)/2之前，先比较a₁和b₁即可，如果两者在精度范围内近似相等，则肯定后面计算出的a₂（即a₁和b₁的算术平均数）也与a₁相等，也就意味着已经达到精度要求无需继续迭代了，这可以省去后面几步执行大数乘法和开方所需消耗的大量运算。

同理，我们再看看更为复杂的Borwein 9次方迭代公式：

t₂ = 1+2r₁

u₂ = (9r₁*(1+r₁+r₁²))^(1/3)

v₂ = t₂² + t₂*u₂ + u₂²

w₂ = 27(1+s₁+s₁²)/v₂

a₂ = w₂*a₁ + 3^(2n-3)*(1-w₂)

s₂ = (1-r₁)³/((t₂ + 2u₂)*v₂)

r₂ = (1-s₂³)^(1/3)

可以看到a₂和a₁的迭代关系，当w₂近似为1的时候，a₂也近似等于a₁，因此我们可以在这步之前检查w₂和1的偏差，如果近似相等则退出循环，无需再进行后面的3步（这3步涉及大量的大数乘法、乘方、开方运算）

增加了8种求PI公式，其Method编号分别为12、13、14、15、16、33、34、35。

其中前5种（12~16）为类BBP公式，即分数级数型，其收敛系数分别为256-1、256-1、4096-1、4096-1、4096-1、65536-1。这些公式均实现了累加求和和P2B3两种方式，后者的速度快得多。注意并非收敛快的公式最终运算速度最高，因为这些公式里往往包含了更多的分项，需要分别计算各分项再求和，所以整体运算速度可能更慢。这点和Machin公式类似（见文档结尾的性能对比表）。

后3种（33~35）均为Borwein迭代公式，使用牛顿迭代法实现，分别为2次收敛、3次收敛、5次收敛公式。类似的，并非迭代收敛速度越大的执行速度越高，例如5次方收敛的每次迭代复杂度比3次方的高很多，因此其执行速度反而更慢。目前为止，Borwein公式中执行速度最快的仍是Method=31的2次迭代和Method=32的4次迭代。

2021年4月22日

1、本轮改进的核心：通过并发来提升性能，以达到千万位以上的计算目标。从几个方面进行了尝试，喜忧参半，不过也在这些尝试摸索中对程序和算法有了更深的认识，最终找到了一线曙光，不过这线曙光在下一版本中才会体现出来。这里先总结一下本轮所做尝试：

（1）SSE2指令：出发点和之前一样，用SSE/AVX优化FFT中的复数乘法运算。为免像上轮FFT缓存那样最终事倍功半，因此还是照例先写测试程序，论证其可行性。

这里没有直接使用汇编指令，而是调用gcc封装的接口<xmmintrin.h>，使用__m128d类型和_mm_mul_pd函数，实现同时对两个双精度浮点数的乘法。考虑到复数乘法如下：

c = a * b = ( a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r )

其中r、i分别代表复数的实部和虚部

先将a、b加载到__m128d变量a1、b1（_mm_load_pd），并将b逆序（交换虚实部次序）加载到到b2（_mm_loadr_pd）：

a1 = [ a.r, a.i ]

b1 = [ b.r, b.i ]

b2 = [ b.i, b.r ]

然后调用2次128位双精度浮点乘（_mm_mul_pd），得到两个128位结果：

c1 = a1 * b1 = [ a.r * b.r, a.i * b.i ]

c2 = a2 * b2 = [ a.r * b.i, a.i * b.r ]

然后将c1、c2分别拆解为double变量（_mm_store_pd）后进行加减运算得到最终结果（本来如能直接调用SSE3的扩展函数_mm_addsub_pd，对c1、c2的高、低部分直接进行加减运算更理想，但遗憾的是测试平台上不支持SSE3）

这个过程也就是把之前的4次double乘，变成了2次128位乘，但因为在执行128位乘法前后的加载、拆解动作，能否有实质性能提升需要实际验证。测试程序分别用传统方式和SSE方式，对同样的两组复数（每组1万个复数）重复计算乘积（1万次），结果非常失望。使用SSE2不仅没有提升，而且耗时比传统方式多1倍。

基于这个测试结果，无法将其引入到最终的代码中，故放弃。

在目前的运行平台上，扩展的SSE3、SSE4、AVX均无法使用，无法判断能否进一步使用这些扩展指令提升性能。

(2) 使用多线程优化大数乘法运算

这个思路是针对所有求PI算法的核心基本运算—大数乘法，用并行计算的思路让其加速，这样可让所有中层函数（如求倒数、平方、阶乘、平方根及n次方根等）和上层算法（各种求PI公式）无需修改都能享受到性能加速。

之前的版本中实现的Karatsuba乘法（wise_div）是最适合使用这种机制的地方。当乘数a、b的长度很大（超过100万），以至于超出FFT精度限制时，需要利用二分法将a、b拆解为较短的长度进行3次FFT乘，然后对结果进行加减合并。如果这3次FFT乘可以3路并行计算，理论上则可将计算性能提高到2~3倍之间。进一步的，如果拆解出的更小片段能由更多的线程并发计算，则性能还可能再提高。

程序实现的基本思路是预派生+任务队列机制：

  1. 使用任务队列（xqueue），存储需要计算的乘法子任务（xtask）
  2. 程序启动时即创建好若干并发线程（xthread），每个线程监听队列中是否有需要计算的子任务，有则取出（peek）进行计算（exec）
  3. 主线程为大数乘法的入口，由其在wise_mul执行第一次拆解，拆解为3个子任务，放入队列（post）
  4. 后续负责计算每个子任务的线程，根据需要可能将该子任务拆解为3个更小的子任务，加入队列
  5. 拆解的线程负责监视并等待3个子任务均计算完毕（wait），然后对3个子任务计算结果进行合并，得到当前任务的计算结果
  6. 为提高CPU利用效率，减少空转，当前线程在等待其子任务计算结果出来的过程中，如果队列中尚有其他子任务（也包括自身拆解出的子任务），则其也将分身去执行1个或多个子任务
  7. 多个线程在存取任务队列时，为防止冲突，使用自己实现的spin-lock轻量锁来进行互斥访问。这里的spin-lock调用gcc内置的atomic函数__sync_bool_compare_and_swap来实现，而不是传统的mutex，不会引起当前线程被调度后切换上下文；同样，在获取锁的重试过程中，调用Intel建议的_mm_pause随机跳过几个指令周期，而不是调用yield()让出CPU，以便更紧凑的执行任务（当然这样类似空循环的操作，会引起CPU占用率较高）

程序设计实现时，在上述基本思路之上，增加了额外的考虑：

  1. 考虑到除了大数乘法运算，可能还有其他可以进行类似分解并发优化的地方，例如进制转换、P2B3模型，代码设计的时候预留了扩展性，可以在xthread基类上派生出执行其他类型任务的计算子类，只需要实现少量的计算和封装函数即可。线程本身的执行调度无需再重写
  2. 实现调试过程中发现了目前的longfloat并非所有成员函数都是可重入的，在多线程并发执行场景下，如果不同线程访问同一个longfloat对象并存在修改可能性，则会导致意想不到的错误。主要的地方是通过dup/swap完成的软拷贝和内存空间重复利用，这导致两个甚至更多个longfloat对象之间彼此关联。解决方法是引入全局的互斥锁（用上面描述的spin-lock实现机制），在涉及到对象之间引用时确保互斥执行。
  3. 考虑到拆解子任务以及队列管理本身也会带来额外开销，线程越多、拆分次数越多，这些额外开销越大。因此一味增加线程数量并不必然导致整体性能增强，甚至可能反而恶化，实验结果也证明了这一点。因此在大数乘法的并发实现时，向队列中加入子任务时（post）引入了控制机制，当子任务的层级（顶层任务层级为0）不为0时，不再追加队列，而是直接执行；同时，当前负责任务拆解的线程，其第一个子任务不追加队列，也是直接执行；并行数量最优为3（含入口主线程）

经过上述的调优，最终效果如下：

  1. 大数位长小于200万个十进制位时，大数乘并发优化效果不明显；位数越大，效果越显著
  2. 位长500万的单次大数乘法，在3路并发时，性能提升90%
  3. 根据不同的PI求解方法，性能提升幅度不同。在利用Borwein 4次迭代（Method=32）求解200万位时，使用3路线程并发（含主线程）的情况下，求解时间从120秒提升到78秒，性能提升幅度约50%；用 Gauss-Legendre迭代公式求解500万位，3路并行乘法将计算耗时从876秒降至570秒，性能提升幅度也是50%；而对于最快的Chudnovsky P2B3计算公式，其计算1000万位时间从180秒降为156秒，提升幅度仅15%。
  4. 总结：大数乘法3线程并行优化比单线程快50%，之所以无法达到原来的2~3倍，是因为1）拆分和并发管理的额外开销；2）线程之间存在数据关联性；3）主线程整体上仍然承担着绝对的计算量，因为大数乘法外的其他运算、以及位长不足以拆分时的乘法运算（百万位以下的乘法运算）均由主线程独立完成，即使在可并行的大数乘法过程中，主线程也承担比其他2个线程略高的计算量。

3  下一步优化的思路基本成型，在更高层级的计算函数（而非基本运算函数）上，拆解为多个可连续执行到底、彼此无干扰的子任务，每个子任务中不仅包含大数乘法，也包含各种其他运算，这样并发效率可大幅提升。仍用上述大数乘法的类似机制进行并发设计和实现。这里考虑了2个应用点：

(1) P2B3模型的并发优化：使用该模型的求PI公式最多，而且其二分递归模型比较适合多路并行

(2) 对不使用P2B3模型、而是逐项累加的，也可以采用类似机制，将其按累加项分割为多个连续子段，每个子段并发各自计算求和，最后进行合并

(3) 牛顿迭代法实现的公式，暂时无法使用上层并发优化手段，只能借助底层的并发大数乘进行加速

(4) 考虑将先前二分递归实现的过程（大数分段乘法、P2B3）改造为非递归模式，这样可一开始就完成所有最细粒度子任务拆解和分配，并将需要的临时变量空间准备就绪，减小在递归调用中频繁队列出入以及空间分配释放的额外开销。目前进制转换就是非递归实现的，可以参考其进行改造

高精度高性能PI值计算程序设计和验证相关推荐

  1. 虹科方案 | 自带灭菌值计算功能—FH值验证干热灭菌强度

    一.干热灭菌 对于能够耐受高温且不容易被蒸汽穿透或者容易被湿热破坏的材料,通常选用干热灭菌法灭菌.比如粉末.油品.油脂.玻璃制品和不锈钢设备的灭菌.由于干热灭菌对微生物的杀灭率低于同一温度下的饱和蒸汽 ...

  2. java中计算器算cos值,Android开发中计算器的sin、cos及tan值计算问题分析

    本文实例讲述了Android开发中计算器的sin.cos及tan值计算问题.分享给大家供大家参考,具体如下: 接到一个需求 :要求计算器sin90=1,拿到知道很疑问 难道不等于一么?测试了四五个手机 ...

  3. android动画sin cos,Android开发中计算器的sin、cos及tan值计算问题分析

    本文实例讲述了Android开发中计算器的sin.cos及tan值计算问题.分享给大家供大家参考,具体如下: 接到一个需求 :要求计算器sin90=1,拿到知道很疑问 难道不等于一么?测试了四五个手机 ...

  4. python顺序结构实验_Python程序设计实验报告二：顺序结构程序设计（验证性实验）...

    安徽工程大学 Python程序设计 实验报告 班级 物流191 姓名姚彩琴学号3190505129 成绩 日期 2020.3.3 指导老师修宇 [实验名称] 实验二 顺序结构程序设计(验证性实验) [ ...

  5. python顺序结构实验设计_实验二 顺序结构程序设计（验证性实验）

    安徽工程大学 Python程序设计实验报告 班级物流192 姓名 徐敏 学号 3190505232 成绩 _____ 日期 2020.3.22 指导老师 修宇 [实验名称] 实验二 顺序结构程序设计( ...

  6. 【数据结构】哈希表、哈希值计算分析

    哈希表.哈希值计算分析 哈希表完整代码 引出哈希表 哈希表(Hash Table) 哈希冲突(Hash Collision) JDK1.8的哈希冲突解决方案 哈希函数 如何生成 key 的哈希值 In ...

  7. adxl345取出值怎么算角度_怎么通过adxl345输出值计算出倾角？

    怎么通过adxl345输出值计算出倾角? 传感器常见问题 / 2012-12-23 此方法是通过勾股定理计算的,请参考以下程序: // 加速的X轴用来算俯仰角;Y轴算横滚角 u16 Gx; u16 G ...

  8. R语言用户自定义函数的语法结构、编写自定义统计值计算函数（使用ifelse结构计算均值和标准差等）、编写自定义日期格式化（format）函数（switch函数使用不同分枝格式化日期数据）、应用自定函数

    R语言用户自定义函数的语法结构.编写自定义统计值计算函数(使用ifelse结构计算均值和标准差等).编写自定义日期格式化(format)函数(switch函数使用不同分枝格式化日期数据).应用自定函数 ...

  9. 项目管理之项目的挣值计算问题

    本文将通过详细讲解历年的四道挣值计算真题来帮助掌握项目管理中的挣值计算问题. 项目管理专业术语与计算公式汇总 真题一:挣值计算(2010.11) 真题二:进度与挣值计算(2012.05) 真题三:进度 ...

最新文章

  1. SQL语句中的select高级用法
  2. 回归分析中的正则化问题
  3. 洛谷 P1340 兽径管理
  4. LVS负载均衡（3）——LVS工作模式与工作原理
  5. Servlet 文件上传
  6. 日常问题——hadoop 任务运行到running job就卡住了 INFO mapreduce.Job: Running job: job_1595222530661_0003
  7. 幻想乡三连A：五颜六色的幻想乡
  8. CSharpGL(53)漫反射辐照度
  9. 基于Python的面部表情识别分析系统
  10. 光谱数据计算色彩指标的软件(功能强大，齐全)
  11. Win11磁盘清理在哪打开？
  12. 王厚祥谈《古诗四帖》基本笔画的书写方法
  13. 数据库系统概论笔记二——画E-R图
  14. 音视频基础知识---音频编码格式
  15. 2022国庆头像制作iAPP安卓源码+附APP成品
  16. 解决 PR 或 AE 启动不了桌面弹出 Crash 文件
  17. 英语的各种 n. adj. vt. vi. 等词性解释
  18. 连接QuickBooks Online实现于IOS App数据同步功能的个人记录
  19. 中英文互译之Excel表格
  20. 国内CDN行业优质服务商

热门文章

  1. 【k8s】Unable to restart cluster, will reset it: apiserver healthz异常
  2. AI绘画工具软件网站合集：这些人工智能绘画生成器效果太赞了
  3. 【Apache Spark 】第 10 章使用 MLlib 进行机器学习
  4. MouseManager
  5. pyqt5中QGraphicsView弹出菜单
  6. gRPC快速入门（三）——Protobuf应用示例
  7. EDA软件—Cadence学习笔记分享（内含安装教程）
  8. 无人驾驶感知篇之融合（一）
  9. matlab repeat until,汇编语言用.REPEAT和.WHILE伪指令实现循环
  10. 运筹学基础【十】 之 盈亏分析模型