bitbank techblog

zkSNARKs をたずねて三千里 第 1 回

はじめに

この連載は、テレビのチャンネルを手で回して変えていた世代の平凡な開発者が、雰囲気で zkSNARKs を作っていく過程を共有するものです。
そのため、正確でない記述が含まれている可能性があることを予めご了承ください。

kamikita

zkSNARKs とは

zkSNARKs は、証明者が なにか について知っているということを、検証者にその事実以外は一切漏らさずに、対話なしで検証者に対して証明する暗号技術です。ゼロ知識証明は、論理を積み上げていくことで証明をする一般的な数学の証明とは違い、確率的な証明をします。

例えば、たくさんある数の中からどの数を選んでも期待通りの結果になるという主張があって、その中から適当な数を選んでもらい、実際に期待通りの結果になることを見せて、それをその証明とするというような感じです。

なにか とは

この連載では、方程式 が成立する変数への値の割り当てを なにか とします。zkSNARKs 自体は、プログラムの実行トレースなど、後述する 算術回路 で表現できるものであればどんなものでも証明できます。

なにか のことを正式には witness と呼びます。

作る zkSNARKs

zkSNARKs は総称で、Groth16PLONK など、その傘の下に入るものがいろいろとありますが、この連載では、ピノキオ という初期型の zkSNARKs を作ります。

作る言語は Rust にしました。

何を証明するのか

この記事では、同じく zkSNARKs の解説記事である Quadratic Arithmetic Programs: from Zero to Hero で使われている

$$
x^3 + x + 5 = 35
$$

を証明の土台として、この 方程式 が成立する正しい witness を証明者が知っていることを証明します。

証明は以下の流れで進んでいきます。

  1. 文字列の方程式 から 方程式の木
  2. 方程式の木 から 算術回路
  3. 算術回路 から Gate
  4. Gate から R1CS
  5. R1CS から QAP
  6. ...

今回は 5 までを実装して、ゼロ知識ではない形で、提供された witness の正しさを検証できるところまでを説明していきます。

作り方についての方針

現場で使われているようなコードは、いろいろと最適化がされていたりして理解が難しい場合が多いですが、この連載では、できるだけ簡単なやり方で効率を考えずに実装をしていくことで、コードを読んで直接 zkSNARKs を理解できるようなものを作っていきたいと考えています。
また、極力既存のライブラリを使わず、自力でコードを書いていきます。

コードの置き場所は ここ です。

有限体Finite field

zkSNARKs では、0 から 素数 − 1 までの整数で構成されている 有限体 上で四則演算を行います。ここでいう 有限体 は位数が素数の 有限体 になります。

また、この上限を決めている素数のことを、位数Order)と呼びます。

例えば、位数 が $7$ の 有限体 は、

0, 1, 2, 3, 4, 5, 6

で構成されます。位数 が $7$ の 有限体 のことを $\mathbb{F}_7$ などと書くようです。他にも $GF(7)$ のような書き方を見かけます。

有限体 上の演算

有限体 上には、足し算と掛け算が定義されています。

足し算

有限体 上の足し算では、有限体 上の要素を足した結果が 位数 を超えた場合は、足した結果を 位数 で割った余りを足し算の結果とすることで、結果が 有限体 上に収まるようにします。

例えば、$\mathbb{F}_7$ で、$3 + 5$ は $8$ になりますが、 $8$ は 有限体 上に存在しないので、$8 / 7$ の余り 1 を結果にします。

このことを、$3 + 5 \equiv 1 \bmod 7$ とも書きます。$\equiv 1 \bmod 7$ で、位数 $7$ の数の世界で余り $1$ のグループに属しているという意味になります。

掛け算

掛け算は、足し算の繰り返しです。例えば、$3 \cdot 3$ は、$3 + 3 + 3$ と同じと考えてください。

引き算と割り算は?

四則演算をするには、加えて引き算と割り算が必要ですが、これらについては、足し算や掛け算と 逆元 を組み合わせることで実現できます。

まず、逆元 の説明に必要な 単位元 について説明します。

単位元

有限体 上の 2 つの演算にはそれぞれ 単位元 が定義されています。単位元 は、演算の相手として使うと何も起こらない 有限体 の要素です。足し算であれば $0$、掛け算であれば $1$ になります。

// 足し算
a + 0 = a  // 単位元は 0

// 掛け算
a * 1 = a  // 単位元は 1
逆元

有限体 上の足し算と掛け算それぞれには、掛け算の $0$ を除いた全ての要素に 逆元 が存在します。ある要素の 逆元 は、その要素と演算すると結果が 単位元 になる 有限体 上の要素です。

// 足し算
a + [a の逆元] = 0

// 掛け算
a * [a の逆元] = 1  // 0 の逆元は存在しない
引き算

ある要素の 逆元 を足すことが、ある要素を引くことになります。

$\mathbb{F}_7$ の場合、$3$ の 逆元有限体 の始まりの要素 $0$ から、$6,5,4$ と 3 つ戻った $4$ となります。確かめてみると $3 + 4 \equiv 0 \bmod 7$ で逆元になっています。

割り算

引き算と同じように考えて、逆元 を掛けることを割り算とします。

例えば、$6 / 2 = 3$ の場合、$6 \cdot 4 = 24 \equiv 3 \bmod 7$ なので、$2$ の 逆元 は $4$ となり、$6 \cdot 4$ が $6 / 2$ に対応する掛け算になります。

有限体 のように 位数 が素数の場合、$0$ 以外の各要素について掛け算の 逆元 が必ず 1 つだけ存在するので、$0$ 以外の要素での割り算が常に成り立ちます。

実装

ここまでに説明してきたことを実装します。

多倍長整数ライブラリ

