NOJ 高精度计算π值

• 针对来者的不同需求，本文将遵循以下的条目给出不同层次的方案（文中出现的“面向**编程”用词仅出于调侃，切勿当真）。
• 声明：本文与笔者皆从编码经验入手，不包含任何主观上的OI成分，若有OIer或ACMer来访的话就当是过眼云烟吧。

Requirement

基于双向链表的存储结构计算并输出指定位数的$\pi$值

Catalogue

1. 面向AC编程（也就图一乐）
2. 面向用例编程（紧急AC入口）
3. 面向实践编程（正文）
4. 面向提高编程（选读）
   a) 结构优化
   b) 内存优化
   c) 性能优化

面向AC编程

---

概述：能跑的BUG都是好程序，能AC的代码都是好答案。——佚名
关键词：打表

NOJ的测试数据贫乏得可怜，如果你只是想在你的NOJ页面上看到“AC”的话，不妨百度/必应/谷歌一下在线PI值计算，并ctrl+c/v一千位写入你的源代码中。下面给出具体代码：

#include <stdio.h>
static const char *PI = "/*将此处替换为你的π值*/";
int main() {
    int n = 0;
    char format_buf[16] = { 0 };
    scanf("%d", &n);
    /**
     * @note 写入格式化字符串，输出字符串长度n+2位，包含“3.”与小数部分的n位。
     */
    sprintf(format_buf, "%%.%ds", n + 2);
    /**
     * @note 当然你简单粗暴地直接输出500位也没有问题，毕竟面向AC编程（bushi）
     */
    printf(format_buf, PI);
    return 0;
}

面向用例编程

---

概述：只要你题干敢蕴含这个条件，我就敢忽略对该条件的断言并且在这个前提下对代码特化。
关键词：特化

问题分解

1. 期望输出是什么？
2. 如何在理论上获取期望输出？
3. 如何存储待输出目标？
4. 如何输出目标？
5. 如何计算目标？
6. 有哪些可用条件？

问题分析

