secp256k1算法详解四（关键点补充说明）

1 magnitude及normalized

由于当前许多项目都用到secp256k1库，比特币作为体量最大的数字货币项目，这里建议直接参考bitcoin-core提供的最新secp256k1源码。仍以field的10x26实现版本为例，相关定义如下：

1. /** This field implementation represents the value as 10 uint32_t limbs in base
2. *  2^26. */
3. typedef struct {
4.    /* A field element f represents the sum(i=0..9, f.n[i] << (i*26)) mod p,
5.     * where p is the field modulus, 2^256 - 2^32 - 977.
6.     *
7.     * The individual limbs f.n[i] can exceed 2^26; the field's magnitude roughly
8.     * corresponds to how much excess is allowed. The value
9.     * sum(i=0..9, f.n[i] << (i*26)) may exceed p, unless the field element is
10.     * normalized. */
11.     uint32_t n[10];
12.     /*
13.      * Magnitude m requires:
14.      *     n[i] <= 2 * m * (2^26 - 1) for i=0..8
15.      *     n[9] <= 2 * m * (2^22 - 1)
16.      *
17.      * Normalized requires:
18.      *     n[i] <= (2^26 - 1) for i=0..8
19.      *     sum(i=0..9, n[i] << (i*26)) < p
20.      *     (together these imply n[9] <= 2^22 - 1)
21.      */
22.     SECP256K1_FE_VERIFY_FIELDS
23. } secp256k1_fe;

复制代码

secp256k1_fe_add
显而易见，对于输出r来说m(r)infinity) {        if (rzr != NULL) {            secp256k1_fe_set_int(rzr, 1);        }        *r = *a;        return;    }    secp256k1_fe_sqr(&z22, &b->z);    secp256k1_fe_sqr(&z12, &a->z);    secp256k1_fe_mul(&u1, &a->x, &z22);    secp256k1_fe_mul(&u2, &b->x, &z12);    secp256k1_fe_mul(&s1, &a->y, &z22); secp256k1_fe_mul(&s1, &s1, &b->z);    secp256k1_fe_mul(&s2, &b->y, &z12); secp256k1_fe_mul(&s2, &s2, &a->z);    secp256k1_fe_negate(&h, &u1, 1); secp256k1_fe_add(&h, &u2);    secp256k1_fe_negate(&i, &s2, 1); secp256k1_fe_add(&i, &s1);    if (secp256k1_fe_normalizes_to_zero_var(&h)) {        if (secp256k1_fe_normalizes_to_zero_var(&i)) {            secp256k1_gej_double_var(r, a, rzr);        } else {            if (rzr != NULL) {                secp256k1_fe_set_int(rzr, 0);            }            secp256k1_gej_set_infinity(r);        }        return;    }    r->infinity = 0;    secp256k1_fe_mul(&t, &h, &b->z);    if (rzr != NULL) {        *rzr = t;    }    secp256k1_fe_mul(&r->z, &a->z, &t);    secp256k1_fe_sqr(&h2, &h);    secp256k1_fe_negate(&h2, &h2, 1);    secp256k1_fe_mul(&h3, &h2, &h);    secp256k1_fe_mul(&t, &u1, &h2);    secp256k1_fe_sqr(&r->x, &i);    secp256k1_fe_add(&r->x, &h3);    secp256k1_fe_add(&r->x, &t);    secp256k1_fe_add(&r->x, &t);    secp256k1_fe_add(&t, &r->x);    secp256k1_fe_mul(&r->y, &t, &i);    secp256k1_fe_mul(&h3, &h3, &s1);    secp256k1_fe_add(&r->y, &h3);    SECP256K1_GEJ_VERIFY(r);}[/code]secp256k1_gej_add_var首先以输入参数a，b中的x，y，z都为magnitude=1（normalized-like）为前提，在此基础上给出按代码顺序的表：
步代码（操作）输入 m输出 m（上限）说明1secp256k1_fe_sqr(&z22, &b->z)b->z:11fe_sqr → 规约为 12secp256k1_fe_sqr(&z12, &a->z)a->z:11 3secp256k1_fe_mul(&u1, &a->x, &z22)1,11fe_mul → 规约为 14secp256k1_fe_mul(&u2, &b->x, &z12)1,11 5secp256k1_fe_mul(&s1, &a->y, &z22)1,11 6secp256k1_fe_mul(&s1, &s1, &b->z)1,11 7secp256k1_fe_mul(&s2, &b->y, &z12)1,11 8secp256k1_fe_mul(&s2, &s2, &a->z)1,11 9secp256k1_fe_negate(&h, &u1, 1)u1:12fe_negate → m_out = 1 + 1 = 210secp256k1_fe_add(&h, &u2)h:2, u2:13fe_add: 2 + 1 = 311secp256k1_fe_negate(&i, &s2, 1)s2:12 12secp256k1_fe_add(&i, &s1)i:2, s1:13 13if (secp256k1_fe_normalizes_to_zero_var(&h)) ...——特殊分支（不进入一般流程）14secp256k1_fe_mul(&t, &h, &b->z)h:3, b.z:11fe_mul → 规约为 115if (rzr != NULL) *rzr = t;t:11 16secp256k1_fe_mul(&r->z, &a->z, &t)a.z:1, t:11 17secp256k1_fe_sqr(&h2, &h)h:31fe_sqr → 规约为 118secp256k1_fe_negate(&h2, &h2, 1)h2:12 19secp256k1_fe_mul(&h3, &h2, &h)h2:2, h:31fe_mul → 规约为 120secp256k1_fe_mul(&t, &u1, &h2)u1:1, h2:21 21secp256k1_fe_sqr(&r->x, &i)i:31 22secp256k1_fe_add(&r->x, &h3)r->x:1, h3:12 23secp256k1_fe_add(&r->x, &t)r->x:2, t:13 24secp256k1_fe_add(&r->x, &t)r->x:3, t:14 25secp256k1_fe_add(&t, &r->x)t:1, r->x:45（t 被覆盖为 t + r->x）26secp256k1_fe_mul(&r->y, &t, &i)t:5, i:31fe_mul → 规约为 127secp256k1_fe_mul(&h3, &h3, &s1)h3:1, s1:11 28secp256k1_fe_add(&r->y, &h3)r->y:1, h3:12 29函数返回r->x:4, r->y:2, r->z:1—这些是按最小合法上界得到的值（未再 normalize 前）由表可知对于secp256k1_gej_add_var函数，虽然运行过程中，零时变量可能会出现相对较高的magnitude值，但最终返回值r->x的magnitude值是满足SECP256K1_GE_X_MAGNITUDE_MAX值的，同理y，z也满足SECP256K1_GEJ_Y_MAGNITUDE_MAX和SECP256K1_GEJ_Z_MAGNITUDE_MAX限制。
3 大数求逆