zkSNARKs は、ペアリングに適した楕円曲線(Pairing-friendly curves)を使用するため、そのような楕円曲線を使って定義される群の 位数 を実装で扱うことができる必要がありますが、セキュリティー上の問題がないレベルの 位数 は $2^{381}$ などとかなり大きくなります。例えば ここ に位数の例が出ています。

Rust の最大の符号なしの数値型は u128 で $2^{128}-1$ までの整数しか扱えないため、zkSNARKs にそのまま使うことはできません。そのため、別途 多倍長整数 ライブラリが必要になりますが、実用的な速さのものを自分で作るのは難しいので、この部分については外部ライブラリの num_bigint を使いました。

有限体 を表す型

有限体 を表す型は Field です。位数 をメンバーに持っています。BigUintnum_bigint ライブラリ に含まれている、符号なし任意精度整数を表す型です。

pub struct Field {
  pub order: Rc<BigUint>,
}

GitHub

有限体 の要素を表す型

有限体 の要素を表す型は FieldElem です。

pub struct FieldElem {
  pub f: Field,
  pub n: BigUint,
}

GitHub

有限体 上の四則演算

四則演算の実装では、使い勝手を上げるために、演算の相手になる rhs (right-hand side の略) を BigUint への変換が定義されている型 ToBigUint としています。

また、変換元の型が BigUint の場合には 位数 以上の大きさの数が入力となる場合があるので、全ての演算の処理の先頭で $\bmod$ をすることで rhs が必ず 有限体 の要素になるようにしています。

足し算

BigUint 上で足し算を行い、結果が 位数 を超えた場合には 位数 を結果から引き、結果を 有限体 の中に収めています。

pub fn plus(&self, rhs: &impl ToBigUint) -> FieldElem {
  let rhs = rhs.to_biguint() % &*self.f.order;
  let mut n = self.n.clone();
  n += &rhs;
  if n >= *self.f.order {
    n -= &(*self.f.order);
  }
  FieldElem { f: self.f.clone(), n }
}

GitHub

掛け算

BigUint 上で掛け算をした後で、位数 で 結果を $\bmod$ しています。

pub fn times(&self, rhs: &impl ToBigUint) -> FieldElem {
  let rhs = rhs.to_biguint() % &*self.f.order;
  let mut n = self.n.clone();
  n *= &rhs.to_biguint();
  n %= &(*self.f.order);
  FieldElem { f: self.f.clone(), n }
}

GitHub

引き算

BigUint での引き算の結果が 有限体 の要素になる場合には、その要素を結果にしています。結果が負の値になる場合には、位数 から負の値を引いたものを結果にしています。

pub fn minus(&self, rhs: &impl ToBigUint) -> FieldElem {
  let rhs = rhs.to_biguint() % &*self.f.order;
  let f = self.f.clone();
  if self.n < rhs {
    let diff = &rhs - &self.n;
    let n = &(*self.f.order) - diff;
    FieldElem { f, n }
  } else {
    let mut n = self.n.clone();
    n -= &rhs;
    FieldElem { f, n }
  }
}

GitHub

割り算

割る数の 逆元 を求めて、逆元 を掛ける演算に変換しています。

pub fn safe_div(&self, rhs: &impl ToBigUint) -> Result<FieldElem, String> {
  let rhs = rhs.to_biguint() % &*self.f.order;
  let inv = self.f.elem(&rhs.to_biguint()).safe_inv()?;
  Ok(self.times(&inv))
}

GitHub

拡張ユークリッドの互除法を使った 逆元 の計算

逆元は、拡張ユークリッドの互除法というアルゴリズムを使って求めています。

このアルゴリズムについては、詳細を Senk さんのブログ記事 から学びました。以下は、記事の内容を自分の言葉で書いたものです。

拡張ユークリッドの互除法

$x \neq 0,\ y \neq 0,\ c=GCD(x,y)$ のときに、
$ax+by=c$ となる整数 $a,\ b$ が存在して、かつ $a,\ b$ を実際に計算することができる。

アルゴリズム使った逆元の求め方

位数 $p$ の 有限体 の、掛け算における $n$ の 逆元 $n^{-1}$ は、

$n^{-1} \cdot n \equiv 1 \bmod p$

と表すことができるので、$1$ を移項して、

$n^{-1} \cdot n - 1 \equiv 0 \bmod p$

$\bmod$ の定義から、上の式は左辺が $p$ で割り切れることを意味するので、ある $m$ が存在して、

$n^{-1} \cdot n - 1 = m \cdot p$

式を変形して、

$n^{-1} \cdot n - m \cdot p = 1$

この式を、$ax+by=c$ に、$a = n^{-1},\ x = n,\ b = ーm,\ y = p,\ c = 1$ を割り当てたものと考えて

n^-1   n  -m   p     1
  \    |   |   |    /
   a * x + b * y = c

アルゴリズムの使用条件が満たされているかを確認すると、

  • 条件 $c=GCD(x,y)$ は $1=GCD(n,p)$ で、$p$ は素数なので OK
  • 条件 $x \neq 0,\ y \neq 0$ は、$n \neq 0,\ p \neq 0$ で、$p$ は素数で、$n$ も掛け算の逆元を求めているため $0$ でないのが前提なので OK

と、満たしているので、この形で拡張ユークリッドの互除法を実行すると、計算結果の $a$ から $n$ の 逆元 を取得することができる。

実装

以下、拡張ユークリッドの互除法で 逆元 を求めるコードです。

