怎樣建設(shè)學(xué)校網(wǎng)站首頁,陽春ycqq人才招聘信息,網(wǎng)站手機(jī)端做排名,畢節(jié)市網(wǎng)站建設(shè)magnitude及normalized由于當(dāng)前許多項(xiàng)目都用到secp256k1庫#xff0c;比特幣作為體量最大的數(shù)字貨幣項(xiàng)目#xff0c;這里建議直接參考bitcoin-core提供的最新secp256k1源碼。仍以field的10x26實(shí)現(xiàn)版本為例#xff0c;相關(guān)定義如下#xff1a;復(fù)制代碼/** This field implem…magnitude及normalized由于當(dāng)前許多項(xiàng)目都用到secp256k1庫比特幣作為體量最大的數(shù)字貨幣項(xiàng)目這里建議直接參考bitcoin-core提供的最新secp256k1源碼。仍以field的10x26實(shí)現(xiàn)版本為例相關(guān)定義如下復(fù)制代碼/** 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;復(fù)制代碼對于magnitude可稱其為“量級”當(dāng)m0時(shí)這時(shí)n[i] 0(i0...9)由此可知此時(shí)secp256k1_fe大數(shù)必為0當(dāng)m1時(shí)n[i] 2*(2^26 - 1)對于i0...8n[9] 2*(2^22 - 1)有些說法是當(dāng)m1時(shí)是將大數(shù)限制到[0,到2p)范圍內(nèi)這是不準(zhǔn)確的此時(shí)secp256k1_fe大數(shù)范圍是[0, 2*(2^256-1)]上限是稍大于2p的magnitude設(shè)計(jì)本質(zhì)是“存儲(chǔ)約束”而非“模p范圍”通過magnitude約束將大數(shù)控制在可高效約簡的范圍內(nèi)magnitude表示大概是2m個(gè)p的量級。對于normalized可稱其為“規(guī)范化歸一化”從注釋中可知每個(gè)n[i]都小于或等于其對應(yīng)的MASK且sum(n[i]) p即規(guī)范化后的大數(shù)是[0, p)之內(nèi)的數(shù)。由此可知規(guī)范化的大數(shù)一定是magnitude為0或1的數(shù)但是magnitude為1的數(shù)不一定是規(guī)范化的大數(shù)。函數(shù)secp256k1_fe_normalize對大數(shù)實(shí)現(xiàn)規(guī)范化操作其本質(zhì)是將大數(shù)看作為特殊的2^26進(jìn)制大數(shù)即N∑ni*wi這里ni 2 * m * (2^26 - 1) for i0..8n9 2 * m * (2^22 - 1)wi2^(i*26)之所以說它特殊是因?yàn)橛捎趍的存在每個(gè)位上數(shù)值是可以大于或等于基數(shù)2^26的正常情況下每個(gè)位上取值小于基數(shù)如10進(jìn)制數(shù)每個(gè)位上數(shù)值都小于基數(shù)10另外最高位取值有特殊限制在secp256k1_fe_normalize函數(shù)中有兩部分操作第一部分先將大數(shù)A規(guī)范為最多為257位的大數(shù)Bm確定時(shí)該數(shù)對應(yīng)確定最大值后續(xù)詳解第二部分將大數(shù)B規(guī)范為小于模數(shù)p的大數(shù)。在secp256k1源碼中其實(shí)定義了各個(gè)參數(shù)大數(shù)中magnitude的最大值在定義VERIFY宏時(shí)會(huì)對參數(shù)合法性進(jìn)行檢查。#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當(dāng)m0時(shí)顯而易見大數(shù)即為0當(dāng)m1且m4時(shí)其實(shí)大數(shù)對應(yīng)最大值分別為pow(2,256)-1)*2*10x1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffepow(2,256)-1)*2*20x3fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffcpow(2,256)-1)*2*30x5fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffapow(2,256)-1)*2*40x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff8當(dāng)m取值為4時(shí)對應(yīng)最大值在調(diào)用規(guī)范化函數(shù)時(shí)第一部分規(guī)范完以后對應(yīng)最大值為0x100003c00000f000000000001e000003c00000f000003c00000f0000f3c00393e為一個(gè)257位的大數(shù)函數(shù)secp256k1_fe_normalize_weak其實(shí)就是secp256k1_fe_normalize第一部分執(zhí)行完以后的輸出該值就是weak函數(shù)在m4下輸出的最大值該值本身為m1的大數(shù)。2 常用函數(shù)magnitude值分析先分析大數(shù)取反原型函數(shù)secp256k1_fe_negate其原型如下復(fù)制代碼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];}復(fù)制代碼這里m是輸入?yún)?shù)a的magnitude當(dāng)m0時(shí)輸出參數(shù)r為2p顯然其對應(yīng)的m1當(dāng)輸入?yún)?shù)a的m1時(shí)其取值范圍為[1, 0x1fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe]則2*(m1)*p-a取值范圍是[0x1fffffffffffffffffffffffffffffffffffffffffffffffffffffffbfffff0be, 0x3fffffffffffffffffffffffffffffffffffffffffffffffffffffffbfffff0bb]所以輸出r的m2所以對于任意輸入m輸出r的mm1。對于secp256k1_fe_add函數(shù)其源碼如下secp256k1_fe_add顯而易見對于輸出r來說m(r)m(r)m(a)r和a都是secp256k1_fe型數(shù)據(jù)。對于secp256k1_fe_mul_int(r, a)a為int型數(shù)據(jù)m(r)m(r)*a。對于secp256k1_fe_mul/secp256k1_fe_sqr函數(shù)在實(shí)現(xiàn)中包含內(nèi)部規(guī)約會(huì)把輸出參數(shù)的“量級”規(guī)約回magnitude 1使輸出可以當(dāng)作normalized-like使用。下面以函數(shù)為例給出各個(gè)步驟的詳細(xì)magnitude值secp256k1_gej_add_var首先以輸入?yún)?shù)ab中的xyz都為magnitude1normalized-like為前提在此基礎(chǔ)上給出按代碼順序的表步 代碼操作 輸入 m 輸出 m上限 說明1 secp256k1_fe_sqr(z22, b-z) b-z:1 1 fe_sqr → 規(guī)約為 12 secp256k1_fe_sqr(z12, a-z) a-z:1 13 secp256k1_fe_mul(u1, a-x, z22) 1,1 1 fe_mul → 規(guī)約為 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)) ... — — 特殊分支不進(jìn)入一般流程14 secp256k1_fe_mul(t, h, b-z) h:3, b.z:1 1 fe_mul → 規(guī)約為 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 → 規(guī)約為 118 secp256k1_fe_negate(h2, h2, 1) h2:1 219 secp256k1_fe_mul(h3, h2, h) h2:2, h:3 1 fe_mul → 規(guī)約為 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 → 規(guī)約為 127 secp256k1_fe_mul(h3, h3, s1) h3:1, s1:1 128 secp256k1_fe_add(r-y, h3) r-y:1, h3:1 229 函數(shù)返回 r-x:4, r-y:2, r-z:1 — 這些是按最小合法上界得到的值未再 normalize 前由表可知對于secp256k1_gej_add_var函數(shù)雖然運(yùn)行過程中零時(shí)變量可能會(huì)出現(xiàn)相對較高的magnitude值但最終返回值r-x的magnitude值是滿足SECP256K1_GE_X_MAGNITUDE_MAX值的同理yz也滿足SECP256K1_GEJ_Y_MAGNITUDE_MAX和SECP256K1_GEJ_Z_MAGNITUDE_MAX限制。3 大數(shù)求逆最新版大數(shù)求逆函數(shù)實(shí)現(xiàn)原型如下復(fù)制代碼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 }復(fù)制代碼3.1 高層概覽1. 主要步驟函數(shù)做的是對域元素x求模逆即求r x-1 mod p這里p 2^256 - 2^32 - 977其實(shí)現(xiàn)流程如下1先把tmpx規(guī)范化保證tmp是庫約定的limb范圍下唯一表示2把規(guī)范化的tmp10x26-bit limbs轉(zhuǎn)換成一種signed30的中間表示若干個(gè)30-bit的有符號limbs以便后面用32-bit算法實(shí)現(xiàn)逆元算法3在signed30表示上運(yùn)行secp256k1_modinv32——這是一個(gè)針對30-bit/32-bit limbs優(yōu)化、并做過constant-time處理的模逆實(shí)現(xiàn)基于safegcd/division-steps改進(jìn)的歐幾里得/半GCD的思想4把得到的signed30結(jié)果再轉(zhuǎn)換回庫的secp256k1_fe10x26表示并寫入r該轉(zhuǎn)換同時(shí)完成必要的約減/規(guī)范化。2. divsteps算法divsteps是一種高效計(jì)算GCDGreatest Common Divisor最大公約數(shù)的方法它通過連續(xù)多次除法步驟減少迭代次數(shù)特別適合硬件實(shí)現(xiàn)和大數(shù)計(jì)算。它是GCD算法的優(yōu)化版本主要特點(diǎn)利用二進(jìn)制表示的優(yōu)勢通過右移操作快速消除因子2一次迭代處理多個(gè) 減法 - 移位 步驟減少循環(huán)次數(shù)使用比較和差值運(yùn)算替代昂貴的除法操作求u和v最大公約數(shù)算法步驟如下1消除因子2統(tǒng)計(jì)并移除兩數(shù)中的所有因子2設(shè)a u / 2^sab v / 2^sb這里sa是u可整除2的最大次數(shù)sb是v可整除2的最大次數(shù)2divsteps迭代比較a和b計(jì)算差值d | a - b |消除d中的所有因子2得到d用min(a, b)和d替換(a, b)重復(fù)直到a b3得出最終結(jié)果gcd(u, v) b * 2^min(sa, sb)算法基于以下假設(shè)如果gcd(u, v)是最大公約數(shù)那么它可以分成兩部分的乘積一部分是2的整數(shù)次冪另一部分是非2的倍數(shù)。步驟1其實(shí)是獲取u、v最大公因數(shù)中能2整數(shù)次冪部分即2min(sa, sb)步驟2獲取非2的倍數(shù)部分。步驟1比較好理解這里分析下步驟2假設(shè)得出最后一步結(jié)論ab時(shí)具體值為ab則a必然滿足a - b b*2x則有a b*(12x)即a是b的倍數(shù)則在上一步必有db*2y a - a由此a a b*2y即a也為b的倍數(shù)由此類推可知最一開始的a和b閉然都是b的倍數(shù)。由以上推理可知最終gcd(u, v) b * 2^min(sa, sb)。以下是算法示例代碼divsteps計(jì)算示例gcd(270, 324)1初始值u270v3242移除因子2sa1sb2 → u135, v813第一次迭代d135-8154移除因子2: sd1 → d27新值: u81, v274第二次迭代d81-2754移除因子2: sd1 → d27新值: u27, v275循環(huán)結(jié)束結(jié)果gcd(270, 324) 27*2^1 54。3.2 逐步詳解步驟1將大數(shù)進(jìn)行規(guī)范化這一步由secp256k1_fe_normalize實(shí)現(xiàn)該函數(shù)已經(jīng)在之前文章中詳細(xì)介紹這里不再進(jìn)行分析。步驟2調(diào)用secp256k1_fe_to_signed30函數(shù)實(shí)現(xiàn)該函數(shù)源碼如下secp256k1_fe_to_signed30該函數(shù)實(shí)現(xiàn)很清晰就是把原有的10個(gè)無符號26bit數(shù)重新按30bit窗口打包到一組32-bit有符號整數(shù)中共9個(gè)元素每個(gè)元素保存30位有效位并且用int32_t保存以允許出現(xiàn)負(fù)數(shù)中間值。步驟3調(diào)用secp256k1_modinv32函數(shù)求模逆該步驟是算法的核心后續(xù)詳細(xì)解析。步驟4調(diào)用secp256k1_fe_from_signed30函數(shù)將30bit格式逆元再轉(zhuǎn)換到標(biāo)準(zhǔn)10x26bit形式。1. 擴(kuò)展歐幾里得算法擴(kuò)展歐幾里得定理對于任意整數(shù)a和b必然存在整數(shù)x和y滿足貝祖等式a·x b·y gcd(a, b)且gcd(a, b)是能表示為a·x b·y形式的最小正整數(shù)。1基礎(chǔ)引理歐幾里得算法的余數(shù)性質(zhì)歐幾里得算法通過反復(fù)求余計(jì)算GCDgcd(a, b) gcd(b, a mod b)其中a mod b a - b*?a/b???為向下取整且0 ≤ a mod b |b|。2數(shù)學(xué)歸納法證明貝祖等式存在性歸納基礎(chǔ)當(dāng)b 0 時(shí)gcd(a, 0) a此時(shí)取 x 1y 0這里y可以取任意值只不過取0后續(xù)計(jì)算最簡單顯然滿足 a·1 0·0 a gcd(a, 0)。歸納步驟假設(shè)對于 (b, a mod b) 存在整數(shù) x 和 y 滿足b·x (a mod b)·y gcd(b, a mod b)根據(jù)余數(shù)定義 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) 也存在解。由數(shù)學(xué)歸納法對任意整數(shù) a, b 均存在這樣的 x, y。3最小性證明設(shè) d gcd(a, b)則 d 整除 a 和 b因此 d 整除 a·x b·y 的所有可能結(jié)果。即任何能表示為 a·x b·y 的整數(shù)都是 d 的倍數(shù)故 d 是其中最小的正整數(shù)。其實(shí)以上結(jié)論是基于一個(gè)前提對于任意兩個(gè)整數(shù)a和b在應(yīng)用歐幾里得算法若干步后都能轉(zhuǎn)化成將gcd(a, b)轉(zhuǎn)換成gcd(b, 0)的形式這個(gè)結(jié)論可以自己查閱相關(guān)整理過程又因?yàn)間cd(a, b) gcd(-a, -b)所以該結(jié)論對于所有非零整數(shù)都成立。接下來看看在算法層面如何推導(dǎo)系數(shù)x和y首先對于gcd(a, b)中的a和b來說都可以看作是各次除法迭代的“余數(shù)r”即gcd(a, b) gcd(b, a%b))如果將a%b看作是newrb看作是ra看作是oldr則有oldr q?r newr這里q是整數(shù)商a//b。接下來因?yàn)榍笙禂?shù)x和y其實(shí)也是一個(gè)迭代推導(dǎo)的過程假設(shè)當(dāng)前有image可以先把最一開始o(jì)ldr和r的由系數(shù)表示出來oldr 1*a 0*br 0*a 1*b即存在初始余數(shù)對(oldr, r) (a, b)初始a的系數(shù)對(oldx, x) (1, 0)初始b的系數(shù)對(oldy, y) (0, 1)滿足方程。再接下來繼續(xù)進(jìn)行迭代q oldr//rnewr oldr - q*r將(1)中oldr和(2)中的r代入該等式并進(jìn)行整理可得image由此可得newx oldx - q*xnewy oldy - q*y即newr仍是ab的線性組合且newxnewy是相應(yīng)系數(shù)。以下是完整算法偽代碼復(fù)制代碼(old_r, r) (a, b)(old_x, x) (1, 0) # 表示 a 的系數(shù)(old_y, y) (0, 1) # 表示 b 的系數(shù)while r ≠ 0:q old_r // r# 更新余數(shù)(old_r, r) (r, old_r - q * r)# 更新系數(shù)(old_x, x) (x, old_x - q * x)(old_y, y) (y, old_y - q * y)# 最終結(jié)果gcd(a,b) old_rx old_x, y old_y復(fù)制代碼以a240b46為例給出算法迭代過程復(fù)制代碼初始值(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)此時(shí)因r0得最終結(jié)果gcd(240, 46)2x-9y47復(fù)制代碼2. 二進(jìn)制GCD算法變種二進(jìn)制GCD算法有以下演進(jìn)時(shí)間線1961年: 原始二進(jìn)制GCD概念出現(xiàn)1967年: Josef Stein正式發(fā)表算法1970年代: Knuth在TAOCP中分析和優(yōu)化1980-90年代: 各種優(yōu)化變種出現(xiàn)包括這個(gè)delta版本2000年代: 被廣泛用于密碼學(xué)庫和大數(shù)運(yùn)算接下來重點(diǎn)分析帶有delta狀態(tài)變量的版本算法如下復(fù)制代碼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ù)制代碼該算法要求第一個(gè)參數(shù)f必須是奇數(shù)第3行并引入了狀態(tài)變量delta用于指導(dǎo)算法選擇不同的約減策略避免陷入低效循環(huán)。該算法證明分兩部分1先證明算法保持 gcd 不變正確性不破壞2然后給出一個(gè)勢函數(shù)potential并用它證明在每次若干步內(nèi)勢函數(shù)會(huì)嚴(yán)格下降從而不可能無限迭代——也就保證終止收斂。同時(shí)指出 delta 在防止振蕩 / 保證下降中的關(guān)鍵作用。第1部分結(jié)論很容易得出在結(jié)合fg奇偶情況下gcd(g, (g - f)//2)gcd(f, (g f)//2)gcd(f, g//2)這三種情況顯然和gcd(f, g)是一樣的。第2部分的證明比較復(fù)雜涉及到勢函數(shù)相關(guān)理論知識盲區(qū)請自行查閱相關(guān)資料這里不再詳細(xì)說明。總之相比如下未引入delta的算法一在某些情況下會(huì)不收斂上述算法是能保證收斂的。gcd_no_delta另外相比如下未引入delta的算法二fg直接判斷屬于“輸入依賴型分支”不同輸入會(huì)導(dǎo)致不同的分支走向容易在硬件層面產(chǎn)生 “分支預(yù)測錯(cuò)誤”尤其在加密算法中輸入的隨機(jī)性會(huì)放大這種問題。而delta引入后分支判斷雖然存在但更多依賴內(nèi)部狀態(tài)變量 delta 的符號delta 0 或 delta ≤ 0而非輸入 f 和 g 的直接比較delta 的更新規(guī)則是固定的1 - delta 或 1 delta其變化模式相對可預(yù)測降低了對輸入隨機(jī)性的敏感度更適合常數(shù)時(shí)間實(shí)現(xiàn)密碼學(xué)算法的關(guān)鍵要求。尤其在大數(shù)運(yùn)算情況下a b 這類比較操作需要完整的減法和符號判斷硬件層面是減法器 符號位檢測而delta的狀態(tài)更新1 - delta 或 1 delta是簡單的算術(shù)操作本質(zhì)是 delta 在0和1之間交替或遞增遞減計(jì)算成本更低。bad_gcd3. 由GCD求模逆對于質(zhì)數(shù)M大于2來說求x0xM的模逆因?yàn)镸是質(zhì)數(shù)所以gcd(M, x)1則存在a*x b*M 1即a*x 1 mod M所以結(jié)合以上1和2小節(jié)的內(nèi)容可以由在gcd算法求解過程中求得x系數(shù)a以下是詳細(xì)算法復(fù)制代碼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復(fù)制代碼該算法基于二進(jìn)制擴(kuò)展GCD算法通過不斷將問題規(guī)模減半來高效計(jì)算模逆元首先明確算法以下關(guān)鍵變量M: 模數(shù)必須是奇數(shù)x: 要求逆元的數(shù)f, g: 跟蹤GCD計(jì)算的兩個(gè)值d, e: 跟蹤系數(shù)維護(hù)關(guān)鍵不變量delta: 控制算法分支的狀態(tài)變量算法關(guān)鍵在于每次迭代都維護(hù)了以下等式λf ≡ d * x (mod M)g ≡ e * x (mod M)這意味著每次迭代d和e始終是f和g在模M下關(guān)于x的系數(shù)這樣只要保證初始時(shí)等式成立以及每次迭代更新fg時(shí)按計(jì)算規(guī)則得到的de和fg仍滿足以上等式這就能保證到最后求得最大公約數(shù)1時(shí)以上等式仍成立及找出了最終的逆元d。以循環(huán)迭代中第1種情況為例進(jìn)行分析根據(jù)規(guī)則更新后的f gg (g - f)//2也就是將g替換為新的f將(g-f)//2替換為新的g那么對應(yīng)的系數(shù)更新規(guī)則必須保證image由于更新后f g由最一開始的等式g ≡ e * x (mod M)可知要保證f ≡ d * x (mod M)成立取d e即可這正是算法中更新d的邏輯。算法中g(shù)的更新邏輯是g (g - f)//2將初始等式λ代入更新邏輯得image這就要求更新后的e滿足image等價(jià)于image而while迭代中第1種情況中的e的更新函數(shù)div2正是實(shí)現(xiàn)了該邏輯在函數(shù)中由于輸入e-d可能為奇數(shù)此時(shí)需要加上奇數(shù)M在模M下不影響最終結(jié)果使得和為偶數(shù)再除于2。同理可以分析while迭代中的另外兩個(gè)分支最終當(dāng)算法終止時(shí)有g(shù)cd(M, x) |f| 1則當(dāng)f 1時(shí)由f≡d*x≡1 mod M可知d即為x逆元當(dāng)f -1時(shí)由d*x≡-1 mod M可知-d即為x逆元正是算法返回值。4. 模逆算法進(jìn)一步優(yōu)化上述每次divstep迭代可以表示為矩陣乘法對向量[f g]和[d e]應(yīng)用如下的轉(zhuǎn)換矩陣1/2*timage其中根據(jù)迭代中分支的不同(u, v, q, r)分別有不同的取值(0, 2, -1, 1)(2, 0, 1, 1)(2, 0, 0, 1)在while每次迭代中都會(huì)用轉(zhuǎn)換矩陣乘于相應(yīng)的向量而在轉(zhuǎn)換矩陣中又包含除于2的操作這樣整個(gè)while循環(huán)會(huì)包含除于2N的操作對于f和g來說算法已經(jīng)保證了每次迭代的除法可以通過移位來實(shí)現(xiàn)而對于d和e來說這些個(gè)除法是非常費(fèi)時(shí)的操作所以可以在每次迭代的矩陣乘法時(shí)不執(zhí)行除于2的操作而把該除法操作留到最后步驟執(zhí)行除于2N操作由此得到以下函數(shù)復(fù)制代碼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)復(fù)制代碼后續(xù)還需進(jìn)行如下計(jì)算image由推導(dǎo)過程可知最后uvqr組成的矩陣乘于初始向量[f g]以后得出的向量每個(gè)元素都是2N的倍數(shù)可由初始值實(shí)際驗(yàn)算下所以有以下的更新函數(shù)復(fù)制代碼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復(fù)制代碼對于d和e來說上述結(jié)論并不成立這里需要一個(gè)和div2函數(shù)類似的在模M下除于2N的函數(shù)借助該函數(shù)可實(shí)現(xiàn)de的更新復(fù)制代碼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)復(fù)制代碼綜合上述所有即可給出執(zhí)行N divsteps版本的modinv函數(shù)復(fù)制代碼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復(fù)制代碼這意味著在實(shí)踐中我們將始終執(zhí)行多個(gè)N個(gè)div步驟。這不是問題因?yàn)橐坏ゞ0進(jìn)一步的divsteps就不再影響f、g、d或e只有δ繼續(xù)增加這就能保證不管輸入是何止算法運(yùn)行步驟一致從而保證了算法時(shí)間一致性。3.3 secp256k1_modinv32源碼分析結(jié)合上述內(nèi)容分析下模逆的核心函數(shù)secp256k1_modinv32其源碼如下復(fù)制代碼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);