最新版大数求逆函数实现原型如下：

1. #define SECP256K1_GE_X_MAGNITUDE_MAX  4
2. #define SECP256K1_GE_Y_MAGNITUDE_MAX  3
3. #define SECP256K1_GEJ_X_MAGNITUDE_MAX 4
4. #define SECP256K1_GEJ_Y_MAGNITUDE_MAX 4
5. #define SECP256K1_GEJ_Z_MAGNITUDE_MAX 1

复制代码

3.1 高层概览

1. 主要步骤

函数做的是对域元素x求模逆，即求r = x-1 mod p，这里p = 2^256 - 2^32 - 977，其实现流程如下：
1）先把tmp（x）规范化，保证tmp是库约定的limb范围下唯一表示；
2）把规范化的tmp（10x26-bit limbs）转换成一种signed30的中间表示（若干个30-bit的有符号limbs），以便后面用32-bit算法实现逆元算法；
3）在signed30表示上运行secp256k1_modinv32——这是一个针对30-bit/32-bit limbs优化、并做过constant-time处理的模逆实现，基于safegcd/division-steps（改进的欧几里得/半GCD）的思想；
4）把得到的signed30结果再转换回库的secp256k1_fe（10x26表示），并写入r（该转换同时完成必要的约减/规范化）。
2. divsteps算法

divsteps是一种高效计算GCD（Greatest Common Divisor最大公约数）的方法，它通过连续多次除法步骤减少迭代次数，特别适合硬件实现和大数计算。它是GCD算法的优化版本，主要特点：

