特別感謝 Eli ben Sasson 一如既往地提供幫助；也特別感謝 Chih-Cheng Liang和Justin Drake 的審閱。

為本系列的第 1 部分和第 2 部分的後續內容，本文將介紹在實際中實現 STARK 的途徑與效果，並使用 python 語言進行實現。STARKs（“可擴充套件的透明知識引數”是一種用於構建關於 f(x)=y的證明的技術。其中， f 可能需要很長時間來計算，但該證明可以非常快速地得到驗證。STARK 是“雙重可擴充套件的”：對於帶有 t 步的計算，其需要大約 O(t * log(t)) 步來生成證明，這可能是最優的；並且其需要 ~O(log2(t)) 步來進行驗證。在這種方式中，哪怕 t 的值很大，只要在一定範圍內，那麼計算過程也比原始計算要快得多。STARK 同樣擁有可進行隱私保護“零知識”屬性。儘管我們將這一屬性應用到本用例中，但是建立可驗證的延遲函式並不需要這個屬性，所以我們不需要擔心。

首先，以下是關於本文的免責宣告：

此程式碼未經過全面稽核，生產用例的可靠性無法保證
這些程式碼遠非最理想的實現（它們都是用 Python 編寫的，你還想怎麼樣）
考慮到特定應用的效率原因，STARK “在現實生活中”（即在 Eli 和 co 的生產實現中實現）傾向於使用二進位制域而不是素域。然而，他們確實在他們的著作中強調，基於素域的方法對於本文描述的 STARK 也是合理並且可行的。
不存在實現 STARK 的“唯一正確道路”。它是一種廣泛的密碼學和數學結構，針對不同應用具有不同的最佳設定。此外，關於減少證明者和驗證者複雜度並提高可靠性的研究仍在繼續。
本文絕對希望你已經瞭解模運算和素域的工作原理，並熟悉多項式、插值和求值的概念。否則的話，請回到本系列的第2部分，以及之前關於二次方程式算術的程式設計的文章。
現在，我們進入正題。

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)]) % modulus
print(“MIMC computed in %.4f sec” % (time.time() - start_time))
return inp
我們之所以選擇 MIMC（參見論文）作為例子，是因為它既（i）易於理解，同時（ii）足夠有趣，並且在現實生活中實在的用處。該函式可被看作下圖的形式：

注意：在許多關於 MIMC 的討論中，你通常會看到人們使用的是 XOR 而不是 +。這是因為 MIMC 通常在二進位制域上完成，而二進位制域的加法就是 XOR。在這裡，我們主要圍繞素域進行。

在本例中，迴圈常量是一個相對較小的列表（例如只包含 64 項），列表中的資料不斷迴圈（也就是說，在 k[64] 之後迴圈回到 k[1]）。

正如我們在這裡所做的那樣，具有非常多輪次的 MIMC 作為可驗證的延遲函式是非常有用的——這是一種難以計算的函式，尤其是無法平行計算，但驗證過程相對容易。 MIMC 本身在某種程度上實現了“零知識”屬性，因為 MIMC 可以“向後”計算（從其相應的“輸出”中恢復“輸入”），但向後計算需要的計算時間比向前計算多 100 倍（並且兩種方向的計算時間都無法通過並行化來顯著加快）。因此，你可以將向後計算函式視為“計算”不可並行化工作量證明的行為，並將前向計算函式計算為“驗證”它的過程。

-我們可以由x -> x^(2p-1)/3 得出 x -> x^3 的倒數。根據費馬小定理，這是正確的。費馬小定理儘管“小”，但毫無疑問，它對數學的重要性要大於更著名的“費馬最後定理”。-

我們在這裡嘗試實現的是，通過使用 STARK 使驗證更有效——相對於驗證者必須自己在前向執行 MIMC，證明者在完成“後向”計算後，將計算“前向”計算的 STARK，並且驗證者只需簡單地驗證 STARK。我們希望計算 STARK 的開銷能夠小於前向執行 MIMC 的速相對於後向的速度差異，因此證明者的時間仍將由最初的“後向”計算而不是（高度可並行化的）STARK 計算主導。無論原始計算的耗時多長，STARK 的驗證都可以相對較快（在我們的 python 實現中，約為 0.05 到 0.3 秒）。

所有計算均以 2^256 – 351 * 2^32 + 1 為模。我們之所以使用這個素域模數是因為它是 2^256 以內最大的素數，它的乘法組包含一個 2^32 階亞組（也就是說，存在數字 g，使得 g 的連續冪模這個素數之後能夠在 2^32 個迴圈以後回到 1） ，其形式為 6k + 5。第一個屬性是必要的，它確保我們的FFT和FRI演算法的有效版本可以發揮作用。第二個屬性確保 MIMC 實際上可以“向後”計算（參見上述關於 x -> x^(2p-1)/3 的使用）。

素域運算

我們首先構建一個可進行素域運算以及在素域上進行多項式運算的方便的類。其程式碼在此。初始的細節如下：

class PrimeField():
def init(self, modulus):
# Quick primality test
assert pow(2, modulus, modulus) == 2
self.modulus = modulus

def add(self, x, y):
    return (x+y) % self.modulus

def sub(self, x, y):
    return (x-y) % self.modulus

def mul(self, x, y):
    return (x*y) % self.modulus

以及用於計算模逆的擴充套件歐幾里德演算法（相當於在素域中計算1 / x）：

