PARI/GP 语言：从入门到实现大素数判定与大数分解

一、PARI/GP简介

PARI/GP 是一种针对数论中的快速计算(大数分解，代数数论，椭圆曲线等) 而设计且应用广泛的计算机代数系统，其具备大量实用的函数来对数学实体进行计算， 诸如矩阵，多项式，幂级数，代数数，以及相当多的超越方程等等。 PARI 也可以作为快速计算的C语言库。这个系统最初是由 Henri Cohen 和他的合作者们开的。

PARI/GP 在早期是一个专门的计算机代数系统，主要针对数论学家，目前已经在许多其他不同的领域得到了很好的应用，例如从拓扑、数值分析直到物理学。

虽然可以进行大量符号运算，但与 Axiom、Magma、Maple、Mathematica、Maxima 或 Reduce 等系统相比，PARI 在此类任务上表现很差，例如多元多项式、形式积分等。另一方面，该系统的三个主要优点是速度快，可以直接使用数学家熟悉的数据类型，以及广泛的代数数论模块(在上述系统中，只有 Magma 提供了类似的功能)。

非数学方面的优点包括可以使用高级脚本语言或 PARI 库编程，PARI 库是一个成熟的系统(其开发始于80年代中期)，用于完成和传播原始的数学研究，同时建立大量用户社区。当然，PARI/GP 是免费软件，受 GNU 通用公共许可证的保护。

我们不妨先对后面可能会出现的名词给出注解。

• PARI 是一个可以快速计算的C语言库；
• gp 是一个容易上手的交互式外壳，它赋予用户调用 PARI 函数的权力；
• GP 是 gp 脚本语言的名字；
• gp2c 是 GP 到 C 的编译器，目前 gp2c 只能运行 GP 语言的一部分。

PARI 有三种不同的用法：

1. 作为 libpari 库，可以从高级语言应用程序调用，例如用 ANSI C 或 C++ 编写；
2. 作为一种复杂的可编程计算器 gp，其语言 GP 包含了 C 等标准语言的大部分控制指令；
3. 编译器 gp2c 将 GP 代码转换为 C 代码，并将其加载到 gp 解释器中。由 gp2c 编译的典型脚本运行速度要快3到10倍。生成的 C 代码可以手工编辑和优化。

---

二、下载 PARI/GP

下载页面：Download PARI/GP
打开下载页面，我们在下述提示中选择下载的软件版本，此处选择的是 Win64 版（2022年12月31日发布）：
下载完成后，一路正常安装即可。如果需要更强大的功能，也可选择下载 PARI/GP 的额外包。

---

三、PARI/GP 编程入门

1. 近似

控制 PARI 的基本原则是，运算和函数首先应该给出尽可能精确的结果；其次，如果它们有任何意义，则是允许的。例如，$1$ 除以 $3$ 并不是 $0.333\cdot \cdot \cdot$，而是有理数 $1/3$ 。要得到浮点实数形式的结果，请这样计算 $1./3$ 或者 $0.+1/3$。

相反，在不精确的对象之间进行运算的结果，尽管本质上不精确，但最终将尽可能精确。例如，考虑两个实数 $x$ 和 $y$ 的相加，结果的准确性是一个先验不可预测的，它取决于 $x$ 和 $y$ 的精度，既取决于它们的大小，也取决于 $x+y$ 的大小。从这些数据中，PARI 计算出结果的正确精度。即便它在 gp 模式的计算中工作，其中有一个默认精度的概念，它的值也只用于将精确类型转换为不精确类型。

特别地，如果一个运算涉及不同精度的对象，一些数字将被 PARI 忽略。例如我们将实数 $x$ 写成 $r+2^eϵ$ ，其中 $r$ 是该实数的有理逼近，$|ϵ|<1$ 。