• 利用二进制表示的优势，通过右移操作快速消除因子2
  
• 一次迭代处理多个 "减法 - 移位" 步骤，减少循环次数
  
• 使用比较和差值运算替代昂贵的除法操作
  

求u和v最大公约数算法步骤如下：
1）消除因子2

• 统计并移除两数中的所有因子2
  
• 设a = u / 2^sa，b = v / 2^sb（这里sa是u可整除2的最大次数，sb是v可整除2的最大次数）
  

2）divsteps迭代

• 比较a和b
  
• 计算差值d = | a - b |
  
• 消除d中的所有因子2得到d'
  
• 用min(a, b)和d'替换(a, b)
  
• 重复直到a == b
  

3）得出最终结果
gcd(u, v) = b * 2^min(sa, sb)
算法基于以下假设：如果gcd(u, v)是最大公约数，那么它可以分成两部分的乘积，一部分是2的整数次幂，另一部分是非2的倍数。步骤1其实是获取u、v最大公因数中能2整数次幂部分，即2min(sa, sb)，步骤2获取非2的倍数部分。步骤1比较好理解，这里分析下步骤2，假设得出最后一步结论a==b时，具体值为a'，b'，则a'必然满足a' - b' = b'*2x，则有a' = b'*(1+2x)，即a'是b'的倍数，则在上一步必有d'=b'*2y = a'' - a'，由此a'' = a' + b'*2y，即a''也为b的倍数，由此类推，可知最一开始的a和b闭然都是b'的倍数。由以上推理可知最终gcd(u, v) = b' * 2^min(sa, sb)。
以下是算法示例代码：

1. pow(2,256)-1)*2*1=0x1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe
2. pow(2,256)-1)*2*2=0x3fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffc
3. pow(2,256)-1)*2*3=0x5fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffa
4. pow(2,256)-1)*2*4=0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff8

复制代码

secp256k1_fe_to_signed30该函数实现很清晰，就是把原有的10个无符号26bit数重新按30bit窗口打包到一组32-bit有符号整数中（共9个元素，每个元素保存30位有效位，并且用int32_t保存以允许出现负数中间值）。
步骤3调用secp256k1_modinv32函数求模逆，该步骤是算法的核心，后续详细解析。
步骤4调用secp256k1_fe_from_signed30函数将30bit格式逆元再转换到标准10x26bit形式。
1. 扩展欧几里得算法

扩展欧几里得定理：对于任意整数a和b，必然存在整数x和y满足贝祖等式：a·x + b·y = gcd(a, b)，且gcd(a, b)是能表示为a·x + b·y形式的最小正整数。
1）基础引理：欧几里得算法的余数性质
欧几里得算法通过反复求余计算GCD：
gcd(a, b) = gcd(b, a mod b)
其中a mod b = a - b*⌊a/b⌋（⌊⌋为向下取整），且0 ≤ a mod b < |b|。
2）数学归纳法证明贝祖等式存在性
归纳基础：
当b = 0 时，gcd(a, 0) = a，此时取 x = 1，y = 0（这里y可以取任意值，只不过取0后续计算最简单），显然满足 a·1 + 0·0 = a = gcd(a, 0)。
归纳步骤：
假设对于 (b, a mod b) 存在整数 x' 和 y' 满足：b·x' + (a mod b)·y' = gcd(b, a mod b)
根据余数定义 a mod b = a - b·⌊a/b⌋，代入上式：b·x' + (a - b·⌊a/b⌋)·y' = gcd(a, b)（因 gcd(a, b) = gcd(b, a mod b)）
整理得：
a·y' + b·(x' - ⌊a/b⌋·y') = gcd(a, b)，令 x = y'，y = x' - ⌊a/b⌋·y'，则有：a·x + b·y = gcd(a, b)
因此，若 (b, a mod b) 存在解，则 (a, b) 也存在解。由数学归纳法，对任意整数 a, b 均存在这样的 x, y。
3）最小性证明
设 d = gcd(a, b)，则 d 整除 a 和 b，因此 d 整除 a·x + b·y 的所有可能结果。即任何能表示为 a·x + b·y 的整数都是 d 的倍数，故 d 是其中最小的正整数。
其实以上结论是基于一个前提：对于任意两个整数a和b，在应用欧几里得算法若干步后，都能转化成将gcd(a, b)转换成gcd(b', 0)的形式，这个结论可以自己查阅相关整理过程，又因为gcd(a, b) = gcd(-a, -b)，所以该结论对于所有非零整数都成立。
接下来看看在算法层面如何推导系数x和y：
首先，对于gcd(a, b)中的a和b来说都可以看作是各次除法迭代的“余数r”，即gcd(a, b) = gcd(b, a%b))，如果将a%b看作是newr，b看作是r，a看作是oldr，则有：oldr = q⋅r + newr，这里q是整数商a//b。
接下来，因为求系数x和y其实也是一个迭代推导的过程，假设当前有：