#Modular inverse using the extended Euclidean algorithm
def inv(self, a):
if a == 0:
return 0
lm, hm = 1, 0
low, high = a % self.modulus, self.modulus
while low > 1:
r = high//low
nm, new = hm-lmr, high-lowr
lm, low, hm, high = nm, new, lm, low
return lm % self.modulus
上述演算法的開銷相對較大。所幸的是，在需要進行眾多模逆計算的特殊情況中，有一個簡單的數學技巧可以幫助我們計算多個逆，我們稱之為蒙哥馬利批量求逆：

-使用蒙哥馬利批量求逆來計算模逆，其輸入為紫色，輸出為綠色，乘法門為黑色，紅色方塊是 唯一的 模逆。-

下述程式碼實現了這個演算法，並附有一些略微醜陋的特殊情況邏輯。如果我們正在求逆的集合中包含零，那麼它會將這些零的逆設定為 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 次項（例如， x^3 + 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§:
y += power_of_x * p_coeff
power_of_x = (power_of_x * x) % self.modulus
return y % self.modulus
思考題：

如果模數為31，那麼f.eval_poly_at([4, 5, 6], 2) 的輸出是多少？

答案是：

6 * 2^2 + 5 * 2 + 4 = 38，38 mod 31 = 7。

還有對多項式進行加、減、乘、除的程式碼；教科書上一般冗長地稱之為加法/減法/乘法/除法。有一個很重要的內容是拉格朗日插值，它將一組 x 和 y 座標作為輸入，並返回通過所有這些點的最小多項式（你可以將其視為多項式求值的逆）：

#構建一個在所有指定 x 座標處返回 0 的多項式
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.modulus for x in root]

def lagrange_interp(self, xs, ys):
# 生成主分子多項式，例如(x – x1) * (x – x2) * … * (x – xn)
root = self.zpoly(xs)

#生成每個值對應的分子多項式，例如，當 x = x2 時，,
#通過用主分子多項式除以對應的 x 座標
# 得到 (x – x1) * (x – x3) * … * (x – xn)
nums = [self.div_polys(root, [-x, 1]) for x in xs]

# 通過求出在每個 x 處的分子多項式來生成分母
denoms = [self.eval_poly_at(nums[i], xs[i]) for i in range(len(xs))]
invdenoms = self.multi_inv(denoms)

# 生成輸出多項式，即每個值對應的分子的總和
# 多項式重新調整為具有正確的 y 值
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.modulus for x in b]

相關數學說明請參閱這篇文章關於“M-of-N”的部分。需要注意的是，我們還有特殊情況方法 lagrange_interp_4 和 lagrange_interp_2 來加速次數小於 2 的拉格朗日插值和次數小於 4 的多項式運算。

快速傅立葉變換

如果你仔細閱讀上述演算法，你可能會注意到拉格朗日插值和多點求值（即求在N個點處次數小於N的多項式的值）都需要執行耗費二次時間。舉個例子，一千個點的拉格朗日插值需要數百萬步才能執行，一百萬個點的拉格朗日插值則需要幾萬億步。這種超低效率的狀況是不可接受的。因此，我們將使用更有效的演算法，即快速傅立葉變換。

FFT 僅需要 O(n * log(n)) 時間（即 1000 個點需要約 10000 步，100 萬個點需要約 2000 萬步），但它的範圍更受限制：其 x 座標必須是滿足 N = 2^k 階的單位根的完整集合。也就是說，如果有 N 個點，則 x 座標必須是某個 p 的連續冪1，p，p^2，p3 …，其中，p^N = 1。該演算法只需要一個小引數調整就可以令人驚訝地用於多點求值或插值運算。

思考題：

找出模337為1的16次單位根，且該單位根的8次冪模 337 不為1。

答案是：

59，146，30，297，278，191，307，40

你可以通過進行諸如 [print(x) for x in range(337) if pow(x, 16, 337) == 1 and pow(x, 8, 337) != 1] 的操作來得到上述答案列表。當然，也有適用於更大模數的更智慧的方法：首先，通過查詢滿足 pow(x, 336 // 2, 337) != 1（這些答案很容易找到，其中一個答案是 5）的值 x 來識別單個模 337 為 1 的原始根（不是完美的正方形） ，然後取它的 (336 / 16) 次冪。

以下是演算法實現（該實現略微簡化，更優化內容請參閱此處的程式碼）：

def fft(vals, modulus, root_of_unity):
if len(vals) == 1:
return vals
L = fft(vals[::2], modulus, pow(root_of_unity, 2, modulus))
R = fft(vals[1::2], modulus, pow(root_of_unity, 2, modulus))
o = [0 for i in vals]
for i, (x, y) in enumerate(zip(L, R)):
y_times_root = y*pow(root_of_unity, i, modulus)
o[i] = (x+y_times_root) % modulus
o[i+len(L)] = (x-y_times_root) % modulus
return o

def inv_fft(vals, modulus, root_of_unity):
f = PrimeField(modulus)
# Inverse FFT
invlen = f.inv(len(vals))
return [(x*invlen) % modulus for x in
fft(vals, modulus, 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 如何基於資料的兩半內容進行兩次 FFT 複製，並進行“粘合”步驟，然後依此類推直到你得到一個元素。

原文連結: https://vitalik.ca/general/2018/07/21/starks_part_3.html
作者: Vitalik
翻譯: 喏貝爾
稿源：公眾號 Unitimes·獨角時代（ID：Uni-times）