// based on extended Euclidean algorithm
pub fn safe_inv(&self) -> Result<FieldElem, String> {
  if self.n == BigUint::zero() {
    return Err("Cannot find inverse of zero".to_string());
  }
  let order = self.f.order.to_bigint().unwrap();
  let v = self.n.to_bigint().unwrap();
  let zero = BigInt::zero();
  let one = BigInt::one();

  // x0*a + y0*b = a
  // x1*a + y1*b = b
  let mut r0 = v.clone();  // initially equals to a
  let mut r1 = order.clone();  // initially equals to b
  let mut x0 = one.clone();
  let mut y0 = zero.clone();
  let mut x1 = zero.clone();
  let mut y1 = one.clone();

  while r1 != zero {
    // a mod b
    // = a - q*b
    // = (x0*a + y0*b) - q*(x1*a + y1*b)
    // = x0*a - q*x1*a + y0*b - q*y1*b
    // = (x0 - x1*q)*a + (y0 - y1*q)*b
    // = r
    let q = &r0 / &r1;
    let r2 = &r0 % &r1;
    // this produces the same result as above r2 using mod
    //let r2 = x2 * order + y2 * v;
    let x2 = &x0 - &x1 * &q;
    let y2 = &y0 - &y1 * &q;

    // do next calculattion based on new and previous equations
    r0 = r1;
    r1 = r2;
    x0 = x1;
    y0 = y1;
    x1 = x2;
    y1 = y2;
  }

  // if the result is not a field element, convert it to a field element
  let mut new_v = x0;
  if new_v < zero.clone() {
    while new_v < zero.clone() {
      new_v += order.clone();
    }
  } else {
    if new_v >= order {
      new_v %= order;
    }
  }
  Ok(FieldElem { f: self.f.clone(), n: new_v.to_biguint().unwrap() })
}

GitHub

証明

方程式 を木で表現する

証明する対象の $x^3 + x + 5 = 35$ という 方程式 が、まず文字列で表現されているとします。この 文字列の方程式 を以下の 方程式の木 に変換します。AddMul に紐づく数字は ID です。

lhs: --------------+
                   v
            ___ Add (4) ___
           /               \
        Mul (2)            Add (3)
       /       \         /        \
    Mul (1)    Var(x)  Var(x)    Num(5)
   /     \
Var(x) Var(x)

rhs: 35

まず、文字列のままでは扱いにくいので、左辺が数式の木、右辺が 有限体 の要素になる形に変換します。以下が 方程式 を表す型です。

pub struct Equation {
  pub lhs: MathExpr,
  pub rhs: FieldElem,
}

GitHub

数式の木 MathExpr は以下のように定義されています。

pub enum MathExpr {
  Equation(Box<MathExpr>, Box<MathExpr>),
  Num(FieldElem),
  Var(String),
  Mul(SignalId, Box<MathExpr>, Box<MathExpr>),
  Add(SignalId, Box<MathExpr>, Box<MathExpr>),
  Div(SignalId, Box<MathExpr>, Box<MathExpr>),
  Sub(SignalId, Box<MathExpr>, Box<MathExpr>),
}

GitHub

左辺は木、右辺は 有限体 の要素 1 つといびつな感じですが、必要最低限の範囲だけをカバーするよう実装したたためこのような形となりました。

方程式 パーサーは文字列を Equation に変換しますが、扱えるのは、

  • 左辺は括弧含む四則演算
  • 値は正か負
  • 右辺は値が一つ

という条件を満たす 方程式 のみです。

方程式 パーサーの詳細についての説明は割愛します。もし興味があれば コード を読んでみてください。

実際に 方程式の木 を作ってみる

以下のコードで、今回の例の 方程式の木 を作ります。

方程式 パーサーは累乗をサポートしていないので、$x^3 + x + 5 = 35$ の $x^3$ の部分を $(x * x * x)$ に置き換えて 方程式 パーサーに入力します。

fn blog_post_1_example_1() {
  let f = &Field::new(&37u8);
  let expr = "(x * x * x) + x + 5 == 35";
  match Parser::parse(f, expr) {
    Ok(eq) => {
      println!("{} -> {:?}", expr, eq);
    },
    Err(_) => panic!(),
  }
}

GitHub

以下が出力です。

(x * x * x) + x + 5 == 35 -> Equation { lhs: Add(4, Mul(2, Var("x"), Mul(1, Var("x"), Var("x"))), Add(3, Var("x"), Num(FieldElem { f: Field { order: 37 }, n: 5 }))), rhs: FieldElem { f: Field { order: 37 }, n: 35 } }

この出力を読みやすくしたのが先ほどの木の図になります。

方程式の木 から Gate の集まりに

次に、方程式 の木 を 算術回路Arithmetic Circuit)を介して、Gate の集まりに変換します。

haccho
算術回路Gate

算術回路Gate を組み合わせて作る木で、Gate は 2 つの入力を 二項演算で処理して 1 つの結果を出力するものです。この記事では足し算と掛け算の Gate を使います。

Gate への入力は、変数、有限体 の要素、一時的な変数のいずれかです。

// 掛け算 Gate と足し算 Gate

   Out          Out
    |            |
   (*)          (+)
  /   \        /   \
In1   In2   In1    In2
// 算術回路の例
- x, y, z は変数
- tn は一時的な変数

        t3
        |
      _(*)_
     |     |
    t1     t2
    /       \
  (+)       (*)
 /   \     /   \
x     4   y     z
算術回路witness の例

算術回路 には、その回路が成り立つ変数と一時的な変数への値の割り当てである witness が存在します。
例えば、上の 算術回路 の例の witness の一例は

x = 2
t1 = 6
y = 3
z = 5
t2 = 15
t3= 90

        90
        |
      _(*)_
     |     |
     6     15
    /       \
  (+)       (*)
 /   \     /   \
2     4   3     5

となります。

方程式 の木から 算術回路

方程式 の木を 算術回路 に変換するには、MulAdd の出力を直接 Gate の入力とはせずに一時的な変数に保存します。こうすることで全ての AddMul を、有限体 の要素、変数または一時的な変数が入出力となる Gate として扱えるようになります。

木のルートは rhs の値ですが、out という変数を置き、そこに rhs の値が割り当てられるようにします。

                  out (=35)
                   |
            ___ Add (4) ___
           /               \
           t2              t3
           |                |
        Mul (2)            Add (3)
       /       \         /        \
      t1      Var(x)  Var(x)    Num(5)
      |
    Mul (1)
   /     \