实际上，取 $s=+1、-1或0$ ，$m$ 为一多精度整数满足 $0\le m<2^l$ ，指数 $e$ 满足 $-2^B\le e<2^B$ (在64位电脑上，B=63)，则
$$|x-r|=|x-ms{2}^e|=2^eϵ<2^{e-l}$$ PARI 的计算方式为，任何小于 $2^e$ 的数都可以被视为精确的零，例如下面的例子：

0.E-28 + 1.E-100 = 0.E-28 ，0.E100 + 1 = 0.E100 。第二项计算应该是 $1<2^e$ ,被视为零了。

另外，我们考虑加法和乘法的区别，若 a = 2^(-100) ，
则 a + 0. = 7.888609052210118054 E-31 ，而乘法会更精确：
a * 1. = 7.8886090522101180541172856528278622967 E-31 。

注意这里符号 E 的意义，0.E10 = 10^(-10) = 0.0000000001，
0.E-1 = 0.E1 = 0.1，1E-2 = 0.01，1E2 = 100 ，2E-2 = 0.02。

2. 运算尽量被允许

第二个原则是 PARI 运算通常是尽量被允许的。例如，取一个向量的指数是没有意义的。然而，经常发生的情况是，我们想要将一个给定的函数应用到一个向量中的所有元素上。使用循环或 apply 内置函数很容易做到这一点，但事实上，PARI 假设这正是您在将标量函数应用到向量时想要做的事情。我们只需求一个向量的指数就可以了，不必写循环。大多数超越函数都以同样的方式进行计算。

其中，第一句代码 1/3 + Mod(1,5) 中的 1/3 代表 3 模 5 的逆元，实际上就是$Z/5Z$ 中的 2 ，后面的 Mod(1,5) 是 $Z/5Z$ 中的 1 ，两者相加则是 $Z/5Z$ 中的 3 。

第二句代码 I + O(5^9) 中，I 代表 $\sqrt{-1}$ ，我们得到一个 $5\text{-adic}$ 数，且该数 $5\text{-adically}$ 地接近 $\sqrt{-1}$ 。

最后一个例子难度可能更大，但是存在从 $Z/15Z$ 和 $Z/10Z$ 到 $Z/5Z$ 的自然映射，并且最终的加法在 $Z/5Z$ 内进行计算。

当然，标准运算符 $+、-、\ast 、/$ 是存在的。我们再次强调，除法是尽可能精确的运算：4 除以 3 得到 $4/3$。除此之外，对整数或多项式的运算，如 $\backslash$ (欧几里得除法)，$\%$(欧几里得余数) 是存在的；对于整数，$/$ 用作除法来计算分数，$\backslash$ 是求整除的商。当指数类型为整型时，还有幂运算符 $\stackrel{ˆ}{}$ ，否则它被认为是一个超越函数。最后的逻辑运算符 $!$ (非算符)、$\&\&$ (且算符)、$||$ (或算符) 也是存在的，给出结果1(真) 或 0(假)。

3. 向量

我们通过下面的方式创建向量：

其中 poltchebi(5) 是创建一个 5 阶的切比雪夫多项式。向量索引如下：

创建取值范围在前 10 个素数内且，满足 $(p,p+2)$ 系孪生素数对的所有素数 $p$ ：

gp > [ p | p <- primes(10), isprime(p+2) ]

输出结果为 [3, 5, 11, 17, 29] 。

创建取值范围在前 10 个素数内且满足模 4 余 1 的素数：

gp > [ p | p <- primes(10), p % 4 == 1 ]

输出结果为 [5, 13, 17, 29]

4. 矩阵

矩阵创建可参考如下两个例子：

矩阵元素的索引类似 Python ：

---

四、大整数素性证明

我们直接引入 isprime() 函数进行素性证明（判定），具体使用如下：

isprime(x, {flag = 0}) 

这个方法能严格地证明或否定一个数字是素数，当 $x$ 确实是一个大素数时，这个过程可能会非常慢。例如，一个 1000 位数的素数在默认算法下需要 15 到 30 分钟。