可以先把最一开始oldr和r的由系数表示出来：oldr = 1*a + 0*b，r = 0*a + 1*b，即存在初始余数对：(oldr, r) = (a, b)，初始a的系数对：(oldx, x) = (1, 0)，初始b的系数对：(oldy, y) = (0, 1)满足方程。
再接下来，继续进行迭代q = oldr//r，newr = oldr - q*r，将(1)中oldr和(2)中的r代入该等式并进行整理可得：

由此可得：newx = oldx - q*x，newy = oldy - q*y，即newr仍是a，b的线性组合，且newx，newy是相应系数。以下是完整算法伪代码：

1. SECP256K1_INLINE static void secp256k1_fe_negate(secp256k1_fe *r, const secp256k1_fe *a, int m) {
2.     r->n[0] = 0x3FFFC2FUL * 2 * (m + 1) - a->n[0];
3.     r->n[1] = 0x3FFFFBFUL * 2 * (m + 1) - a->n[1];
4.     r->n[2] = 0x3FFFFFFUL * 2 * (m + 1) - a->n[2];
5.     r->n[3] = 0x3FFFFFFUL * 2 * (m + 1) - a->n[3];
6.     r->n[4] = 0x3FFFFFFUL * 2 * (m + 1) - a->n[4];
7.     r->n[5] = 0x3FFFFFFUL * 2 * (m + 1) - a->n[5];
8.     r->n[6] = 0x3FFFFFFUL * 2 * (m + 1) - a->n[6];
9.     r->n[7] = 0x3FFFFFFUL * 2 * (m + 1) - a->n[7];
10.     r->n[8] = 0x3FFFFFFUL * 2 * (m + 1) - a->n[8];
11.     r->n[9] = 0x03FFFFFUL * 2 * (m + 1) - a->n[9];
12. }

复制代码

以a=240，b=46为例，给出算法迭代过程：

1. SECP256K1_INLINE static void secp256k1_fe_add(secp256k1_fe *r, const secp256k1_fe *a) {
2.     SECP256K1_FE_VERIFY(r);
3.     SECP256K1_FE_VERIFY(a);
4.     VERIFY_CHECK(r->magnitude + a->magnitude <= 32);
6.     secp256k1_fe_impl_add(r, a);
7.     r->magnitude += a->magnitude;
8.     r->normalized = 0;
10.     SECP256K1_FE_VERIFY(r);
11. }

复制代码

2. 二进制GCD算法变种

二进制GCD算法有以下演进时间线：

