1 关键结构体
1.1 secp256k1_fe
secp256k1库域元素field element，其具体定义如下
/** This field implementation represents the value as 10 uint32_t limbs in base * 2^26. */ typedef struct { /* A field element f represents the sum(i=0..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 field's magnitude roughly * corresponds to how much excess is allowed. The value * sum(i=0..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 i=0..8 * n[9] <= 2 * m * (2^22 - 1) * * Normalized requires: * n[i] <= (2^26 - 1) for i=0..8 * sum(i=0..9, n[i] << (i*26)) < p * (together these imply n[9] <= 2^22 - 1) */ SECP256K1_FE_VERIFY_FIELDS } secp256k1_fe;
由定义可知，secp256k1_fe由10个uint32_t变量进行数据存储，9 个 26 位字（n[0]到n[8]）和1 个 22 位字（n[9]），总位数为：9*26+22=256位，每个字大小为：n[i] ≤ 2^26 - 1（i=0..8），n[9] ≤ 2^22 - 1。在定义VERIFY时，m（magnitude）是一个缩放因子（scaling factor），用于控制每个字的最大允许值。具体来说：
对于前 9 个字（i=0..8）：
n[i] ≤ 2 × m × (2^26 - 1)
这意味着每个字的最大值被放大了 2m 倍。
对于最后一个字（i=9）：
n[9] ≤ 2 × m × (2^22 - 1)
同样被放大了 2m 倍。
在椭圆曲线运算（如模乘、点加）中，中间结果可能会超出单个字的范围。通过允许每个字暂时存储更大的值（通过m缩放），可以避免频繁的归一化（normalization）操作，从而提高效率。这里的归一化，是将每个字的值压缩到标准范围（如 n[i] ≤ 2^26 - 1）的过程。当m=1时：
前9个字的最大值为2 × (2^26 - 1) = 134,217,726（二进制下 27 位）。
最后一个字的最大值为2 × (2^22 - 1) = 8,388,606（二进制下 24 位）。
此外根据椭圆曲线方程，x坐标三次方运算可能导致中间结果快速增长，需要更大的临时存储空间，所以x能需要更大缩放希数，如下定义
/** Maximum allowed magnitudes for group element coordinates * in affine (x, y) and jacobian (x, y, z) representation. */ #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
1.2 secp256k1_fe_storage
typedef struct {
uint32_t n[8];
} secp256k1_fe_storage;
正如其名称显示，该数据结构用8个uint32_t表示大数，最大支持256位大数。
1.3 secp256k1_ge & secp256k1_gej
群元素group element，椭圆曲线上点的仿射坐标（affine coordinates），x和y坐标都是secp256k1_fe类型大数，infinity变量指明点是否为无穷远点。
/** A group element in affine coordinates on the secp256k1 curve, * or occasionally on an isomorphic curve of the form y^2 = x^3 + 7*t^6. * Note: For exhaustive test mode, secp256k1 is replaced by a small subgroup of a different curve. */ typedef struct { secp256k1_fe x; secp256k1_fe y; int infinity; /* whether this represents the point at infinity */ } secp256k1_ge;
椭圆曲线上点的Jacobian射影坐标，射影坐标(x, y, z)对应的仿射坐标为(x/z^2, y/z^3)，同样infinity变量指明点是否为无穷远点。
/** A group element of the secp256k1 curve, in jacobian coordinates. * Note: For exhastive test mode, secp256k1 is replaced by a small subgroup of a different curve. */ typedef struct { secp256k1_fe x; /* actual X: x/z^2 */ secp256k1_fe y; /* actual Y: y/z^3 */ secp256k1_fe z; int infinity; /* whether this represents the point at infinity */ } secp256k1_gej;
1.4 secp256k1_ge_storage
以存储形式表示的group element，该形式坐标无法表示无穷远点。
typedef struct { secp256k1_fe_storage x; secp256k1_fe_storage y; } secp256k1_ge_storage;
2 关键理论
2.1 Schoolbook 算法
Schoolbook 算法（竖式乘法算法）是最基础的多精度整数乘法方法，直接模拟人类手算乘法的过程。它易于理解和实现，适用于中小规模的数值计算。
一 核心原理
Schoolbook算法的核心思想是将多位数乘法拆解为多个单字乘法和累加操作，类似于小学竖式计算：
逐位相乘：用乘数的每一位分别去乘被乘数的每一位，得到部分积。
错位累加：将每个部分积按位对齐后累加，得到最终结果。
二 算法步骤（以十进制为例）
假设计算两个多位数A*B，其中：
A = an-1an-2...a1a0（共n位）
B = bm-1bm-2...b1b0（共m位）
具体步骤：
1 初始化结果数组：创建一个长度为n+m的数组C，用于存储中间结果和最终乘积。
2 逐位相乘并累加：
1 for i in range(m): # 遍历乘数的每一位 2 carry = 0 3 for j in range(n): # 遍历被乘数的每一位 4 # 计算当前位的乘积和进位 5 product = a[j] × b[i] + C[i+j] + carry 6 C[i+j] = product % 10 # 当前位的结果 7 carry = product // 10 # 进位 8 C[i+n] = carry # 处理最高位的进位
3 结果处理：
数组C即为最终乘积，可能包含前导零，需按需处理。
三 计算机实现细节（多字乘法）
在计算机中，多精度整数通常以字数组（Word Array）形式存储，每个字（Word）是一个固定位数的整数（如 32 位或 64 位）。此时 Schoolbook 算法的实现需注意：
1 单字乘法的处理：
两个w位的字相乘可能产生2w位的结果，需用两个字表示。例如，32位 × 32位 → 64位结果，需拆分为高位和低位两部分。
2 进位管理：
每次乘法产生的进位需传递到高位，累加中间结果时需处理多字进位。
1 A = [0xFF, 0xFF] 2 B = [0xFF, 0xFF, 0xFF] 3 C = [0, 0, 0, 0, 0] # 结果数组 4 w = 8 5 m = 3 6 n = 2 7 for i in range(m): 8 carry = 0 9 10 for j in range(n): 11 product = A[j] * B[i] 12 product_high = product >> w 13 product_low = product & ((1 << w) - 1) 14 15 # 累加低位部分到当前位置 16 sum_low = product_low + C[i+j] 17 carry = sum_low >> w # 进位传递到下一位 18 C[i+j] = sum_low & ((1 << w) - 1) 19 20 # 累加高位部分及进位到下一个位置 21 C[i+j+1] += product_high+carry 22 23 print("result\n0x", end="") 24 for i in range(m+n): 25 print("%02x"%C[m+n-1-i], end="") 26 print("")
程序中第11~13行将乘积分成高低字节两部分，第16行将低字节部分和当前位置现有值相加，17行获取相加结果的进位，18行将相加结果的低字节部分保留在当前位置，21行将之前乘积的高字节部分和产生的进位一起累加到下一个位置，最后输出最终结果。
2.2 二次剩余及相关理论
一 欧拉函数
对于正整数n，欧拉函数φ(n)表示小于n且与n互质的数的个数。如果n是质数，则φ(n)=n-1。如果n可以分解成两个互质的整数之积，则其欧拉函数等于其互质因子欧拉函数之积，即φ(n)=φ(p1p2)=φ(p1)φ(p2)。
二 欧拉定理
费马小定理：若p是质数，且整数a满足gcd(a, p) = 1，则ap-1 ≡ 1 mod p。其经典证明如下：
1 构造互质的剩余希
考虑小于p的正整数集合S = { 1, 2, 3, ..., p-1 }，由于p是质数，且gcd(a, p) = 1，则集合S中每个元素与p互质，且a与p也互质。
构造新集合T = { a*1 mod p, a*2 mod p, a*(p-1) mod p}，即对S中每个元素乘以a后取模p。
2 证明T是S的一个排列
需证明T中的元素与S中的元素完全相同（仅顺序可能不同）：
互质性：由于gcd(a, p) = 1且gcd(k, p) = 1（k ∈ S），根据互质性质，gcd(a*k, p) = 1，因此a*k mod p仍属于S（即T ⊆ S）。
唯一性：若a*k1 ≡ a*k2 mod p，则a*(k1 - k2) ≡ 0 mod p，由于gcd(a, p) = 1，可推出k1 ≡ k2 mod p，又因为k1，k2 < p，故k1 = k2，即T中元素互不相同。
因此，T是S的一个排列（元素完全一致，顺序不同）。
3 两边乘积的同余关系
由于T是S的排列，两集合中的所有元素的乘积模p相等：
(a*1)*(a*2)*...*(a*(p-1)) ≡ 1*2*...*(p-1) mod p
左边可化简为：
ap-1*(1*2*...*(p-1) ≡ (1*2*...*(p-1) mod p，令W = 1*2*...*(p-1)，则上式变为：
ap-1*W ≡ W mod p
4 消去W得到结论
由于p是质数，W = (p-1)!与p互质，因此可以同余式两边同时除于W（或乘以W的逆元），即得ap-1 ≡ W mod p。
欧拉定理是费马小定理的推广，使用于更一般的情形。其定义：若两个正整数a和n互质（即gcd(a, n) = 1），则aφ(n) ≡ 1 mod n，φ(n)是欧拉函数。其证明过程类似。
三 二次剩余
设p是一个正整数（通常为质数），且整数a与p互质（即gcd(a, p) = 1），如果存在整数x使得：x2 ≡ a mod p，则称a是模p的二次剩余；否则，称a是模p的二次非剩余。简单来说，二次剩余是“模意义下能开平方的数”。例如：
模7时，12 ≡ 1，22 ≡ 4，32 ≡ 2，42 ≡ 2，52 ≡ 4，62 ≡ 1，因此1、2、4是模7的二次剩余，3、5、6是模7的二次非剩余。
二次剩余有如下关键性质
1 数量对称性
若p是奇质数，则模p的二次剩余和二次非剩余各有(p-1)/2个。
2 欧拉判别法
对于奇质数p和与p互质的整数a，a是模p的二次剩余的充要条件是：a(p-1)/2 ≡ 1 mod p，若a是二次非剩余，则：a(p-1)/2 ≡ -1 mod p。
3 勒让德符号
为简化二次剩余的表示，定义勒让德符号
其核心性质是：
后续的secp256k1_fe_sqrt函数实现用到了这里的二次剩余理论，在这里先提前说一下，因为a是模p的二次剩余的充要条件是a(p-1)/2 ≡ 1 mod p，则有a*a(p-1)/2 ≡ a mod p，即a(p+1)/2 ≡ a mod p，如果存在k使得p=4*k+3，那么(p+1)/4正好是整数，这时有(a(p+1)/4)2 ≡ a mod p，也就是说如果a是模p的二次剩余，则a(p+1)/4必是a的开方之一（模p操作下），这正是secp256k1_fe_sqrt函数的理论基础。
3 关键函数解析
3.1 secp256k1_fe_normalize
该函数实现大数的归一化操作，其实现如下
1 static void secp256k1_fe_normalize(secp256k1_fe *r) { 2 uint32_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4], 3 t5 = r->n[5], t6 = r->n[6], t7 = r->n[7], t8 = r->n[8], t9 = r->n[9]; 4 5 /* Reduce t9 at the start so there will be at most a single carry from the first pass */ 6 uint32_t m; 7 uint32_t x = t9 >> 22; t9 &= 0x03FFFFFUL; 8 9 /* The first pass ensures the magnitude is 1, ... */ 10 t0 += x * 0x3D1UL; t1 += (x << 6); 11 t1 += (t0 >> 26); t0 &= 0x3FFFFFFUL; 12 t2 += (t1 >> 26); t1 &= 0x3FFFFFFUL; 13 t3 += (t2 >> 26); t2 &= 0x3FFFFFFUL; m = t2; 14 t4 += (t3 >> 26); t3 &= 0x3FFFFFFUL; m &= t3; 15 t5 += (t4 >> 26); t4 &= 0x3FFFFFFUL; m &= t4; 16 t6 += (t5 >> 26); t5 &= 0x3FFFFFFUL; m &= t5; 17 t7 += (t6 >> 26); t6 &= 0x3FFFFFFUL; m &= t6; 18 t8 += (t7 >> 26); t7 &= 0x3FFFFFFUL; m &= t7; 19 t9 += (t8 >> 26); t8 &= 0x3FFFFFFUL; m &= t8; 20 21 /* ... except for a possible carry at bit 22 of t9 (i.e. bit 256 of the field element) */ 22 23 /* At most a single final reduction is needed; check if the value is >= the field characteristic */ 24 x = (t9 >> 22) | ((t9 == 0x03FFFFFUL) & (m == 0x3FFFFFFUL) 25 & ((t1 + 0x40UL + ((t0 + 0x3D1UL) >> 26)) > 0x3FFFFFFUL)); 26 27 /* Apply the final reduction (for constant-time behaviour, we do it always) */ 28 t0 += x * 0x3D1UL; t1 += (x << 6); 29 t1 += (t0 >> 26); t0 &= 0x3FFFFFFUL; 30 t2 += (t1 >> 26); t1 &= 0x3FFFFFFUL; 31 t3 += (t2 >> 26); t2 &= 0x3FFFFFFUL; 32 t4 += (t3 >> 26); t3 &= 0x3FFFFFFUL; 33 t5 += (t4 >> 26); t4 &= 0x3FFFFFFUL; 34 t6 += (t5 >> 26); t5 &= 0x3FFFFFFUL; 35 t7 += (t6 >> 26); t6 &= 0x3FFFFFFUL; 36 t8 += (t7 >> 26); t7 &= 0x3FFFFFFUL; 37 t9 += (t8 >> 26); t8 &= 0x3FFFFFFUL; 38 39 /* If t9 didn't carry to bit 22 already, then it should have after any final reduction */ 40 41 /* Mask off the possible multiple of 2^256 from the final reduction */ 42 t9 &= 0x03FFFFFUL; 43 44 r->n[0] = t0; r->n[1] = t1; r->n[2] = t2; r->n[3] = t3; r->n[4] = t4; 45 r->n[5] = t5; r->n[6] = t6; r->n[7] = t7; r->n[8] = t8; r->n[9] = t9; 46 }
域的模数p = 2256 - 232 - 977，域中大数r被分为10个limbs，其标准表示为每个limb是26-bit数（除了最后一个是22-bit数），即：
如之前所描述，由于magnitude的存在，在运算过程中，会使某些limbs超过26-bit，导致进位错误。归一化目标：1 确保最终值 < p; 2 所有limbs小于226（除了t9，小于222）。
在函数中共进行了两次模约减操作，首轮模约减主要目的是去除高位t9的22-bit之后多余部分：
uint32_t x = t9 >> 22; // 取出超过 256 位的部分 t9 &= 0x03FFFFFUL; // 保留低 22 位
这里的x值为r/2256向下取整，由SECP256K1_GE_X_MAGNITUDE_MAX等宏定义可知x最大为4-bit数，因为2256 ≡ 232 + 977 mod p，所以超出2256部分即为 x*2256 ≡ x*(232 + 977) mod p ≡ x*(1 << 32) + x*977 mod p ≡ x*(1 << 6)*226 + x*977 mod p，由以上推理可知，x*2256在模p操作下，可分解为两部分，其中x*977（即x*0x3D1）可加到t0部分，而x*232可以将乘积的高于26-bit部分x<<6加到t1，由于x最大为4-bit数，所以x<<6最大为10-bit数，也就是x*2256高位在模p意义下，只能循环传导到t0和t1，对应代码第10行。之后11~19行，依次按标准模式计算t0~t8，本级保留低26-bits，进位加到下一级，最终所有进位传导到t9。运行进行到这里会有两种情况：1, t9 >> 22为1，即 t9有进位，说明大数归一化值还大于2256，需要进一步归约；2, t9 == 222 - 1（即最大值），t2~t8也都为0x3FFFFFF最大值，因为2256 - p = 0x1000003d1 = 0x40*226 + 0x3d1，所以如果t1 + 0x40 + ((t0 + 0x3D1UL) >> 26大于0x3FFFFFF，则必然产生向t2的进位，从而直接导致最终t9超过最大值，也即和超过2256，假设大数归一化值是A，可用不等式A+(2256 - p) ≥ 2256，即A ≥ p，当A等于p时，其加法可以通过下图清楚的看到各进位情况。简单的说就是这两种情况就是第一次归一化后的值大于2256或者大于或等于p，对应代码24~25行，这两种情况下就仍需进一步归一化操作（可以进一步验证在magnitude不大于4的情况下，第一次归一化操作后进位x最大为1）。
类似于第一次归一化，代码28~37行，将可能为1的进位x依次传导到t0~t9中，然后第42行将已经“补偿”到低位的最高位用且操作消除掉，最后44~45将结果赋值给远大数变量。了解了该函数的原理，函数secp256k1_fe_normalize_weak就更好理解了，它只是进行了一次归一化操作，也就是说最终结果t9可能还在23bit处包含1个最高进位。
3.2 secp256k1_fe_normalizes_to_zero
static int secp256k1_fe_normalizes_to_zero(secp256k1_fe *r) { uint32_t t0 = r->n[0], t1 = r->n[1], t2 = r->n[2], t3 = r->n[3], t4 = r->n[4], t5 = r->n[5], t6 = r->n[6], t7 = r->n[7], t8 = r->n[8], t9 = r->n[9]; /* z0 tracks a possible raw value of 0, z1 tracks a possible raw value of P */ uint32_t z0, z1; /* Reduce t9 at the start so there will be at most a single carry from the first pass */ uint32_t x = t9 >> 22; t9 &= 0x03FFFFFUL; /* The first pass ensures the magnitude is 1, ... */ t0 += x * 0x3D1UL; t1 += (x << 6); t1 += (t0 >> 26); t0 &= 0x3FFFFFFUL; z0 = t0; z1 = t0 ^ 0x3D0UL; t2 += (t1 >> 26); t1 &= 0x3FFFFFFUL; z0 |= t1; z1 &= t1 ^ 0x40UL; t3 += (t2 >> 26); t2 &= 0x3FFFFFFUL; z0 |= t2; z1 &= t2; t4 += (t3 >> 26); t3 &= 0x3FFFFFFUL; z0 |= t3; z1 &= t3; t5 += (t4 >> 26); t4 &= 0x3FFFFFFUL; z0 |= t4; z1 &= t4; t6 += (t5 >> 26); t5 &= 0x3FFFFFFUL; z0 |= t5; z1 &= t5; t7 += (t6 >> 26); t6 &= 0x3FFFFFFUL; z0 |= t6; z1 &= t6; t8 += (t7 >> 26); t7 &= 0x3FFFFFFUL; z0 |= t7; z1 &= t7; t9 += (t8 >> 26); t8 &= 0x3FFFFFFUL; z0 |= t8; z1 &= t8; z0 |= t9; z1 &= t9 ^ 0x3C00000UL; /* ... except for a possible carry at bit 22 of t9 (i.e. bit 256 of the field element) */ return (z0 == 0) | (z1 == 0x3FFFFFFUL); }
该函数在不完全归一化（normalize）整个字段元素的情况下，判断其值是否为零或等于模数，代码中z0用于检测是否所有位都为 0，z1用于检测是否所有位构成了模数p。类似secp256k1_fe_normalizes_to_zero_var函数是secp256k1_fe_normalizes_to_zero的变体（var）版，其功能与前者一致：判断一个字段元素是否为0或模数p（即归一化后是否为 0），但使用了更快路径和懒惰归一化以提高性能，该函数的关键代码如下：
t0 = r->n[0]; t9 = r->n[9]; x = t9 >> 22; t0 += x * 0x3D1UL; z0 = t0 & 0x3FFFFFFUL; z1 = z0 ^ 0x3D0UL; if ((z0 != 0UL) & (z1 != 0x3FFFFFFUL)) { return 0; }
若t9有进位（超过 22 位），说明有“≥2256”成分，我们提取出这部分（乘以 0x3D1 并加回 t0），这是归一化关键。然后立即对 t0 做一次判断：z0 == 0 表示可能是 0，z1 == 0x3FFFFFFUL 表示可能是 P（因为 t0 应该是 0x3D0），如果两个都不满足，直接返回 0，这是个非常快的判断分支。其他部分代码就和secp256k1_fe_normalizes_to_zero一致了。
3.3 secp256k1_fe_mul_inner
该函数是大数乘法的内部实现，libsecp256k1的作者Pieter Wuille的设计理念可以总结为："所有公共API都要求输入是完全归一化的，但内部函数可以处理非归一化形式。中间值的大小必须受控，但可以比26位大（比如在该函数内部被限制在30bit），为的是避免频繁的carry"。以下是完整函数
1 SECP256K1_INLINE static void secp256k1_fe_mul_inner(uint32_t *r, const uint32_t *a, const uint32_t * SECP256K1_RESTRICT b) { 2 uint64_t c, d; 3 uint64_t u0, u1, u2, u3, u4, u5, u6, u7, u8; 4 uint32_t t9, t1, t0, t2, t3, t4, t5, t6, t7; 5 const uint32_t M = 0x3FFFFFFUL, R0 = 0x3D10UL, R1 = 0x400UL; 6 7 VERIFY_BITS(a[0], 30); 8 VERIFY_BITS(a[1], 30); 9 VERIFY_BITS(a[2], 30); 10 VERIFY_BITS(a[3], 30); 11 VERIFY_BITS(a[4], 30); 12 VERIFY_BITS(a[5], 30); 13 VERIFY_BITS(a[6], 30); 14 VERIFY_BITS(a[7], 30); 15 VERIFY_BITS(a[8], 30); 16 VERIFY_BITS(a[9], 26); 17 VERIFY_BITS(b[0], 30); 18 VERIFY_BITS(b[1], 30); 19 VERIFY_BITS(b[2], 30); 20 VERIFY_BITS(b[3], 30); 21 VERIFY_BITS(b[4], 30); 22 VERIFY_BITS(b[5], 30); 23 VERIFY_BITS(b[6], 30); 24 VERIFY_BITS(b[7], 30); 25 VERIFY_BITS(b[8], 30); 26 VERIFY_BITS(b[9], 26); 27 28 /** [... a b c] is a shorthand for ... + a<<52 + b<<26 + c<<0 mod n. 29 * for 0 <= x <= 9, px is a shorthand for sum(a[i]*b[x-i], i=0..x). 30 * for 9 <= x <= 18, px is a shorthand for sum(a[i]*b[x-i], i=(x-9)..9) 31 * Note that [x 0 0 0 0 0 0 0 0 0 0] = [x*R1 x*R0]. 32 */ 33 34 d = (uint64_t)a[0] * b[9] 35 + (uint64_t)a[1] * b[8] 36 + (uint64_t)a[2] * b[7] 37 + (uint64_t)a[3] * b[6] 38 + (uint64_t)a[4] * b[5] 39 + (uint64_t)a[5] * b[4] 40 + (uint64_t)a[6] * b[3] 41 + (uint64_t)a[7] * b[2] 42 + (uint64_t)a[8] * b[1] 43 + (uint64_t)a[9] * b[0]; 44 /* VERIFY_BITS(d, 64); */ 45 /* [d 0 0 0 0 0 0 0 0 0] = [p9 0 0 0 0 0 0 0 0 0] */ 46 t9 = d & M; d >>= 26; 47 VERIFY_BITS(t9, 26); 48 VERIFY_BITS(d, 38); 49 /* [d t9 0 0 0 0 0 0 0 0 0] = [p9 0 0 0 0 0 0 0 0 0] */ 50 51 c = (uint64_t)a[0] * b[0]; 52 VERIFY_BITS(c, 60); 53 /* [d t9 0 0 0 0 0 0 0 0 c] = [p9 0 0 0 0 0 0 0 0 p0] */ 54 d += (uint64_t)a[1] * b[9] 55 + (uint64_t)a[2] * b[8] 56 + (uint64_t)a[3] * b[7] 57 + (uint64_t)a[4] * b[6] 58 + (uint64_t)a[5] * b[5] 59 + (uint64_t)a[6] * b[4] 60 + (uint64_t)a[7] * b[3] 61 + (uint64_t)a[8] * b[2] 62 + (uint64_t)a[9] * b[1]; 63 VERIFY_BITS(d, 63); 64 /* [d t9 0 0 0 0 0 0 0 0 c] = [p10 p9 0 0 0 0 0 0 0 0 p0] */ 65 u0 = d & M; d >>= 26; c += u0 * R0; 66 VERIFY_BITS(u0, 26); 67 VERIFY_BITS(d, 37); 68 VERIFY_BITS(c, 61); 69 /* [d u0 t9 0 0 0 0 0 0 0 0 c-u0*R0] = [p10 p9 0 0 0 0 0 0 0 0 p0] */ 70 t0 = c & M; c >>= 26; c += u0 * R1; 71 VERIFY_BITS(t0, 26); 72 VERIFY_BITS(c, 37); 73 /* [d u0 t9 0 0 0 0 0 0 0 c-u0*R1 t0-u0*R0] = [p10 p9 0 0 0 0 0 0 0 0 p0] */ 74 /* [d 0 t9 0 0 0 0 0 0 0 c t0] = [p10 p9 0 0 0 0 0 0 0 0 p0] */ 75 76 c += (uint64_t)a[0] * b[1] 77 + (uint64_t)a[1] * b[0]; 78 VERIFY_BITS(c, 62); 79 /* [d 0 t9 0 0 0 0 0 0 0 c t0] = [p10 p9 0 0 0 0 0 0 0 p1 p0] */ 80 d += (uint64_t)a[2] * b[9] 81 + (uint64_t)a[3] * b[8] 82 + (uint64_t)a[4] * b[7] 83 + (uint64_t)a[5] * b[6] 84 + (uint64_t)a[6] * b[5] 85 + (uint64_t)a[7] * b[4] 86 + (uint64_t)a[8] * b[3] 87 + (uint64_t)a[9] * b[2]; 88 VERIFY_BITS(d, 63); 89 /* [d 0 t9 0 0 0 0 0 0 0 c t0] = [p11 p10 p9 0 0 0 0 0 0 0 p1 p0] */ 90 u1 = d & M; d >>= 26; c += u1 * R0; 91 VERIFY_BITS(u1, 26); 92 VERIFY_BITS(d, 37); 93 VERIFY_BITS(c, 63); 94 /* [d u1 0 t9 0 0 0 0 0 0 0 c-u1*R0 t0] = [p11 p10 p9 0 0 0 0 0 0 0 p1 p0] */ 95 t1 = c & M; c >>= 26; c += u1 * R1; 96 VERIFY_BITS(t1, 26); 97 VERIFY_BITS(c, 38); 98 /* [d u1 0 t9 0 0 0 0 0 0 c-u1*R1 t1-u1*R0 t0] = [p11 p10 p9 0 0 0 0 0 0 0 p1 p0] */ 99 /* [d 0 0 t9 0 0 0 0 0 0 c t1 t0] = [p11 p10 p9 0 0 0 0 0 0 0 p1 p0] */ 100 101 c += (uint64_t)a[0] * b[2] 102 + (uint64_t)a[1] * b[1] 103 + (uint64_t)a[2] * b[0]; 104 VERIFY_BITS(c, 62); 105 /* [d 0 0 t9 0 0 0 0 0 0 c t1 t0] = [p11 p10 p9 0 0 0 0 0 0 p2 p1 p0] */ 106 d += (uint64_t)a[3] * b[9] 107 + (uint64_t)a[4] * b[8] 108 + (uint64_t)a[5] * b[7] 109 + (uint64_t)a[6] * b[6] 110 + (uint64_t)a[7] * b[5] 111 + (uint64_t)a[8] * b[4] 112 + (uint64_t)a[9] * b[3]; 113 VERIFY_BITS(d, 63); 114 /* [d 0 0 t9 0 0 0 0 0 0 c t1 t0] = [p12 p11 p10 p9 0 0 0 0 0 0 p2 p1 p0] */ 115 u2 = d & M; d >>= 26; c += u2 * R0; 116 VERIFY_BITS(u2, 26); 117 VERIFY_BITS(d, 37); 118 VERIFY_BITS(c, 63); 119 /* [d u2 0 0 t9 0 0 0 0 0 0 c-u2*R0 t1 t0] = [p12 p11 p10 p9 0 0 0 0 0 0 p2 p1 p0] */ 120 t2 = c & M; c >>= 26; c += u2 * R1; 121 VERIFY_BITS(t2, 26); 122 VERIFY_BITS(c, 38); 123 /* [d u2 0 0 t9 0 0 0 0 0 c-u2*R1 t2-u2*R0 t1 t0] = [p12 p11 p10 p9 0 0 0 0 0 0 p2 p1 p0] */ 124 /* [d 0 0 0 t9 0 0 0 0 0 c t2 t1 t0] = [p12 p11 p10 p9 0 0 0 0 0 0 p2 p1 p0] */ 125 126 c += (uint64_t)a[0] * b[3] 127 + (uint64_t)a[1] * b[2] 128 + (uint64_t)a[2] * b[1] 129 + (uint64_t)a[3] * b[0]; 130 VERIFY_BITS(c, 63); 131 /* [d 0 0 0 t9 0 0 0 0 0 c t2 t1 t0] = [p12 p11 p10 p9 0 0 0 0 0 p3 p2 p1 p0] */ 132 d += (uint64_t)a[4] * b[9] 133 + (uint64_t)a[5] * b[8] 134 + (uint64_t)a[6] * b[7] 135 + (uint64_t)a[7] * b[6] 136 + (uint64_t)a[8] * b[5] 137 + (uint64_t)a[9] * b[4]; 138 VERIFY_BITS(d, 63); 139 /* [d 0 0 0 t9 0 0 0 0 0 c t2 t1 t0] = [p13 p12 p11 p10 p9 0 0 0 0 0 p3 p2 p1 p0] */ 140 u3 = d & M; d >>= 26; c += u3 * R0; 141 VERIFY_BITS(u3, 26); 142 VERIFY_BITS(d, 37); 143 /* VERIFY_BITS(c, 64); */ 144 /* [d u3 0 0 0 t9 0 0 0 0 0 c-u3*R0 t2 t1 t0] = [p13 p12 p11 p10 p9 0 0 0 0 0 p3 p2 p1 p0] */ 145 t3 = c & M; c >>= 26; c += u3 * R1; 146 VERIFY_BITS(t3, 26); 147 VERIFY_BITS(c, 39); 148 /* [d u3 0 0 0 t9 0 0 0 0 c-u3*R1 t3-u3*R0 t2 t1 t0] = [p13 p12 p11 p10 p9 0 0 0 0 0 p3 p2 p1 p0] */ 149 /* [d 0 0 0 0 t9 0 0 0 0 c t3 t2 t1 t0] = [p13 p12 p11 p10 p9 0 0 0 0 0 p3 p2 p1 p0] */ 150 151 c += (uint64_t)a[0] * b[4] 152 + (uint64_t)a[1] * b[3] 153 + (uint64_t)a[2] * b[2] 154 + (uint64_t)a[3] * b[1] 155 + (uint64_t)a[4] * b[0]; 156 VERIFY_BITS(c, 63); 157 /* [d 0 0 0 0 t9 0 0 0 0 c t3 t2 t1 t0] = [p13 p12 p11 p10 p9 0 0 0 0 p4 p3 p2 p1 p0] */ 158 d += (uint64_t)a[5] * b[9] 159 + (uint64_t)a[6] * b[8] 160 + (uint64_t)a[7] * b[7] 161 + (uint64_t)a[8] * b[6] 162 + (uint64_t)a[9] * b[5]; 163 VERIFY_BITS(d, 62); 164 /* [d 0 0 0 0 t9 0 0 0 0 c t3 t2 t1 t0] = [p14 p13 p12 p11 p10 p9 0 0 0 0 p4 p3 p2 p1 p0] */ 165 u4 = d & M; d >>= 26; c += u4 * R0; 166 VERIFY_BITS(u4, 26); 167 VERIFY_BITS(d, 36); 168 /* VERIFY_BITS(c, 64); */ 169 /* [d u4 0 0 0 0 t9 0 0 0 0 c-u4*R0 t3 t2 t1 t0] = [p14 p13 p12 p11 p10 p9 0 0 0 0 p4 p3 p2 p1 p0] */ 170 t4 = c & M; c >>= 26; c += u4 * R1; 171 VERIFY_BITS(t4, 26); 172 VERIFY_BITS(c, 39); 173 /* [d u4 0 0 0 0 t9 0 0 0 c-u4*R1 t4-u4*R0 t3 t2 t1 t0] = [p14 p13 p12 p11 p10 p9 0 0 0 0 p4 p3 p2 p1 p0] */ 174 /* [d 0 0 0 0 0 t9 0 0 0 c t4 t3 t2 t1 t0] = [p14 p13 p12 p11 p10 p9 0 0 0 0 p4 p3 p2 p1 p0] */ 175 176 c += (uint64_t)a[0] * b[5] 177 + (uint64_t)a[1] * b[4] 178 + (uint64_t)a[2] * b[3] 179 + (uint64_t)a[3] * b[2] 180 + (uint64_t)a[4] * b[1] 181 + (uint64_t)a[5] * b[0]; 182 VERIFY_BITS(c, 63); 183 /* [d 0 0 0 0 0 t9 0 0 0 c t4 t3 t2 t1 t0] = [p14 p13 p12 p11 p10 p9 0 0 0 p5 p4 p3 p2 p1 p0] */ 184 d += (uint64_t)a[6] * b[9] 185 + (uint64_t)a[7] * b[8] 186 + (uint64_t)a[8] * b[7] 187 + (uint64_t)a[9] * b[6]; 188 VERIFY_BITS(d, 62); 189 /* [d 0 0 0 0 0 t9 0 0 0 c t4 t3 t2 t1 t0] = [p15 p14 p13 p12 p11 p10 p9 0 0 0 p5 p4 p3 p2 p1 p0] */ 190 u5 = d & M; d >>= 26; c += u5 * R0; 191 VERIFY_BITS(u5, 26); 192 VERIFY_BITS(d, 36); 193 /* VERIFY_BITS(c, 64); */ 194 /* [d u5 0 0 0 0 0 t9 0 0 0 c-u5*R0 t4 t3 t2 t1 t0] = [p15 p14 p13 p12 p11 p10 p9 0 0 0 p5 p4 p3 p2 p1 p0] */ 195 t5 = c & M; c >>= 26; c += u5 * R1; 196 VERIFY_BITS(t5, 26); 197 VERIFY_BITS(c, 39); 198 /* [d u5 0 0 0 0 0 t9 0 0 c-u5*R1 t5-u5*R0 t4 t3 t2 t1 t0] = [p15 p14 p13 p12 p11 p10 p9 0 0 0 p5 p4 p3 p2 p1 p0] */ 199 /* [d 0 0 0 0 0 0 t9 0 0 c t5 t4 t3 t2 t1 t0] = [p15 p14 p13 p12 p11 p10 p9 0 0 0 p5 p4 p3 p2 p1 p0] */ 200 201 c += (uint64_t)a[0] * b[6] 202 + (uint64_t)a[1] * b[5] 203 + (uint64_t)a[2] * b[4] 204 + (uint64_t)a[3] * b[3] 205 + (uint64_t)a[4] * b[2] 206 + (uint64_t)a[5] * b[1] 207 + (uint64_t)a[6] * b[0]; 208 VERIFY_BITS(c, 63); 209 /* [d 0 0 0 0 0 0 t9 0 0 c t5 t4 t3 t2 t1 t0] = [p15 p14 p13 p12 p11 p10 p9 0 0 p6 p5 p4 p3 p2 p1 p0] */ 210 d += (uint64_t)a[7] * b[9] 211 + (uint64_t)a[8] * b[8] 212 + (uint64_t)a[9] * b[7]; 213 VERIFY_BITS(d, 61); 214 /* [d 0 0 0 0 0 0 t9 0 0 c t5 t4 t3 t2 t1 t0] = [p16 p15 p14 p13 p12 p11 p10 p9 0 0 p6 p5 p4 p3 p2 p1 p0] */ 215 u6 = d & M; d >>= 26; c += u6 * R0; 216 VERIFY_BITS(u6, 26); 217 VERIFY_BITS(d, 35); 218 /* VERIFY_BITS(c, 64); */ 219 /* [d u6 0 0 0 0 0 0 t9 0 0 c-u6*R0 t5 t4 t3 t2 t1 t0] = [p16 p15 p14 p13 p12 p11 p10 p9 0 0 p6 p5 p4 p3 p2 p1 p0] */ 220 t6 = c & M; c >>= 26; c += u6 * R1; 221 VERIFY_BITS(t6, 26); 222 VERIFY_BITS(c, 39); 223 /* [d u6 0 0 0 0 0 0 t9 0 c-u6*R1 t6-u6*R0 t5 t4 t3 t2 t1 t0] = [p16 p15 p14 p13 p12 p11 p10 p9 0 0 p6 p5 p4 p3 p2 p1 p0] */ 224 /* [d 0 0 0 0 0 0 0 t9 0 c t6 t5 t4 t3 t2 t1 t0] = [p16 p15 p14 p13 p12 p11 p10 p9 0 0 p6 p5 p4 p3 p2 p1 p0] */ 225 226 c += (uint64_t)a[0] * b[7] 227 + (uint64_t)a[1] * b[6] 228 + (uint64_t)a[2] * b[5] 229 + (uint64_t)a[3] * b[4] 230 + (uint64_t)a[4] * b[3] 231 + (uint64_t)a[5] * b[2] 232 + (uint64_t)a[6] * b[1] 233 + (uint64_t)a[7] * b[0]; 234 /* VERIFY_BITS(c, 64); */ 235 VERIFY_CHECK(c <= 0x8000007C00000007ULL); 236 /* [d 0 0 0 0 0 0 0 t9 0 c t6 t5 t4 t3 t2 t1 t0] = [p16 p15 p14 p13 p12 p11 p10 p9 0 p7 p6 p5 p4 p3 p2 p1 p0] */ 237 d += (uint64_t)a[8] * b[9] 238 + (uint64_t)a[9] * b[8]; 239 VERIFY_BITS(d, 58); 240 /* [d 0 0 0 0 0 0 0 t9 0 c t6 t5 t4 t3 t2 t1 t0] = [p17 p16 p15 p14 p13 p12 p11 p10 p9 0 p7 p6 p5 p4 p3 p2 p1 p0] */ 241 u7 = d & M; d >>= 26; c += u7 * R0; 242 VERIFY_BITS(u7, 26); 243 VERIFY_BITS(d, 32); 244 /* VERIFY_BITS(c, 64); */ 245 VERIFY_CHECK(c <= 0x800001703FFFC2F7ULL); 246 /* [d u7 0 0 0 0 0 0 0 t9 0 c-u7*R0 t6 t5 t4 t3 t2 t1 t0] = [p17 p16 p15 p14 p13 p12 p11 p10 p9 0 p7 p6 p5 p4 p3 p2 p1 p0] */ 247 t7 = c & M; c >>= 26; c += u7 * R1; 248 VERIFY_BITS(t7, 26); 249 VERIFY_BITS(c, 38); 250 /* [d u7 0 0 0 0 0 0 0 t9 c-u7*R1 t7-u7*R0 t6 t5 t4 t3 t2 t1 t0] = [p17 p16 p15 p14 p13 p12 p11 p10 p9 0 p7 p6 p5 p4 p3 p2 p1 p0] */ 251 /* [d 0 0 0 0 0 0 0 0 t9 c t7 t6 t5 t4 t3 t2 t1 t0] = [p17 p16 p15 p14 p13 p12 p11 p10 p9 0 p7 p6 p5 p4 p3 p2 p1 p0] */ 252 253 c += (uint64_t)a[0] * b[8] 254 + (uint64_t)a[1] * b[7] 255 + (uint64_t)a[2] * b[6] 256 + (uint64_t)a[3] * b[5] 257 + (uint64_t)a[4] * b[4] 258 + (uint64_t)a[5] * b[3] 259 + (uint64_t)a[6] * b[2] 260 + (uint64_t)a[7] * b[1] 261 + (uint64_t)a[8] * b[0]; 262 /* VERIFY_BITS(c, 64); */ 263 VERIFY_CHECK(c <= 0x9000007B80000008ULL); 264 /* [d 0 0 0 0 0 0 0 0 t9 c t7 t6 t5 t4 t3 t2 t1 t0] = [p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 265 d += (uint64_t)a[9] * b[9]; 266 VERIFY_BITS(d, 57); 267 /* [d 0 0 0 0 0 0 0 0 t9 c t7 t6 t5 t4 t3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 268 u8 = d & M; d >>= 26; c += u8 * R0; 269 VERIFY_BITS(u8, 26); 270 VERIFY_BITS(d, 31); 271 /* VERIFY_BITS(c, 64); */ 272 VERIFY_CHECK(c <= 0x9000016FBFFFC2F8ULL); 273 /* [d u8 0 0 0 0 0 0 0 0 t9 c-u8*R0 t7 t6 t5 t4 t3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 274 275 r[3] = t3; 276 VERIFY_BITS(r[3], 26); 277 /* [d u8 0 0 0 0 0 0 0 0 t9 c-u8*R0 t7 t6 t5 t4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 278 r[4] = t4; 279 VERIFY_BITS(r[4], 26); 280 /* [d u8 0 0 0 0 0 0 0 0 t9 c-u8*R0 t7 t6 t5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 281 r[5] = t5; 282 VERIFY_BITS(r[5], 26); 283 /* [d u8 0 0 0 0 0 0 0 0 t9 c-u8*R0 t7 t6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 284 r[6] = t6; 285 VERIFY_BITS(r[6], 26); 286 /* [d u8 0 0 0 0 0 0 0 0 t9 c-u8*R0 t7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 287 r[7] = t7; 288 VERIFY_BITS(r[7], 26); 289 /* [d u8 0 0 0 0 0 0 0 0 t9 c-u8*R0 r7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 290 291 r[8] = c & M; c >>= 26; c += u8 * R1; 292 VERIFY_BITS(r[8], 26); 293 VERIFY_BITS(c, 39); 294 /* [d u8 0 0 0 0 0 0 0 0 t9+c-u8*R1 r8-u8*R0 r7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 295 /* [d 0 0 0 0 0 0 0 0 0 t9+c r8 r7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 296 c += d * R0 + t9; 297 VERIFY_BITS(c, 45); 298 /* [d 0 0 0 0 0 0 0 0 0 c-d*R0 r8 r7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 299 r[9] = c & (M >> 4); c >>= 22; c += d * (R1 << 4); 300 VERIFY_BITS(r[9], 22); 301 VERIFY_BITS(c, 46); 302 /* [d 0 0 0 0 0 0 0 0 r9+((c-d*R1<<4)<<22)-d*R0 r8 r7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 303 /* [d 0 0 0 0 0 0 0 -d*R1 r9+(c<<22)-d*R0 r8 r7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 304 /* [r9+(c<<22) r8 r7 r6 r5 r4 r3 t2 t1 t0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 305 306 d = c * (R0 >> 4) + t0; 307 VERIFY_BITS(d, 56); 308 /* [r9+(c<<22) r8 r7 r6 r5 r4 r3 t2 t1 d-c*R0>>4] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 309 r[0] = d & M; d >>= 26; 310 VERIFY_BITS(r[0], 26); 311 VERIFY_BITS(d, 30); 312 /* [r9+(c<<22) r8 r7 r6 r5 r4 r3 t2 t1+d r0-c*R0>>4] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 313 d += c * (R1 >> 4) + t1; 314 VERIFY_BITS(d, 53); 315 VERIFY_CHECK(d <= 0x10000003FFFFBFULL); 316 /* [r9+(c<<22) r8 r7 r6 r5 r4 r3 t2 d-c*R1>>4 r0-c*R0>>4] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 317 /* [r9 r8 r7 r6 r5 r4 r3 t2 d r0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 318 r[1] = d & M; d >>= 26; 319 VERIFY_BITS(r[1], 26); 320 VERIFY_BITS(d, 27); 321 VERIFY_CHECK(d <= 0x4000000ULL); 322 /* [r9 r8 r7 r6 r5 r4 r3 t2+d r1 r0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 323 d += t2; 324 VERIFY_BITS(d, 27); 325 /* [r9 r8 r7 r6 r5 r4 r3 d r1 r0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 326 r[2] = d; 327 VERIFY_BITS(r[2], 27); 328 /* [r9 r8 r7 r6 r5 r4 r3 r2 r1 r0] = [p18 p17 p16 p15 p14 p13 p12 p11 p10 p9 p8 p7 p6 p5 p4 p3 p2 p1 p0] */ 329 }
函数中第5行对M，R0，R1变量进行了初始化，M取值为了进行异或操作，使得每个limb保留26bit，R0和R1取值和secp256k1_fe_normalize中分析类似，只不过是将大数超过2260部分归一化到低位limb，即 2260 ≡ 0x1000003d10 mod p ≡ (1 << 36) + 0x3d10 mod p ≡ 0x400*226 + 0x3d10 mod p，归一化最低limb部分对应参数R0=0x3d10，次低limb部分对应参数R1=0x400。接下来7~26行是对输入参数的合法性检查，本节开头已经对此做出说明。之后第28~31行的注释说明了数组形式大数的具体含义，以及代码注释中px的含义，此外[x 0 0 0 0 0 0 0 0 0 0] = [x*R1 x*R0]，等式表示大数中超过2260部分x，可以通过之前提到的归一化操作，借助R0和R1归一化到最低两个位置，即x*2260 ≡ x*0x400*226 + x*0x3d10 mod p ≡ x*R1 + x*R0 mod p。之后34~43行先求了p9位置的累加和d，46行将累加和d的低26bit值t9留在了p9位置，而后对累加和进行右移位操作，保留高位作为进位放在了p10位置。51行求得了p0位置的值，第54~62行求p10位置累加和，并将该累加和加到了之前的进位d上，第65行，u0是累加和的低26bit值对应位置是p10，继续进行右移位操作保留高位值d（对应位置p11），该行最后一条语句提前为u0向低位归一化做准备（69行注释说明了大数乘积当前各位置对应的实际取值，因为u0还没有完成归一化，所以最低p0位置应为新t加上u0*R0之前的值），70行完成了高位u0向最低两位的归一化操作，此时c为向p1位置的进位，t0为p0位置的值（注释73是归一化之前的对应值，74行是u0完成归一化后各位置对应值，此时已经消去u0）。
接下来，76~77行继续求p1位置对应的大数乘法累加和，并将其加到之前的进位t上，80~87求p11位置对应的大数乘法累加和，并将其加到上一次的进位值d上，第90行保留d低26bit作为p11位置上实际值，并右移保留高位作为向p12位置的进位，同样改行最后一条语句为位置p11向次低位p1做归一化做准备，95行完成归一化操作并移除u1，并获取p1位置值t1。再之后，第101~124行获取p2位置值t2，第126~149行获取p3位置值t3，第151~174行获取p4位置值t4，第176~199行获取p5位置值t5，第201~224行获取p6位置值t6，第226~251行获取p7位置值t7，第253~295行获取p8位置t8，在此期间又将t3~t7存储到了最终结果r[3]~r[7]中（应该是程序做过极限测试，后续的归一化操作已经不会影响到该部分最终的取值）。
第296行将来自p8位置的进位和p18位置d对应归一化低位值以及最一开始获得的t9进行相加获取p9位置的值c，之后，第299行只保留低22bit值赋值给r[9]（最终归一化结果最高位置只保留22bit）,右移22bit保留高位c，最后将p18位置d对应归一化高位值和c进行相加（在此之前都是按照260位大数进行归一化操作，而进行到这一步，c已经相当于超出大数256位之后对应的值），接下来，将按照256位大数进行归一化操作，第306行首先c到最低位部分归一化操作，将对应值和t0相加，第309行将和d低26bit保存到r[0]，右移保留高位d，第313行将对应次低位部分归一化值和t1以及上一步的d进行累加，并更新为d，第318~326行按相同方式依次更新t[1]和t[2]，进行到这里所有归一化进位已经处理完毕算法结束。而大数平方函数secp256k1_fe_sqr_inner的处理过程和该函数基本一致，不再进行详细分析。
3.4 secp256k1_fe_sqrt
之前的二次剩余理论部分已经得出a如果是模p的二次剩余，则a(p+1)/4必是a的开平方，先直接给出函数实现
1 /* 2 * (p+1)/4 = 3 * 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 4 * 111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 5 * 111111111111110111111111111111111111100001100 6 * 可以将该二进制按连续1进行分块，分三块，连续1的个数分别是2，22，223 7 */ 8 static int secp256k1_fe_sqrt(secp256k1_fe *r, const secp256k1_fe *a) { 9 /** Given that p is congruent to 3 mod 4, we can compute the square root of 10 * a mod p as the (p+1)/4'th power of a. 11 * 12 * As (p+1)/4 is an even number, it will have the same result for a and for 13 * (-a). Only one of these two numbers actually has a square root however, 14 * so we test at the end by squaring and comparing to the input. 15 * Also because (p+1)/4 is an even number, the computed square root is 16 * itself always a square (a ** ((p+1)/4) is the square of a ** ((p+1)/8)). 17 */ 18 secp256k1_fe x2, x3, x6, x9, x11, x22, x44, x88, x176, x220, x223, t1; 19 int j; 20 21 /** The binary representation of (p + 1)/4 has 3 blocks of 1s, with lengths in 22 * { 2, 22, 223 }. Use an addition chain to calculate 2^n - 1 for each block: 23 * 1, [2], 3, 6, 9, 11, [22], 44, 88, 176, 220, [223] 24 */ 25 secp256k1_fe_sqr(&x2, a); 26 secp256k1_fe_mul(&x2, &x2, a); // a^3 0b11 27 28 secp256k1_fe_sqr(&x3, &x2); // a^6 0b110 29 secp256k1_fe_mul(&x3, &x3, a); // x3 a^7 0b111 30 31 x6 = x3; 32 for (j=0; j<3; j++) { // 0b111000 33 secp256k1_fe_sqr(&x6, &x6); 34 } 35 secp256k1_fe_mul(&x6, &x6, &x3); // x6 a^63 0b111111 36 37 x9 = x6; 38 for (j=0; j<3; j++) { // 0b111111000 39 secp256k1_fe_sqr(&x9, &x9); 40 } 41 secp256k1_fe_mul(&x9, &x9, &x3); // x9 a^511 0b111111111 42 43 x11 = x9; 44 for (j=0; j<2; j++) { // 0b11111111100 45 secp256k1_fe_sqr(&x11, &x11); 46 } 47 secp256k1_fe_mul(&x11, &x11, &x2); // x11 a^2047 0b11111111111 48 49 x22 = x11; 50 for (j=0; j<11; j++) { // 0b1111111111100000000000 51 secp256k1_fe_sqr(&x22, &x22); 52 } 53 secp256k1_fe_mul(&x22, &x22, &x11); // x22 0b1111111111111111111111 54 55 x44 = x22; 56 for (j=0; j<22; j++) { // 0b11111111111111111111110000000000000000000000 57 secp256k1_fe_sqr(&x44, &x44); 58 } 59 secp256k1_fe_mul(&x44, &x44, &x22); // 0b11111111111111111111111111111111111111111111 60 61 x88 = x44; 62 for (j=0; j<44; j++) { 63 secp256k1_fe_sqr(&x88, &x88); 64 } 65 secp256k1_fe_mul(&x88, &x88, &x44); // 0b1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 66 67 x176 = x88; 68 for (j=0; j<88; j++) { 69 secp256k1_fe_sqr(&x176, &x176); 70 } 71 secp256k1_fe_mul(&x176, &x176, &x88); 72 73 x220 = x176; 74 for (j=0; j<44; j++) { 75 secp256k1_fe_sqr(&x220, &x220); 76 } 77 secp256k1_fe_mul(&x220, &x220, &x44); 78 79 x223 = x220; 80 for (j=0; j<3; j++) { 81 secp256k1_fe_sqr(&x223, &x223); 82 } 83 secp256k1_fe_mul(&x223, &x223, &x3); 84 85 /* The final result is then assembled using a sliding window over the blocks. */ 86 t1 = x223; 87 for (j=0; j<23; j++) { // 左移23bit，给01111111111111111111111留出位置 88 secp256k1_fe_sqr(&t1, &t1); 89 } 90 secp256k1_fe_mul(&t1, &t1, &x22); 91 for (j=0; j<6; j++) { // 000011 length is 6 92 secp256k1_fe_sqr(&t1, &t1); 93 } 94 secp256k1_fe_mul(&t1, &t1, &x2); 95 secp256k1_fe_sqr(&t1, &t1); // for the last two 0 96 secp256k1_fe_sqr(r, &t1); 97 98 /* Check that a square root was actually calculated */ 99 secp256k1_fe_sqr(&t1, r); 100 return secp256k1_fe_equal(&t1, a); 101 }
如2~6行注释中做出的解释，(p+1)/4二进制共有3段连续的二进制1，算法正是围绕着该特点进行实现的。
第25行求x2=a2，26行求得x2=a3=a0b11，第28~29行求得x3=a7=a0b111，第31~35行求得x6=a63=a0b111111，之后类似，直到第83行求得x223=a0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffff，至此(p+1)/4最高位部分连续223个二级制1部分的平方已经计算完毕，之后第86行~90行，首先通过平方操作实现“左移23位操作”，将接下来的22个连续二级制1部分内容乘到结果中（这里面包含一个0位）。然后第91行~94行，将最后2个连续1乘到结果内，之后95~96行实现最后“左移2位操作”。最后99~100行结合二次剩余理论，通过开方的平方是否和原a相等，判断a是否真的“可开方“。
3.5 secp256k1_fe_inv
该函数的实现和secp256k1_fe_sqrt类似，不再进行详细分析，只给出完整的函数实现。
1 /* 2 * 在有限域中求数a的模逆，由欧拉公式可知a^fai(p) = 1 mod p，即a^(p-1) = 1 mod p 3 * 所以a^(p-2)即为a的逆元,0b11111111111111111111111111111111111111111111111111111111111111111 4 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 5 1111111111111111111111111111111111111111111111111111111111011111111111111111111110000101101 6 */ 7 static void secp256k1_fe_inv(secp256k1_fe *r, const secp256k1_fe *a) { 8 secp256k1_fe x2, x3, x6, x9, x11, x22, x44, x88, x176, x220, x223, t1; 9 int j; 10 11 /** The binary representation of (p - 2) has 5 blocks of 1s, with lengths in 12 * { 1, 2, 1, 22, 223 }. Use an addition chain to calculate 2^n - 1 for each block: 13 * [1], [2], 3, 6, 9, 11, [22], 44, 88, 176, 220, [223] 14 */ 15 secp256k1_fe_sqr(&x2, a); 16 secp256k1_fe_mul(&x2, &x2, a); 17 18 secp256k1_fe_sqr(&x3, &x2); 19 secp256k1_fe_mul(&x3, &x3, a); 20 21 x6 = x3; 22 for (j=0; j<3; j++) { 23 secp256k1_fe_sqr(&x6, &x6); 24 } 25 secp256k1_fe_mul(&x6, &x6, &x3); 26 27 x9 = x6; 28 for (j=0; j<3; j++) { 29 secp256k1_fe_sqr(&x9, &x9); 30 } 31 secp256k1_fe_mul(&x9, &x9, &x3); 32 33 x11 = x9; 34 for (j=0; j<2; j++) { 35 secp256k1_fe_sqr(&x11, &x11); 36 } 37 secp256k1_fe_mul(&x11, &x11, &x2); 38 39 x22 = x11; 40 for (j=0; j<11; j++) { 41 secp256k1_fe_sqr(&x22, &x22); 42 } 43 secp256k1_fe_mul(&x22, &x22, &x11); 44 45 x44 = x22; 46 for (j=0; j<22; j++) { 47 secp256k1_fe_sqr(&x44, &x44); 48 } 49 secp256k1_fe_mul(&x44, &x44, &x22); 50 51 x88 = x44; 52 for (j=0; j<44; j++) { 53 secp256k1_fe_sqr(&x88, &x88); 54 } 55 secp256k1_fe_mul(&x88, &x88, &x44); 56 57 x176 = x88; 58 for (j=0; j<88; j++) { 59 secp256k1_fe_sqr(&x176, &x176); 60 } 61 secp256k1_fe_mul(&x176, &x176, &x88); 62 63 x220 = x176; 64 for (j=0; j<44; j++) { 65 secp256k1_fe_sqr(&x220, &x220); 66 } 67 secp256k1_fe_mul(&x220, &x220, &x44); 68 69 x223 = x220; 70 for (j=0; j<3; j++) { 71 secp256k1_fe_sqr(&x223, &x223); 72 } 73 secp256k1_fe_mul(&x223, &x223, &x3); 74 75 /* The final result is then assembled using a sliding window over the blocks. */ 76 t1 = x223; 77 for (j=0; j<23; j++) { 78 secp256k1_fe_sqr(&t1, &t1); 79 } 80 secp256k1_fe_mul(&t1, &t1, &x22); 81 for (j=0; j<5; j++) { 82 secp256k1_fe_sqr(&t1, &t1); 83 } 84 secp256k1_fe_mul(&t1, &t1, a); 85 for (j=0; j<3; j++) { 86 secp256k1_fe_sqr(&t1, &t1); 87 } 88 secp256k1_fe_mul(&t1, &t1, &x2); 89 for (j=0; j<2; j++) { 90 secp256k1_fe_sqr(&t1, &t1); 91 } 92 secp256k1_fe_mul(r, a, &t1); 93 }
原创作者: zhaoweiwei 转载于: https://www.cnblogs.com/zhaoweiwei/p/18945891/secp256k1_2
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
