Python NZMATH 模块的数论计算探索之一

前言

在上一篇文章中，笔者给出了 NZMATH 模块的介绍和安装，以及进行 ECPP 大整数的素性判定。本文我们将继续探索 NZMATH 的数论计算功能。

---

一、Arith1 方法 —— 算术函数杂项汇总

1.1 开根取整

对自然数 $N$ 的根号取下整：

from nzmath import *

arith1.floorsqrt(22)

以上输出结为 $4$ 。事实上，$\sqrt{22}\approx 4.6904$ ，取下整即朝着 $0$ 的方向取整，结果的确为 $4$ 。注意，此方法只能针对正整数进行运算。

对自然数 $N$ 开高次根号并取下整：

arith1.floorpowerroot(1000, 4)

上述代码输出结果 $5$ ，系对 $1000$ 开 $4$ 次根号并取下整。

1.2 Legendre (Jacobi) 符号

Legendre 符号对于数论爱好者都并不陌生，即下述定义（$p$ 为奇素数） $$\left(\frac{a}{p}\right)=\{\begin{array}{rlrl}1& ,& & p\nmid a且存在整数x使得x^2\equiv a\mathrm{mod}p\\ -1& ,& & 不存在整数x使得x^2\equiv a\mathrm{mod}p\\ 0& & & p\mid a\end{array}$$ 下面求解勒让德符号 $\left(\frac{399541}{999907}\right)$ ，其中分子分母皆为素数。

arith1.legendre(399541,999907) 

输出结果为 $1$ ，即 $399541$ 是模 $999907$ 的二次剩余。这表明二次同余方程 $x^2\equiv 399541\mathrm{mod}999907$ 有解。那么这个解的一个特殊情形又是什么呢 ？我们使用下述语句来求 $x$：

arith1.modsqrt(399541, 999907)

输出结果为 $108094$ ，可以验证，$x\equiv 108094\mathrm{mod}999907$ 是 $x^2\equiv 399541\mathrm{mod}999907$ 的解。它还有另一个伴随解 $x\equiv -108094\equiv 891813\mathrm{mod}999907$ ，如果不计此伴随解，则二次同余式 $x^2\equiv 399541\mathrm{mod}999907$ 仅有唯一的一组解。

下面继续求解勒让德符号 $\left(\frac{343917}{999907}\right)$ ，其中的分子为合数。

arith1.legendre(343917,999907)

输出结果为 $1$ ，表明二次同余方程 $x^2\equiv 343917\mathrm{mod}999907$ 有解。

arith1.modsqrt(343917, 999907) 

通过上述代码求其解为 $x\equiv 468521\mathrm{mod}999907$ ，该解也是唯一的。

Jacobi 符号是勒让德符号的推广，由勒让德符号定义如下：$$\left(\frac{a}{n}\right)={\left(\frac{a}{p_1}\right)}^{{\alpha }_1}{\left(\frac{a}{p_2}\right)}^{{\alpha }_2}\cdots {\left(\frac{a}{p_k}\right)}^{{\alpha }_k}$$ 其中 $n=p_1^{{\alpha }_1}{p}_2^{{\alpha }_2}\cdots p_k^{{\alpha }_k}$ 为大于 $1$ 的奇数。

现在我们求解雅可比符号 $\left(\frac{844143}{985067}\right)$ ，其中分子分母皆为合数：

arith1.legendre(844143,985067)

输出结果为 $1$ ，对于雅可比符号，此时还不足以说明二次同余方程 $x^2\equiv 844143\mathrm{mod}985067$ 有解。我们考虑下面的运算：

arith1.modsqrt(844143, 985067)

上述代码输出结果为 $639990$ ，但实际上，$639990$ 并不是二次同余方程 $x^2\equiv 844143\mathrm{mod}985067$ 的解。出现这个问题的原因是，arith1.modsqrt(a, p) 中的第二个参数 $p$ 必须是素数，否则会出现如上错误的结果。 那么该同余方程是否有解 ？我们通过下述代码进行搜索：

rst = []
for i in range(1,985068):
    if (i**2 - 844143) % 985067 == 0:
        rst.append(i)
        
rst

结果输出的是一个空列表，说明二次同余方程 $x^2\equiv 844143\mathrm{mod}985067$ 无解。因此雅可比符号为 $1$ ，也不能说明二次同余方程有解。但雅可比符号为 $-1$ ，则一定有二次同余方程无解。例如：

arith1.legendre(308081,688781)

上述雅可比符号 $\left(\frac{308081}{688781}\right)$ 的计算结果为 $-1$ ，表明 $x^2\equiv 308081\mathrm{mod}688781$ 无解。

最后举一个雅可比符号为 $1$ ，且有解的例子 $\left(\frac{308071}{688781}\right)$ ，（其中分子分母皆为合数）：

arith1.legendre(308071,688781)

计算输出雅可比符号的值为 $1$ ，而且此时二次同余方程 $x^2\equiv 308071\mathrm{mod}688781$ 有多组解 $$x\equiv 16977,162705,526076,671804\mathrm{mod}688781$$ 即便去掉伴随解，此时仍有不唯一的两组解 $$x\equiv 16977,162705\mathrm{mod}688781\text{\hspace{1em}(1)}$$ 值得注意的是 $n=23\cdot 29947$ ，则有 $$\left(\frac{308071}{688781}\right)=\left(\frac{308071}{23}\right)\left(\frac{308071}{29947}\right)$$ 由程序计算知，$$\left(\frac{308071}{23}\right)=\left(\frac{308071}{29947}\right)=1$$ 因此二次同余方程 $x^2\equiv 308071\mathrm{mod}29947$ 有唯一解 $$x\equiv 16977\mathrm{mod}29947\text{\hspace{1em}(2)}$$ 事实上，满足 $(1)$ 式中 $x\equiv 16977\mathrm{mod}688781$ 的 $x$，也一定满足 $(2)$ 式。故同余方程 $x^2\equiv 308071\mathrm{mod}688781$ 与 $x^2\equiv 308071\mathrm{mod}29947$ 各自都有一部分解是相同的，但两者解集不完全相同。

以下我们验证雅可比符号的二次互反律，其中要求 $\left(\frac{a}{n}\right)$ 的 $a$ 和 $n$ 为大于 $1$ 的奇数：

a = 308071
n = 688781
A = arith1.legendre(a,n)
B = (-1)**((a-1)*(n-1)//4)*arith1.legendre(n,a)
A == B

输出结果为 True ，这样就验证了二次互反律。

1.3 模 $q$ 的乘积逆元

对于整数 $n$ 和 $q$ , $(n,q)=1$ ，若有 $m$ 使得 $n\cdot m\equiv 1\mathrm{mod}q$ ，则称 $m$ 为模 $q$ 下 $n$ 的逆元。我们可以很容易地进行求解，比如 $111111\cdot m\equiv 1\mathrm{mod}993341$ ，来求解逆元 $m$：

arith1.inverse(111111, 993341)

输出结果为 $m=408347$ 。

1.4 中国剩余定理的求解

中国剩余定理为求解如下的一次同余方程组：$$\{\begin{array}{rl}& x\equiv a_1\mathrm{mod}{m}_1\\ & x\equiv a_2\mathrm{mod}{m}_2\\ & \cdots \\ & x\equiv a_n\mathrm{mod}{m}_n\end{array}$$ 其中要求 $m_1,m_2,...,m_n$ 两两互素。例如求解如下的方程组：$$\{\begin{array}{rl}& x\equiv 2\mathrm{mod}250777\\ & x\equiv 3\mathrm{mod}555677\\ & x\equiv 5\mathrm{mod}794111\end{array}$$ 该方程组的求解代码为：

# 列表中的第一个元素是余数，第二个是模数
arith1.CRT([[2, 250777], [3, 555677], [5, 794111]])

输出的结果为 $79143032176845751$ ，而实际的结果为 $$x\equiv 79143032176845751\mathrm{mod}110660170719250219$$ 其中模数为 $110660170719250219=m_1\cdot m_2\cdot m_3$ 。

1.5 判断平方数

print(arith1.issquare(1234321))
print(arith1.issquare(7654321))

输出结果 $1111$ 和 $0$ ，前者表明 $1234321=1111^2$ ，后者表明 $7654321$ 不是一个平方数。

二、Arygcd 方法 —— 三类最大公约数的计算

2.1 整数的最大公因子

arygcd.binarygcd(113,1017)

输出结果为 $113$ 。

2.2 高斯整数的最大公因子

我们也可以求两个高斯整数的最大公因子：

arygcd.arygcd_i(1,13,13,11)

输出结果为 $(1,1)$ ，即高斯整数 $1+13i$ 与 $13+11i$ 的最大公因子为 $1+i$ 。

2.3 艾森斯坦整数的最大公因子

进一步，开可以计算两个艾森斯坦整数的最大公因子。艾森斯坦整数的定义如下：$$z=a+b\omega ,\omega =e^{\frac{2\pi i}{3}}$$ 以下给出计算代码：

arygcd.arygcd_w(2, 13, 33, 15)

输出结果为 $(4,5)$ ，即 $2+13\omega$ 与 $33+15\omega$ 的最大公因子为 $4+5\omega$ 。

至此，我们已经初步体会了 NAMATH 在数论上的编程实践，也欢迎读者浏览本文的后续内容。