Var(x) Var(x)
算術回路 から Gate の集まりに

次に、算術回路 を Gate の集まりに分解して、各 Gate を $a\ \ (op)\ \ b = c$ の形で表します。

例えば、

     t2
     |
  Mul (2)       ---> t1 * x = t2
 /       \
t1      Var(x)

という感じです。

割り算と引き算

方程式 には割り算と引き算が含まれています。これらについてはどう扱えば良いでしょうか?

まず割り算ですが、移項して掛け算に書き換えるのが一般的なようです。

例えば、こんな 算術回路 があったとします。

      e
      |
     (/)
    /   \
   c     d
   |
  (+)
 /   \
a     b

まず 算術回路Gate の集まりに分解して、$a\ \ (op)\ \ b = c$ の形で表します。

a + b = c
c / d = e

次に、割り算の両辺に $d$ を掛けて左辺と右辺を入れ替えることで、割り算を掛け算に変換します。

a + b = c
e * d = c   // == c / d = e

こうすると、算術回路 の木として考えた場合には前後のつながりがおかしくなるのですが、
算術回路 から Gate の集まりに分解した後で Gate を結合して 算術回路 に戻すことはないので、Gate レベルで式が成立していれば良いようです。

逆元 を使って、$c \cdot d^{-1} = c$ とすることもできると思いますが、割り算の 逆元 は求めるコストが高いため式自体を書き変えるのだと思います。

引き算についても、同じ様に足し算の式に変換することも 逆元 を足すことも可能ですが、一般的には足し算の 逆元 を求めるコストが低いためか 逆元 を使うようです。

このあたりは、この分野に詳しい Barry Whitehat さんに教えてもらいました。初歩的な質問に対して大変親切に答えていただき非常に助かりました。zk 界隈には人が少ないためか、初学者に対して親切にしてくれる人が多い気がします。

実装

まず、方程式 の木から 算術回路 への変換です。

算術回路 のノードの入出力に使われるのは、

  • 有限体 の要素
  • 変数
  • 一時的な変数

なので、これらの要素だけで回路が構成されるように変換するのが一番素直だと思うのですが、
今回の実装では、この段階で次の段階の R1CS での処理の準備も併せて行っています。

例えば Gate の型は、以下のようになっていますが、

pub struct Gate {
  pub a: Term,
  pub b: Term,
  pub c: Term,
}

GitHub

$a,\ b,\ c$ の型である Term には、有限体 の要素、変数、一時的な変数の他に、次に説明する R1CS で使う Sum, One, Out という選択肢が含まれています。

Sum は足し算、One は $1$、Out算術回路 の最終出力を表しています。

pub enum Term {
  Num(FieldElem),
  One,
  Out,
  Sum(Box<Term>, Box<Term>),  // Sum will not contain Out and Sum itself
  TmpVar(SignalId),
  Var(String),
}

GitHub

変換のコードは ここ にあります。

次に、演算の種類ごとに Gate への変換方法を説明します。

掛け算

直接 Gate 化しています。

MathExpr::Mul(signal_id, left, right) => {
  let a = Gate::traverse_lhs(f, left, gates);
  let b = Gate::traverse_lhs(f, right, gates);
  let c = Term::TmpVar(*signal_id);
  gates.push(Gate { a, b, c: c.clone() });
  c
}

GitHub

割り算

掛け算に変換しています。

MathExpr::Div(signal_id, left, right) => {
  let a = Gate::traverse_lhs(f, left, gates);
  let b = Gate::traverse_lhs(f, right, gates);
  let c = Term::TmpVar(*signal_id);
  // a / b = c
  // -> b * c = a
  gates.push(Gate { a: b, b: c.clone(), c: a.clone() });
  // send c to next as the original division does
  c
},

GitHub

足し算

足し算を括弧で括ったものにに $1$ を掛ける掛け算に変換しています。この形を取る理由については R1CS のところで説明します。

MathExpr::Add(signal_id, left, right) => {
  let a = Gate::traverse_lhs(f, left, gates);
  let b = Gate::traverse_lhs(f, right, gates);
  let c = Term::TmpVar(*signal_id);
  // a + b = c
  // -> (a + b) * 1 = c
  let sum = Term::Sum(Box::new(a), Box::new(b));
  gates.push(Gate { a: sum, b: Term::One, c: c.clone() });
  c
},

GitHub

引き算

移項で足し算化をして、足し算として処理をしています。

MathExpr::Sub(signal_id, left, right) => {
  let a = Gate::traverse_lhs(f, left, gates);
  let b = Gate::traverse_lhs(f, right, gates);
  let c = Term::TmpVar(*signal_id);
  // a - b = c
  // -> b + c = a
  // -> (b + c) * 1 = a
  let sum = Term::Sum(Box::new(b), Box::new(c.clone()));
  gates.push(Gate { a: sum, b: Term::One, c: a.clone() });
  c
},

GitHub

heiwaen
実際に Gate の集まりに変換してみる

このコードで、方程式の木 から Gate の集まりへの変換を実行できます。

fn blog_post_1_example_1() {
  let f = &Field::new(&37u8);
  let expr = "(x * x * x) + x + 5 == 35";
  let eq = Parser::parse(f, expr).unwrap();
  let gates = &Gate::build(f, &eq);
  println!("{:?}", gates);
}

GitHub

見やすく整形した実行結果は、

[
  "x" * "x" = t1,
  "x" * t1 = t2,
  ("x" + 5) * 1 = t3,
  (t2 + t3) * 1 = out
]

となります。$t_n$ は一時的な変数です。

Gate の集まりから R1CS

制約

Gate制約constraint)とも呼ばれます。ここからは Gate制約 と呼んでいきます。

R1CS 形式