Q：期望输出是什么？
A：一个$n$位小数精度的$\pi$值。
Q：如何在理论上获取期望输出？
A：期望输出为有限位$\pi$，故实际问题为如何计算$\pi$值；存在多种迭代、级数和计算方法，本文考虑从下式计算$\pi$值：
$$\frac{\pi }{2}=\sum _{k=0}^{\infty }\frac{k!}{(2k+1)!!},wheren!!=\prod _{k=0}^{\lceil \frac{n}{2}\rceil -1}(n-2k)$$
Q：如何存储待输出目标？
A：根据需求，确定使用双向链表存$\pi$值；若令链表节点的存储数据为$t$位十进制位，且$\pi$值对应的十进制数从高位到低位在链表中有序存储，则对应精度的$\pi$值表示则为一个$\lceil \frac{n+1}{t}\rceil$个节点的双向链表；其中，小数点的位置必须被明确定义。
Q：如何输出目标？
A：顺序输出节点数值，遇到小数点的位置就输出小数点。注意$t$位补齐以及整数部分的前缀零剔除（若整数部分为孤立零则需要保留）和小数部分的后缀零剔除。
Q：如何计算目标？
A：已经有数学表达式，就可以将问题归结为如何用代码表示出该式子。首先对计算式进行简易的转化得到：
$$\{\begin{array}{l}{R}_0=2\\ \\ R_n=\frac{n}{2n+1}{R}_{n-1},n\ge 1\\ \\ \pi ={\sum }_{k=0}^{\infty }{R}_k\end{array}$$根据上式，$\pi$的计算过程即可以转化为有限次$R_k$迭代和$\pi$的求和，其中$R_k$的计算涉及大数的乘法以及除法，$\pi$的求和涉及大数的加法，故只需实现大数初始化、加法、乘法、除法四项操作。
Q：有哪些可用条件？
A：(1) 输出的结果$\pi$是个定值，而整数部分只有$3$一位数字，故可以选择只存储并计算小数部分，而实现该操作只需稍微调整计算公式的初值，使计算过程中产生的值始终为小数。(2) 最终输出的$\pi$值位数是确定的，也即可以舍弃链表的动态插删功能，选择一次性构造$n$个节点，让双向链表的流于形式。(3) 选择的公式每迭代$2\sim 5$次可以产生一个有效位，而式子中最大的操作数为$2n+1$，也即$int$类型便能够支持目标迭代$4$万次（根据乘法操作不溢出的情形粗略判断）(4) 所有出现的过程数皆为正数，故不需要考虑数的符号。

下面给出具体代码：

#include <stdio.h>
#include <malloc.h>
#include <stddef.h>

#define N 4     //! 每个节点存储四位十进制
#define R 10000 //! 模数10^4

typedef unsigned int int_t;

typedef struct node_s {
    int_t value;
    struct node_s *prev, *next;
} node_t;

typedef struct number_s {
    node_t *head, *tail;
} number_t;

int main() {
    int n = 0, t, k;
    scanf("%d", &n);

    /**
     * @note 令pi的初值为R0+R1+R2，r的初值为R2，
     *  则迭代计算步骤为：
     *  r = r * k
     *  r = r / (2 * k + 1)
     *  pi = pi + r
     */
    number_t pi, r;
    node_t *p, *q;

    /**
     * @brief 初始化链表
     * @note 创建ceil(n/N)+1个节点，
     *  多申请一个存储节点防止最后节点的精度丢失
     */
    t = (n + N - 3) / N + 1;
    p = pi.head = (node_t*)malloc(sizeof(node_t));
    q = r.head  = (node_t*)malloc(sizeof(node_t));
    p->value = 9333;
    q->value = 2666;
    p->prev = NULL;
    q->prev = NULL;
    while (--t > 0) {
        p->next = (node_t*)malloc(sizeof(node_t));
        q->next = (node_t*)malloc(sizeof(node_t));
        p->next->prev = p;
        q->next->prev = q;
        p = p->next;
        q = q->next;
        p->value = 3333;
        q->value = 6666;
    }
    p->next = q->next = NULL;
    pi.tail = p;
    r.tail  = q;

    /**
     * @brief π值计算
     * @note 固定迭代2500次，
     * 	需依据所需精度适当调整
     */
    t = 2;
    while (++t < 2500) {
        int_t carrier   = 0; //! 进位数
        int_t remainder = 0; //! 余数

        //! r = r * k
        k = t;
        q = r.tail;
        while (q) {
            q->value = carrier + q->value * k;
            carrier = q->value / R;
            q->value %= R;
            q = q->prev;
        }

        //! r = r / (2 * k + 1)
        k = 2 * t + 1;
        q = r.head;
        while (q) {
            q->value += remainder * R;
            remainder = q->value % k;
            q->value /= k;
            q = q->next;
        }

        //! pi = pi + r
        carrier = 0;
        p = pi.tail;
        q = r.tail;
        while (p) {
            p->value += q->value + carrier;
            carrier = p->value / R;
            p->value %= R;
            p = p->prev;
            q = q->prev;
        }
    }

    /**
     * @brief 输出π值
     */
    putchar('3');
    putchar('.');
    p = pi.head;
    t = (n - 1) / N;
    while (t-- > 0) {
        printf("%04d", p->value);
        p = p->next;
    }
    t = N - (n - 1) % N - 1;
    while (t-- > 0) {
        p->value /= 10;
    }
    printf("%d", p->value);

    /**
     * @brief 释放内存
     */
    p = pi.head;
    q = r.head;
    while (p->next) {
        p = p->next;
        q = q->next;
        free(p->prev);
        free(q->prev);
    }
    free(p);
    free(q);

    return 0;
}

关于上述代码的补充解释

• 关于初始化操作：在问题分析中已经确定只保存小数部分，则计算过程产生的任意值都应该满足有效值为小数范围。若令r的初值为$R_0$，则在乘的过程中将产生大于$1$的数。而$R_2$是第一个在第一步乘操作不会进位到整数部分的项，故可以令r的迭代从该项开始。同时，应该使pi的初值为 { R k ∣ k ∈ N } \{R_k|k\in{N}\} {Rk​∣k∈N}前三项之和的小数部分。
• 关于乘除顺序问题：因为有限精度的表示方法，不论以何种次序计算，结果值总会与理论真实值存在误差。然而，在保证乘法操作不产生大于等于$1$的数后，除法操作总是能得到期望真实值。故而最终期望真实值的误差产生于除法期望真实值与理论真实值的误差对乘法操作的干扰（我的天啊，这段话简直说得乱七八糟，不过你大概知道是什么意思就可以（其实不知道也没关系因为这不重要））。
  ① 记 x 为 上 一 轮 迭 代 除 法 操 作 的 真 实 值 ， T 为 任 意 正 整 数 真 实 值 ： R = T ( [ x ] + { x } ) 2 T + 1 计 算 值 ： R ′ = T [ x ] 2 T + 1 ε 1 = [ R ] − [ R ′ ] = [ [ T [ x ] 2 T + 1 ] + { T [ x ] 2 T + 1 } + T { x } 2 T + 1 ] − [ T [ x ] 2 T + 1 ] = [ { T [ x ] 2 T + 1 } + T { x } 2 T + 1 ] ⇒ ε 1 ∈ { 0 , 1 } ② 记 y 为 当 前 轮 的 输 入 值 （ 有 y = [ y ] 成 立 ） ， T 为 任 意 正 整 数 真 实 值 ： R = T y 2 T + 1 计 算 值 ： R ′ = T ⋅ [ y 2 T + 1 ] ε 2 = [ R ] − [ R ′ ] = [ T ⋅ ( [ y 2 T + 1 ] + { y 2 T + 1 } ) ] − [ T ⋅ [ y 2 T + 1 ] ] = [ T ⋅ { y 2 T + 1 } ] ⇒ ε 2 ∈ { x ∣ x ∈ N ∧ x ≤ T } \begin{array}{l} ①记x为上一轮迭代除法操作的真实值，T为任意正整数 \\ \\ 真实值：R=\frac{T([x]+\{x\})}{2T+1} \\ \\ 计算值：R'=\frac{T[x]}{2T+1} \\ \\ \varepsilon_1=[R]-[R'] =[[\frac{T[x]}{2T+1}]+\{\frac{T[x]}{2T+1}\}+\frac{T\{x\}}{2T+1}]-[\frac{T[x]}{2T+1}] =[\{\frac{T[x]}{2T+1}\}+\frac{T\{x\}}{2T+1}] \\ \\ \Rightarrow\varepsilon_1\in\{0,1\} \\ \\ ②记y为当前轮的输入值（有y=[y]成立），T为任意正整数 \\ \\ 真实值：R=\frac{Ty}{2T+1} \\ \\ 计算值：R'=T\cdot[\frac{y}{2T+1}] \\ \\ \varepsilon_2=[R]-[R'] =[T\cdot([\frac{y}{2T+1}]+\{\frac{y}{2T+1}\})]-[T\cdot[\frac{y}{2T+1}]] =[T\cdot\{\frac{y}{2T+1}\}] \\ \\ \Rightarrow\varepsilon_2\in\{x|x\in{N}\land{x}\leq{T}\} \end{array} ①记x为上一轮迭代除法操作的真实值，T为任意正整数真实值：R=2T+1T([x]+{x})​计算值：R′=2T+1T[x]​ε1​=[R]−[R′]=[[2T+1T[x]​]+{2T+1T[x]​}+2T+1T{x}​]−[2T+1T[x]​]=[{2T+1T[x]​}+2T+1T{x}​]⇒ε1​∈{0,1}②记y为当前轮的输入值（有y=[y]成立），T为任意正整数真实值：R=2T+1Ty​计算值：R′=T⋅[2T+1y​]ε2​=[R]−[R′]=[T⋅([2T+1y​]+{2T+1y​})]−[T⋅[2T+1y​]]=[T⋅{2T+1y​}]⇒ε2​∈{x∣x∈N∧x≤T}​
  其中①②分别对应先乘后除、先除后乘的误差分析，显然法①带来的误差更小。（证不出来，谁能算摊销误差谁算去吧，我直接摆烂了；别问为什么①更好，问就是显然+实践经验，①的误差基本维持在一个节点之内，而②当$T$增大时误差直接爆炸）
