钙沿雌闲1 magnitude及normalized由于当前许多项目都用到secp256k1库比特币作为体量最大的数字货币项目这里建议直接参考bitcoin-core提供的最新secp256k1源码。仍以field的10x26实现版本为例相关定义如下复制代码/** This field implementation represents the value as 10 uint32_t limbs in base* 2^26. */typedef struct {/* A field element f represents the sum(i0..9, f.n[i] (i*26)) mod p,* where p is the field modulus, 2^256 - 2^32 - 977.** The individual limbs f.n[i] can exceed 2^26; the fields magnitude roughly* corresponds to how much excess is allowed. The value* sum(i0..9, f.n[i] (i*26)) may exceed p, unless the field element is* normalized. */uint32_t n[10];/** Magnitude m requires:* n[i] 2 * m * (2^26 - 1) for i0..8* n[9] 2 * m * (2^22 - 1)** Normalized requires:* n[i] (2^26 - 1) for i0..8* sum(i0..9, n[i] (i*26)) p* (together these imply n[9] 2^22 - 1)*/SECP256K1_FE_VERIFY_FIELDS} secp256k1_fe;复制代码对于magnitude可称其为“量级”当m0时这时n[i] 0(i0...9)由此可知此时secp256k1_fe大数必为0当m1时n[i] 2*(2^26 - 1)对于i0...8n[9] 2*(2^22 - 1)有些说法是当m1时是将大数限制到[0,到2p)范围内这是不准确的此时secp256k1_fe大数范围是[0, 2*(2^256-1)]上限是稍大于2p的magnitude设计本质是“存储约束”而非“模p范围”通过magnitude约束将大数控制在可高效约简的范围内magnitude表示大概是2m个p的量级。对于normalized可称其为“规范化归一化”从注释中可知每个n[i]都小于或等于其对应的MASK且sum(n[i]) p即规范化后的大数是[0, p)之内的数。由此可知规范化的大数一定是magnitude为0或1的数但是magnitude为1的数不一定是规范化的大数。函数secp256k1_fe_normalize对大数实现规范化操作其本质是将大数看作为特殊的2^26进制大数即N∑ni*wi这里ni 2 * m * (2^26 - 1) for i0..8n9 2 * m * (2^22 - 1)wi2^(i*26)之所以说它特殊是因为由于m的存在每个位上数值是可以大于或等于基数2^26的正常情况下每个位上取值小于基数如10进制数每个位上数值都小于基数10另外最高位取值有特殊限制在secp256k1_fe_normalize函数中有两部分操作第一部分先将大数A规范为最多为257位的大数Bm确定时该数对应确定最大值后续详解第二部分将大数B规范为小于模数p的大数。在secp256k1源码中其实定义了各个参数大数中magnitude的最大值在定义VERIFY宏时会对参数合法性进行检查。#define SECP256K1_GE_X_MAGNITUDE_MAX 4#define SECP256K1_GE_Y_MAGNITUDE_MAX 3#define SECP256K1_GEJ_X_MAGNITUDE_MAX 4#define SECP256K1_GEJ_Y_MAGNITUDE_MAX 4#define SECP256K1_GEJ_Z_MAGNITUDE_MAX 1由以上定义可知在secp256k1算法中0 m 4当m0时显而易见大数即为0当m1且m4时其实大数对应最大值分别为pow(2,256)-1)*2*10x1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffepow(2,256)-1)*2*20x3fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffcpow(2,256)-1)*2*30x5fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffapow(2,256)-1)*2*40x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff8当m取值为4时对应最大值在调用规范化函数时第一部分规范完以后对应最大值为0x100003c00000f000000000001e000003c00000f000003c00000f0000f3c00393e为一个257位的大数函数secp256k1_fe_normalize_weak其实就是secp256k1_fe_normalize第一部分执行完以后的输出该值就是weak函数在m4下输出的最大值该值本身为m1的大数。2 常用函数magnitude值分析先分析大数取反原型函数secp256k1_fe_negate其原型如下复制代码SECP256K1_INLINE static void secp256k1_fe_negate(secp256k1_fe *r, const secp256k1_fe *a, int m) {r-n[0] 0x3FFFC2FUL * 2 * (m 1) - a-n[0];r-n[1] 0x3FFFFBFUL * 2 * (m 1) - a-n[1];r-n[2] 0x3FFFFFFUL * 2 * (m 1) - a-n[2];r-n[3] 0x3FFFFFFUL * 2 * (m 1) - a-n[3];r-n[4] 0x3FFFFFFUL * 2 * (m 1) - a-n[4];r-n[5] 0x3FFFFFFUL * 2 * (m 1) - a-n[5];r-n[6] 0x3FFFFFFUL * 2 * (m 1) - a-n[6];r-n[7] 0x3FFFFFFUL * 2 * (m 1) - a-n[7];r-n[8] 0x3FFFFFFUL * 2 * (m 1) - a-n[8];r-n[9] 0x03FFFFFUL * 2 * (m 1) - a-n[9];}复制代码这里m是输入参数a的magnitude当m0时输出参数r为2p显然其对应的m1当输入参数a的m1时其取值范围为[1, 0x1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe]则2*(m1)*p-a取值范围是[0x1fffffffffffffffffffffffffffffffffffffffffffffffffffffffbfffff0be, 0x3fffffffffffffffffffffffffffffffffffffffffffffffffffffffbfffff0bb]所以输出r的m2所以对于任意输入m输出r的mm1。对于secp256k1_fe_add函数其源码如下secp256k1_fe_add显而易见对于输出r来说m(r)m(r)m(a)r和a都是secp256k1_fe型数据。对于secp256k1_fe_mul_int(r, a)a为int型数据m(r)m(r)*a。对于secp256k1_fe_mul/secp256k1_fe_sqr函数在实现中包含内部规约会把输出参数的“量级”规约回magnitude 1使输出可以当作normalized-like使用。下面以函数为例给出各个步骤的详细magnitude值secp256k1_gej_add_var首先以输入参数ab中的xyz都为magnitude1normalized-like为前提在此基础上给出按代码顺序的表步 代码操作 输入 m 输出 m上限 说明1 secp256k1_fe_sqr(z22, b-z) b-z:1 1 fe_sqr → 规约为 12 secp256k1_fe_sqr(z12, a-z) a-z:1 13 secp256k1_fe_mul(u1, a-x, z22) 1,1 1 fe_mul → 规约为 14 secp256k1_fe_mul(u2, b-x, z12) 1,1 15 secp256k1_fe_mul(s1, a-y, z22) 1,1 16 secp256k1_fe_mul(s1, s1, b-z) 1,1 17 secp256k1_fe_mul(s2, b-y, z12) 1,1 18 secp256k1_fe_mul(s2, s2, a-z) 1,1 19 secp256k1_fe_negate(h, u1, 1) u1:1 2 fe_negate → m_out 1 1 210 secp256k1_fe_add(h, u2) h:2, u2:1 3 fe_add: 2 1 311 secp256k1_fe_negate(i, s2, 1) s2:1 212 secp256k1_fe_add(i, s1) i:2, s1:1 313 if (secp256k1_fe_normalizes_to_zero_var(h)) ... — — 特殊分支不进入一般流程14 secp256k1_fe_mul(t, h, b-z) h:3, b.z:1 1 fe_mul → 规约为 115 if (rzr ! NULL) *rzr t; t:1 116 secp256k1_fe_mul(r-z, a-z, t) a.z:1, t:1 117 secp256k1_fe_sqr(h2, h) h:3 1 fe_sqr → 规约为 118 secp256k1_fe_negate(h2, h2, 1) h2:1 219 secp256k1_fe_mul(h3, h2, h) h2:2, h:3 1 fe_mul → 规约为 120 secp256k1_fe_mul(t, u1, h2) u1:1, h2:2 121 secp256k1_fe_sqr(r-x, i) i:3 122 secp256k1_fe_add(r-x, h3) r-x:1, h3:1 223 secp256k1_fe_add(r-x, t) r-x:2, t:1 324 secp256k1_fe_add(r-x, t) r-x:3, t:1 425 secp256k1_fe_add(t, r-x) t:1, r-x:4 5 t 被覆盖为 t r-x26 secp256k1_fe_mul(r-y, t, i) t:5, i:3 1 fe_mul → 规约为 127 secp256k1_fe_mul(h3, h3, s1) h3:1, s1:1 128 secp256k1_fe_add(r-y, h3) r-y:1, h3:1 229 函数返回 r-x:4, r-y:2, r-z:1 — 这些是按最小合法上界得到的值未再 normalize 前由表可知对于secp256k1_gej_add_var函数虽然运行过程中零时变量可能会出现相对较高的magnitude值但最终返回值r-x的magnitude值是满足SECP256K1_GE_X_MAGNITUDE_MAX值的同理yz也满足SECP256K1_GEJ_Y_MAGNITUDE_MAX和SECP256K1_GEJ_Z_MAGNITUDE_MAX限制。3 大数求逆最新版大数求逆函数实现原型如下复制代码1 static void secp256k1_fe_impl_inv(secp256k1_fe *r, const secp256k1_fe *x) {2 secp256k1_fe tmp *x;3 secp256k1_modinv32_signed30 s;45 secp256k1_fe_normalize(tmp);6 secp256k1_fe_to_signed30(s, tmp);7 secp256k1_modinv32(s, secp256k1_const_modinfo_fe);8 secp256k1_fe_from_signed30(r, s);9 }复制代码3.1 高层概览1. 主要步骤函数做的是对域元素x求模逆即求r x-1 mod p这里p 2^256 - 2^32 - 977其实现流程如下1先把tmpx规范化保证tmp是库约定的limb范围下唯一表示2把规范化的tmp10x26-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_fe10x26表示并写入r该转换同时完成必要的约减/规范化。2. divsteps算法divsteps是一种高效计算GCDGreatest Common Divisor最大公约数的方法它通过连续多次除法步骤减少迭代次数特别适合硬件实现和大数计算。它是GCD算法的优化版本主要特点利用二进制表示的优势通过右移操作快速消除因子2一次迭代处理多个 减法 - 移位 步骤减少循环次数使用比较和差值运算替代昂贵的除法操作求u和v最大公约数算法步骤如下1消除因子2统计并移除两数中的所有因子2设a u / 2^sab v / 2^sb这里sa是u可整除2的最大次数sb是v可整除2的最大次数2divsteps迭代比较a和b计算差值d | a - b |消除d中的所有因子2得到d用min(a, b)和d替换(a, b)重复直到a b3得出最终结果gcd(u, v) b * 2^min(sa, sb)算法基于以下假设如果gcd(u, v)是最大公约数那么它可以分成两部分的乘积一部分是2的整数次幂另一部分是非2的倍数。步骤1其实是获取u、v最大公因数中能2整数次幂部分即2min(sa, sb)步骤2获取非2的倍数部分。步骤1比较好理解这里分析下步骤2假设得出最后一步结论ab时具体值为ab则a必然满足a - b b*2x则有a b*(12x)即a是b的倍数则在上一步必有db*2y a - a由此a a b*2y即a也为b的倍数由此类推可知最一开始的a和b闭然都是b的倍数。由以上推理可知最终gcd(u, v) b * 2^min(sa, sb)。以下是算法示例代码divsteps计算示例gcd(270, 324)1初始值u270v3242移除因子2sa1sb2 → u135, v813第一次迭代d135-8154移除因子2: sd1 → d27新值: u81, v274第二次迭代d81-2754移除因子2: sd1 → d27新值: u27, v275循环结束结果gcd(270, 324) 27*2^1 54。3.2 逐步详解步骤1将大数进行规范化这一步由secp256k1_fe_normalize实现该函数已经在之前文章中详细介绍这里不再进行分析。步骤2调用secp256k1_fe_to_signed30函数实现该函数源码如下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基础引理欧几里得算法的余数性质欧几里得算法通过反复求余计算GCDgcd(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 1y 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 yy 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看作是newrb看作是ra看作是oldr则有oldr q?r newr这里q是整数商a//b。接下来因为求系数x和y其实也是一个迭代推导的过程假设当前有image可以先把最一开始oldr和r的由系数表示出来oldr 1*a 0*br 0*a 1*b即存在初始余数对(oldr, r) (a, b)初始a的系数对(oldx, x) (1, 0)初始b的系数对(oldy, y) (0, 1)满足方程。再接下来继续进行迭代q oldr//rnewr oldr - q*r将(1)中oldr和(2)中的r代入该等式并进行整理可得image由此可得newx oldx - q*xnewy oldy - q*y即newr仍是ab的线性组合且newxnewy是相应系数。以下是完整算法伪代码复制代码(old_r, r) (a, b)(old_x, x) (1, 0) # 表示 a 的系数(old_y, y) (0, 1) # 表示 b 的系数while r ≠ 0:q old_r // r# 更新余数(old_r, r) (r, old_r - q * r)# 更新系数(old_x, x) (x, old_x - q * x)(old_y, y) (y, old_y - q * y)# 最终结果gcd(a,b) old_rx old_x, y old_y复制代码以a240b46为例给出算法迭代过程复制代码初始值(old_r, r) (240, 46), (old_x, x) (1, 0), (old_y, y) (0, 1)第1步240/465...10, q 5, (old_r, r) (46, 10), (old_x, x) (0, 1), (old_y, y) (1, -5)第2步46/104...6, q 4, (old_r, r) (10, 6), (old_x, x) (1, -4), (old_y, y) (-5, 21)第3步10/61...4, q 1, (old_r, r) (6, 4), (old_x, x) (-4, 5), (old_y, y) (21, -26)第4步6/41...2, q 1, (old_r, r) (4, 2), (old_x, x) (5, -9), (old_y, y) (-26, 47)第5步4/22...0, q 2, (old_r, r) (2, 0), (old_x, x) (-9, 23), (old_y, y) (47, -120)此时因r0得最终结果gcd(240, 46)2x-9y47复制代码2. 二进制GCD算法变种二进制GCD算法有以下演进时间线1961年: 原始二进制GCD概念出现1967年: Josef Stein正式发表算法1970年代: Knuth在TAOCP中分析和优化1980-90年代: 各种优化变种出现包括这个delta版本2000年代: 被广泛用于密码学库和大数运算接下来重点分析带有delta状态变量的版本算法如下复制代码1 def gcd(f, g):2 Compute the GCD of an odd integer f and another integer g.3 assert f 1 # require f to be odd4 delta 1 # additional state variable5 while g ! 0:6 assert f 1 # f will be odd in every iteration7 if delta 0 and g 1:8 delta, f, g 1 - delta, g, (g - f) // 29 elif g 1:10 delta, f, g 1 delta, f, (g f) // 211 else:12 delta, f, g 1 delta, f, (g ) // 213 return abs(f)复制代码该算法要求第一个参数f必须是奇数第3行并引入了状态变量delta用于指导算法选择不同的约减策略避免陷入低效循环。该算法证明分两部分1先证明算法保持 gcd 不变正确性不破坏2然后给出一个势函数potential并用它证明在每次若干步内势函数会严格下降从而不可能无限迭代——也就保证终止收敛。同时指出 delta 在防止振荡 / 保证下降中的关键作用。第1部分结论很容易得出在结合fg奇偶情况下gcd(g, (g - f)//2)gcd(f, (g f)//2)gcd(f, g//2)这三种情况显然和gcd(f, g)是一样的。第2部分的证明比较复杂涉及到势函数相关理论知识盲区请自行查阅相关资料这里不再详细说明。总之相比如下未引入delta的算法一在某些情况下会不收敛上述算法是能保证收敛的。gcd_no_delta另外相比如下未引入delta的算法二fg直接判断属于“输入依赖型分支”不同输入会导致不同的分支走向容易在硬件层面产生 “分支预测错误”尤其在加密算法中输入的随机性会放大这种问题。而delta引入后分支判断虽然存在但更多依赖内部状态变量 delta 的符号delta 0 或 delta ≤ 0而非输入 f 和 g 的直接比较delta 的更新规则是固定的1 - delta 或 1 delta其变化模式相对可预测降低了对输入随机性的敏感度更适合常数时间实现密码学算法的关键要求。尤其在大数运算情况下a b 这类比较操作需要完整的减法和符号判断硬件层面是减法器 符号位检测而delta的状态更新1 - delta 或 1 delta是简单的算术操作本质是 delta 在0和1之间交替或递增递减计算成本更低。bad_gcd3. 由GCD求模逆对于质数M大于2来说求x0复制代码def div2(M, x):Helper routine to compute x/2 mod M (where M is odd).assert M 1if x 1: # If x is odd, make it even by adding M.x M# x must be even now, so a clean division by 2 is possible.return x // 2def modinv(M, x):Compute the inverse of x mod M (given that it exists, and M is odd).assert M 1delta, f, g, d, e 1, M, x, 0, 1while g ! 0:# Note that while division by two for f and g is only ever done on even inputs, this is# not true for d and e, so we need the div2 helper function.if delta 0 and g 1:delta, f, g, d, e 1 - delta, g, (g - f) // 2, e, div2(M, e - d)elif g 1:delta, f, g, d, e 1 delta, f, (g f) // 2, d, div2(M, e d)else:delta, f, g, d, e 1 delta, f, (g ) // 2, d, div2(M, e )# Verify that the invariants df/x mod M, eg/x mod M are maintained.assert f % M (d * x) % Massert g % M (e * x) % Massert f 1 or f -1 # |f| is the GCD, it must be 1# Because of invariant d f/x (mod M), 1/x d/f (mod M). As |f|1, d/f d*f.return (d * f) % M复制代码该算法基于二进制扩展GCD算法通过不断将问题规模减半来高效计算模逆元首先明确算法以下关键变量M: 模数必须是奇数x: 要求逆元的数f, g: 跟踪GCD计算的两个值d, e: 跟踪系数维护关键不变量delta: 控制算法分支的状态变量算法关键在于每次迭代都维护了以下等式λf ≡ d * x (mod M)g ≡ e * x (mod M)这意味着每次迭代d和e始终是f和g在模M下关于x的系数这样只要保证初始时等式成立以及每次迭代更新fg时按计算规则得到的de和fg仍满足以上等式这就能保证到最后求得最大公约数1时以上等式仍成立及找出了最终的逆元d。以循环迭代中第1种情况为例进行分析根据规则更新后的f gg (g - f)//2也就是将g替换为新的f将(g-f)//2替换为新的g那么对应的系数更新规则必须保证image由于更新后f g由最一开始的等式g ≡ e * x (mod M)可知要保证f ≡ d * x (mod M)成立取d e即可这正是算法中更新d的逻辑。算法中g的更新逻辑是g (g - f)//2将初始等式λ代入更新逻辑得image这就要求更新后的e满足image等价于image而while迭代中第1种情况中的e的更新函数div2正是实现了该逻辑在函数中由于输入e-d可能为奇数此时需要加上奇数M在模M下不影响最终结果使得和为偶数再除于2。同理可以分析while迭代中的另外两个分支最终当算法终止时有gcd(M, x) |f| 1则当f 1时由f≡d*x≡1 mod M可知d即为x逆元当f -1时由d*x≡-1 mod M可知-d即为x逆元正是算法返回值。4. 模逆算法进一步优化上述每次divstep迭代可以表示为矩阵乘法对向量[f g]和[d e]应用如下的转换矩阵1/2*t