1. static void secp256k1_gej_add_var(secp256k1_gej *r, const secp256k1_gej *a, const secp256k1_gej *b, secp256k1_fe *rzr) {
2.     /* 12 mul, 4 sqr, 11 add/negate/normalizes_to_zero (ignoring special cases) */
3.     secp256k1_fe z22, z12, u1, u2, s1, s2, h, i, h2, h3, t;
4.     SECP256K1_GEJ_VERIFY(a);
5.     SECP256K1_GEJ_VERIFY(b);
7.     if (a->infinity) {
8.         VERIFY_CHECK(rzr == NULL);
9.         *r = *b;
10.         return;
11.     }
12.     if (b->infinity) {
13.         if (rzr != NULL) {
14.             secp256k1_fe_set_int(rzr, 1);
15.         }
16.         *r = *a;
17.         return;
18.     }
20.     secp256k1_fe_sqr(&z22, &b->z);
21.     secp256k1_fe_sqr(&z12, &a->z);
22.     secp256k1_fe_mul(&u1, &a->x, &z22);
23.     secp256k1_fe_mul(&u2, &b->x, &z12);
24.     secp256k1_fe_mul(&s1, &a->y, &z22); secp256k1_fe_mul(&s1, &s1, &b->z);
25.     secp256k1_fe_mul(&s2, &b->y, &z12); secp256k1_fe_mul(&s2, &s2, &a->z);
26.     secp256k1_fe_negate(&h, &u1, 1); secp256k1_fe_add(&h, &u2);
27.     secp256k1_fe_negate(&i, &s2, 1); secp256k1_fe_add(&i, &s1);
28.     if (secp256k1_fe_normalizes_to_zero_var(&h)) {
29.         if (secp256k1_fe_normalizes_to_zero_var(&i)) {
30.             secp256k1_gej_double_var(r, a, rzr);
31.         } else {
32.             if (rzr != NULL) {
33.                 secp256k1_fe_set_int(rzr, 0);
34.             }
35.             secp256k1_gej_set_infinity(r);
36.         }
37.         return;
38.     }
40.     r->infinity = 0;
41.     secp256k1_fe_mul(&t, &h, &b->z);
42.     if (rzr != NULL) {
43.         *rzr = t;
44.     }
45.     secp256k1_fe_mul(&r->z, &a->z, &t);
47.     secp256k1_fe_sqr(&h2, &h);
48.     secp256k1_fe_negate(&h2, &h2, 1);
49.     secp256k1_fe_mul(&h3, &h2, &h);
50.     secp256k1_fe_mul(&t, &u1, &h2);
52.     secp256k1_fe_sqr(&r->x, &i);
53.     secp256k1_fe_add(&r->x, &h3);
54.     secp256k1_fe_add(&r->x, &t);
55.     secp256k1_fe_add(&r->x, &t);
57.     secp256k1_fe_add(&t, &r->x);
58.     secp256k1_fe_mul(&r->y, &t, &i);
59.     secp256k1_fe_mul(&h3, &h3, &s1);
60.     secp256k1_fe_add(&r->y, &h3);
62.     SECP256K1_GEJ_VERIFY(r);
63. }

复制代码

接下来重点分析带有delta状态变量的版本，算法如下：

1. 1 static void secp256k1_fe_impl_inv(secp256k1_fe *r, const secp256k1_fe *x) {
2. 2     secp256k1_fe tmp = *x;
3. 3     secp256k1_modinv32_signed30 s;
4. 4
5. 5     secp256k1_fe_normalize(&tmp);
6. 6     secp256k1_fe_to_signed30(&s, &tmp);
7. 7     secp256k1_modinv32(&s, &secp256k1_const_modinfo_fe);
8. 8     secp256k1_fe_from_signed30(r, &s);
9. 9 }

复制代码