• 为什么不使用别的收敛更快的级数或高斯-勒让德这类迭代算法？

1. 在本节“面向用例编程”中，你实现的是最拉胯的大数支持，它甚至连减法都不支持！在如此简陋的大数支持外加必须使用双链表的约束下，用高斯-勒让德算法迭代实在是太难为人了。如果你硬要这么做的话，那么请不要在此节逗留，您需要的是“面向SCORE编程”而不是混个AC。
2. 更快的收敛速度需要付出使用复杂表达式的代价，那么问题来了，你能确保你提供的大数-$int$乘除操作能够撑到收敛到给定精度的时候而不发生溢出吗？或者说，你能确保在给定的操作数及操作符不会使这固定节点数的“大数”莫名其妙损失大量精度吗？
3. 本节/本文选用的级数收敛虽然慢是慢了点，但是它更简洁、更稳，这意味着你跟容易和它玩到一块而不是抗起那天杀的调试器来打个天昏地暗，拼个你死我活——有可能世界都毁灭了，你还在那悲催地调试并思考“这段为什么会错”以及“这段凭什么能跑”之类的问题。

面向实践编程

---

概述：相比于高度特异化的编码方式，更普遍的选择是使用相对标准的一套工具，进行细微的改动使其能够适配于目标问题。此节所谓“面向实践”，即指在遵循“各司其职”的规范下对个功能模块进行独立，封装并提高代码的可复用性。
关键词：一般化实现

