可逆需求的提出
在 VtuberMusic 网站前端开发的过程中，有一个需求是需要实现随机播放的功能。这本来是一个很平常的功能，一开始的做法是，随机播放模式下，上一首和下一首都是随机挑选一首歌播放。但是有个新需求又被提了出来：希望在随机播放模式下也能正常返回上一首。我们在网易云音乐上做了实验，发现网易云的随机播放支持不断地回退，能够准确溯源之前的若干首歌，甚至回退之后还能再准确地“前进”。我意识到也许可以实现一个“可逆”的随机数生成器，它不仅可以根据当前状态生成“下一个”随机数，还可以回退状态，生成“上一个”随机数，这样就可以完全解决这个问题。而使用常见的线性同余算法来生成伪随机数的话，这是非常容易做到的。于是我把这个任务接了下来，并且在后来把这个模块抽取出来进行了一下简单封装，开源了一份 JS 版本的代码：github 仓库。
模块的具体使用方法在仓库中可以直接看到，这里主要是介绍一下这个“可逆” PRNG 到底是如何实现的。
线性同余算法
简介
线性同余 (Linear Congruence) 算法是一个常见的伪随机数产生算法，简单来说，它用以下方式产生伪随机数序列 $\lbrace I_n \rbrace$：
$$I_{n+1}=(aI_n+c) \bmod m$$
例如，取 $a=51,\ c=3,\ m=100$，再随意地取 $I_0=1$，我们可以得到一串“随机数”：$54,57,10,13,66,...$
线性同余算法被广泛地采用。例如，C++11 标准中，<random> 库内置了三种随机数引擎的实现，其中之一就是线性同余，另外两个则是梅森缠绕器 ( mt19937 ) 算法和带进位减（一种延迟斐波那契）算法。线性同余的优点是计算简单，缺点则是随机数的质量不够好（具体详见后面的小节），但是对于大多数要求不高的随机数使用场合，使用线性同余就足够了。
Hull-Dobell 定理
线性同余算法虽然简单，但是在参数 $(a,c,m)$ 的选择上颇有讲究，可以说直接关系到“随机数”的质量。不加证明地给出 Hull-Dobell 定理：
当满足：
- $c$ 与 $m$ 互质
- $(a-1)$ 可以被 $m$ 的所有质因数整除
- 如果 $4 \mid m$，则必须 $4 \mid (a-1)$
那么线性同余生成的伪随机序列 $\lbrace I_n \rbrace$ 达到最大周期 $m$。也就是说，此序列在一个周期中取遍 $0$ 到 $(m-1)$ 的所有取值。
原根
Hull-Dobell 定理保证了伪随机序列的周期是足够大的，而周期无疑是衡量伪随机数性质的一个重要指标。不过，我们不一定会选择满足 Hull-Dobell 定理的参数。另外一种参数的选取基于以下方法：
定义：如果 $a$ 模 $n$ 的阶是 $\varphi(n)$，则称 $a$ 为模 $n$ 的一个原根。
其中，$\varphi(n)$ 是欧拉函数，表示小于 $n$ 的正整数中与 $n$ 互质的个数，特别地，$\varphi(1)=1$。
至于阶，可以简单这样理解：阶是满足 $a^i \equiv 1 \pmod n$ 的所有正整数 $i$ 中最小的 $i$。
我们考虑，如果取模数 $m$ 为一个质数，那么就有 $\varphi(m)=m-1$，根据定义，这是显然的。我们再取 $a$ 为 $m$ 的一个原根，取 $c=0$，我们可以保证这样构造的伪随机序列的周期是 $(m-1)$，一个周期中取遍 $1$ 到 $(m-1)$ 的所有可能取值。因为根据上述内容，对于任意的正整数 $k\in[1, m-1]$，要令 $a^ik \equiv k \pmod m$，$i$ 至少为 $(m-1)$。
基于这种方法选择的参数也能保证周期足够大，但是要注意初值 $k$ 不能为 $0$，同时它也无法生成随机数 $0$。
更多性质
上述两种方法都保证了伪随机序列的周期足够大，并且也保证了伪随机数的分布服从离散域上的均匀分布。原因是显然的，因为每个周期中每个可能的取值正好被取了一遍。
另外介绍一个衡量随机数质量的指标：Knuth's Spectral Test。它要求随机数序列在更高维的空间中也尽量均匀。
注：下面这一段属于个人猜测，可能不一定是真实的做法。
举个例子，我们有随机数序列 $0,1,2,3,...$ 满足 $I_k = k \bmod n$。从周期上来看，只要 $n$ 足够大，周期就能足够长。从一维的角度来看，它确实也保证每个周期取遍 $0$ 到 $(n-1)$ 之间的所有整数。但是如果我们连续取出两个数 $(a,b)$ 作为一组，表示二维平面上的一个点的话，不难发现我们获得的所有点并没有在平面上均匀分布，而是都落在若干条直线上。因此，这样的伪随机序列就“不够好”。
摘抄我阅读的知乎回答（详见文末参考）中的一段话：
要在伪随机数生成器界混，仅仅入门是不够的。从工程的角度来讲，的值要（在合理的范围内）足够小，以避免溢出的问题。从安全（实用）性的角度来讲，还要满足良好的随机性，这一点可以通过 Knuth's Spectral Test 来评估，要通过 2, 3, 4, 5 以及 6 维的 Spectral Test 才行。
因此 $(a,c,m)$ 的参数选择是非常讲究的。关于一些可用的参数，可以选择去此知乎回答找，也可以看后面的“JS 代码细节”一节，里面有介绍几组常用参数。
可逆的基本实现
由于线性同余算法生成随机数是利用“当前”随机数去推算“下一个”随机数，因此我们如果能找到逆向推算的方法，用“下一个”随机数反推出“当前”随机数，那么自然就实现了可逆的效果。幸运的是，这是可以做到的。
定义：如果满足 $aa_{inv} \equiv 1 \pmod m$，就称 $a_{inv}$ 为 $a$ 在模 $m$ 意义下的乘法逆元，也可形式地记作 $a^{-1} \pmod m$。
逆元又可以被称为“数论倒数”，因为它具有和倒数非常类似的性质。利用逆元 $a_{inv}$，我们来尝试推出 $I_n$ 与 $I_{n+1}$ 的关系。
$$I_{n+1}=(aI_n+c)\bmod m$$
不失一般性，不妨令 $c\in[0, m-1]$。于是可以两边消去 $c$：
$$(I_{n+1}+m-c)\bmod m=aI_n\bmod m$$
两边乘上 $a_{inv}$：
$$a_{inv}(I_{n+1}+m-c)\bmod m=I_n \bmod m$$
从而：
$$I_n=a_{inv}(I_{n+1}+m-c)\bmod m$$
至此，我们完成了利用 $I_{n+1}$ 推算 $I_n$ 的任务。用代码来描述的话就是：
prev = a_inv * (cur + m - c) % m
于是只要实现一个基于线性同余的简单 PRNG，然后保存参数 a 在模 m 意义下的逆元 a_inv，就能够实现可逆随机数生成器了。另外必须指出，这样的设计允许你设置初始的状态，或者说初始的“随机数”。
把这样一个随机数引擎给出的随机数变换到任意整数区间 $[\mathrm{min},\mathrm{max}]$ 上也是容易的，只需要把输出模上区间长度，再加上区间的左端点即可，可以认为这样得到的分布仍然“近似均匀”。同时，如果需要在这个区间上指定初值，我们会发现有许多不同的初始状态都可以对应这个区间上的初值，此时可以随机选择一个，保证多次设置后的随机序列是尽量不同的。具体可参见 github 仓库。
虽然逆元有如此美妙的性质，但我们还有两个问题亟待解决：
- 逆元一定存在吗？
- 逆元怎么求？
可逆的数学保证
数论：裴蜀定理
要解答逆元存在性的问题，我们引入数论领域的经典定理：裴蜀定理。
裴蜀定理 (Bézout's identity)：
对于整数 $a,b$，存在整数 $x, y$ 使得：
$$\gcd(a,b) = ax+by$$
引理 1：
不定方程 $ax+by=c$ 在 $\gcd(a,b) \nmid c$ 时无解。
证明：$\gcd(a,b)$ 整除等式左边但不整除右边，显然无解。
推论 1：
整数 $a,b$ 互质的充分必要条件是存在整数 $x, y$ 使得：
$$ax+by=1$$
证明：必要性由裴蜀定理直接给出。充分性：若存在这样的 $x,y$，由引理 1 知 $\gcd(a,b) \mid 1$，从而 $\gcd(a,b)=1$。
推论 2：
整数 $a,n$ 互质的充分必要条件是存在整数 $x$ 使得：
$$ax \equiv 1 \pmod n$$
证明：由推论 1，显然。
显然，这里的 $x$ 就是 $a$ 在模 $n$ 意义下的逆元。所以结论很简单：如果 $a,m$ 互质，那么逆元就存在。
接下来我们来初等地证明一下裴蜀定理。
已知整数 $a,b$，记集合 $S=\lbrace ax+by\mid x,y\in \mathbb{Z} \ \land \ ax+by>0 \rbrace$
记整数 $d=ax_0+by_0$ 为集合 $S$ 中的最小元素。如果能证明 $d=\gcd(a,b)$，那么就证明了裴蜀定理，因为 $d\in S$。
设 $a=qd+r$，其中 $q$ 和 $r$ 是 $a$ 模 $d$ 的商和余数，因此有 $0\leq r<d$。代入 $d=ax_0+by_0$，有：
$$r=a-q(ax_0+by_0)=a(1-qx_0)-b(qy_0)$$
这意味着 $r$ 必须为 $0$。否则就有 $r\in S$ 并且 $r<d$，这与我们假定的 $d$ 为 $S$ 中最小元素相矛盾。
从而我们有 $d\mid a$。同理有 $d\mid b$。从而 $d\mid \gcd(a,b)$。又注意到 $d=ax_0+by_0$，因此 $\gcd(a,b)\mid d$。
综上可得，$d=\gcd(a,b)$。于是裴蜀定理得证。
算法：扩展欧几里得算法
要解答逆元如何求的问题，我们引入算法领域经典的扩展欧几里得算法。
首先考虑欧几里得算法，它是经典的求最大公约数 (GCD) 的算法，也称辗转相除法，无须多言，贴上 C++ 代码：
int gcd (int a, int b) { return b == 0 ? a : gcd(b, a % b); }
而扩展欧几里得算法就是在递归进行上述计算的同时，计算出一组使得 $\gcd(a,b)=ax+by$ 的解 $(x,y)$。
先考虑递归基。上述算法将终止于 $a_n=\gcd(a,b),\ b_n=0$ 处，其中 $n$ 表示递归的层数。此时，$\gcd(a,b)=a_nx_n+b_ny_n$ 的一组解为 $x_n=1,\ y_n=0$。
考虑递归的退栈过程。假设我们求出了第 $n$ 层时的一组解 $(x_n,y_n)$，我们注意到 $(n-1)$ 层和 $n$ 层之间有如下关系：
$$a_n=b_{n-1},\ b_n=a_{n-1}-qb_{n-1}$$
其中，$q$ 等于 $a_{n-1}$ 整除 $b_{n-1}$ 的商，也即 a / b。
代入等式 $\gcd(a,b)=a_nx_n+b_ny_n$，我们有：
$$\gcd(a,b)=b_{n-1}x_n+(a_{n-1}-qb_{n-1})y_n=a_{n-1}y_n+b_{n-1}(x_n-qy_n)$$
与 $(n-1)$ 层的解 $\gcd(a,b)=a_{n-1}x_{n-1}+b_{n-1}y_{n-1}$ 比较系数得：$x_{n-1}=y_n,\ y_{n-1}=x_n-qy_n$
转换为代码就是：
int tmp = y; y = x - (a / b) * y; x = tmp;
从而得到扩展欧几里得算法的 C++ 代码：
int exgcd (int a, int b, int &x, int &y) { if (b == 0) { x = 1, y = 0; return a; } int r = exgcd(b, a % b, x, y); int tmp = y; y = x - (a / b) * y; x = tmp; return r; }
计算得到的 x, y 就是原不定方程的一组解。对于求逆元问题，我们只需要调用 exgcd(a, n, x, y)，得到的 x 就满足条件：a * x % n == 1。但是我们想要把 x 变换到 [0, n-1] 上来。考虑到 x 可能为负，因此一种安全的转换方法是：
x = (x % n + n) % n;
至此，我们就求得了 $a$ 在模 $n$ 意义下的逆元。
值得指出的是，就算 $a$ 和 $n$ 不互质，此算法仍然会给出一个解，但是显然它不满足逆元的条件。因此最好还需要检验一下所得到的解，也就是：
assert(a * x % n == 1);
JS 代码细节
JS 整数计算精确性问题
需要指出的是，因为 JS 的一些特性，我们并不能随心所欲地选择模数 $m$。这是因为 JS 的数值都是 Number 类型，并不严格区分整数和浮点数，导致浮点数的运算误差被引入整数。JS 中有一个参数 Number.MAX_SAFE_INTEGER，表示能够安全进行运算的最大整数。举个例子，C++11 标准中的线性同余随机数引擎的两组预定义参数分别为 $(16807, 0, 2147483647)$ 和 $(48271, 0, 2147483647)$，而这两组参数无法在 JS 实现中使用，因为在做乘法等运算的时候会溢出最大安全整数。
作为妥协，我在 JS 实现中选择了常用的一组 "Magic Number" $(9301, 49297, 233280)$ 作为默认参数。同时，在求出逆元 inv_a 后检查是否满足 a * inv_a % m == 1 来检查 a 和 m 的互质性，以及检查是否溢出最大安全整数。另外，JS 新标准中引入了新的数据类型大整数类 BigInt，利用它就可以摆脱整数运算的误差，从而可以设置任意的模数 $m$。但是由于是新标准，目前还没有得到广泛支持，因此使用时需要注意兼容性。我另外实现了一个使用了 BigInt 的版本，在此版本中的默认参数为 $(48271, 0, 2147483647)$。具体可以参见 github 仓库。
实际上，使用 $(9301, 49297, 233280)$ 这组参数是有潜在问题的，原因在于 233280 不是一个质数，如果把生成的随机数模上 233280 的一个因子，会使得这个 PRNG 的周期变成这个因子，比如模上 3 以后，PRNG 的周期就变成了 3。考虑到使用场景是随机播放，这显然是不好的。又因为在随机播放任务中没必要考虑更高维度的分布均匀性，所以实际上选择一个质数作为 $m$，再选择它的一个原根（这是为了到达 $(m-1)$ 的最大周期）作为 $a$，同时选择 $c=0$ 即可。
JS 实现扩展欧几里得
JavaScript 是引用语义的，因此在 JS 中实现扩展欧几里得算法需要特别注意这一点。直接上代码（只展示没有使用 BigInt 的版本，使用了 BigInt 的版本区别也不大）：
// Extended Euclidean Algorithm. function exgcd (a, b, vec2) { if (b == 0) { vec2.x = 1; vec2.y = 0; return a; } var r = exgcd(b, a % b, vec2); var tmp = Number(vec2.y); // copy vec2.y = vec2.x - Math.floor(a / b) * vec2.y; vec2.x = tmp; return r; } // Find the inverses of a (mod n). // That is, find x s.t. a * x % n = 1. function getInverses (a, n) { var vec2 = { x: 0, y: 0 }; exgcd(a, n, vec2); // We don't need the gcd of a and n. return (vec2.x % n + n) % n; }
使用了一个对象 vec2 作为 x 和 y 传递的载体，同时用 Math.floor 实现类似整除的效果。在 JS 编程过程中涉及到引用语义和值语义区分的地方一定要非常小心仔细地对待，而且最好熟知引用语义的原理，否则仍然非常容易出错。
参考
特别感谢 @lkmcfj 大佬的全程帮助和解答！