快速Euclidean算法

算法手册：Clegg L E, Mac Gabhann F. Molecular mechanism matters: Benefits of mechanistic computational models for drug development[J]. Pharmacological Research, 2015, 99: 149-154.

扩展欧几里得算法

给定$a,b$，计算$(x,y,d)$，使满足：$ax+by=d,d=gcd(a,b)$

易知：
$$\begin{array}{r}a\cdot 1+b\cdot 0=a\\ a\cdot 0+b\cdot 1=b\end{array}$$
以它为初值，令：
$$\begin{array}{rl}a\cdot x_1+b\cdot y_1=& u\cdots \cdots (\ast )\\ a\cdot x_2+b\cdot y_2=& v\cdots \cdots (\ast \ast )\end{array}$$
如果$u\ne 0$，那么执行欧几里得算法：

• 计算

$$q=\lfloor \frac{v}{u}\rfloor$$

• 令

$$\begin{array}{rl}(\ast )\to & (\ast \ast )\\ (\ast \ast )-q(\ast )\to & (\ast )\end{array}$$

• 即

$$\begin{array}{rl}(x_2,y_2,v)\leftarrow & (x_1,y_1,u)\\ (x_1,y_1,u)\leftarrow & (x_2-q{x}_1,y_2-q{y}_1,v-qu)\end{array}$$

注意数据依赖，实例化时要使用临时变量暂存。

直到$u=0$结束，有$a{x}_2+b{y}_2=v$，得出结果$(x=x_2,y=y_2,d=v)$

def exgcd(a,b):
	'xa + yb = d, (x,y,d)'

	if a>b:
		x1,y1,d1 = 1,0,a
		x2,y2,d2 = 0,1,b
	else:
		x1,y1,d1 = 0,1,b
		x2,y2,d2 = 1,0,a

	while d2 != 0:
		t = d1//d2
		x1,y1,d1 = x1-t*x2, y1-t*y2, d1-t*d2
		x1,y1,d1,x2,y2,d2 = x2,y2,d2,x1,y1,d1

	return x1,y1,d1

二进制欧几里得算法

给定$a,b$，计算$d$，使满足：$d=gcd(a,b)$

假设$a,b$都是偶数，那么
$$\gcd (a,b)=2\gcd (\frac{a}{2},\frac{b}{2})$$
因此，计算$e=2^k$，使得$a^{\prime }=\frac{a}{e},b^{\prime }=\frac{b}{e}$不全为偶数

那么
$$\begin{array}{rl}\gcd (a,b)=& e\cdot \gcd (a^{\prime },b^{\prime })\\ 2\nmid & gcd(a^{\prime },b^{\prime })\end{array}$$
我们只需计算$\gcd (a^{\prime },b^{\prime })$。令$u=a^{\prime },v=b^{\prime }$

若$u$是偶数且$v$是奇数，那么
$$gcd(u,v)=gcd(\frac{u}{2},v)$$
当$u\ne 0$，执行欧几里得算法：

• 若$u$是偶数，那么反复执行$u\leftarrow u/2$
• 若$v$是偶数，那么反复执行$v\leftarrow v/2$
• 当$u\ge v$，令$u\leftarrow u-v$；否则$v\leftarrow v-u$

当$u=0$时结束，此时$v=gcd(a^{\prime },b^{\prime })$，得出结果$d=e\cdot v$

注意，反复执行除二操作，可以用一次右移代替。

二进制求逆算法

对于奇素数$p$，给定$a\in Z_p^{\ast }$，计算$x=a^{-1}\mathrm{mod}p$

易知：
$$\begin{array}{r}a\cdot 1+p\cdot 0=a\\ a\cdot 0+p\cdot 1=p\end{array}$$
以它为初值，令：
$$\begin{array}{rl}a\cdot x_1+p\cdot y_1=& u\cdots \cdots (\ast )\\ a\cdot x_2+p\cdot y_2=& v\cdots \cdots (\ast \ast )\end{array}$$
我们只计算$x_1,x_2$，不计算$y_1,y_2$ (也不容易计算)

易知，$x\in Z_p^{\ast }$

如果$u\ne 1\wedge v\ne 1$，那么执行欧几里得算法：

• 若$u$是偶数，那么反复执行
  • $u\leftarrow u/2$
  • 若$x_1$是偶数，令$x_1\leftarrow x_1/2$；否则$x_1\leftarrow (x_1+p)/2$
• 若$v$是偶数，那么反复执行
  • $v\leftarrow v/2$
  • 若$x_2$是偶数，令$x_2\leftarrow x_2/2$；否则$x_2\leftarrow (x_2+p)/2$
• 当$u\ge v$，令$u\leftarrow u-v,x_1=x_1-x_2$；否则$v\leftarrow v-u,x_2=x_2-x_1$

当$u=1$时，$a\cdot x_1\equiv u\mathrm{mod}p$，得出结果$x=x_1$

当$v=1$时，$a\cdot x_2\equiv v\mathrm{mod}p$，得出结果$x=x_2$

如果初值设为$x_1=b\in Z_p$而非$x_1=1$，那么上述算法的结果为$b{a}^{-1}\mathrm{mod}p$，即：计算除法不再需要额外的乘法！