·PS·由于源代码的长度与复杂度的原因，不对其全部或大块地展示。阅读全部代码请点击[此处]

问题引入

在“面向用例编程”中，双链表被抛弃了动态性，变成了功能与数组相差无几但是性能却落后于数组的一个摸鱼结构。另一点是，对于上文的代码，哪怕目标问题发生极其细微的改动，都可能引起其贯穿始终的细节变更，复用性差，调试调优也较为困难。对于实际生产，这样的现象显然是yue~的，于是引入健康的编码风格、代码框架便成了燃眉之急（夸张了，只要够自信，随你怎么写，甚至反向缩进都行（假的，这种人建议刀了））。

问题罗列

• 名不副实的双链表
• 又臭又长的初始化操作
• 被丑陋地暴露出来的四则运算细节
• 难以使用的输出代码
• 弱自适应的迭代终点
• 不（必然）安全的常量定义
• 大量无效的空运算

解决办法

【Pr.1: 双链表应用】

链表最大的优点：$O(1)$的插删操作时间复杂度
废弃上文使用固定节点数的做法，当数需要被拓展时插入新的节点，当无效节点（非有效数字的前缀、后缀零）出现时删除该节点。
链表的插删操作将发生在所有对目标数据发生写的方法中，故而所有存在改写节点信息的方法都应该在内部完成对双链表的维护。
附言：都决定实现个正经的双链表了，那就别再持有“只需要保存小数部分”这个执念了！同时支持整数部分和小数部分其实并非困难的事情，只要能够定下一个定点（运算的基准点），那么不论是只保存小数部分还是保存全部，它们的本质其实都是一样的：向定点两段拓展的整数序列。显然，将小数点往左的第一个整数单元作为这个定点是不错的选择。

以下是一个可供参考的双链表结构：

typedef struct node_s {
	IntType value;
	struct node_s *prev, *next;
} node_t;
typedef struct number_s {
	note_t *base; //! 整数部分低位第一个单元
	node_t *head; //! 最高位节点
	node_t *tail; //! 最低位节点
} number_t;

【Pr.2: 初始化】

初始化一个大数就得手动写一个没什么技术含量的循环这真的是很令人烦躁的事。当然要求初始化能够实现无理数赋值，这个条件又显得有些苛刻。初始化至少是易于使用的，对于一个大数类，至少得存在一种初始化方式使其能够被赋为任意精度允许范围内的值。
提供一个字符串作为参数来初始化这是很人性化的操作，而提供低成本的整数序列也是不错的选择。
至于提及的无理数初始化，这其实是一个伪命题。由于使用有限精度，无理数其实也就退化成了一个有限长度的小数，这时候依旧可以使用上述的初始化方式。但不可否认，当该长度足够长时，初始化为无理数总是一个耗费时间与精力的事情。
而初始化一个有理数——即便它是无限长度的——就没有那么复杂。若已知一个有理数的分数形式或有限步四则运算表达式，我相信你更愿意计算出它的初值而不是手动赋给它——这其中造成的性能损失相较于节省下来的劳力损耗往往是可以接受的。