ここまでに作られた 制約 の集まりを、以下の R1CS 形式に変換していきます。

$$
\mathbf{a} \circ \mathbf{w} * \mathbf{b} \circ \mathbf{w} - \mathbf{c} \circ \mathbf{w} = 0
$$

各記号は、

  • $\mathbf{w}$ は変数に 有限体 の要素が割り当てられた witness ベクター
  • $\mathbf{a},\ \mathbf{b},\ \mathbf{c}$ は $w$ の要素を選択するベクター
  • $\circ$ は ベクターの内積を作る演算

と考えてください。

ベクターの内積

ベクターは配列で、ベクターの内積は同じ長さの 2 つのベクターの各要素を掛け合わせてできるベクターのことです。内積の計算をコードで書いてみると、

function innerProduct(a, b) {
  let c = []
  for (let i=0; i<a.length; ++i) {
    c.push(a[i] * b[i])
  }
  return c
}

const v1 = [1, 2, 3]
const v2 = [2, 3, 4]

const ip = innerProduct(v1, v2)
console.assert(JSON.stringify(ip) === "[2,6,12]")

こんな感じです。

制約R1CS の項の対応

制約 $a \cdot b = c$ を移項して、$a \cdot b - c = 0$ としたものと、R1CS 形式の $ \mathbf{a} \circ \mathbf{w} * \mathbf{b} \circ \mathbf{w} - \mathbf{c} \circ \mathbf{w} = 0 $ は、

  • $ \mathbf{a} \circ \mathbf{w} $ が $a$
  • $ \mathbf{b} \circ \mathbf{w} $ が $b$
  • $ \mathbf{c} \circ \mathbf{w} $ が $c$

というように対応しています。

制約 の項をベクターの内積で置き換える

witness ベクターの作成

まず、先頭に $1$、その後に、Gate の集まりに現れる変数、一時的変数と out に割り当てられた 有限体 の要素を配置した witness ベクター(以降 $w$ と呼びます)を作ります。

今回の例の $w$ は以下のようになります。

     1 x  t1  t2 t3  out
w = [1, 3, 9, 27, 8, 35]
実際に witness ベクターをつくってみる

今回の例の、変数を置き換える前の状態の witness ベクターをこのコードで作ることができます。

fn blog_post_1_example_1() {
  let f = &Field::new(&37u8);
  let expr = "(x * x * x) + x + 5 == 35";
  let eq = Parser::parse(f, expr).unwrap();
  let gates = &Gate::build(f, &eq);
  let r1cs_tmpl = R1CSTmpl::from_gates(f, gates);

  println!("{:?}", r1cs_tmpl.witness);
}

GitHub

出力は、以下の通りです。

[1, "x", t1, t2, t3, out]
witness 要素を抽出するベクターの作成

次に、$w$ から必要な要素を抽出するベクターを作ります。

変数の抽出

今回の例の 制約 の集まりは以下の通りでした。

[
  "x" * "x" = t1,
  "x" * t1 = t2,
  ("x" + 5) * 1 = t3,
  (t2 + t3) * 1 = out
]

例えば、変数のみで構成されている 2 行目

"x" * t1 = t2,

の各項を抽出してみます。

$x$ は、$w$ の左から 2 つ目なので、

[0, 1, 0, 0, 0, 0]

というベクターを作って $w$ との内積を取ると、

$ [0, 1, 0, 0, 0, 0] \circ [1, x, t1, t2, t3, out] = 0 \cdot 1 + 1 \cdot x + 0 \cdot t1 + 0 \cdot t2 + 0 \cdot t3 + 0 \cdot out = x$

と、$x$ が抽出できます。$t1, t2$ についても同じやり方で抽出可能です。

有限体 の要素の抽出

有限体 の要素の抽出は少しやり方違います。

仮に

1 * 5 = t1

という制約があったとすると、$1$ は $w$ の 1 つ目の $1$ を選択するベクター $[1, 0, 0, 0, 0, 0]$ で抽出します。

$5$ については、$1$ を $5$ 倍するベクター $[5, 0, 0, 0, 0, 0]$ で抽出することができます。

足し算 Gate の扱い

3 行目の

("x" + 5) * 1 = t3,

のように括弧が含まれる場合は、括弧の部分の抽出に内積が要素ごとの足し算であることを利用します。

具体的には $1$ を $5$ 倍したものと $x$ を選択するベクター $[5, 1, 0, 0, 0, 0]$ との内積をとると、

$$
[5, 1, 0, 0, 0, 0] \circ w = 5 * 1 + x * 1 + t1 * 0 + t2 * 0 + t3 * 0 + out * 0 = 5 + x
$$

となり、期待通りの結果になります。

実装

実装は、R1CSTmplR1CS に分かれていて、どちらも Constraint という型を使っています。

Constraint

制約 を表す型です。$a,\ b,\ c$ のベクターのほとんどの要素は $0$ なので、無駄を省くために、SparseVec という $0$ の要素を持たない型にデーターを保存しています。

pub struct Constraint {
  pub a: SparseVec,
  pub b: SparseVec,
  pub c: SparseVec,
}

GitHub

R1CSTmpl

R1CSTmpl は、Gate の集まりから作られる、変数のままの状態の witness と全ての 制約 を持つ型です。
indices は、要素名から、要素の配置位置のインデックスを取得するためのディクショナリです。

pub struct R1CSTmpl<'a> {
  pub f: &'a Field,
  pub constraints: Vec<Constraint>,
  pub witness: Vec<Term>,
  pub indices: HashMap<Term, usize>,  // Term's index in witness vector
}

GitHub

R1CS

R1CS は、R1CSTmplwitness への 有限体 の要素の割り当てから作られる、有限体 の要素が割り当てられた witness ベクターと、全ての 制約 を持つ型です。

pub struct R1CS {
  pub constraints: Vec<Constraint>,
  pub witness: SparseVec,
}