使用 ispseudoprime() 可以快速检查一个合数。使用 primecert() 来获得一个原始证明，而不是一个是或否的答案。

isprime() 函数接受向量或矩阵作为参数，如果 flag=0 ，则使用了如下方法的组合：

• Baillie-Pomerance-Selfridge-Wagstaff 综合检验，能快速检验一个合数；
• Selfridge “p-1” 测试，如果 $x-1$ 足够光滑；
• Adleman-Pomerance-Rumely-Cohen-Lenstra (APRCL)，用于一般中型的整数 $x$ (小于 1500 bits)；
• Atkin-Morain 椭圆曲线素性证明（ECPP），对于较大 $x$ 的时候。

如果 flag=1，则使用 Selfridge-Pocklington-Lehmer “p−1” 测试，这需要部分分解各种辅助整数，而且可能会非常慢。

如果 flag=2，则只使用 APRCL。

如果 flag=3，则只使用 ECPP。

例如，我们判定整数 $300000000000000000053$ 的素性，之前通过 NZMATH 中的 ECPP 耗时 $63.667$ 秒才将其判定为素数。而通过 GP 语言中的 isprime() ，我们可以瞬间判定其为素数。

当判定一个 93 位整数的素性时 ：
$$\begin{array}{rl}& 153211620887015423991278431667808361439217294295901387715\\ & 486473457925534859044796980526236853\end{array}$$ 前面通过 NZMATH 中的 ECPP 耗时 $381.815$ 秒才将其判定为素数！而通过 GP 语言中的 isprime() ，我们可以瞬间判定其为素数。

而对于一个 117 位的整数：
$$\begin{array}{rl}& 6518718497852438209581196842331429821705369381473495010499\\ & 55875423544371888519632103159753836519854149224951224466197\end{array}$$ 我们也几乎在 1 秒钟内完成其为素数的判定证明。

至于一个更大的素数，比如 $10^{451}+451$ ，如果直接使用命令：

gp > isprime(10^451+451)

则会报下述错误：

我们如何解决 “the PARI stack overflows !” 的错误呢？上图中 hint 给出的提示实际上并不好办，因为按要求设置了 gprc 配置文件后，还是不能解决这个问题。下面提供一个更容易实现的方法，就是添加一行语句：

default(parisize, 10000000000);

则此时 “new stack size” 就变成了 9536.743 Mbytes (约为 9.313 GB)，比默认的 7.629 Mbytes 要大多了。

此后再通过 isprime(10^451+451) 判断，就能在大约 6 秒内完成是素数的判断。

对于最后这个 1025 位的大整数：
$$\begin{array}{rl}& 4573150652070000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000000000000000000000\\ & 0000000000000000000000000000000000001\end{array}$$ 我们同样通过 isprime() 也可以瞬间判定上述大整数的确为素数。实际上，所费时间竟不到 1 秒钟！这解决了本人前文大整数的椭圆曲线素性证明——从ECPP简介到Python实现中，在最后提出的 “提高 ECPP 代码的运行效率” 的问题。

上述所有计算的运行代码如下：

---

五、大整数的素因子分解

我们继续分解上次在大整数分解与两类特殊的素数及其猜测文章中未分解彻底的 $(2^{353}+1)/3$ ，这是一个 106 位的合数，通过在线的椭圆曲线和二次数域筛法经过长达 7 个多小时计算，仍未分解出来。而经过运用 GP 语言中的 factor() 方法，即可实现更快速的大整数分解。

经过几个小时的计算(具体时间忘记了，但远未达到 6小时)，最后 GP 给我们分解出了 $(2^{353}+1)/3$ 的两个素因子，即：
$$\begin{array}{rl}\frac{2^{353}+1}{3}=& 3803909572078746837295094051706948091\times \\ & 1607818533384485707707842837146335251\\ & 451162017762519557029955613946641\end{array}$$ 在上式的两个素因子中，前者是一个 37 位的素数，后者是一个 70 位的素数。

