STARKs Part III
This article is the Japanese translation of
[STARKs, Part III: Into the Weeds](https://vitalik.ca/general/2018/07/21/starks_part_3.html) originally written on 2018 Jul 21.
Eli Ben-Sasson さんに、やさしく補助していたき、いつも通り多大なる感謝をもうしあげます。ならびに、Chih-Cheng Liang さんと Justin Drake さんの査読にも感謝いたします。また Ben Fisch さんの reverse MIMC 技術の提案にも感謝いたします。
警告：数学 と python コードをたくさん含みます。
[前回](/ethereum/412-STARK2.md/)と[前々回](/ethereum/411-STARK1.md/)の投稿に続き、この記事は実際に STARK を実装するとどうなるかどうか、python を使って実装していきましょう。
STARKs とは、Scalabel Transparent Argument of Knowledge の略称で、f(x) = y を満たすことの証明をつくる技術です。
ここで、f を計算することは、潜在的に長い時間を要しやすいのですが、検証は素早く終えることがきます。
STARK は『二重にスケーラブル』です。
というのは、t ステップの計算の証明をひとつ作るのに、およそ $O(t \cdot log(t))$ ステップかかり、ほぼ最適化されています。
かつ、検証にはおよそ $~ O(log(log(t)))$ ステップを要し、t の値が適度に大きくなろうと、元の計算と比べるととても速いことがわかります。
STARKs には『ゼロ知識』というプライバシー保護の性質もあり、
使用例などはこの記事に書くつもりではあるが、
検証可能な遅延関数をつくることには、この性質は必要なく、
そのことについて心配する必要はありません。
最初にいくつか注意点を挙げておきます。
- このコードは細部まで監査されたものではなく、プロダクションでの使用におけるプロトコルの健全性は保証されていません。
- このコードの最適化は部分的なものです。（ python で書かれています。何を期待できましょうか。）
- 実際のSTARK（つまり Eli や co さんによる実装）は、効率性の観点から有限二進数体（簡単なものだと例えば、 $\mathbb{F}_2[X]/(x^3+x+1)$ は位数8の体となります。）上でプロトコルを展開していますが、彼らは、ここで用いる素数位数有限体によるアプローチをしても法則が破られないことを強調しています。$\mathbb{F}_2[X]/(f)$ が体となるには、f が既約多項式でないといけませんが、有限体上の既約性判定は、例えば、[Berlekamp Algorithm](https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1930-02.pdf) などが使えます。）
- STARK を実装する『真の方法』はありません。STARK は暗号学的かつ数学的構築のある幅広いカテゴリを指しており、異なるアプリケーション間では違う設定が採用されます。また証明者と検証者の複雑度を減らしたり健全性を証明するといった研究目的としても使われており、一つの実装はまだ存在していません。
- この記事は、合同算術と素数位数有限体の仕組みに関する前提知識と、ラグランジュ補間や評価に関する多項式についての基本知識を要します。これらについて、知らなければ[前回の記事]()や以前の[二次算術プログラム]()の記事を参照してください。
それでは、はじめましょう！
# MIMC
我々の STARK では次の関数をつかうことにしましょう。
~~~python
def mimc(inp, steps, k):
start_time = time.time()
for i in range(steps-1):
inp = (inp**3 + k[i % len(k)]) % modulus
print("MIMC computed in %.4f sec" % (time.time() - start_time))
return inp
~~~
この MiMC 関数（[論文](https://eprint.iacr.org/2016/492.pdf)）を選んだ理由は、シンプルで理解しやすく、かつ実際に使えるほどには興味深いものであることの二つです。この関数は視覚化すると次のようになります。
我々の例では、round定数 `k` は比較的小さな（例えば 64 つの）リストとし、何度も巡回して使用されます。
つまり `k[63]` の次は `k[0]` を使います。
ここでやるように、とても大きな数のラウンドを用いた MIMC は、検証可能な遅延関数として役に立ちます。
それは、計算するのが難しく、特に並列計算するのが難しいような関数であるが、検証は容易いような関数を指します。
MIMC は、出力から入力を復元するような、逆向き計算が可能ですが、それは前向きの計算のおよそ100倍の時間を要します。
そのような理由から、MIMC はある程度、検証可能な遅延関数としての性質を満たします。
逆向き計算を、並列化不能な PoW の採掘計算として、前向き計算を PoW の検証作業のようなものだとして捉えることが可能です。
より $x \mapsto x^3$ の逆数は $x \mapsto x^{\frac{2p−1}{3}}$ です
](/image/vb/stark3_2.png)
ここでの課題は、STARK を使った検証の効率化です。
検証者が、前向き計算のMIMC をする代わりに、
証明者が、逆向き計算の完了後に、前向き計算における STARK を計算します。
検証者は、このSTARK を検証するだけ、という形にしてみましょう。
嬉しいことに、STARK の計算は、逆向き計算と前向き計算の計算速度差よりも小さい計算量ですみます。
そのため、証明者は逆向きの計算に多くの時間を取れらますが、STARK の計算に費やす時間はほとんどありません。
元の計算がどれほど長くとも、STARK の検証は、比較的に速く終わります（我々の python 実装では、~0.05-0.3秒ほどです。）
全ての計算は、$2^{256} - 351 \cdot 2^{32} + 1$ で除した有限体で行われます。
この素数位数有限体を使うのは、$2^{256}$ 以下の有限体で、位数 $2^{32}$ の乗法部分群を含んでいる最大の有限体であり、
かつ6で除した余りが5である（つまりある k があって 6k+5 と表される）という理由からです。
一つ目は、我々のFFT と FRI アルゴリズムを効率よく使うための条件であり、
二つ目は、逆計算を可能にするための制約です。つまりは、$x \mapsto x^{\frac{2p−1}{3}}$ において、$\frac{2p−1}{3} = \frac{12k+10-1}{3} = 4k+3$
とすることが可能であるように選んでいます。
# Prime field operations
有限体上の計算に便利なクラスを一つ定義するとこから始めましょう。
もちろんその上の多項式計算にも使えます。コードは次の通りです。
~~~python
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
~~~
合同算術上の逆数を求める関数も定義します。
~~~py
# 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-lm*r, high-low*r
lm, low, hm, high = nm, new, lm, low
return lm % self.modulus
~~~
このアルゴリズムは比較的高価です
幸いなことに多くの逆数計算を一度に計算する数学的なトリックがあります、
それは [Montgomery batch inversion](https://books.google.com/books?id=kGu4lTznRdgC&pg=PA54&lpg=PA54&dq=montgomery+batch+inversion&source=bl&ots=tPJcPPOrCe&sig=Z3p_6YYwYloRU-f1K-nnv2D8lGw&hl=en&sa=X&ved=0ahUKEwjO8sumgJjcAhUDd6wKHWGNA9cQ6AEIRDAE#v=onepage&q=montgomery%20batch%20inversion&f=false) と呼ばれるものです。
以下のコードはこのアルゴリズムを実装したものです。
逆数を求める集合に零が含まれている場合も加味したため、少しだけ読みにくいコードとなっています。
~~~py
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
~~~
このアルゴリズムは後々、評価した多項式の割り算を行う時に重要となります。
さてある多項式演算を考えましょう。
我々は多項式を配列として、扱います。
つまり、 $x^3+2x+1$ は、`[1,2,0,1]` と次数の低い係数から表現することとします。
~~~py
# 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.modulus
return y % self.modulus
~~~
多項式の加減乗除の実装は、課題とします。
非自明な実装のひとつは、ラグランジュ補間です。
x座標とy座標の組の集合を入力として取り、それら全ての点を通る最小多項式を返します。
（これを多項式評価の逆プロセスと捉えることができます）
~~~py
# 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.modulus 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.modulus for x in b]
~~~
ここでの数学については、[こちらの記事](https://blog.ethereum.org/2014/08/16/secret-sharing-erasure-coding-guide-aspiring-dropbox-decentralizer/)の M of N の節をみてください。
我々は、かなり頻繁に生じる２次未満多項式と４次未満多項式のラグランジュ補間の速度をあげるために、
`lagrange_interp_4` と `lagrange_interp_2` という二つのメソッドも持っていることに注意してください。
# Fast Fourier Transforms
上のアルゴリズムを注意深く観察すると、
ラグランジュ補間と多項式多点評価（つまりは、N次未満多項式をN点で評価すること）の両者共に実行に、二乗時間かかることがわかります。
つまり、1000個の点に対するラグランジュ補間であれば、実行に数百万ステップ費やし、もし1000000個の点であれば、数兆ステップかかることとなります。
効率性の観点でこれは許容しがいため、より効率的なアルゴリズムつまりは、
高速フーリエ変換 FFT を使用してみましょう。
FFT は $O(n⋅log(n))$ しかかかりません（つまり、1000個の点に対し、高々10000ステップで、1000000個の点では、高々2000万ステップです）。
しかし、スコープ内では、さらに制限されます。
x 軸は、ある位数 $N=2^k$ に対して、[１の $N$ 乗根]となっていなければなりません。
つまりは、 $N$ 個の点があり、x軸は、ある $p^N = 1$ を満たす $p$ に対して乗数の列 $1,p,p^2,p^3,...$ となっていなければなりません。
このアルゴリズムは、驚くことに、パラメータを少し調整するだけで、多点評価や補間に利用することができます。
~~~py
課題
mod 337 における 1の16乗根であって、1の8乗根でないものを見つけよ。
答え
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 となる値 x を探して、
位数 337 の有限体上で、1 の原始根を一つ探します。
これは簡単に見つかり、一つの答えは 5 となります
そして、その (336/16) th powerを取ります
~~~
以下は FFT のアルゴリズムです（少し簡略化した形です；もう少し最適化したものはこちらのコードを参照してください。
~~~py
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` を使うと、期待通りの答えが得られます。例えば
~~~py
>>> 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の仕組みや理由、多項式計算一般について、より直感的に理解するために、[こちら](https://web.cecs.pdx.edu/~maier/cs584/Lectures/lect07b-11-MG.pdf)をお勧めします。
また、DFTとFFTの比較については、[このスレッド](https://dsp.stackexchange.com/questions/41558/what-are-some-of-the-differences-between-dft-and-fft-that-make-fft-so-fast?rq=1)で詳しく説明しています。
ただし、フーリエ変換に関するほとんどの文献は、実数および複素数上のフーリエ変換について述べているので、注意が必要です。
素数有限体上の話ではありません。
もし、これが難しすぎて理解したくないと思うのであれば、奇妙な呪文か何かのように扱ってください。
何度かコードを実行して動作することを確認したならば、視界がクリアになることでしょう。
# Thank Goodness It's FRI-day
(that's "Fast Reed-Solomon Interactive Oracle Proofs of Proximity")
注意：今こそ、[第2部](/ethereum/412-STARK2.md/)を読み直す良い機会かもしれません。
さて、ここからは低次の証明を作るための[コード](https://github.com/ethereum/research/blob/master/mimc_stark/fri.py)に入ります。
低次証明とは、与えられた値の集合のうち、少なくともある高い割合（例えば80％）が、
与えられた値の数よりもはるかに低い次数を持つ特定の多項式の評価を表していることを（確率的に）証明することである。
直感的には、「多項式を表すと主張するある Merkle ルートが、少しの誤差はあるが、実際に多項式を表している」という証明だと思えばよい。
入力として、次の4つを取ります。
- 低次多項式を評価した値の集合
- 累乗根。多項式が評価されるx座標はこの根の累乗により連続する列となる
- 多項式次数が厳密に N 未満である証明をするための値 N
- 合同算術の除数
我々のアプローチは再帰的なものであり、2つのケースがある。
まず、次数が十分に低い場合は、値のリスト全体を証明として提供するだけです。
これは「基底の場合」です。
基底の場合の検証は簡単で、FFTあるいはラグランジュ補間などでそれらの値を表す多項式を補間し、その次数が N 未満であることを検証すればよいのです。
そうでない場合は、次数がある設定された最小値よりも高ければ、[第2部](/ethereum/412-STARK2.md/)の下部で説明した垂直・対角線のトリックを行います。
まず始めに、`values` （つまりは）を Merkle 木に入れ、 $x$ 座標としては Merkle ルートを使って得られる疑似乱数 `special_x` を選択します。その後「列」を計算します。
~~~py
# Calculate the set of x coordinates
xs = get_power_cycle(root_of_unity, modulus)
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))
~~~
これは、数行のコードに多くのことを詰め込んでいます。
大まかなアイデアは、多項式を再解釈することです。
$P(x)$ は多項式 $Q(x,y)$ であり、 $P(x)=Q(x,x^4)$ となります。
$P$ が次数 $\lt N$ ならば $P'(y)=Q(specal\_x,y)$ は次数 $\lt N^4$ となります。
Qを実際に係数形式で計算するのは面倒なので（まだ比較的厄介で高価なFFTが必要！）代わりに別のトリックを使用します。
$Q(x,x^4)$ の引数となる任意の与えられた $x^4$ に対して、 $x$ の値を4つ与えてやることができます。
`x`、`modulus - x` と、今の合同算術体系内に二つ存在する $\sqrt{-1}$ をかけた `sqrt(-1) * x` です。ですので、今既に 4 つの $Q(?,x^4)$ の値を持っていることになり、 $R(x)=Q(x,x^4)$ の補間を行うことができます。そして、 $R(special\_x) = Q(specal_x,x^4) = P'(x^4)$ を計算すれば良いのです。 $x^4$ の取りうる値は $\frac{N}{4}$ 個あり、またそれらすべてを計算するのは容易いことです。
我々の証明では、$x^4$ の値をマークルルートのハッシュをシードとして利用した疑似乱数による探索により、例えば 40 個選び出して検証します。各探索に対して、 $Q(?,x^4)$ の五つの値のマークルブランチを提供します。
~~~py
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)])
~~~
検証者の仕事は、これら5つの値が実際にある一つの４次未満多項式に乗っているかどうかを検証します。
そこから、 列に対して、FRI を適用し、それを繰り返し、列が実際に $\frac{N}{4}$ 次未満であることを検証します。FRI に対してはそれで全てです。
課題として、エラーがある多項式評価の低次証明を作ってみて、いくつのエラーで検証機を通せるか試してみるのもいいでしょう（ヒント：prove_low_degree 関数を修正する必要があります；デフォルトの証明機では、1つのエラーでも膨らんで検証が失敗します）。
# The STARK
リマインダー：今こそ、[第1部](/ethereum/411-STARK1.md/)を見直す良い機会かもしれません。
さて、ここからは、これらのピースをまとめる実際の肉付けに入ります。
`def mk_mimc_proof(inp, steps, round_constants)`(コードは[こちら](https://github.com/ethereum/research/blob/master/mimc_stark/mimc_stark.py))は、
ある数のステップを入力として、MIMC関数を実行した結果の証明を生成します。
まず、いくつかの前提条件を与えます。
~~~py
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
~~~
この extension factor というのは、我々がどこまで計算トレース（MIMC 関数の実行による中間的な値の集合）を伸ばしていくのかという量です。
われわれは、$k>32$ となる k に対して、1 の $2^k$ 乗根を持たないため、
extension factor を掛けた step 数は高々 $2^32$ であることを必要とします。
最初の計算は、計算トレースを生成することです。
つまりは、入力から、計算のすべての中間的値を、出力まで集めてきたものです。
~~~py
# 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)]) % modulus)
output = computational_trace[-1]
~~~
そして、その計算トレースを多項式に変換します。トレースにおける $g_1^{steps} = 1$ となるような累乗根 $g_1$ の後継冪数の上でのトレースにおいて、その後継となる数を順に並べていきます。
次に $g_2^{steps\cdot 8} = 1$ をみたす累乗根 $g_2$ の後継冪数からなるより大きな集合のもとで多項式を評価します。（つまりは $g_2^8 = g_1$ となります。）
~~~py
computational_trace_polynomial = inv_fft(computational_trace, modulus, subroot)
p_evaluations = fft(computational_trace_polynomial, modulus, root_of_unity)
~~~
MIMC のラウンド定数を多項式に変換することができます。
これらラウンド定数のループは頻繁に、我々のテストでは64ステップ毎に、回されるため、64次多項式を構成することがわかります。
このことにより、その式や、その式を用いた別の式を評価するのが非常に容易になります。
~~~py
skips2 = steps // len(round_constants)
constants_mini_polynomial =
fft(round_constants, modulus, 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, modulus, f.exp(root_of_unity, skips2))
~~~
実行ステップ数が8192、ラウンド定数が64であるとします。
我々は FFT を使って、ラウンド定数である $g_2 = (g_1)^{128}$ から生成される 64 点を多項式へ変換する計算をします。
そして、ラウンド定数の間を $0$ で埋めて、 $g_1$ で生成される8192点を表す多項式へ拡張する計算をします。
$(g_1)^{128}$ は、64 ステップでループするので、$g_1$ もまた同様にループすることがわかります。
この拡張は 512 ステップを計算するだけです。なぜらなば、512 ステップで同じことを繰り返すことがわかっているからです。
さて、では パート１でのフィボナッチ数の例でみた $C(P(x))$ を計算しようと思いますが、
今回は、 $C(P(x))$ の代わりに $C(P(x), P(g_1 \cdot x),K(x))$ を使います。
~~~py
# 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)])
% modulus for i in range(precision)]
print('Computed C(P, K) polynomial')
~~~
ここで注意したいのが、もはや我々は多項式を係数形式のデータとして扱っていないということです。
我々は、多項式を累乗根の冪の列の上の評価の列という形式のデータとして扱っています。
`c_of_p` は、 $Q(x) = C(P(x), P(g_1\cdot x), K(x)) = P(g_1 \cdot x) − P(x)^3 − K(x)$ となるように意図しています。
この計算のゴールは、我々が計算トレースとして輪状に配置した任意の $x$ に対して、（最後の次のステップはないので最後だけは除きますが）
トレースの次の値は、前の値の３乗 + ラウンド定数 に等しくなるということです。
パート１で見たフィボナッチ数の例では座標 k の次は k+１でしたが、それとは異なり、ここでは累乗根 $g_1$ の冪を継続数としてとることで計算トレースを並べています。
ですので、ある一つの計算ステップの x 座標は $x = (g_1)^i$ の形で表現され、次のステップの座標は、$(g_1)^{i+1} = (g_1)^i \cdot g_1 = x \cdot g_1$ となります。
従って、最後を除く全ての $(g_1)$ の冪に対して、 $P(x \cdot g_1) = P(x)^3+K(x)$ あるいは $P(x \cdot g_1)−P(x)^3−K(x) = Q(x) = 0$ が成立します。
このように $Q(x)$ は最後を除くすべての累乗根 $g_1$ の冪の後続数の列に対して、 $0$ と等しくなるでしょう。
「Q(x) がこれら全ての x 座標で零と等しいならば、Q(x) は Z(x) で割り切れる」という代数学の定理があります。
ここで、$Z(x) = (x−x_1) \cdot (x−x_2) \cdot ... \cdot (x−x_n)$ です。
我々が検査したいすべての点において $Q(x)$ が零と等しいことをみることはとても大変です。
（そのような証明検証は元の計算をするよりも計算時間が長くなります）
その代わりに、我々は Q(x) が Z(x) の倍数多項式であることを証明するのに、間接的な（つまり確率的な）アプローチを使います。
さて、どのような方法なのでしょうか。もちろん、 $D(x) = \frac{Q(x)}{Z(x)}$ を提供し、FRI を使って、それが実際に余りなしで割り切れる商であることを証明すればよいのです。
$Z(x) = \frac{x^{steps} - 1}{x - (g_1)^{steps-1}}$ として表現できるため、
$Z(x)$ で割る操作はとても簡単なものになります。
~~~py
# 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 % modulus for cp, zd, zni in zip(c_of_p_evaluations, z_den_evaluations, z_num_inv)]
print('Computed D polynomial')
~~~
ポイントとしては、多項式 $Z$ の分母分子を直接先に評価計算をしており、
`multi_inv` を使用し、まとめて逆数を効率よく計算しているというところです。
その後各点毎に $Q(x)$ を掛けて $D(x)$ を求めています。
元の計算トレースの一部分である、低次元拡張の部分に沿った部分においては、 $Z(x) = 0$ であることに注意してください。
つまり、逆数をとると、計算が壊れてしまいます。
これはあまり嬉しくない事項ですが、
単に我々が、これらの穴を乱数による検証サンプルとして選ばないように、FRI アルゴリズムを改良すればよいだけです。
ですので、決してこのことが問題となることはありません。
$Z(x)$ の表現がとてもコンパクトであるので、我々は更なる恩恵を受けます。
検証者はどのような $x$ に対しても前処理などもなく素早く $Z(x)$ を計算できます。
証明者がステップ数のサイズの多項式を扱うのは、まあそれは仕方のないことでしょうが、
しかし検証者に対して同じことを強いるのはよくありません。検証者の仕事は、本質をつきかつ簡素なものでなくてはなりません。
（高速に処理する必要があります）
ある少しの選ばれた点の上で、$D(x) \cdot Z(x) = Q(x)$ を確率的に検証することで、『遷移条件』を検証することができます。
遷移条件とは、つまり、各計算ステップが、直前の計算ステップの有効な結果のものであるという制約です。
しかし、それだけでなく、『境界条件』もまた満たす必要があります。
境界条件とは、つまり、入力と出力が証明者の主張するものと一致しているかということです。
証明者に対して、計算トレースの最初と最後の値 $P(1), D(1), P(g^{steps-1}), D(g^{steps-1})$ だけを要求するだけでは、
検証はあまりにも、脆いものとなってしまいます。その他の点で、データが多項式上に乗っているということを示す証拠がどこにもありません。
ですので、多項式の割り算のトリックと似た種類のテクニックを使い、これを解決します。
~~~py
# 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) % modulus for p, i, invq in zip(p_evaluations, i_evaluations, inv_z2_evaluations)]
print('Computed B polynomial')
~~~
ここでの議論はつぎの通りです。
証明者は $P(1) = input$ であることと、 $P(g^{-1}) = output$ であることを証明したいわけです。
$I(x)$ をこれら二点を補間することで得られる多項式とします。つまりは、$(1,input)$ と $(steps-1, output)$ の2点を通る直線です。
そうすると $P(x) - I(x)$ は、これら2点上では $0$ となることでしょう。
つまり、$P(x) - I(x)$ は $Z_2(x) := (x-1)(x-(steps-1))$ の倍数多項式となっているはずです。
つまり、その商 $\frac{P(x) - I(x)}{(x-1)(x-(steps-1))}$ を与えれば良いのです！
~~~
課題
Suppose you wanted to also prove that
the value in the computational trace after the 703rd computational step is equal to 8018284612598740.
How would you modify the above algorithm to do that?
Mouseover below for answer Set I ( x ) to be the interpolant of
(1, input) , (g^703, 8018284612598740), (last_step, output),
and make a proof by providing the quotient
B(x) = (P(x) − I(x)) / ((x − 1)⋅(x − g^703)⋅(x − last_step))
~~~
それでは $P$ , $D$ , $B$ をつなげあわせて、マークルルートを作成しましょう。
~~~py
# 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')
~~~
今、我々は $P$ , $D$ , $B$ が実際に正しい最大次数を持つ多項式であるかどうか証明する必要があります。
しかし、FRI の証明はとても大きくかつ高価であるため使用したくありません。
代わりに、疑似乱数による $P$ , $D$ , $B$ の線形結合を計算します。
（乱数のシードとして $P$ , $D$ , $B$ のマークルルートを使用します）
そしてこれを行った上で FRI 証明に移ります。
~~~py
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 % modulus)
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) % modulus
for i in range(precision)]
~~~
これらすべての多項式が正しい低次数を持っていない限りは、
ランダムに選択された線形結合が正しい値を返すということはほぼ不可能であり、
そのためこれで十分となります。
我々は、多項式 $D$ の次数が $2 \cdot steps$ 未満であることと、
多項式 $P$ と $D$ の次数が $steps$ 未満であることを証明したい状況であり、
実際、 $P \cdot x^{steps}$, $B$ , $B^{steps}$ と $D$ のランダムな線型結合 をつくり、
この結合の次数が $2 \cdot steps$ 未満であることを検査します。
さて、全ての多項式に対して、スポット検査を行います。
いくつかランダムにインデックスを生成し、これらのインデックスの上での多項式評価であるマークルブランチを提供します。
~~~py
# 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)
~~~
`get_pseudorandom_indices` 関数は `[0 ... precision-1]` の間の乱数を返し、
パラメータ `exclude_multiples_of` は、ある与えられたパラメータ(ここでは `extension_factor`)の倍数ではないという制約を乱数に与えます。
これにより、もとの計算トレース上からサンプルしていないということを確約することができます。
元の計算トレース上では、間違った答えを得る可能性が高くなります。
おより、~250-500 キロバイトからなる証明は、マークルルートの集合と、スポット検査のブランチと、ランダムな線型結合からなる低次数証明ひとつ、の3点で構成されます。
~~~py
o = [mtree[1],
l_mtree[1],
branches,
prove_low_degree(l_evaluations, root_of_unity, steps * 2, modulus, exclude_multiples_of=extension_factor)]
~~~
証明の大部分は、実践的には、マークルブランチと、FRI 証明であり、とてもたくさんのブランチから構成されています。
そして最終的に、検証者がおいしくいただくことができる『お肉』にあたる検証コードは次の通りです。
~~~py
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) % modulus == 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)) % modulus == 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) % modulus == 0
~~~
証明者が提供したマークル木の証明の位置 $x$ に対して、
検証者は、 マークル木の証明であることをまず検査し、
次に $C(P(x), P(g_1 ⋅ x), K(x)) = Z(x) ⋅ D(x)$ と
$B(x) ⋅ Z_2(x) + I(x) = P(x)$ を検査します。
（リマインダ : 元の計算トレース上にない $x$ に対しては、Z(x) は零ではなく、
$C(P(x), P(g_1 ⋅ x), K(x))$ が零となるような可能性はほぼありません。）
検証者は線型結合の正しさも検査し、
FRI 証明の検証をするために `verify_low_degree_proof(l_root, root_of_unity, fri_proof, steps * 2, modulus, exclude_multiples_of=extension_factor)` を呼びます。
これで全てです！
余談として、
FRI 証明の検査や クロス多項式のスポット検査がどれほど必要であるかについてを証明する健全性解析は、実に技巧的です。
少なくともあなたが、より頭のおかしい最適化を行わない限りにおいては、ここにコードとして全て示してあります。
このコードを私が走らせたとき、300-400x 倍の証明のオーバーヘッドがかかりました。
（例えば、計算するのに 0.2 秒かかる MIMC の計算は、その証明を計算するのに 60 秒かかります）
これが示すように、4-core のノートパソコンで前方方向への MIMC 計算の STARK 計算は、
後方方向への MIMC 計算と比べて実際に非常に早いということがわかります。
また、これらはどちらも python による相対的に見て不十分な実装であり、最適化された場合の比率は異なる可能性はあります。
MIMC 上の STARK 証明のオーバーヘッドは、特筆すべきほど小さいということも指摘しておくべきでしょう。
というのは、MIMC はほぼ完璧に『算術化可能』だからです。
その数学的な形式というのはとてもシンプルです。
平均的な計算でありば、より算術的に複雑な操作（例えば大小関係）をもち、
オーバーヘッドは非常に高くなる可能性があります。
10000-50000x くらいにはなるでしょうか。