GitHub

その他の詳細

その他の詳細については R1CSTmpl のコードR1CS のコード を読んでみてください。特に複雑な部分はありません。

takano
実際に 制約 の集まりから R1CS へ変換してみる

このコードで、制約 の集まりから R1CS への変換を実行できます。

fn blog_post_1_example_2() {
  let f = &Field::new(&37u8);
  let expr = "(x * x * x) + x + 5 == 35";
  let eq = Parser::parse(f, expr).unwrap();
  let gates = &Gate::build(f, &eq);
  let r1cs_tmpl = R1CSTmpl::from_gates(f, gates);

  println!("w = {:?}", r1cs_tmpl.witness);
  println!("{:?}", r1cs_tmpl.constraints);
}

GitHub

見やすく整形した実行結果は、以下の通りです。

w = [1, "x", t1, t2, t3, out]
[
  [0,1,0,0,0,0] . w * [0,1,0,0,0,0] . w - [0,0,1,0,0,0] . w = 0,
  [0,1,0,0,0,0] . w * [0,0,1,0,0,0] . w - [0,0,0,1,0,0] . w = 0,
  [5,1,0,0,0,0] . w * [1,0,0,0,0,0] . w - [0,0,0,0,1,0] . w = 0,
  [0,0,0,1,1,0] . w * [1,0,0,0,0,0] . w - [0,0,0,0,0,1] . w = 0
]

R1CS から QAP

次は、内積ベースの R1CS を多項式ベースの QAP に変換します。

R1CS の行列化

ここ によると、より正式には 制約 の一つ一つを R1CS と呼ぶのではなく、項ごとに 制約 のベクターをまとめて 3 つの行列にしたものの全体を R1CS と呼ぶようです。
行列は同じ長さのベクターを並べて 1 つにまとめたものです。

ここで、R1CS のベクターの集まりを行列に変換します。

    | 0 1 0 0 0 0 |
a = | 0 1 0 0 0 0 |
    | 5 1 0 0 0 0 |
    | 0 0 0 1 1 0 |

    | 0 1 0 0 0 0 |
b = | 0 0 1 0 0 0 |
    | 1 0 0 0 0 0 |
    | 1 0 0 0 0 0 |

    | 0 0 1 0 0 0 |
c = | 0 0 0 1 0 0 |
    | 0 0 0 0 1 0 |
    | 0 0 0 0 0 1 |

R1CS 行列の転置

行列の転置とは、行列の左下を固定して、横軸が縦軸になるように、左下に向けて右上が左下にくるように裏返すような操作です。

|0 0 1|  --転置--> |0 2|
|2 0 0|                       |0 0|
                      |1 0|

コードで書くと、

function transpose(m) {
  let n = []
  for (let row=0; row<m.length; ++row) {
    let cols = []
    for (let col=0; col<m[row].length; ++col) {
      cols.push(m[col][row])
    }
    n.push(cols)
  }
  return n
}

m = [[1,2], [3,4]]
m = transpose(m)
console.assert(JSON.stringify(m) === JSON.stringify([[1,3],[2,4]]))

となります。

これを、R1CS の行列に適用していきます。こういうイメージです。

     __ witness __
        | 0 0 1 0 0 0
制約 | 1 0 0 0 1 0
      ...

           |
          転置
           |
           V

         _制約_
        | 0 1
        | 0 0
witness | 1 0
        | 0 0
        | 1 0
        | 0 0
具体例

今回の例の $a$ の行列を転置してみます。

わかりやすくなるように、縦軸の見出しには制約番号を、横軸の見出しには対応する witness の要素を置きました。元の行列は

[a の行列]

           1  x  t1 t2 t3 out
制約 1 | 0  1  0  0  0  0 |
制約 2 | 0  1  0  0  0  0 |
制約 3 | 5  1  0  0  0  0 |
制約 4 | 0  0  0  1  1  0 |

で、この行列の転置は、

[a の転置行列]

      制 制 制  制
      約 約 約  約
      1  2  3  4
     ____________
 1  | 0  0  5  0 |
 x  | 1  1  1  0 |
t1  | 0  0  0  0 |
t2  | 0  0  0  1 |
t3  | 0  0  0  1 |
out | 0  0  0  0 |

となります。

各行について 多項式 をつくる

転置後の行列の各行は witness 要素の各 制約 での値を表しています。

この各行を、ラグランジュ補完で $x=n$ で評価すると 制約 n の値になる多項式に変換します。

ラグランジュ補完

ラグランジュ補完は、点の集まりから、すべての点を通る曲線の式を導出する方法です。

行列の行を点の集まりに

まず、行列の各行を以下のルールで点の集まりに変換します。

  • $x$ 座標は 制約 $n$ の $n$
  • $y$ 座標は、対象となる witness の要素の行の 制約 $n$ の値
  • $y=0$ となる点については除外する

例えば、今回の例の witness 要素 $x$ の行

     制 制 制  制
     約 約 約  約
     1  2  3  4
    ____________
   | ...        |
x  | 1  1  1  0 |

を変換すると、以下の 3 つの点の集まりになります。

$$
(1,\ 1), (2,\ 1), (3,\ 1)
$$

ラグランジュ補完で 多項式 をつくる
手順
  1. それぞれの点について、その点の $x$ 座標で評価すると $y$ 座標になり、その他の $x$ 座標で評価すると $0$ になる多項式を作る
  2. 全ての多項式を足し上げて 1 つの多項式にする
具体例

$x$ の行に対応する点の集まり $(1,\ 1), (2,\ 1), (3,\ 1)$ から多項式を作ってみます。

点 $(1,\ 1)$ に対応する多項式

$(1,\ 1)$ に対応する多項式は、$x=1$ で評価すると $1$ になり、$x=2,\ 3,\ 4$ で評価すると $0$ になる必要があります。