---

六、关于Wagstaff 素数

事实上，我们将形如 $\frac{2^p+1}{3}$ 的素数称为 Wagstaff 素数（瓦格斯塔夫素数），其中 $p$ 是奇素数。瓦格斯塔夫素数出现在新梅森猜想中，并在密码学中得到应用。

新梅森猜想是指，对于任何奇自然数 $p$，如果以下任何两个条件成立，那么第三个条件也成立：

• $p=2k\pm 1$ 或 $p=4k\pm 3$ 对于某个自然数 $k$ 成立。
• $2^p-1$ 是素数（梅森素数）
• $(2^p+1)/3$ 是素数（瓦格斯塔夫素数）

如果 $p$ 是奇合数，则 $2^p-1$ 和 $(2^p+1)/3$ 都是合数。因此，只需要测试奇素数 $p$ 来验证猜想的真实性。

截止目前2023年1月，以上所有三个条件都成立的已知数字只有 9 个：$3、5、7、13、17、19、31、61、127$ 。

Renaud Lifchitz 通过系统地测试所有已知其中一个条件成立的素数，证明了新梅森猜想对于所有小于等于 $32582656$ 的整数都是真的。即要么三个都满足，要么只满足其中一项。

截至2023年1月，产生瓦格斯塔夫素数的指数 $p$ 为（已经被证明是素数的有33个）：

3， 5， 7， 11， 13， 17， 19， 23， 31， 43， 61， 79， 101， 127， 167， 191， 199， 313， 347，701， 1709， 2617， 3539， 5807， 10501， 10691， 11279， 12391， 14479， 42737， 83339， 95369， 117239 。

同样截至2023年1月，产生瓦格斯塔夫可能素数 (PRP) 的已知指数 $p$ 为（概率素数共有11个）：

127031、138937、141079、267017、269987、374321、986191、4031399、…、13347311、13372531、15135397。

其中，$4031399$ 发现于2010年2月，它有 $1213572$ 位数。$13347311$ 和 $13372531$ 发现于2013年9月，它们分别为具有 $4017941$ 位和 $4025533$ 位的概率素数。目前尚不清楚 $4031399$ 和 $13347311$ 之间是否有其它指数产生瓦格斯塔夫可能素数。2021年6月，$15135397$ 被发现，这是一个具有 $4556209$ 位数的素数，该数在目前的 PRP 记录中排在第三位。

特别地，Lucas–Lehmer–Riesel 测试可以检验瓦格斯塔夫概率素数，如果 $p$ 是瓦格斯塔夫素数的指数，则
$$25^{2^{p-1}}\equiv 25(mod{2}^p+1)$$ 另一方面，我们有瓦格斯塔夫素数的推广形式：
$$Q(b,n)=\frac{b^n+1}{b+1}$$ 其中基底 $b\ge 2$ 。如果 $n$ 为奇数，则我们有：
$$Q(b,n)=\frac{(-b{)}^n-1}{(-b)-1}=R_n(-b)$$ 其中 $R_n(-b)$ 表示以 $-b$ 为基底的 Repunit 数。

对于 $b=10$ ，此时的广义瓦格斯塔夫素数 $Q(10,n)$ 有如下六个：
$$\begin{array}{rl}& 9091、909091、909090909090909091、\\ & 909090909090909090909090909091、\\ & 9090909090909090909090909090909090909090909090909091、\\ & 9090909090909090909090909090909090909090909090909090\\ & 90909090909091\end{array}$$ 它们依次为 $4$ 位、$6$ 位、$18$ 位、$30$ 位、$52$ 位、$66$ 位的素数（位数都是偶数）。

进一步，当 $n=5,7,19,31,53,67,293,641,2137,3011,268207,1600787$ 时，$Q(10,n)$ 都是素数（或概率素数）。

我们也可以猜测，广义瓦格斯塔夫素数有无穷多个，这又是一个啃不动的硬骨头。