以下是可供参考的初始化形式：

/**
 * @brief 从字符串创建大数对象
 * @param precision: 最大小数位精度
 * @param num: 大数字符串
 */
number_t* create(int precision, const char *num);

/**
 * @brief 根据字符串构造大数对象
 * @param precision: 最大小数位精度
 * @param num: 大数字符串
 */
Number::Number(int precision, const std::string_view &num);

/**
 * @brief 根据整数序列构造大数对象
 * @param precision: 最大小数位精度
 * @param ints: 整数部分序列
 * @param decs: 小数部分序列
 */
Number::Number(int precision, const std::initializer_list<IntType> &ints = {},
	const std:::initializer_list<IntType> &decs = {});

【Pr.3: 函数封装】

封装是一件很简单的事，就算仅仅是把原本在你main函数的代码丢进一个名叫demo的函数里，这也算是封装——它完成了应该干的事，用户可以轻易明白名字叫demo的函数执行了一段示例代码，并且他们不必关注其内部细节（如果demo函数确实干了demo该干的事的话）。
然而，封装并不是如此简单粗暴的事，高质量的封装则更加不是。不过好在，本文并不侧重什么是以及如何做到高质量的封装，而只是怂恿你尽量封装并且能够选择在比较合适的时机封装。
Pr.2的初始化已经展现了封装的一角，而为了使主程序逻辑尽可能清晰，至少该为大数的加、乘、除、输出操作都封装一层。

以下是可供参考的封装形式：

//! C形式
number_t* num_add(number_t *lhs, const number_t *rhs); //! 大数相加
number_t* num_mul(number_t *lhs, IntType rhs); //! 大数-整形乘法
number_t* num_div(number_t *lhs, IntType rhs); //! 大数-整形除法
void num_print(const number_t *num);           //! 大数输出

//! C++形式
//! 以成员函数/友元函数定义
Number& Number::operator+=(const Number &rhs); //! 大数相加
Number& Number::operator*=(IntType rhs);       //! 大数-整形乘法
Number& Number::operator/=(IntType rhs);       //! 大数-整形除法
std::ostream operator<<(std::ostream &os, const Number &num); //! 大数输出

/* 以上加、乘、除操作的结果皆写入左操作数 */

【Pr.4: 迭代终点】

在第二节提供的代码中，算法迭代次数固定。如此操作，一方面当需求精度较小时，会存在大量无效的迭代过程，而当需求精度较大时，迭代次数又会不足，从而使得算法无法给出正确结果。如何针对目标问题确定迭代终点，这是Pr.4的第一个问题。当某个大数达到某个常数，又或许作为累加器的大数的变化量足够小，迭代就应该终止。这个终点是由具体问题确定的，无法预先给出通解。
另一个可以预先确定的迭代终点是，除法终止的时机。显然当数被除尽时，除法就可以终止。但若除不尽呢？我们可以提供给大数最大小数精度和当前小数精度两个参量，让其描述小数位的情况。这样一来，即便除法无法穷尽，当当前小数精度达到最大小数精度时，其也应该完成值的更新并立刻终止除操作。
可以发现，当该参量被引入后，本文目标问题的迭代终止条件也很容易确定：令pi和r的最大小数精度都匹配为需求精度n，由于r的初始小数精度不为$0$且r随着迭代进行逐渐减小，故而当r的当前小数精度回归$0$时（此时r的值即为$0$，也即变化量被置零），也就是$\pi$值计算迭代的终点。
同Pr.1，任何改变小数精度的操作都应该在方法内部完成对当前小数精度的维护。

【Pr.5: 常量检查】