そのために $x=2,\ 3,\ 4$ で評価した場合に $0$ になるように、分子に $(x-2)(x-3)(x-4)$ を置きます。こうすると、例えば $x=2$ であれば $(x-2)=(2-2)=0$ となり、$x=2,\ 3,\ 4$ で式全体が $0$ になります。

また、$x=1$ の場合の分子 $(1-2)(1-3)(1-4)$ を分母に置き、$x=1$ で全体が $1$ なるようにします。

多項式全体はこうなります。

$$
\frac{(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)}
$$

点 $(2,\ 1), (3,\ 1)$ に対応する多項式

点 $(2,\ 1), (3,\ 1)$ についてもやり方は同じです。多項式はこうなります。

  • 点 $(2,\ 1)$
    $$
    \frac{(x-1)(x-3)(x-4)}{(2-1)(2-3)(2-4)}
    $$

  • 点 $(3,\ 1)$
    $$
    \frac{(x-1)(x-2)(x-4)}{(3-1)(3-2)(3-4)}
    $$

全ての点に対応する多項式

この 3 つの多項式を全て足し上げると、witness の $x$ 要素の多項式になります。
$$
\frac{(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)} +
\frac{(x-1)(x-3)(x-4)}{(2-1)(2-3)(2-4)} +
\frac{(x-1)(x-2)(x-4)}{(3-1)(3-2)(3-4)}
$$

$y$ 座標が $1$ でない点の扱い

今回の例の witness の $1$ の要素の、制約 $3$ のように、$y$ 座標が $1$ でない場合があります。

      制 制 制  制
      約 約 約  約
      1  2  3  4
     ____________
 1  | 0  0  5  0 |
     ...

この場合は、多項式全体に $y$ 座標の値を掛けると希望する多項式になります。以下のような感じです。

$$
\frac{(x-1)(x-2)(x-4)}{(3-1)(3-2)(3-4)} \cdot 5
$$

実装

このコード で、R1CS を転置後、$a,\ b,\ c$ の各 witness 要素に対応する行について build_polynomial_for_target_values を呼び、多項式への変換を行った上で式を展開して、次数ごとに項をまとめています。

多項式の変換に必要な、多項式の足し算、引き算、掛け算は、Polynomial 型に実装されています。実装方法は以下の通りです。

足し算と引き算

同じ次数の係数を足し引きしています。どちらか片方にしかない次数の場合は、その片方の項をそのまま結果に含めています。

掛け算

2 つ項の係数を掛け合わたものを計算結果の係数に、2 つの項の次数の合計を計算結果の次数にしています。どの項同士を掛け合わせるのかについては分配律を使って決めています。

この実装はとても効率が悪いのですが、フーリエ変換を使うことで高速化ができるようです。この 動画 が参考になりました。

演算における変数の扱い

以上の演算は変数に 有限体 の要素を割り当てずに行っているため、演算結果に変数が含まれます。演算後に変数への 有限体 要素の割り当てを決めることで、多項式を 有限体 の要素に評価できるようになっています。

実際に R1CS と $w$ から 多項式に変換してみる

以下のコードで、R1CS と $w$ から多項式への変換を実行できます。

fn blog_post_1_example_1() {
  let f = &Field::new(&37u8);
  let expr = "(x * x * x) + x + 5 == 35";
  let eq = Parser::parse(f, expr).unwrap();
  let gates = &Gate::build(f, &eq);
  let r1cs_tmpl = R1CSTmpl::from_gates(f, gates);

  // build witness
  /*
    x = 3
    t1 = x(3) * x(3) = 9
    t2 = t1(9) * x(3) = 27
    t3 = x(3) + 5 = 8
    out = t2(27) + t2(8) = 35
  */
  let witness = {
    use crate::snarks::term::Term::*;
    HashMap::<Term, FieldElem>::from([
      (Term::var("x"), f.elem(&3u8)),
      (TmpVar(1), f.elem(&9u8)),
      (TmpVar(2), f.elem(&27u8)),
      (TmpVar(3), f.elem(&8u8)),
      (Out, eq.rhs),
    ])
  };

  let r1cs = R1CS::from_tmpl(f, &r1cs_tmpl, &witness).unwrap();
  let qap = QAP::build(f, &r1cs, &ApplyWitness::Beginning);

  println!("a:\n{}", qap.a_polys.pretty_print());
  println!("b:\n{}", qap.b_polys.pretty_print());
  println!("c:\n{}", qap.c_polys.pretty_print());
}

GitHub

以下は、実行結果に witness 要素の見出しをつけて、整形して読みやすくしたものです。

a:
  1 [20,2,36,16]
  x [6,13,3,18]
 t1 [0,0,0,0]
 t2 [10,31,10,23]
 t3 [29,27,29,26]
out [0,0,0,0]

b:
  1 [3,1,21,12]
  x [12,24,23,18]
 t1 [20,30,1,23]
 t2 [0,0,0,0]
 t3 [0,0,0,0]
out [0,0,0,0]

c:
  1 [0,0,0,0]
  x [0,0,0,0]
 t1 [36,35,32,17]
 t2 [23,16,3,32]
 t3 [32,18,28,33]
out [2,21,2,12]

多項式が、左から右に向けて次数が大きくなる順序で係数が並べられたベクターとして出力されています。

例えば $a$ の witness 要素 $1$ のベクター [20,2,36,16] は、

$$
16x^3 + 36x^2 + 2x + 20
$$

という多項式を表しています。

misen

QAP での個別の 制約 成立チェック

多項式と $w$ から 制約 を復元する

多項式を特定の $x$ で評価してできるベクターと $w$ との内積を取ると、対応する 制約 の要素を復元することができます。

例えば $a$ で $x=1$ の場合のベクターと $w$ の内積を取ると、

         +--- 全ての witness 要素に対応する多項式を x=1 で評価して出来たベクター
  w      v
