2021SC@SDUSC
本章综述
大整数ZZ类主要实现了任意长度大整数表示、最大公因数、Jacobi符号和素性检验。笔者将通过逐个分析ZZ.cpp源代码中函数的形式来一步步向读者展示NTL是如何实现上述功能的。
计算最大公因数（gcd）
（1）数学基础：（广义）欧几里得除法
知识储备（定理，公立，公式）
·如果 b|a ,则(a,b) = b;
·如果a,b为两整数，则(a,b) = (b,a)
·如果 p为素数，a为整数，且p ∤ a，则a和p互素
证明：设（a,p） = d ,则有d|p，且d|a 。因为p是素数，所以d = 1或者d = p
对于d = p ,和p ∤ a矛盾，所以d = 1，即，（a，p） = 1.结论成立
·设b为任意正整数，（0，b） = b
·设a,b,c≠0且为整数。若c|a,c|b,则c|(sa+tb)。（若c|a,c|b，则c整除a,b的任意线性组合）
证明：s*a+t*b = s*θ*c+t*β*
c = c*(s* θ
+t* β
)
·设a,b,c为三个不全为零的整数，如果a = q*b+c，其中q为整数，则（a,b） = (b,c)
证明：d = (a,b) ，=> d|a，d|b =>d|(a+(-q)b) =>d|c ,d|b =>d为b,c的公因数 =>d≤d1
d1 = (b,c),同理得：d1≤d。=>d1=d =>(a,b) = (b,c)
广义欧几里得除法：
欧几里得算法又称辗转相除法，是指用于计算两个非负整数a，b的最大公约数。应用领域有数学和计算机两个方面。计算公式gcd(a,b) = gcd(b,a mod b)。欧几里得算法是用来求两个正整数最大公约数的算法。古希腊数学家欧几里得在其著作《The Elements》中最早描述了这种算法,所以被命名为欧几里得算法。扩展欧几里得算法可用于RSA加密等领域。
假如需要求 1997 和 615 两个正整数的最大公约数,用欧几里得算法，是这样进行的：
1997 / 615 = 3 (余 152)
615 / 152 = 4(余7)
152 / 7 = 21(余5)
7 / 5 = 1 (余2)
5 / 2 = 2 (余1)
2 / 1 = 2 (余0)
至此，最大公约数为1
以除数和余数反复做除法运算，当余数为 0 时，取当前算式除数为最大公约数，所以就得出了 1997 和 615 的最大公约数 1。
（2）代码分析：
广义欧几里得除法求最大公因数
long GCD(long a, long b)
{
long u, v, t, x;
if (a < 0) {
if (a < -NTL_MAX_LONG) ResourceError("GCD: integer overflow");
a = -a; //判断输入的长整型是否溢出
}
if (b < 0) {
if (b < -NTL_MAX_LONG) ResourceError("GCD: integer overflow");
b = -b; //判断输入的长整型是否溢出
}
if (b==0)
x = a;
else {
u = a;
v = b;
do {
t = u % v; //a = qb+r欧几里得核心算法（详见上知识储备-广义欧几里得除法）
u = v;
v = t;
} while (v != 0);
x = u;
}
return x;
}贝祖公式（广义欧几里得的逆方法）
利用广义欧几里得的算法步骤一步步回溯，就可以找到这样一组（x,y）（详细代码见@元解~殇怀）
·贝祖公式：
求s,t的一种方法也是简单方法：是利用广义欧几里得除法先得到（a,b) 即a,b的最大公因数，然后回代得到整数s,t。
2.2.4关于模运算
数学基础：乘法逆元求解
如下例：
乘法逆元：已知A和N互素，则存在一个整数B，使得A*B=1(mod N)
利用广义欧几里得除法（辗转相除法）对余数N进行辗转相除，然后将商逆序排列（如上图第一行黑色字），然后利用贝祖公式的变形求出上图第三行的相关数据（其中第三行第一的数字永远为1，第二个数字和第一行的第一个数字一样），如：5 = 1+1*4 ；9 = 4+5*1 ；14 = 5+9*1 ；23 = 9+14*1 ；37 = 14+23*1 ……
最后求出550即为550关于模1769的乘法逆元，即（550*550）mod 1769 = 1。
long InvModStatus(long& x, long a, long n) //求a关于n的乘法逆元，并将其赋给X。同时返回a和ns是否互素。
{
long d, s, t;
XGCD(d, s, t, a, n);
if (d != 1) {
x = d;
return 1;
}
else {
if (s < 0)
x = s + n; //规定乘法逆元要是正的，如果求出来不是正数，则要加上一倍的n转成正数输出
else
x = s;
return 0;
}
}
long InvMod(long a, long n)
{
long d, s, t;
XGCD(d, s, t, a, n);
if (d != 1) {
InvModError("InvMod: inverse undefined");
}
if (s < 0)
return s + n;
else
return s;
}
long PowerMod(long a, long ee, long n)
{
long x, y;
unsigned long e;
if (ee < 0)
e = - ((unsigned long) ee);
else
e = ee;
x = 1;
y = a;
while (e) {
if (e & 1) x = MulMod(x, y, n); //算法加速，一次循环，两次计算乘积，通过右移指令和按位取与指令做到快速的循环迭代，减少循环次数，加快运算速度。
//利用公式：((a{x} mod n)*(a{y} mod n))mod n = a{x+y} mod n ---证明见后
y = MulMod(y, y, n);
e = e >> 1;
}
if (ee < 0) x = InvMod(x, n);
return x;
}证明：((a{x} mod n)*(a{y} mod n))mod n = a{x+y} mod n
原式左侧 = ((a{x}-k1*n)*(a{y}-k2*n))mod n
= (a{x+y}-n(k2* a{x}+k1* a{y})+n{2}*k1*k2)mod n
= (a{x+y}) mod n