被涯旅购上述每次divstep迭代可以表示为矩阵乘法对向量[f g]和[d e]应用如下的转换矩阵1/2*timage其中根据迭代中分支的不同(u, v, q, r)分别有不同的取值(0, 2, -1, 1)(2, 0, 1, 1)(2, 0, 0, 1)在while每次迭代中都会用转换矩阵乘于相应的向量而在转换矩阵中又包含除于2的操作这样整个while循环会包含除于2N的操作对于f和g来说算法已经保证了每次迭代的除法可以通过移位来实现而对于d和e来说这些个除法是非常费时的操作所以可以在每次迭代的矩阵乘法时不执行除于2的操作而把该除法操作留到最后步骤执行除于2N操作由此得到以下函数复制代码def divsteps_n_matrix(delta, f, g):Compute delta and transition matrix t after N divsteps (multiplied by 2^N).u, v, q, r 1, 0, 0, 1 # start with identity matrixfor _ in range(N):if delta 0 and g 1:delta, f, g, u, v, q, r 1 - delta, g, (g - f) // 2, 2*q, 2*r, q-u, r-velif g 1:delta, f, g, u, v, q, r 1 delta, f, (g f) // 2, 2*u, 2*v, qu, rvelse:delta, f, g, u, v, q, r 1 delta, f, (g ) // 2, 2*u, 2*v, q , rreturn delta, (u, v, q, r)复制代码后续还需进行如下计算image由推导过程可知最后uvqr组成的矩阵乘于初始向量[f g]以后得出的向量每个元素都是2N的倍数可由初始值实际验算下所以有以下的更新函数复制代码def update_fg(f, g, t):Multiply matrix t/2^N with [f, g].u, v, q, r tcf, cg u*f v*g, q*f r*g# (t / 2^N) should cleanly apply to [f,g] so the result of t*[f,g] should have N zero# bottom bits.assert cf % 2**N 0assert cg % 2**N 0return cf N, cg N复制代码对于d和e来说上述结论并不成立这里需要一个和div2函数类似的在模M下除于2N的函数借助该函数可实现de的更新复制代码def div2n(M, Mi, x):Compute x/2^N mod M, given Mi 1/M mod 2^N.assert (M * Mi) % 2**N 1# Find a factor m such that m*M has the same bottom N bits as x. We want:# (m * M) mod 2^N x mod 2^N# m mod 2^N (x / M) mod 2^N# m mod 2^N (x * Mi) mod 2^Nm (Mi * x) % 2**N# Subtract that multiple from x, cancelling its bottom N bits.x - m * M# Now a clean division by 2^N is possible.assert x % 2**N 0return (x N) % Mdef update_de(d, e, t, M, Mi):Multiply matrix t/2^N with [d, e], modulo M.u, v, q, r tcd, ce u*d v*e, q*d r*ereturn div2n(M, Mi, cd), div2n(M, Mi, ce)复制代码综合上述所有即可给出执行N divsteps版本的modinv函数复制代码def modinv(M, Mi, x):Compute the modular inverse of x mod M, given Mi1/M mod 2^N.assert M 1delta, f, g, d, e 1, M, x, 0, 1while g ! 0:# Compute the delta and transition matrix t for the next N divsteps (this only needs# (N1)-bit signed integer arithmetic).delta, t divsteps_n_matrix(delta, f % 2**N, g % 2**N)# Apply the transition matrix t to [f, g]:f, g update_fg(f, g, t)# Apply the transition matrix t to [d, e]:d, e update_de(d, e, t, M, Mi)return (d * f) % M复制代码这意味着在实践中我们将始终执行多个N个div步骤。这不是问题因为一旦g0进一步的divsteps就不再影响f、g、d或e只有δ继续增加这就能保证不管输入是何止算法运行步骤一致从而保证了算法时间一致性。3.3 secp256k1_modinv32源码分析结合上述内容分析下模逆的核心函数secp256k1_modinv32其源码如下复制代码1 /* Compute the inverse of x modulo modinfo-modulus, and replace x with it (constant time in x). */2 static void secp256k1_modinv32(secp256k1_modinv32_signed30 *x, const secp256k1_modinv32_modinfo *modinfo) {3 /* Start with d0, e1, fmodulus, gx, zeta-1. */4 secp256k1_modinv32_signed30 d {{0}};5 secp256k1_modinv32_signed30 e {{1}};6 secp256k1_modinv32_signed30 f modinfo-modulus;7 secp256k1_modinv32_signed30 g *x;8 int i;9 int32_t zeta -1; /* zeta -(delta1/2); delta is initially 1/2. */1011 /* Do 20 iterations of 30 divsteps each 600 divsteps. 590 suffices for 256-bit inputs. */12 for (i 0; i 20; i) {13 /* Compute transition matrix and new zeta after 30 divsteps. */14 secp256k1_modinv32_trans2x2 t;15 zeta secp256k1_modinv32_divsteps_30(zeta, f.v[0], g.v[0], t);16 /* Update d,e using that transition matrix. */17 secp256k1_modinv32_update_de_30(d, e, t, modinfo);18 /* Update f,g using that transition matrix. */19 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(f, 9, modinfo-modulus, -1) 0); /* f -modulus */20 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(f, 9, modinfo-modulus, 1) 0); /* f modulus */21 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(g, 9, modinfo-modulus, -1) 0); /* g -modulus */22 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(g, 9, modinfo-modulus, 1) 0); /* g modulus */2324 secp256k1_modinv32_update_fg_30(f, g, t);2526 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(f, 9, modinfo-modulus, -1) 0); /* f -modulus */27 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(f, 9, modinfo-modulus, 1) 0); /* f modulus */28 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(g, 9, modinfo-modulus, -1) 0); /* g -modulus */29 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(g, 9, modinfo-modulus, 1) 0); /* g modulus */30 }3132 /* At this point sufficient iterations have been performed that g must have reached 033 * and (if g was not originally 0) f must now equal /- GCD of the initial f, g34 * values i.e. /- 1, and d now contains /- the modular inverse. */3536 /* g 0 */37 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(g, 9, SECP256K1_SIGNED30_ONE, 0) 0);38 /* |f| 1, or (x 0 and d 0 and f modulus) */39 VERIFY_CHECK(secp256k1_modinv32_mul_cmp_30(f, 9, SECP256K1_SIGNED30_ONE, -1) 0 ||40 secp256k1_modinv32_mul_cmp_30(f, 9, SECP256K1_SIGNED30_ONE, 1) 0 ||41 (secp256k1_modinv32_mul_cmp_30(x, 9, SECP256K1_SIGNED30_ONE, 0) 0 42 secp256k1_modinv32_mul_cmp_30(d, 9, SECP256K1_SIGNED30_ONE, 0) 0 43 secp256k1_modinv32_mul_cmp_30(f, 9, modinfo-modulus, 1) 0));4445 /* Optionally negate d, normalize to [0,modulus), and return it. */46 secp256k1_modinv32_normalize_30(d, f.v[8], modinfo);47 *x d;48 }复制代码函数整体流程和上一节的介绍基本一致不再进行详细分析这里重点分析下函数中容易让人困惑的几个点。1. 参数zeta在代码一开始引入映射zeta -(delta1/2)则在delta取值1/2时zeta取值为-1以此为起点进行迭代在zeta引入后有以下迭代分支zeta0 且 g1zeta-zeta-2对应之前delta0且g1分支(zeta0 且 g1)或g为偶数zetazeta-1对应之前delta另外另种分支这样进行精心设计的好处是能用能用移位异或等操作完全替换比较操作在密码学领域消除分支判断实现运行恒定时间。secp256k1_modinv32_divsteps_30代码中zeta核心更新逻辑如下复制代码1 /* Compute conditional masks for (zeta 0) and for (g 1). */2 c1 zeta 31;3 mask1 c1;4 c2 g 1;5 mask2 -c2;6 mask1 mask2;7 /* Conditionally change zeta into -zeta-2 or zeta-1. */8 zeta (zeta ^ mask1) - 1;复制代码第2行处代码用于获取zeta符号当zeta为非负数时c10当zeta为负数时c10xFFFFFFFF所以对于mask1来说要么取值为0要么为0xFFFFFFFF。第3行当g为奇数时c21对应mask2为0xFFFFFFFF当g为偶数时c20对应mask2为0。第6行表明当g 奇数 且 zeta 0时mask1取值为0xFFFFFFFF这时第8行zeta(zeta ^ -1) - 1 (~zeta) - 1 -zeta -1 - 1 -zeta-2即对应之前说的分支1而其他情况下zeta zeta -1。2. 迭代次数在Bernstein–Yang的论文 “Fast constant-time gcd computation and modular inversion” 中给出对于n位整数输入需要的最大divsteps数量大约是2nC其中C是一个很小的常数。另外文中也说了在1.45n divsteps已经足够保证收敛但是现实里通常会再留安全裕度。在secp256k1的实现里选择600 divsteps作为固定上限方便分组且有足够裕度实际是只需要590 divsteps就能保证覆盖256位输入的最坏情况这里590的数字来自对 safegcd 算法的 worst-case 分析参考 secp256k1 开发者文档和论文中的 bound即在 256 位范围内任意输入的 GCD 计算最多 590 步就一定能结束。在库实现中最终进行了20x30分组这样的好处是批处理 (batching)每次不做 1 个 divstep而是“打包”多个这里是 30 个用矩阵乘法批量更新系数这样减少分支、提高效率。20 组 × 30 步 600 步实现里用固定循环展开保证 constant-time无论输入是什么都走一样的循环次数。所以这两个数字30 和 20并不是数学限制而是工程上选择的「批大小」和「循环次数」。