对Miracl大数运算库中函数xgcd()的理解与思考
背景
在完成C语言实现中国剩余定理的实践课题的代码编写过程中，xgcd()函数的理解与使用方法至关重要。
MIracl用户手册中是这样介绍 xgcd() 的：
int xgcd (big x, big y, big xd, big yd, big z)
Calculates extended Greatest Common Divisor of two big numbers. Can be used to calculate modular
inverses. Note that this routine is much slower than a mad() operation on numbers of similar size.Parameters:
←x
→y
→xd
→yd
→z = gcd(x, y) = (x * xd) + (y * yd)Returns:
GCD as integer, if possible, otherwise MR_TOOBIG
Precondition:
If xd and yd are not distinct, only xd is returned. The GCD is only returned if z distinct from both xd
and ydExample:
xgcd(x, p, x, x, x,); // x = 1/x mod p (p is prime)
从中，我们可以看出，xgcd()本身实现了扩展的Euclid算法，同时，也可被用作求模逆。下面对两种功能分别作以说明。
xgcd() 实现扩展的Euclid算法
扩展的Euclid算法，是指在通过Euclid算法及辗转相除法求出整数x,y的最大公约数 gcd(x,y) 后，对上述一系列带余除法的过程量进行整理从而得到 a*x + b*y = gcd(x,y) 这一等式的过程，即获得最大公约数 gcd(x,y) 关于整数 x 和 y 的整系数线性表出。
值得注意的是，当 x, y > 0 时，因为 gcd(x,y) <= x, y ，所以等式 a*x + b*y = gcd(x,y) 中的 a, b 两个系数必异号，及一正一负，这一点之后会有体现。
当使用 xgcd() 求 x, y 经扩展的Euclid算法后的结果时，输入参数 x, y 需预先赋值，而参数 xd, yd, z 则只需填入用来接收函数返回值的变量，返回值分别是上述等式 a*x + b*y = gcd(x,y) 中的 a, b 和 gcd(x,y)，这一用法并不难理解。
然而，在测试过程中，却出现了问题。
输入下列测试数据得到相应运算结果：
//测试A
//Input
x = 13
y = 17
xgcd(x, y, xd, yd, z)
//Output
xd = 4
yd = 3
z = 1
结果中， z = 1 不难理解， 因为 13 与 17 互素，因此 z = gcd(13,17) = 1 。然而 xd 和 yd 却让人迷惑。
根据上文对于扩展Euclid算法的描述，xd 和 yd 的应分别为等式 a*x + b*y = gcd(x,y) 中的 a 和 b，然而带入发现，13*4 + 17*3 != 1，经过观察可以看出，存在等式 13*4 - 17*3 = 1，即似乎存在等式 x*xd - y*yd = z 。
鉴于这一发现在用户手册中并未说明，为了验证其真实性，我们将 x, y 对调，得到如下输出结果：
//测试B
//Input
x = 17
y = 13
xgcd(x, y, xd, yd, z)
//Output
xd = 10
yd = 13
z = 1
可以看出，存在等式 17*10 - 13*13 = 1，证实了先前的猜想：即函数 xgcd(x,y,xd,yd,z) 运算后得到的结果满足 xd*x - yd*y = z 。
这一发现不由得让人有些匪夷所思，这一重要特性为什么没有在用户手册中说明呢？考虑到Miracl使用手册在先前的使用过程中确实存在诸多表述不清等现象，我特地翻看了用来编译 Miracl.lib 的C源码文件中的 mrxgcd.c ，终于在其中的函数注释中得到了答案：
int xgcd(_MIPD_ big x,big y,big xd,big yd,big z)
{ /* greatest common divisor by Euclids method *
* extended to also calculate xd and yd where *
* z = x.xd + y.yd = gcd(x,y) *
* if xd, yd not distinct, only xd calculated *
* z only returned if distinct from xd and yd *
* xd will always be positive, yd negative */
}
xd will always be positive, yd negative
注释中清楚地说明了：系数 xd 的符号为正，而系数 yd 的符号为负。
xgcd() 实现模逆运算
对于等式 a*x = 1 (mod m) ，我们称整数 a 为整数 x 模 m 的逆元，其中要求 gcd(a,m) = 1 ，即 a 与 m 互素。
值得注意的是，尽管在用户手册中写明了 xgcd() 可被用于实现模逆，然而具体的参数含义以及与扩展Euclid算法这一用法的区别却并未明确指出，只在最后的 Example 栏中写了如下样例：
xgcd(x, p, x, x, x,); // x = 1/x mod p (p is prime)
因此，根据初步猜测，我做了一下尝试：
//测试C
//Input
x = 13
y = 17
xgcd(x, y, z, z, z)
//Output
z = 4
//测试D
//Input
x = 17
y = 13
xgcd(x, y, z, z, z)
//Output
z = 10
根据测试结果可以看出，输出 z 确实为 x 在模 y 下的逆， 即 z*x = 1 (mod y) 。
至此，我们基本确定，按照 xgcd(x, y, z, z, z) 的形式输入参数，返回值 z 即为 x 在模 y 下的逆。
联系与思考
经过上述两种用法的分析与总结，我们有这样的思考：
模逆运算与扩展Euclid算法共用一个函数，而根据所学知识，由扩展Euclid算法得到的等式 a*x + b*y = gcd(x,y) 中，若 gcd(x,y) = 1，即 x, y 互素，那么对等式两边同时模 y ，则能得到 a*x = 1 (mod y)，a 就为模 y 下 x 的逆，而反之，如果对等式两边模 x ，b 却并非是模 x 下 y 的逆，因为根据上文可知，xgcd() 的返回值 b 实际上自带负号。从测试 A、B、C、D 四次的测试结果中我们可以看到，A 中 xd 与 C 中的 z 相同，均为 4 ，也就是 13 模 17 的逆，而A中的 yd 并非 17 模 13 的逆，B中的 xd 才是。因此，求谁的逆，就需要把谁赋值给 x ，如此得出的 xd 才为模 y 的逆。
由于在 xgcd() 扩展的Euclid算法的用法中 xd 的值，与模逆运算用法中 z 的值相等，这不由得让我们好奇 xgcd() 函数的内部实现，即如何通过输入参数的不同来提供不同的函数功能的。
用户手册有这样一段模糊的表述：
If xd and yd are not distinct, only xd is returned. The GCD is only returned if z distinct from both xd
and yd
意思是”如果 xd 与 yd 相同，则只返回 xd，如果 z 与 xd 和 yd 均不相同则返回最大公约数 gcd”。
同时我们在 mrxgcd.c 中的函数注释中看到一段更为清晰的表述：
* if xd, yd not distinct, only xd calculated *
* z only returned if distinct from xd and yd *
进一步说明了“当 xd 与 yd 相同时，只计算 xd ，当 z 与 xd 和 yd 均不相同时才会被返回”
最后，我们再结合函数源码（片段）：
if (xd!=yd)
{
negify(x,mr_mip->w2);
mad(_MIPP_ mr_mip->w2,mr_mip->w3,mr_mip->w1,y,mr_mip->w4,mr_mip->w4);
copy(mr_mip->w4,yd);
}
copy(mr_mip->w3,xd);
if (z!=xd && z!=yd) copy(mr_mip->w1,z);
MR_OUT
return (size(mr_mip->w1));
我们可以得到如下结论：
xgcd() 函数首先对输入的参数 x, y 做Euclid算法，得出 gcd(x,y) 并保留辗转相除过程中的所有过程信息，之后进行判断，根据传入参数 xd 和 yd 是否相同来选择性地为 yd 赋值；接着在外部对 xd 进行赋值，再而判断参数 z 是否与 xd 和 yd 均不相同，若不同则再为 z 赋值。
至此，我们不难理解 xgcd() 如何同时实现两种功能了：进行Euclid算法后，函数在必然给参数 xd 赋值的前提下，根据后三个参数的不同，有选择性地为 yd、z 赋值，从而实现:
- xd、yd、z 均不同时，函数 xgcd() 给引用 xd、yd、z 赋值，xd、yd、z 满足 x*xd - y*yd = z = gcd(x,y)。
- xd 与 yd 相同，z 与 xd 不同时，函数 xgcd() 给引用 xd 和 z 赋值，xd、z 满足 x*xd - y*yd = z = gcd(x,y)。
- xd 与 yd、z 相同时，函数 xgcd()给引用 xd 赋值，xd 满足 x*xd - y*yd = z = gcd(x,y)，也即满足 x*xd = 1 (mod y)。