该算法要求第一个参数f必须是奇数（第3行），并引入了状态变量delta，用于指导算法选择不同的约减策略，避免陷入低效循环。该算法证明分两部分：
1）先证明算法保持 gcd 不变（正确性不破坏）；
2）然后给出一个势函数（potential），并用它证明在每次若干步内势函数会严格下降，从而不可能无限迭代——也就保证终止（收敛）。同时指出 delta 在防止振荡 / 保证下降中的关键作用。
第1部分结论很容易得出，在结合f，g奇偶情况下gcd(g, (g - f)//2)，gcd(f, (g + f)//2)，gcd(f, g//2)这三种情况显然和gcd(f, g)是一样的。第2部分的证明比较复杂，涉及到势函数相关理论，知识盲区，请自行查阅相关资料，这里不再详细说明。
总之，相比如下未引入delta的算法一（在某些情况下会不收敛），上述算法是能保证收敛的。

1. #include <stdint.h>
2. #include <stdlib.h>
4. // 计算二进制表示中后缀0的个数
5. static int count_trailing_zeros(uint64_t x) {
6.     if (x == 0) return 64;
7.     int cnt = 0;
8.     while ((x & 1) == 0) {
9.         cnt++;
10.         x >>= 1;
11.     }
12.     return cnt;
13. }
15. uint64_t gcd_divsteps(uint64_t u, uint64_t v) {
16.     if (u == 0) return v;
17.     if (v == 0) return u;
19.     // 步骤1: 移除所有因子2
20.     int sa = count_trailing_zeros(u);
21.     int sb = count_trailing_zeros(v);
22.     u >>= sa;
23.     v >>= sb;
25.     // 步骤2: divsteps迭代
26.     while (u != v) {
27.         if (u < v) { uint64_t t = u; u = v; v = t; }  // 确保u >= v
28.         
29.         uint64_t d = u - v;
30.         int sd = count_trailing_zeros(d);
31.         d >>= sd;  // 一次消除所有因子2
32.         
33.         u = d;
34.     }
36.     // 步骤3: 恢复共同的因子2
37.     return u << (sa < sb ? sa : sb);
38. }

复制代码

gcd_no_delta另外相比，如下未引入delta的算法二，f，g直接判断属于“输入依赖型分支”，不同输入会导致不同的分支走向，容易在硬件层面产生 “分支预测错误”（尤其在加密算法中，输入的随机性会放大这种问题）。而delta引入后，分支判断虽然存在，但更多依赖内部状态变量 delta 的符号（delta > 0 或 delta ≤ 0），而非输入 f 和 g 的直接比较，delta 的更新规则是固定的（1 - delta 或 1 + delta），其变化模式相对可预测，降低了对输入随机性的敏感度，更适合常数时间实现（密码学算法的关键要求）。尤其在大数运算情况下，a > b 这类比较操作需要完整的减法和符号判断（硬件层面是减法器 + 符号位检测），而delta的状态更新（1 - delta 或 1 + delta）是简单的算术操作（本质是 delta 在0和1之间交替，或递增递减），计算成本更低。

1. static void secp256k1_fe_to_signed30(secp256k1_modinv32_signed30 *r, const secp256k1_fe *a) {
2.     const uint32_t M30 = UINT32_MAX >> 2;
3.     const uint64_t a0 = a->n[0], a1 = a->n[1], a2 = a->n[2], a3 = a->n[3], a4 = a->n[4],
4.                    a5 = a->n[5], a6 = a->n[6], a7 = a->n[7], a8 = a->n[8], a9 = a->n[9];
6.     r->v[0] = (a0       | a1 << 26) & M30;
7.     r->v[1] = (a1 >>  4 | a2 << 22) & M30;
8.     r->v[2] = (a2 >>  8 | a3 << 18) & M30;
9.     r->v[3] = (a3 >> 12 | a4 << 14) & M30;
10.     r->v[4] = (a4 >> 16 | a5 << 10) & M30;
11.     r->v[5] = (a5 >> 20 | a6 <<  6) & M30;
12.     r->v[6] = (a6 >> 24 | a7 <<  2
13.                         | a8 << 28) & M30;
14.     r->v[7] = (a8 >>  2 | a9 << 24) & M30;
15.     r->v[8] =  a9 >>  6;
16. }

复制代码

bad_gcd3. 由GCD求模逆


对于质数M（大于2）来说求x（0 0 and g & 1:            delta, f, g, u, v, q, r = 1 - delta, g, (g - f) // 2, 2*q, 2*r, q-u, r-v        elif g & 1:            delta, f, g, u, v, q, r = 1 + delta, f, (g + f) // 2, 2*u, 2*v, q+u, r+v        else:            delta, f, g, u, v, q, r = 1 + delta, f, (g    ) // 2, 2*u, 2*v, q  , r    return delta, (u, v, q, r)[/code] 
后续还需进行如下计算：

由推导过程可知，最后u，v，q，r组成的矩阵乘于初始向量[f g]以后得出的向量每个元素都是2N的倍数（可由初始值实际验算下），所以有以下的更新函数：

1. (old_r, r) = (a, b)
2. (old_x, x) = (1, 0)  # 表示 a 的系数
3. (old_y, y) = (0, 1)  # 表示 b 的系数
5. while r ≠ 0:
6.     q = old_r // r
7.      
8.     # 更新余数
9.     (old_r, r) = (r, old_r - q * r)
10.    
11.     # 更新系数
12.     (old_x, x) = (x, old_x - q * x)
13.     (old_y, y) = (y, old_y - q * y)
14.    
15. # 最终结果   
16. gcd(a,b) = old_r
17. x = old_x, y = old_y

复制代码

对于d和e来说上述结论并不成立，这里需要一个和div2函数类似的在模M下除于2N的函数，借助该函数可实现d，e的更新：

1. 初始值：(old_r, r) = (240, 46), (old_x, x) = (1, 0), (old_y, y) = (0, 1)
2. 第1步：240/46=5...10, q = 5, (old_r, r) = (46, 10), (old_x, x) = (0, 1), (old_y, y) = (1, -5)
3. 第2步：46/10=4...6, q = 4, (old_r, r) = (10, 6), (old_x, x) = (1, -4), (old_y, y) = (-5, 21)
4. 第3步：10/6=1...4, q = 1, (old_r, r) = (6, 4), (old_x, x) = (-4, 5), (old_y, y) = (21, -26)
5. 第4步：6/4=1...2, q = 1, (old_r, r) = (4, 2), (old_x, x) = (5, -9), (old_y, y) = (-26, 47)
6. 第5步：4/2=2...0, q = 2, (old_r, r) = (2, 0), (old_x, x) = (-9, 23), (old_y, y) = (47, -120)
7. 此时因r=0，得最终结果gcd(240, 46)=2，x=-9，y=47

复制代码

综合上述所有，即可给出执行N divsteps版本的modinv函数：

1. 1961年: 原始二进制GCD概念出现
2. 1967年: Josef Stein正式发表算法
3. 1970年代: Knuth在TAOCP中分析和优化
4. 1980-90年代: 各种优化变种出现，包括这个delta版本
5. 2000年代: 被广泛用于密码学库和大数运算

复制代码

这意味着在实践中，我们将始终执行多个N个div步骤。这不是问题，因为一旦g=0，进一步的divsteps就不再影响f、g、d或e（只有δ继续增加），这就能保证不管输入是何止，算法运行步骤一致，从而保证了算法时间一致性。
3.3 secp256k1_modinv32源码分析

结合上述内容分析下模逆的核心函数secp256k1_modinv32，其源码如下：

1. 1 def gcd(f, g):
2. 2     """Compute the GCD of an odd integer f and another integer g."""
3. 3     assert f & 1  # require f to be odd
4. 4     delta = 1     # additional state variable
5. 5     while g != 0:
6. 6         assert f & 1  # f will be odd in every iteration
7. 7         if delta > 0 and g & 1:
8. 8             delta, f, g = 1 - delta, g, (g - f) // 2
9. 9         elif g & 1:
10. 10             delta, f, g = 1 + delta, f, (g + f) // 2
11. 11         else:
12. 12             delta, f, g = 1 + delta, f, (g    ) // 2
13. 13     return abs(f)

复制代码

函数整体流程和上一节的介绍基本一致，不再进行详细分析，这里重点分析下函数中容易让人困惑的几个点。
1. 参数zeta

在代码一开始引入映射zeta = -(delta+1/2)，则在delta取值1/2时，zeta取值为-1，以此为起点进行迭代，在zeta引入后，有以下迭代分支：
zeta0且g&1分支；
(zeta>=0 且 g&1)或g为偶数：zeta=zeta-1，对应之前delta另外另种分支；
这样进行精心设计的好处是能用能用移位异或等操作完全替换比较操作（在密码学领域，消除分支判断，实现运行恒定时间）。

1. def gcd_no_delta(f, g):
2.     assert f & 1
3.     while g != 0:
4.         assert f & 1
5.         if g & 1:
6.             g = (g - f) // 2
7.         else:
8.             g = g // 2
9.     return abs(f)

复制代码

第2行处代码用于获取zeta符号，当zeta为非负数时c1=0，当zeta为负数时c1=0xFFFFFFFF，所以对于mask1来说，要么取值为0，要么为0xFFFFFFFF。第3行当g为奇数时，c2=1，对应mask2为0xFFFFFFFF，当g为偶数时，c2=0，对应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数量大约是2n+C，其中C是一个很小的常数。另外文中也说了在