上节代码中给出的宏值N与R并不总是那么可信，一个有效的N值必须满足对于任意已提供操作，其方法过程中不会产生值溢出。这表明N的大小是受到IntType表示范围的约束并且是存在明确上界的。
在本文中，大数的四则运算依赖竖式计算方法逐节点计算完成，每个节点单元存储的数据单元取值在$[0,R)$，若假定IntType型的右操作数同样在该取值当中，则两数相乘的最大值为$(R-1{)}^2$，产生的进位为$\lfloor \frac{(R-1{)}^2}{R}\rfloor$，将其作用于下一单元，则乘法产生的最大值即为$(R-1{)}^2+\lfloor \frac{(R-1{)}^2}{R}\rfloor$；约定$R\ge 10$，则最大值为$R^2-R-1$（$R=10^N$）。
即：一个合法的N应当至少满足$R^2-R-1\le MaximumofIntType$。

以下是基于模板的编译期静态检查策略（详细代码请查看源码）：

template <typename T, size_t N>
using Require = typename std::enable_if_t<
	//! 保证类型是无符号整形（本节依旧不引入数符）
	std::is_unsigned<T>::value
	//! 确保10^N是10的幂次（10^N计算过程本身也会溢出）
	&& ispow10_t<pow10_t<N>::value>::value
	//! 确保10^(2N)是10的幂次
	&& ispow10_t<pow10_squre_t<N>::value>::value
	//! 满足R²-R-1<=Maximum of IntType
	&& pow10_squre_t<N>::value - pow10_t<N>::value - 1 < std::numeric_limits<T>::max(),
	int>;

面向提高编程

---

概述：别拦着我，我还能再优化！
关键词：优化

结构优化

结构优化不是静态、一成不变的工作。良好的程序结构往往是从粗糙的大框架开始，而在编码过程中进行局部调整、调优，以使得接口对于编码者更友好，应用层对于用户更友好。
对于中等及以上规模并且大概率需要后期维护的程序，哪怕你有极大的自信心，上手即编码都是绝不可取的。另一方面，花费大量的时间在设计上，没有穷尽地吹毛求疵，过度设计也是不可取的。结构的设计终究是面向应用的，而非理想主义者的玩具。
故而，本文对于结构的优化不会单独讨论，而是选择穿插在下文中分散叙述。

·PS·如果有余裕的话，可以选择阅读文末GitHub地址上的源码，希冀其存在对你有所帮助的地方。

内存优化

考虑以下操作：

int *ptrs[65536] { nullptr };
for (auto &ptr : ptrs) {
	ptr = new int;
}
readonly_func(ptrs);
for (int i = 32768; i < 65536; ++i) {
	delete ptrs[i];
	ptrs[i] = nullptr;
	readonly_func(ptrs);
}
for (int i = 0; i < 32768; ++i) {
	delete ptrs[i];
	ptrs[i] = nullptr;
	ptrs[i + 32768] = new int;
	readonly_func(ptrs);
}
for (auto &ptr : ptrs) {
	if (ptr != nullptr) {
		delete ptr;
	}
}

