二进制扩展欧几里得算法(binary extended euclidean algorithm)

最近实现了下椭圆曲线Secp251k1，频繁用到的一个算法是求模逆元，也就是 $ax = 1 \mod p$ 这里,p是素数,我使用的前面所写的算法,把尾递归改成迭代。发现速度仍然比tomlib慢了10倍，内存也是大了很多，矩阵乘法局部变量太多，我猜测是这样。这里,看了他们源码,实现的是一个叫做 binary extended euclidean algorithm 的算法，在扩展欧几里得的基础上做了优化，不用每次做商求余，只需要每次右移。而移位操作的代价是很小的。

前文中扩展欧几里得算法中所描述的,最终是要找到 $ax + by = 1$的等式。

二进制扩展欧几里得算法用到了下面几个事实 ,gcd 最大公约数

a b 都是偶数,那么 gcd(a,b) = 2 gcd(a/2,b/2),
a 是偶数,b是奇数,那么 gcd(a,b) = gcd(a/2,b)
a b 都是奇数,那么 gcd(a,b) = gcd((a-b)/2,b)

根据这两个，我们可以不断的把 $\begin{bmatrix} a\\ p \end{bmatrix} = A \begin{bmatrix} c\\ x \end{bmatrix}$ 这里,c 是最大公约数,本例中,p是素数,c = 1,而扩展欧几里得算法可以表示为一系列的线性变换。这里以，求 12，13 的过程为例,所有运算都是在素数域13下,mod 13, 求2的逆元很简单,就是(13 + 1) / 2 = 7

a	p	步骤	状态	备注
12	13	初始状态	\| 1 ,0 \| \| 0 ,1 \|
6	13	a = a / 2	\| 7 ,0 \| \| 0 ,1 \|	矩阵列第一列同时乘以2的逆元 ,这里是(13 + 1) / 2 = 7
3	13	a = a / 2	\| 10 ,0\| \| 0 ,1 \|	矩阵列第一列同时乘以2的逆元7, 7 * 7 % 13 = 10
3	10	p = p - a	\| 10 ,4 \| \| 0 ,1 \|	第二列减去第一列, 0 - 10 % 13 = 3
3	5	p = p - a	\| 10 ,8 \| \| 0 ,7 \|	第2列除以2, 3 * 7 % 13 = 8
3	2	p = p - a	\| 10 ,11 \| \| 0 ,7 \|	第2列除以2, (8 -10) % 13 = 11
2	3	交换 a p	\| 11 ,10 \| \| 7 ,0 \|	交换 12 列
1	3	a = a/2	\| 12 ,10 \| \| 10 ,0 \|	11 * 7 % 13 = 12

那么有

\[\begin{bmatrix}1\\3\end{bmatrix} = \begin{bmatrix}12&10\\10&0\end{bmatrix} \cdot \begin{bmatrix} 12\\13 \end{bmatrix} \mod 13\]

所以 12 模13 逆元是 12,也就是他本身。

可以看到，矩阵$A$只有第一行起到了作用,实现过程中,我们可以忽略第二行，只需要两个变量就能表示中间状态。js 代码如下

避免使用乘法

经过测试之后,发现如果所有a/2 都写成乘以2的逆元,那么会运用到大量乘法操作,即使采用分治法,速度也不令人满意。下面就运用一个技巧快速计算$a/2 \mod p$

如果a是偶数,那么直接取a的一半 $a/2$
如果a是奇数,$a = 2n + 1$ , 那么 $a/2 = (2n + 1) * i = 2 * i * n + i = n + i \mod p$ 这样,就不需要乘法操作,直接移位和加法就能快速的完成运算。

i 是 2 的逆元, i = (p+1)/2, p 是奇数

 

function binaryExGcd(a,b){


    function div2(a,p,invers2){
        /***
         *   a/2 mod p ,
         *   如果 a 是偶数,那么直接 a/2
         *   如果 a 是奇数, 那么 a = 2n + 1
         *   a/2  =  a * I  =  2 * I * n + I = n + I 
         *   这里 I 是2 的逆元, I = (p + 1) / 2
         */
        if((a & 1n) == 0){
            return a>>1n;
        }else{
            const half = a >> 1n;
            return half + invers2
        }
    }


    if(a >= b){
        console.error("a must less than b")
        return [0,0]
    }
    if((b & 1n) == 0n){
        console.error(b, "is not odd number")
        return;
    }

    var invers2 = (b + 1n) >> 1n;
    var A = a , B = b 
    var u = 1n, v = 0n

    var gcd = 1n;
    while(true){
        console.log(A,B,u,v)
        if(A == 1n){
            break;
        }
        if(B == 0n){
            gcd *= A;
            return [gcd,u];
        }

        if (A > B){
            var c = A 
            A = B
            B = c 
            
            c = u 
            u = v 
            v = c 

            continue;
        }

        let isAOdd = (A & 1n) == 1n
        let isBOdd = (B & 1n) == 1n

        if(!isBOdd && ! isAOdd ){
            gcd <<= 1n;
        }
        
        if (isAOdd && isBOdd){
            B =  (B - A) >> 1n;
            v = (v - u) 
            v = div2(v,b,invers2);
            
            if(v < 0){
                v += b
            }
            
           
            


        }
        else if(!isAOdd){
            A >>= 1n;
            u = div2(u,b,invers2);
            if(u < 0){
                u += b
            }

        }
        else if(!isBOdd){
            B >>= 1n;
            v = div2(v,b,invers2);

            if(v < 0){
                v += b
            }
        }
    }
    u %= b 
    if(u < 0n){
        u += b
    }
    return [gcd,u];
}
 
var count  = 1;
while (count --)
{
    var a = BigInt(Math.ceil(Math.random() * 100000000))
    var b =  BigInt(Math.ceil(Math.random() * 100000000000))
    if( (b & 1n) == 0){
        b += 1n;
    }

    a = 311n 
    b = 997n

    var cc = binaryExGcd(a,b);
    var gcd = cc[0]
    var c = cc[1]
    if(gcd){
        console.log("result:",c ,gcd, (a * c )% b)
    }
    else{
    }
    
    if(gcd && (a * c )% b != gcd){
        console.log("error:",a,b,c ,gcd, (a * c )% b)
        break;
    }
}

module.exports = binaryExGcd

TODO

上面只处理了，p是奇数的情况，无法处理偶数的情况，（偶数，2没有逆元）

L 记事本