數環演算法
① 解九連環一共用多少下
解開九連環共需要341步,只要上或下一個環,就算一步,不是在框架上滑動。
解一連環需要1步,解二連環需要2步,由此可知,解三連環需要5步,解四連環需要10步,解五連環需要21步,解六連環需要42步,解七連環需要85步,解八需要170步,解九連環需要341步,解十連環需要682步……以後的類推。
② 九連環資料,越多越好!
九連環是一種流傳於山西省的漢族民間的智力玩具。它用九個圓環相連成串,以解開為勝。據明代楊慎《丹鉛總錄》記載,曾以玉石為材料製成兩個互貫的圓環,「兩環互相貫為一,得其關捩,解之為二,又合而為一」。後來,以銅或鐵代替玉石,成為婦女兒童的玩具。它在中國差不多有二千年的歷史,卓文君在給司馬相如的信中有「九連環從中折斷」的句子。清代,《紅樓夢》中也有林黛玉巧解九連環的記載。周邦彥也留下關於九連環的名句「縱妙手、能解連環。」
起源
九連環
西漢才女,辭賦家司馬相如之妻卓文君曾提及九連環:……,七弦琴無心彈,八行書無可傳,九連環從中折斷,十里長亭望眼欲穿;百思想,千懷念,萬般無奈把郎怨。……
卓文君生於西漢,諸葛亮生於東漢末年,其時漢室江山已分崩離析。二人相差幾百年。也就是說,在諸葛亮之前幾百年的西漢,九連環已經存在。故「九連環由諸葛亮發明」之說並不正確,可能系後世誤傳。也有人認為卓文君作詞的故事似元朝杜撰,因為詞風明顯不是漢朝時所有。
2003年3月8日,中國甘肅省嘉峪關市的王仲斌以3分57秒成功解出九連環,進入吉尼斯世界紀錄大全。
2012年10月25日CCTV新聞頻道報道,江西理工大學學生楊咸陽創造最快拆解九連環的記錄,時間為161秒(蒙眼)[1]
歷史
傳說九連環源於中國古代漢族民間,一說發明於戰國時代,另一說發明於三國時期,但能確認就是九連環的記載是明代楊慎(1488-1559,號升庵)的《丹鉛總錄》(見《升庵集》卷六十八),並不早於歐洲。
在中國,戰國時代名家惠施曾著立《連環可解》的立論。惠施所說連環是指《戰國策》卷第十三中提到的玉連環,南宋鮑彪注稱這種玉連環是「兩環相貫」,顯然不是這里所說的九連環。據說三國時期,諸葛亮常帶兵打仗,為排遣妻子寂寞而發明。於明代普及,明代中期時,流傳更是極廣。清代上至士大夫,下至販夫走卒,個個愛玩「九連環」。《紅樓夢》中曾有描寫在深閨中玩九連環的細節。
在西方,16世紀前,歐洲有了九連環的記載。1550年,巴黎刊行的數學文獻,清楚地討論過這「中國難題」。著名義大利數學家卡當的著作中將之稱為「中國九連環」。1685年,英國數學家瓦里斯對此作了詳細的數學說明。19世紀,格羅斯用二進位數給了它一個十分優美的解答。
九連環無論在任何時候,都有這聰明的象徵。在古代,對於人們來說,九連環不算是一種玩具,而是代表智慧的象徵。電視劇看多的人應該會有這樣一個印象:出使天朝的外邦,有些比較囂張的,都會拿出九連環來刁難朝中大臣。在大家都束手無策,外邦使官洋洋得意的時候,總會有一個比較聰明的人出來解出九連環,挽回天朝的顏面。因此,九連環總是會被賦予聰明,有智慧的帽子。
構造
九
各式各樣的九連環 (14張)
連環流行極廣,形式多樣,規格不一。其製作,用金屬絲製成圓形小環九枚,九環相連,套在條形橫板或各式框架上,其框柄有劍形、如意形、蝴蝶形、梅花形等,各環均以銅桿與之相接。玩時,依法使九環全部聯貫子銅圈上,或經過穿套全部解下。其解法多樣,可分可合,變化多端。得法者需經過81次上下才能將相連的九個環套入一柱,再用256次才能將九個環全部解下。此外,也可套成花籃、綉球、宮燈等狀。
同時,九連環也是按照一種順序來解的。解九連環需要相當一段時間,這也可以訓練人的耐心。不僅 如此,九連環還可以根據需要自行增加環數提高難度,但環數增加將使解開步驟呈幾何級數遞增,且本質上並沒有改變解環方法,因此通常所見仍是九環為主。
拆解原理
編輯
鎏金九連環
解開九連環共需要341步,只要上或下一個環,就算一步,不是在框架上滑動。希望大家能夠通過獨立思考,解決這個問題。九連環的解下和套上是一對逆過程。解法跟計算機的格雷碼是同一原理。
九連環的每個環互相制約,只有第一環能夠自由上下。要想下/上第n個環,就必須滿足兩個條件(第一個環除外)。一、第n-1個環在架上;二、第n-1個環前面的環全部不在架上。玩九連環就是要努力滿足上面的兩個條件。解下九連環本質上要從後面的環開始下,而先下前面的環,是為了下後面的環,前面的環還要裝上,不算是真正地取下來。
我們先從最簡單的一連環開始。解一連環需要1步:一下。解二連環需要2步:二下,一下。那解三連環呢?需要5步:一下,三下,一上,二下,一下。也就是解一個連環,再把最後一個環解下,再上一個一環,再解一個二連環。那解一個四連環,需要10步:二下,一下,四下,一上,二上,一下,三下,一上,二下,一下。也就是解一個二連環,再解最後一個環,再上一個二連環,再解一個三連環。
也就是說,解N連環,就是先解一個N-2連環,再解最後一個環,再上N-2連環,再解N-1連環。
解一連環需要1步,解二連環需要2步,由此可知,解三連環需要5步,解四連環需要10步,解五連環需要21步,解六連環需要42步,解七連環需要85步,解八需要170步,解九連環需要341步,解十連環需要682步……以後的類推。
破解方法
編輯
基本方法
把框架和九個圓環分開,如左手持框架柄,右手握環,從右到左編號為1-9將環套入框架為「上」,取出為「下」。
九連環拆解共256步:
下9:
下1(結果98765432在上):下1
下3(結果987654在上):下3上1下12
下5(結果9876在上):下5上12下1上3上1下12下4上12下1下3上1下12
下7(結果98在上):下7上12下1上3上1下12上4上12下1下3上1下12上5上12下1上3上1下12下4上12下1下3上1下12下6上12下1上3上1下12上4上12下1下3上1下12下5上12下1上3上1下12下4上12下1下3上1下12
下9(結果8在上):下9;
九連環的解法
下8:
上2(結果82在上):上12下1
上3(結果83在上):上3上1下12
上4(結果84在上):上4上12下1下3上1下12
上5(結果85在上):上5上12下1上3上1下12下4上12下1下3上1下12
上6(結果86在上):上6上12下1上3上1下12上4上 12下1下3上1下12下5上12下1上3上1下12下4上12下1下3上1下12
上7(結果87在上):上7上12下1上3上1下12上4上12下1下3上1下12上5上12下1上3上1下 12下4上12下1下3上1下12下6上12下1上3上1下12上4上12下1下3上1下12下5上12下1上3上1下12下4上12下1下3上1下12
下8(結果7在上):下8;
下7:
上2(結果72在上):上12下1
上3(結果73在上):上3上1下12
上4(結果74在上):上4上12下1下3上1下12
上5(結果75在上):上5上12下1上3上1下12下4上12下1下3上1下12
上6(結果76在上):上6上12下1上3上1下12上4上12下1下3上1下12下5上12下1上3上1下12下4上12下1下3上1下12
下7(結果6在上):下7;
下6:
上2(結果62在上):上12下1
上3(結果63在上):上3上1下12
上4(結果64在上):上4上12下1下3上1下12
上5(結果65在上):上5上12下1上3上1下12下4上12下1下3上1下12
下6(結果5在上):下6;
下5:
上2(結果52在上):上12下1
上3(結果53在上):上3上1下12
上4(結果54在上):上4上12下1下3上1下12
下5(結果4在上):下5;
下4:
上2(結果42在上):上12下1
上3(結果43在上):上3上1下12
下4(結果3在上):下4;
下3:
上2(結果32在上):上12下1
下3(結果2在上):下3;
下12:
下12(結果拆解完成):上1下12。
九連環安裝共341步:
上98:
上2(結果2在上):上12下1
上3(結果3在上):上3上1下12
上4(結果4在上):上4上12下1下3上1下12
上5(結果5在上):上5上12下1上3上1下12下4上12下1下3上1下12
上6(結果6在上):上6上12下1上3上1下12上4上12下1下3上1下12下5上12 下1上3上1下12下4上12下1下3上1下12
上7(結果7在上):上7上12下1上3上1下12上4上12下1下3上1下12上5上12下1上3上1下12下4上12下1下3上1下12下6上12下1上3上1下12上4上12下1下3上1下12下5上12下1上3上1下12下4上12下1下3上1下12
上8(結果8在上):上8上12下1上3上1下12上4上12下1下3上1下12上5上12下1上3上1下12下4上12下1下3上1下12上6上12下1上3上1下12上4上12下1下3上1下12下5上12下1上3上1下12 下4上12下1下3上1下12下7上12下1上3上1下12上4上12下1下3上1下12上5上12下1上3上1下12下4上12下1下3上1下12下6上12下1上3上1下12上4上12下1下3上1下12下5上12下1上3上1下12下4上12下1下3上1下12
上9(結果98在上):上9
上76:
九連環的解法
上2(結果982在上):上12下1
上3(結果983在上):上3上1下12
上4(結果984在上):上上4上12下1下3上1下12
上5(結果985在上):上5上12下1上3上1下12下4上12下1下3上1下12
上6(結果986在上):上6上12下1上3上1下12上4上12 下1下3上1下12下5上12下1上3上1下12下4上12下1下3上1下12
上7(結果9876在上):上7
上54:
上2(結果98762在上):上12下1
上3(結果98763在上):上3上1下12
上4(結果98764在上):上4上12下1下3上1下12
上5(結果987654在上):上5
上32:
上2(結果9876542在上):上12下1
上3(結果9876532在上):上3
上1:
上1(結果安裝完成):上1。
遞歸破解
用「遞歸」描述九連環的拆裝方法比較容易理解和記憶。所謂遞歸就是第n個步驟的解決辦法可以用已知的n-1步(或更早)的辦法來解決。對於九連環來說,就是拆下第n個環的方法可以用拆下第n-1個環的方法來描述。把拆下第n個環的問題轉化成為如何拆下第n-1個環的問題,也就是我們會才第n-1個環就會拆第n個環。以下是具體的拆裝方法描述:
n 拆下第1個環的方法:(D1)
1.把第1個環推出橫桿,從橫桿上面穿下去。
n 裝上第1個環的方法:(U1 )
1.把第1個環從橫桿下面穿上去,拉到外面後套進橫桿。
n 拆下第2個環的方法:(D2 )
1. 把第1個換裝上;(U1)
2. 把第2個和第1個環一起推出橫桿,把第二個環從橫桿上面穿下去;(卸下第2環)
3. 把第1個換再拆下。(D1)
n 裝上第2個環的方法:(U2)
1. 把第1個換裝上;(U1)
2. 把第2個環從橫桿下面穿上去,拉到前面後套進橫桿;(裝上第2環)
3. 把第1個換再拆下。(D1)
n 拆下第n環的方法:( Dn )
1. 把第n-1環裝上去;(Un-1)
2.把第n和n-1環一起推出橫桿,把第n環從橫桿上面穿下去;(卸下第n環)
3. 再把第n-1環卸下。(Dn-1)
n 裝上第n個環的方法:(Un)
1. 把第n-1環裝上去;(Un-1)
2. 把第n個環從橫桿下面穿上去,拉到前面後套進橫桿;(裝上第n環)
3. 把第n-1個換再拆下。(Dn-1)
n 為了加快速度,可以把第n+1個環和第n個環一起卸下去:(Dn.n+1)
1. 把第n+1和n環一起推出橫桿,把n+1環從橫桿上面穿下去;(卸下第n+1環)
2. 把第n-1環裝上去;(Un-1)
3、把第n和n-1環一起推出橫桿,把第n環從橫桿上面穿下去;(卸下第n環)
4. 再把第n-1環卸下。(Dn-1)
舉例說明:
在第1個和第2個環都卸下去的情況下,如何把第3、4環卸下(n=3,D3.4):
1. 把第4和3環一起推出橫桿,把4環從橫桿上面穿下去;(卸下第4環)
2. 把第2環裝上去;(按方法U2)
3. 把第3和2環一起推出橫桿,把第3環從橫桿上面穿下去;(卸下第3環)
4. 再把第2環卸下。(按照方法D2)
步驟計算
編輯
九連環是中國漢族民間玩具。規定環在桿上用1表示,環在下面用0表示。規定最左邊的環是可以任意上下的那一環,輸出數據中最右邊必須是1,也就是說,010100要寫成0101。
現今是X連環,由於「輸出數據中最右邊必須是1」,所以X可以理解為無限大,右邊多餘的0在輸出時都省略。初始化各環都是0,以下是前9步的情況:
1.1
2.11
3.01
4.011
5.111
6.101
7.001
8.0011
9.1011
問在X連環裝上過程中,第n步完成後,具體情況是怎麼樣的。
答案:將n轉化為二進制,求其格雷碼。將二進制的格雷碼逆序輸出,即得具體情況。
注意:這個演算法揭示了傳統的九連環與現代的格雷碼的重要關系!
程序實現
一種解法
方法一:首先你不斷地數1,2,1,2,1,2,1……
數1的時候上或下第1個環
數2的時候先看看從第一個環數起第一個在框上的環是第幾個,就上/下它的下一個環。如:第1個在框上,則只需上/下第2個環;
第1至第5個都在框下,第6個在框上就上/下第7個環。
一直,堅持數完341個數就解出來了。
方法二:解九連環,首先要知道1、2環怎麼解開,會拆1、2環(如何拆卸請見第一步),就能拆後面的環了。
拆卸方法
第一環:把第一環從左邊拿起,從上面放下。(有圖解,每環步驟均相同)(綠紅方塊幫助辨別方向,箭頭指引移動方向) 重要步驟!
第二環:和第一環步驟一樣(第一環不能先下)。
拆卸提示(順序)
卸下第一環
卸下第三環
卸下第二環
卸下第五環
用前面的方法裝上第三環,此時可卸下第四環。
卸下第三環和更左邊的環。
卸下第七環
裝上第五環,卸下第六環。
裝上第四環,卸下第五環。
......卸完為止。
卸下第九環,離成功更進一步。
想辦法卸下第八環。
想辦法卸下更前面的環。
大功告成(其實有一百多個步驟)
③ 如何用python編寫一個素數環
此文主要目的,是向大家展示如何才能用python語言,來部署STARK演算法。
STARKs(可擴容的透明知識論證)是創建一種證明的技術,這項證明中f(x)=y,其中f可能要花很長的時間來進行計算,但是這個證明可以被很快驗證。STARK是「雙重擴容」:對於一個需要t步驟的計算,這會花費大約O(t * log(t))步驟才能完成這個證明,這可能是最優的情況,而且這需要通過~O(log2(t))個步驟才能驗證,對於中等大小的T值,它比原始計算快得多。STARKs也擁有隱私保護的「零知識證明」的特性,雖然我們將這類使用案例應用到其中,從而完成可驗證的延遲功能,不需要這類性質,所以我們不用擔心。
首先,先請幾項說明:
這個代碼還沒有完全審核;在實際使用案例中的情況,還不能保證
這部分代碼是還沒有達到理想狀態(是用Python語言寫的)
STARKs 的「真實情況」 傾向於使用二進制欄位而不是素數域的特定應用程序效率的原因;但是,他們確實也表現出,這里寫出的代碼是合法並且可用的。
沒有一個真實的方法來使用STARK。它是一個非常寬泛的加密和數學架構,同時為不同的應用有不同的設置,以及連續的研究來減少證明者和驗證者的復雜性,同時提高可用性。
此文希望大家能夠知道,模運算和素數域是如何運行的,
並且和多項式概念,插值和估值進行結合。
現在,讓我們一起來了解吧!
MIMC
下面是STARK的功能展示:
def mimc(inp, steps, round_constants): start_time = time.time() for i in range(steps-1): inp = (inp**3 + round_constants[i % len(round_constants)]) % molus print("MIMC computed in %.4f sec" % (time.time() - start_time)) return inp
我們選擇MIMC作為案例,因為它(i)很容易理解,(ii)在真實世界使用的很多。函數功能見下圖:
注意:在很多關於MIMC的討論中,你可以典型地看出使用了XOR,而不是+;這是因為MIMC可以在二進制情況下使用,其中添加是XOR;這里我們會在素數領域進行。
在我們的案例中,常數相對而言會是比較小的列表(例如,64位),這會一直連續地進行周期循環(也就說,在k[64]之後)。MIMC自身可以獲得這個特性,因為MIMC可以向後進行計算(從相應的輸出獲得輸入),但是往後計算需要比向前計算多花費100倍的時間(並且沒有方向可以同步進行)。所以你可以將往後計算的功能想像成計算不能同步的工作量證明,並且往前方向計算的功能可以作為驗證的過程。
x -> x(2p-1)/3 是x -> x3 的反函數;根據費馬小定理,這是真實的,盡管這個定理沒有費馬大定理出名,但是依然對數學的貢獻很大。
我們嘗試使用STARK來進行更加有效的驗證,而不是讓驗證者必須在向前方向運行MIMC,在完成向後計算之後,證明者可以在向前方向進行STARK計算,並且驗證者可以很簡單地驗證STARK。我們希望計算STARK可以比MIMC向前和向後之間的運行速度差別要小,所以證明者的時間仍然是有初始的向後計算來主導的。而並不是STARK計算。STARK的認證會相對較快(在python語言演算法中,可以是0.05-0.3秒),不論初始的計算時間有多長。
所有的計算會在2256 – 351 * 232 + 1個模內完成;我們使用素數模,因為它是小於2256 最大的素數,其中乘法群包含了232 個子集(也就是說,有這樣一個數g,從而在完全232次循環之後,G素數環的連續冪模繞回到1),而且是按照6k+5的形式。首個特性是保證FFT和FRI演算法的有效版本,其次是保證MIMC實際上可以向後計算(請見上面提到的x -> x(2p-1)/3 使用方法)。
素域操作
我們通過建立方便的等級來進行素域的操作,同時也有多項式的操作。代碼如下,收首先是小數位數:
class PrimeField(): def __init__(self, molus): # Quick primality test assert pow(2, molus, molus) == 2 self.molus = molus def add(self, x, y): return (x+y) % self.molus def sub(self, x, y): return (x-y) % self.molus def mul(self, x, y): return (x*y) % self.molus
並且使用擴展歐幾里得演算法,來計算模塊逆轉(這和在素域中計算1/x相同):
# Molar inverse using the extended Euclidean algorithm def inv(self, a): if a == 0: return 0 lm, hm = 1, 0 low, high = a % self.molus, self.molus while low > 1: r = high//low nm, new = hm-lm*r, high-low*r lm, low, hm, high = nm, new, lm, low return lm % self.molus
上面的演算法是相對昂貴的;幸運地是,對於特定的案例,我們需要做很多的模逆計算,有一個數學方法可以讓我們來計算很多逆運算,被稱為蒙哥馬利批量求逆:
使用蒙哥馬利批量求逆來計算模逆,其輸入為紫色,輸出為綠色,乘法門為黑色,紅色方塊是唯一的模逆。
下面的代碼是演算法的體現,其中包含一些特別的邏輯。如果我們正在求逆的集合中包含零,那麼它會將這些零的逆設置為 0 並繼續前進。
def multi_inv(self, values): partials = [1] for i in range(len(values)): partials.append(self.mul(partials[-1], values[i] or 1)) inv = self.inv(partials[-1]) outputs = [0] * len(values) for i in range(len(values), 0, -1): outputs[i-1] = self.mul(partials[i-1], inv) if values[i-1] else 0 inv = self.mul(inv, values[i-1] or 1) return outputs
這部分演算法接下來會驗證稱為非常重要的東西,特別是當我們開始和不同階的多項式進行計算的時候。
現在我們來看看一些多項式計算。我們把多項式當做一個數據集,其中的i是第i階(例如,x3 + 2x + 1變成[1, 2, 0, 1])。下面就是在一個點進行多項式估算的方法:
# Evaluate a polynomial at a point def eval_poly_at(self, p, x): y = 0 power_of_x = 1 for i, p_coeff in enumerate(p): y += power_of_x * p_coeff power_of_x = (power_of_x * x) % self.molus return y % self.molus
困難和挑戰
f.eval_poly_at([4, 5, 6], 2)的輸出是多少?模是31嗎?
下面的解釋就是答案
.其實也有代碼是多項式加法,減法,乘法和除法;這是很長的加減乘除運算。有一個很重要的內容是拉格朗日插值,它將一組 x 和 y 坐標作為輸入,並返回通過所有這些點的最小多項式(你可以將其視為多項式求值的逆):
# Build a polynomial that returns 0 at all specified xs def zpoly(self, xs): root = [1] for x in xs: root.insert(0, 0) for j in range(len(root)-1): root[j] -= root[j+1] * x return [x % self.molus for x in root] def lagrange_interp(self, xs, ys): # Generate master numerator polynomial, eg. (x - x1) * (x - x2) * ... * (x - xn) root = self.zpoly(xs) # Generate per-value numerator polynomials, eg. for x=x2, # (x - x1) * (x - x3) * ... * (x - xn), by dividing the master # polynomial back by each x coordinate nums = [self.div_polys(root, [-x, 1]) for x in xs] # Generate denominators by evaluating numerator polys at each x denoms = [self.eval_poly_at(nums[i], xs[i]) for i in range(len(xs))] invdenoms = self.multi_inv(denoms) # Generate output polynomial, which is the sum of the per-value numerator # polynomials rescaled to have the right y values b = [0 for y in ys] for i in range(len(xs)): yslice = self.mul(ys[i], invdenoms[i]) for j in range(len(ys)): if nums[i][j] and ys[i]: b[j] += nums[i][j] * yslice return [x % self.molus for x in b]
相關數學知識請參見此文的M-N部分。需要注意,我們也會有特別的方法lagrange_interp_4和lagrange_interp_2來加速次數小於 2 的拉格朗日插值和次數小於 4 的多項式運算。
快速傅立葉變換
如果你仔細閱讀上面的演算法,你也許會發現拉格朗日插值和多點求值(即求在N個點處次數小於N的多項式的值)都需要耗費2次時間,例如對於1000個點求拉格朗日插值,需要幾百萬個步驟,而且100萬個點的拉格朗日插值需要萬億個步驟。這是不可接受的低效率,所以我們需要使用更加有效的演算法,快速傅立葉變換。
FFT只需要花費O(n * log(n))的時間(也就是說,1000個點的計算需要10,000步,100萬個點的計算需要2000步),雖然它的范圍更受限制;x坐標必須是單位根部的完全集合,必須滿足N = 2k 階。也就是說,如果有N個點,那麼x坐標必須某個P值的連續冪,1, p, p2, p3…,其中pN = 1。這個演算法能夠用來進行多點計算和插值計算,而且只需要調整一個小參數。
下面就是演算法詳情(這是個簡單的表達方式;更詳細內容可以參閱此處代碼)
def fft(vals, molus, root_of_unity): if len(vals) == 1: return vals L = fft(vals[::2], molus, pow(root_of_unity, 2, molus)) R = fft(vals[1::2], molus, pow(root_of_unity, 2, molus)) o = [0 for i in vals] for i, (x, y) in enumerate(zip(L, R)): y_times_root = y*pow(root_of_unity, i, molus) o[i] = (x+y_times_root) % molus o[i+len(L)] = (x-y_times_root) % molus return o def inv_fft(vals, molus, root_of_unity): f = PrimeField(molus) # Inverse FFT invlen = f.inv(len(vals)) return [(x*invlen) % molus for x in fft(vals, molus, f.inv(root_of_unity))]
你可以自己通過一些輸入來運行代碼,並且看看是否能得到想要的結果,當你使用eval_poly_at的時候,給出你期望得到的答案。例如:
>>> fft.fft([3,1,4,1,5,9,2,6], 337, 85, inv=True) [46, 169, 29, 149, 126, 262, 140, 93] >>> f = poly_utils.PrimeField(337) >>> [f.eval_poly_at([46, 169, 29, 149, 126, 262, 140, 93], f.exp(85, i)) for i in range(8)] [3, 1, 4, 1, 5, 9, 2, 6]
傅里葉變換會把[x[0] …. x[n-1]]作為輸入,並且它的目標是輸出x[0] + x[1] + … + x[n-1]作為首個元素,x[0] + x[1] * 2 + … + x[n-1] * w**(n-1)作為第二個元素,等等;快速傅里葉變換可以通過把數據分為兩半,來完成這個,在兩邊都進行FFT,然後將結果結合在一起。
上圖就是信息如何進行FFT運算的解釋。請注意FFT是如何進行兩次數據復制,並且進行粘合,直到你得到一個元素。
現在,我們把所有部分組合起來,看看整件事情是如何:def mk_mimc_proof(inp, steps, round_constants),它生成運行 MIMC 函數的執行結果的證明,其中給定的輸入為步驟數。首先,是一些 assert 函數:
# Calculate the set of x coordinates xs = get_power_cycle(root_of_unity, molus) column = [] for i in range(len(xs)//4): x_poly = f.lagrange_interp_4( [xs[i+len(xs)*j//4] for j in range(4)], [values[i+len(values)*j//4] for j in range(4)], ) column.append(f.eval_poly_at(x_poly, special_x))
擴展因子是我們將要拉伸的計算軌跡(執行 MIMC 函數的「中間值」的集合)。
m2 = merkelize(column) # Pseudo-randomly select y indices to sample # (m2[1] is the Merkle root of the column) ys = get_pseudorandom_indices(m2[1], len(column), 40) # Compute the Merkle branches for the values in the polynomial and the column branches = [] for y in ys: branches.append([mk_branch(m2, y)] + [mk_branch(m, y + (len(xs) // 4) * j) for j in range(4)])
我們需要步數乘以擴展因子最多為 2^32,因為當 k > 32 時,我們沒有 2^k 次的單位根。
computational_trace_polynomial = inv_fft(computational_trace, molus, subroot) p_evaluations = fft(computational_trace_polynomial, molus, root_of_unity)
我們首個計算會是得出計算軌跡;也就是說,所有的計算中間值,從輸入到輸出。
assert steps <= 2**32 // extension_factor assert is_a_power_of_2(steps) and is_a_power_of_2(len(round_constants)) assert len(round_constants) < steps
然後,我們會從將計算軌跡轉換為多項式,在單位根 g (其中,g^steps = 1)的連續冪的軌跡上「放下」連續值,然後我們對更大的集合——即單位根 g2 的連續冪,其中 g2^steps * 8 = 1(注意 g2^8 = g)的多項式求值。
# Generate the computational trace computational_trace = [inp] for i in range(steps-1): computational_trace.append((computational_trace[-1]**3 + round_constants[i % len(round_constants)]) % molus) output = computational_trace[-1]
黑色: g1 的冪。紫色: g2 的冪。橙色:1。你可以將連續的單位根看作一個按這種方式排列的圓圈。我們沿著 g1的冪「放置」計算軌跡,然後擴展它來計算在中間值處(即 g2 的冪)的相同多項式的值。
我們可以將MIMC的循環常數轉換為多項式。因為這些循環常數鏈是非常通常發生地(在我們的測試中,每64個步驟都會進行),最終證明他們形成了64階的多項式,而且外面可以很容易計算出它的表達式,以及擴展式:
skips2 = steps // len(round_constants) constants_mini_polynomial = fft(round_constants, molus, f.exp(subroot, skips2), inv=True) constants_polynomial = [0 if i % skips2 else constants_mini_polynomial[i//skips2] for i in range(steps)] constants_mini_extension = fft(constants_mini_polynomial, molus, f.exp(root_of_unity, skips2))
假設其中有8192個步驟,並且有64個循環常數。這是我們想要做的:我們正在進行FFT,從而計算循環常數來作為g1128 的功能。然後我們在之間加入很多零,來完成g1本身的功能。因為g1128 大約每64步進行循環,我們知道g1這個功能也會同樣。我們只計算這個擴展中的512個步驟,因為我們知道這個擴展會在每512步之後重復。現在,我們按照斐波那契案例中那樣,計算C(P(x)),除了這次是計算,需要注意,我們不在計算使用系數形式的多項式;而是根據高次單位根的連續冪來對多項式進行求值。
c_of_p需要滿足Q(x) = C(P(x), P(g1*x),K(x)) = P(g1*x) – P(x)**3 – K(x);目標是對於任何我們放入計算軌道的x(除了最後一步,因為在最後一步之後,就沒有步驟),計算軌跡中的下個數值就和之前的相等,再加上循環常量。與第1部分中的斐波那契示例不同,其中如果某個計算步驟是在k向量,下個就會是k+1向量,我們把低次單位根( g1 )的連續冪放下計算軌跡,所以如果某個計算步驟是在x = g1i ,下個步驟就會在g1i+1 = g1i * g1 = x * g1。因此,對於低階單位根( g1 )的每一個冪,我們希望最終會是P(x*g1) = P(x)**3 + K(x),或者P(x*g1) – P(x)**3 – K(x) = Q(x) = 0。因此,Q(x) 會在低次單位根 g 的所有連續冪上等於零(除了最後一個)。
# Create the composed polynomial such that # C(P(x), P(g1*x), K(x)) = P(g1*x) - P(x)**3 - K(x) c_of_p_evaluations = [(p_evaluations[(i+extension_factor)%precision] - f.exp(p_evaluations[i], 3) - constants_mini_extension[i % len(constants_mini_extension)]) % molus for i in range(precision)] print('Computed C(P, K) polynomial')
有個代數定理證明,如果Q(x)在所有這些x坐標,都等於零,那麼最小多項式的乘積就會在所有這些x坐標等於零:Z(x) = (x – x_1) * (x – x_2) * … * (x – x_n)。通過證明在任何單個的坐標,Q(x)是等於零,我們想要證明這個很難,因為驗證這樣的證明比運行原始計算需要耗費更長的時間,我們會使用一個間接的方式來證明Q(x)是Z(x)的乘積。並且我們會怎麼做呢?通過證明D(x) = Q(x) / Z(x),並且使用FRI來證明它其實是個多項式,而不是個分數。
我們選擇低次單位根和高次單位根的特定排列,因為事實證明,計算Z(x),而且除以Z(x)也十分簡單:Z 的表達式是兩項的一部分。
需要注意地是,直接計算Z的分子和分母,然後使用批量模逆的方法將除以Z轉換為乘法,隨後通過 Z(X) 的逆來逐點乘以 Q(x) 的值。需要注意,對於低次單位根的冪,除了最後一個,都可以得到Z(x) = 0,所以這個計算包含其逆計算就會中斷。這是非常不幸的,雖然我們會通過簡單地修改隨機檢查和FRI演算法來堵住這個漏洞,所以就算我們計算錯誤,也沒關系。
因為Z(x)可以簡潔地表達,我們也可以獲得另個好處:驗證者對於任何特別的x,可以快速計算Z(x),而且還不需要任何提前計算。對於證明者來說,我們可以接受證明者必須處理大小等於步數的多項式,但我們不想讓驗證者做同樣的事情,因為我們希望驗證過程足夠簡潔。
# Compute D(x) = Q(x) / Z(x) # Z(x) = (x^steps - 1) / (x - x_atlast_step) z_num_evaluations = [xs[(i * steps) % precision] - 1 for i in range(precision)] z_num_inv = f.multi_inv(z_num_evaluations) z_den_evaluations = [xs[i] - last_step_position for i in range(precision)] d_evaluations = [cp * zd * zni % molus for cp, zd, zni in zip(c_of_p_evaluations, z_den_evaluations, z_num_inv)] print('Computed D polynomial')
在幾個隨機點上,進行概念檢測D(x) * Z(x) = Q(x),從而可以驗證轉賬約束,每個計算步驟是之前步驟的有效結果。但是我們也想驗證邊界約束,其中計算的輸入和輸出就是證明者所說的那樣。只是要求證明者提供P(1), D(1), P(last_step)還有D(last_step)的數值,這些都是很脆弱的;沒有證明,那些數值都是在同個多項式。所以,我們使用類似的多項式除法技巧:
# Compute interpolant of ((1, input), (x_atlast_step, output)) interpolant = f.lagrange_interp_2([1, last_step_position], [inp, output]) i_evaluations = [f.eval_poly_at(interpolant, x) for x in xs] zeropoly2 = f.mul_polys([-1, 1], [-last_step_position, 1]) inv_z2_evaluations = f.multi_inv([f.eval_poly_at(quotient, x) for x in xs]) # B = (P - I) / Z2 b_evaluations = [((p - i) * invq) % molus for p, i, invq in zip(p_evaluations, i_evaluations, inv_z2_evaluations)] print('Computed B polynomial')
那麼,我們的論證如下。證明者想要證明P(1) == input和P(last_step) == output。如果我們將I(x)作為插值,那麼就是穿越(1, input)和(last_step, output)亮點的線,於是P(x) – I(x)就會在這亮點上等於零。因此,它會證明P(x) – I(x)是P(x) – I(x)的乘積,並且我們通過提高商數來證明這點。
紫色:計算軌跡多項式 (P) 。綠色:插值 (I)(注意插值是如何構造的,其在 x = 1 處等於輸入(應該是計算軌跡的第一步),在 x=g^(steps-1) 處等於輸出(應該是計算軌跡的最後一步)。紅色:P-I。黃色:在x = 1和 x=g^(steps-1)(即 Z2)處等於 0 的最小多項式。粉紅色:(P – I) / Z2。
現在,我們來看看將P,D和B的默克爾根部組合在一起。
現在,我們需要證明P,D和B其實都是多項式,並且是最大的正確階數。但是FRI證明是很大且昂貴的,而且我們不想有三個FRI證明,所以,我們計算 P,D 和 B 的偽隨機線性組合,並且基於它來進行FRI證明:
# Compute their Merkle roots mtree = merkelize([pval.to_bytes(32, 'big') + dval.to_bytes(32, 'big') + bval.to_bytes(32, 'big') for pval, dval, bval in zip(p_evaluations, d_evaluations, b_evaluations)]) print('Computed hash root')
除非所有這三個多項式有正確的低階,不然幾乎不可能有隨機選擇的線性組合,所以這很足夠。
我們想要證明D的階數小於2 * steps,而且P 和 B 的次數小於steps,所以我們其實使用了隨機的P, P * xsteps, B, Bsteps 和 D的隨機組合,並且可以看出這部分組合是小於2 * steps。
現在,我們來檢查下所有的多項式組合。我們先獲得很多隨機的索引,然後在這些索引上為默克爾樹枝提供多項式:
k1 = int.from_bytes(blake(mtree[1] + b'\x01'), 'big') k2 = int.from_bytes(blake(mtree[1] + b'\x02'), 'big') k3 = int.from_bytes(blake(mtree[1] + b'\x03'), 'big') k4 = int.from_bytes(blake(mtree[1] + b'\x04'), 'big') # Compute the linear combination. We don't even bother calculating it # in coefficient form; we just compute the evaluations root_of_unity_to_the_steps = f.exp(root_of_unity, steps) powers = [1] for i in range(1, precision): powers.append(powers[-1] * root_of_unity_to_the_steps % molus) l_evaluations = [(d_evaluations[i] + p_evaluations[i] * k1 + p_evaluations[i] * k2 * powers[i] + b_evaluations[i] * k3 + b_evaluations[i] * powers[i] * k4) % molus for i in range(precision)]
get_pseudorandom_indices函數會回復[0…precision-1]范圍中的隨機索引,而且exclude_multiples_of參數並不會給出特定參數倍數的值。這就保證了,我們不會沿著原始計算軌跡進行采樣,否則就會獲得錯誤的答案。
證明是由一組默克爾根、經過抽查的分支以及隨機線性組合的低次證明組成:
# Do some spot checks of the Merkle tree at pseudo-random coordinates, excluding # multiples of `extension_factor` branches = [] samples = spot_check_security_factor positions = get_pseudorandom_indices(l_mtree[1], precision, samples, exclude_multiples_of=extension_factor) for pos in positions: branches.append(mk_branch(mtree, pos)) branches.append(mk_branch(mtree, (pos + skips) % precision)) branches.append(mk_branch(l_mtree, pos)) print('Computed %d spot checks' % samples)
整個證明最長的部分是默克爾樹分支,還有FRI證明,這是有更多分支來組成的。這是驗證者的實質結果:
o = [mtree[1], l_mtree[1], branches, prove_low_degree(l_evaluations, root_of_unity, steps * 2, molus, exclude_multiples_of=extension_factor)]
在每個位置,證明者需要提供一個默克爾證明,從而讓驗證者能夠檢查這個默克爾證明,並且檢查C(P(x), P(g1*x), K(x)) = Z(x) * D(x)以及B(x) * Z2(x) + I(x) = P(x)(提醒:對於不在初始計算軌道上的x,Z(x)不會是零,所以C(P(x), P(g1*x), K(x)也不會是零)。驗證者也會檢查線性組合是正確的,然後調用。
for i, pos in enumerate(positions): x = f.exp(G2, pos) x_to_the_steps = f.exp(x, steps) mbranch1 = verify_branch(m_root, pos, branches[i*3]) mbranch2 = verify_branch(m_root, (pos+skips)%precision, branches[i*3+1]) l_of_x = verify_branch(l_root, pos, branches[i*3 + 2], output_as_int=True) p_of_x = int.from_bytes(mbranch1[:32], 'big') p_of_g1x = int.from_bytes(mbranch2[:32], 'big') d_of_x = int.from_bytes(mbranch1[32:64], 'big') b_of_x = int.from_bytes(mbranch1[64:], 'big') zvalue = f.div(f.exp(x, steps) - 1, x - last_step_position) k_of_x = f.eval_poly_at(constants_mini_polynomial, f.exp(x, skips2)) # Check transition constraints Q(x) = Z(x) * D(x) assert (p_of_g1x - p_of_x ** 3 - k_of_x - zvalue * d_of_x) % molus == 0 # Check boundary constraints B(x) * Z2(x) + I(x) = P(x) interpolant = f.lagrange_interp_2([1, last_step_position], [inp, output]) zeropoly2 = f.mul_polys([-1, 1], [-last_step_position, 1]) assert (p_of_x - b_of_x * f.eval_poly_at(zeropoly2, x) - f.eval_poly_at(interpolant, x)) % molus == 0 # Check correctness of the linear combination assert (l_of_x - d_of_x - k1 * p_of_x - k2 * p_of_x * x_to_the_steps - k3 * b_of_x - k4 * b_of_x * x_to_the_steps) % molus == 0
其實還沒有完成成功;證明對跨多項式檢查和 FRI 所需的抽查次數的可靠性分析是非常棘手的。但是這些就是所有代碼,至少你不用擔心進行瘋狂的優化。當我運行以上代碼的時候,我們會獲得STARK證明,會有300-400倍的證明成本例如,一個需要 0.2 秒的 MIMC 計算需要 60 秒來證明)。這就使得4核機器計算MIMC中的 STARK,實際上可以比後向計算 MIMC 更快。也就是說,在python語言,這會相對低效的實現,並且這也會證明運行時間比例會不同。同時,也值得指出,MIMC 的 STARK 證明成本非常低,因為MIMC幾乎是完美地可計算,它的數學形式很簡單。對於平均計算,會包含更少的清晰計算(例如,檢查一個數是大於還是小於另一個),其計算成本可能會更高,會有大約10000-50000倍。