在以上代码中，内存被密集而反复地申请和释放。对于每个被申请的内存单元，只执行了三次只读操作（假定readonly_func是对有效指针的遍历操作）。
假设readonly_func是程序的主功能，那么相比于readonly_func的调用，内存操作的代价是否可以接受？
或者说，一次内存操作的代价换三次readonly_func的执行，是否值得？
在多数情况下，语言标准库实现的内存申请释放操作都是摊销的常数时间复杂度，它们是高效且安全的。但是高效的根本是需要维护的，而安全则需要以性能为代价。维护代价是否足够大，安全是否并不重要，从这两个问题出发，大抵就能确定内存上的优化是否必要。
上例中，半数内存被释放之后紧跟着等量的内存被申请，假若将本该被释放的内存移用给被申请的部分，那么移用的代价是否要比new/delete的代价更小？这个结果依旧是不确定的，它取决于当前时刻堆内存的状态。在大多数情况下，连续小单元内存的申请释放操作是$O(1)$的，因此，这样的改进对于性能的提升可能并不明显。而当申请的内存大小不定，申请时机不规律，由于内存碎片的产生，维护内存堆的代价将被计入，此时选择自行实现一个内存池作为new/delete的替代则可以获得较为明显的性能提升。(本段仅为模糊的说法，非绝对，具体情况具体判断）
在本文中，目标问题为$\pi$值计算。鉴于乘除操作将带来前后缀零的擦除以及迭代计算的存在，考虑重新实现new/delete，如下：

Node* Number::allocate(IntType x, Node *prev, Node *next) {
	如果闲置节点链表非空:
		取出该节点
	否则:
		申请新节点
	初始化节点并返回
}
void Number::deallocate(Node *node) {
	断开节点的前驱后继并将该节点插入闲置节点链表
}

如上，被擦除的节点将会被保留，以便新节点的快速分配。但是相对的，若大数总是向双端增长，那么allocate就会退化为new的操作；若是数字足够大且零擦除的频率远大于申请新节点的速率，那么大量内存的保留就会造成内存的浪费，甚至于拖慢其它模块的内存申请操作。没有策略是普适的，编码者能做的只是选择最适合目标的那个策略。
至于将内存申请释放独立为一个方法而非在使用处进行代码替换，这是结构上的设计。若需要变更内存分配的方式，则只需要修改此函数而不需要改动其他源码。

性能优化

引子：编程是一项实际的工作，而非完全由理论支配。语言底层，编译器优化，硬件架构等诸多外界因素共同影响着程序的运行效率。理论复杂度优秀的算法，可能是对CPU不友好的；而过度手动优化，又大概率是对编译器不友好的。故而编程者要想编写出高效的程序，必须在理论支持的基础上，综合各类现实因素，压榨硬件等环境配件。

简单来说，程序就是一串指令，每一条指令都会消耗一个常数的时钟周期。故而提升性能的关键就是设计办法使得程序消耗的总始终周期数更少。达成这个目标大体有以下两点：

• 优化逻辑结构，替换算法，减少指令的总数目
• 使用消耗的时钟周期更少的等价的指令序列

本文的算法已确定，不再做大刀阔斧的更改。该节针对CPU的结构与技术，在消耗的时钟周期更少上做两方面的讨论。

① 分支预测

简介：CPU流水线技术下，对于分支逻辑CPU将预测分支并预先加载该分支的指令，预测成功则继续执行，否则丢弃已加载的指令并执行另一分支
性能提升点：预先加载指令而不必等待条件判断完成
惩罚：预测失败造成的预指令清除

分支预测分为静态预测和动态预测，静态预测在编译期由编译器完成，动态预测在运行期由预测的历史上下文动态调整该选择哪一分支。

• 静态预测：使用C的likely/unlikely宏或C++20的相应关键字指定某一个分支更容易被执行。
• 动态预测
  [1] 使分支条件呈现阶段性，分支预测将更容易成功
  [2] 减少分支语句的密集调用，减少分支预测的发生
  [3] 尽可能合并分支条件，提高分支预测的成功率
  [4] 当指令执行代价较大时，根据情况适当引入分支逻辑
  [5] 当指令较为简单时，如果可能的话使用三目或位运算等替代（不推荐）

对于目标问题，主要采取本节中动态预测[4]的策略，具体操作为：对于乘除操作，引入分支条件，当节点值为0时，取消本次运算。
大数乘除中的0并无规律可言，故而该分支预测的成功率是难以保证的。之所以采取这样的操作，则是假定在整个操作过程中，取消除法与取余操作的提升大于清除指令的总代价。

② 快速缓存

简介：CPU与内存的交互是慢的，现代CPU中的多级缓存架构就是为了解决这样的问题。在该架构下，从内存中读取数据会同时读取相邻的数据载入缓存行。由于空间和时间的局部性原理，相邻的数据大概率会在接下来的时间被访问，此时则可以从高速缓存行直接读取该数据而不必再与内存交互。
性能：现代CPU架构一般为三级缓存，越接近CPU的缓存响应速度越快。其中最快的一级缓存性能可以达到内存交互的两百倍。

一级缓存的大小一般为32字节或64字节。本文已有的代码实现中，每个节点只存储一个$IntType$，完全没有利用到高速缓存的优势。为了使尽可能多的有效数据被载入高速缓存，令每个节点存储4个$IntType$，并重写相关的初始化、内存分配、运算操作，完成基于快速缓存的性能优化。

·PS·对运算操作的改写连带了unroll操作，该操作在CPU并行计算方面也将一定程度上提升程序性能。

性能分析

测试目标

1. Plain：原始版本
2. Freespace：引入闲置节点链表
3. BranchMul：对乘法操作启用零值判断分支
4. BranchDiv：对除法操作启用零值判断分支
5. Freespace_BranchMul
6. Freespace_BranchDiv
7. Freespace_BranchMul_BranchDiv
8. Freespace_Cache：闲置节点链表+缓存优化

Performance of 1~4

n/ms	Plain	Freespace	BranchMul	BranchDiv
32	0.0377	0.0103	0.0194	0.0101
Rate	100%	27.35%	51.65%	26.68%
128	0.2553	0.2025	0.1904	0.1940
Rate	100%	79.32%	74.58%	75.98%
512	1.7027	1.6552	1.7334	1.7977
Rate	100%	97.21%	101.80%	105.58%
4096	89.5917	90.2696	84.6139	85.6809
Rate	100%	100.75%	94.44%	95.63%
16384	1524.8486	1530.8360	1471.8772	1476.7833
Rate	100%	100.39%	96.52%	96.84%
65536	25156.7890	25319.5595	24423.8712	24565.4975
Rate	100%	100.64%	97.08%	97.64%

Performance of 5~8

n/ms	Freespace_BranchMul	Freespace_BranchDiv	Freespace_BranchMul_BranchDiv	Freespace_Cache
32	0.0100	0.0100	0.0100	0.0322
Rate	26.70%	26.68%	26.70%	85.45%
128	0.1773	0.1821	0.1798	0.1512
Rate	69.43%	71.31%	70.42%	59.21%
512	1.7823	1.8040	1.8037	1.6016
Rate	104.67%	105.94	105.93	94.06%
4096	85.4299	85.2916	85.1736	79.3627
Rate	95.35%	95.20%	95.06%	88.58%
16384	1471.9906	1469.7545	1476.4153	1249.2727
Rate	96.53%	96.38%	96.82%	81.92%
65536	24541.7651	24527.6621	24558.4590	19244.1821
Rate	97.55%	97.49%	97.62%	76.49%

·PS·n为输入规模，以上数据取自1024次迭代的算术均值。

简单观察可以发现，提供的若干优化版本整体上较原始版本有了提升，但是在不同的数据规模下表现各有差异。输入规模较大时，该优化的效果并不明显，甚至于对于特定规模，无优化的版本反而占优势。
对于此现象，首先得清楚一点：最终的程序是经手编译器的代码。
编译器若干年的发展，各种优化手段都比较成熟。或许我们花大力气做的优化，人家编译器早已经默认在做了。而手动的优化在一定程度上破坏了源程序的结构，反而不利于编译器做进一步的优化。这也能解释为什么我们所谓的无优化版本存在比优化版本更快的现象。
故而，如果所谓的优化仅仅是将预先计算常量，简化分支，乘除取余化为位运算，分支替换为三目这类小trick的话，则完全不需要去废这个心神。打开编译器的二级优化，这些操作编译器就基本帮你完成了。
真正需要编码者动脑筋做的优化应该是选择更适配目标问题与目标机器的算法与数据结构，以及如何设计更利于开发与维护的结构、大框架。
另外由上图可见，对于缓存的优化是极为成功的。显然，编译器所能做的优化大抵是局部的等价替代，而做不到这种较为离谱的改变。而充分利用高速缓存，实际上也可以说是选择了一种更实适配目标机器的数据结构（压榨机器性能）。

附：为了更好的体会编译器的优化有多强大，下图放出-O3选项下n=65536规模下各目标的运行耗时。

·PS·-O3选项给出的优化可能导致意料之外的篡写，为了防止您的程序产生不明不白的BUG，请慎用该选项。

---

·PS·所有代码已上传[Github]https://github.com/ZYMelaii/NWPU-NOJ