|  1|   |0|
|  x|   |1|
| t1| o |0| = 1*0 + t1*1 + x*0 + t2*0 + t3*0 + out*0 = t1
| t2|   |0|
| t3|   |0|
|out|   |0|

$t1$ が復元されます。

制約 の個別成立確認

次に 制約 全体を復元する式 $a(x) \cdot b(x) - c(x) = 0$ を定義します。

$w$ の要素は実際には変数ではなく、割り当てられた 有限体 の要素になります。

a(x) =
  1 * aの多項式 1 (x) +
  x * aの多項式 2 (x) +
 t1 * aの多項式 3 (x) +
 t2 * aの多項式 4 (x) +
 t3 * aの多項式 5 (x) +
out * aの多項式 6 (x)

b(x) =
  1 * bの多項式 1 (x) +
  x * bの多項式 2 (x) +
 t1 * bの多項式 3 (x) +
 t2 * bの多項式 4 (x) +
 t3 * bの多項式 5 (x) +
out * bの多項式 6 (x)

c(x) =
  1 * cの多項式 1 (x) +
  x * cの多項式 2 (x) +
 t1 * cの多項式 3 (x) +
 t2 * cの多項式 4 (x) +
 t3 * cの多項式 5 (x) +
out * cの多項式 6 (x)

a(x) * b(x) - c(x) = 0

この式を全ての 制約 に対応する $x$ の値について評価して、全てが成立すれば $w$ は正しいことになります。

制約 の一括成立確認

まず Minimum Polynomial について説明します。

Minimum Polynomial

Minimum Polynomial は、全ての 制約 に対応する $x$ の値 $n$ について $(x-n)$ という項を作り、それら全てを掛け合わせた多項式です。ここからはその多項式を $Z(x)$ と呼びます。

例えば、$x=1,\ 2,\ 3,\ 4$ の場合は、$Z(x)$ は $(x-1)(x-2)(x-3)(x-4)$ になります。

確認の方法

正しい $w$ から $a(x) * b(x) - c(x)$ が作られた場合には、全ての 制約 に対応する $x$ について評価結果が $0$ になりますが、このような性質を持つ多項式は、$Z(x)$ の倍数になるため、$w$ が正しい場合には、ある多項式 $H(x)$ が存在して下の式が成り立ちます。

$$
a(x) * b(x) - c(x) = Z(x) * H(x)
$$

これを利用して、$a(x) * b(x) - c(x)$ を $Z(x)$ で割り、余りがでるかを確認することで、全ての 制約 が正しいかどうかをまとめて確認できます。

$w$ は余りが出なければ正しく、余りが出るようであれば不正です。

確認の実装

多項式の割り算と $Z(x)$ を作るコードを実装して、$a(x) * b(x) - c(x)$ を $Z(x)$ で割り余りがでるかどうかを確認します。

割り算

割り算については、高校数学で学ぶらしい簡単なアルゴリズムを実装しました。
詳細については ここ などの説明を読んでみてください。コードは ここ にあります。

$Z(x)$ を作るコード

この部分でわかりにくいところはないと思います。コードは ここ です。

実際に 制約 の一括確認をしてみる

以下のコードで、witness が 正しいケースと正しくないケースで、制約 の一括確認を実行できます。

fn blog_post_1_example_2() {
  let f = &Field::new(&37u8);
  let expr = "(x * x * x) + x + 5 == 35";
  let eq = Parser::parse(f, expr).unwrap();
  let gates = &Gate::build(f, &eq);
  let r1cs_tmpl = R1CSTmpl::from_gates(f, gates);

  // build witness
  let good_witness = {
    use crate::snarks::term::Term::*;
    HashMap::<Term, FieldElem>::from([
      (Term::var("x"), f.elem(&3u8)),
      (TmpVar(1), f.elem(&9u8)),
      (TmpVar(2), f.elem(&27u8)),
      (TmpVar(3), f.elem(&8u8)),
      (Out, f.elem(&35u8)),
    ])
  };
  let bad_witness = {
    use crate::snarks::term::Term::*;
    HashMap::<Term, FieldElem>::from([
      (Term::var("x"), f.elem(&4u8)),  // replaced 3 with 4
      (TmpVar(1), f.elem(&9u8)),
      (TmpVar(2), f.elem(&27u8)),
      (TmpVar(3), f.elem(&8u8)),
      (Out, f.elem(&35u8)),
    ])
  };

  for test_case in vec![("good", good_witness), ("bad", bad_witness)] {
    let (name, witness) = test_case;
    let r1cs = R1CS::from_tmpl(f, &r1cs_tmpl, &witness).unwrap();

    let qap = QAP::build(f, &r1cs, &ApplyWitness::Beginning);
    let (a, b, c) = qap.get_flattened_polys();
    let t = a * &b - &c;

    let num_constraints = f.elem(&gates.len());
    let z = QAP::build_z(f, &num_constraints);

    let is_witness_valid = match t.divide_by(&z) {
      DivResult::Quotient(_) => true,
      DivResult::QuotientRemainder(_) => false,
    };
    println!("is {} witness valid? -> {}", name, is_witness_valid);
  }
}

GitHub

おわりに

まだ ピノキオ の全工程の途中でゼロ知識では全くない状態ですが、まだここまでしかできていないので、今回はここで終わりです。

この後は、Trusted Setup で作られた誰も知らない値を楕円曲線上の点に移したものを複数使い、計算内容を秘匿化しながら、ある点で $a(x) * b(x) - c(x)$ が Minimum Polynomial で割り切れるのかを確認するような流れになるようです。

この計算にはペアリングという手法を使った掛け算が必要です。ペアリングと聞いても何のことか想像がつかないと思いますが、自分もそうでした。実際のところまだよくわからないのですが、理解が進んだ時点でまたブログ記事を書こうと思います。ではその時まで。

platform
Author image
About exfinen
expand_less