有限体の複数の元の逆数を一度に求める高速な方法

初めに

有限体の要素の逆数（逆元）操作は重たい処理です。
ここでは複数（$n$ 個）の要素の逆数をまとめて計算することで高速化する方法を紹介します。

2個の場合

$a$, $b$ を有限体の要素とします。求めたいものは $1/a$, $1/b$ です。
ここで、$1/a=b/ab$, $1/b=a/ab$ なので、まず $ab$ の逆数 $1/ab$ を求めてから $b$ 倍, $a$ 倍すると $1/a$, $1/b$ が求まります。
演算コストは1回の乗算を $M$, 1回の逆数計算を $I$ とすると $3M+I$ です。
BLS12-381の曲線では $I$ が $M$ より100倍前後重たいので、2回逆数を求める $2I$ に比べて約半分のコストとなります。

3個の場合

$a$, $b$, $c$ を有限体の要素とします。2個の場合と同様に考えると $abc$ を求めてその逆数 $1/abc$ を計算し、$bc$, $ac$, $ab$ を掛けて $1/a$, $1/b$, $1/c$ を求めればよいです。演算コストは $8M+I$.

4個以上の場合

$n=2$ や $n=3$ が簡単だったので一般化も容易そうに思いますが、$n=4$ で少し戸惑います。$a$, $b$, $c$, $d$ に対して $1/abcd$ を計算すると、$bcd$, $acd$, $abc$ が必要になります。
計算量を減らすために共通部分を探すと、$bc$ を求めたら、$(bc)d$, $a(bc)$ は1回ずつの乗算、$acd$ は共通項が無いので2回の乗算です。このやり方では $n$ が増えると大変そうです。
$n$ が大きくなっても同じ形の処理ができるように方針を少し変えましょう。
$ab$, $abc$, $abcd$ を計算してから $1/abcd$ が求まったところで

• $abc \times 1/abcd=1/d$,
• $ab \times (d \times 1/abcd)=1/c$,
• $a \times (cd \times 1/abcd)=1/b$,
• $bcd \times 1/abcd = 1/a$

とします。これなら1個あたり3回の乗算ですみます。

一時バッファを用意する

より一般的に $n$ 個の有限体の要素 $x_0$, $x_1, \dots, x_{n-1}$ が与えられたとします。
$n$ 個のバッファ $T$ を用意し、$x_0$, $x_0 x_1$, $x_0 x_1 x_2, \dots, x_0 \cdot \dots \cdot x_{n-1}$ と先頭から順次掛けたものを保持します。

$i$	0	1	2	...	$n-2$	$n-1$
$T[i]$	$x_0$	$x_0 x_1$	$x_0 x_1 x_2$	...	$x_0 \cdot \dots \cdot x_{n-2}$	$x_0 \cdot \dots \cdot x_{n-1}$

このテーブル作成のコストは $(n-1)M$ です。
そして $T[n-1]$ の逆数 $t=1/(x_0 \cdot \dots \cdot x_{n-1})$ を計算し、一つ前の $T[n-2]=x_0 \cdot \dots \cdot x_{n-2}$ と掛けると $1/x_{n-1}$ が求まります。
$t$ に $x_{n-1}$ を掛けると $t ← t x_{n-1} = 1/(x_0 \cdot \dots \cdot x_{n-2})$. これに $T[n-3]=x_0 \cdot \dots \cdot x_{n-3}$ を掛けて $1/x_{n-2}$.
以下同様に $t$ を更新しつつ、後ろから前にテーブルの値を掛けていくと全ての逆数が求まります。
演算コストは大体 $3nM + I$ です。$n$ が十分大きいときは一つあたりの逆数のコストは $3M$ に近づきます。

0を含むときの注意点

$x_0, \dots, x_{n-1}$ の中に0が含まれていると $T[n-1]=0$ となり、全ての逆数の値が0になってしまいます。従ってこのアルゴリズムでは $x_i=0$ ならスキップしなければなりません。
また $x_i=1$ のときは逆数も $x_i=1$ なので、スキップすると少しだけ演算コストが減ります。

C++による実装例

mclで使われているC++でのコードを抜粋します。$n$ 個の要素を持つ配列型の const Tin& x を入力として、それらの逆数を Tout& y に出力します。引数の T *t はテーブル作成のためのワーク領域です。
T::mul(T& z, const T& x, const T& y) は $z=x \times y$ を求める関数です。
前半で x[i] が0でも1でもないときだけ掛け算しながらテーブル T に配置し、後半で t[pos-1] の逆数を求めて、後ろから掛け算しながら逆数を求めています。

template<class Tout, class Tin, class T>
size_t invVecWork(Tout& y, Tin& x, size_t n, T *t)
{
  size_t pos = 0;
  for (size_t i = 0; i < n; i++) {
    if (!(x[i].isZero() || x[i].isOne())) {
      if (pos == 0) {
        t[pos] = x[i];
      } else {
        T::mul(t[pos], t[pos - 1], x[i]);
      }
      pos++;
    }
  }
  const size_t retNum = pos;
  T inv;
  if (pos > 0) {
    T::inv(inv, t[pos - 1]);
    pos--;
  }
  for (size_t i = 0; i < n; i++) {
    const size_t idx = n - 1 - i;
    if (x[idx].isZero() || x[idx].isOne()) {
      y[idx] = x[idx];
    } else {
      if (pos > 0) {
        T::mul(y[idx], inv, t[pos - 1]);
        inv *= x[idx];
        pos--;
      } else {
        y[idx] = inv;
      }
    }
  }
  return retNum;
}

楕円曲線の複数の点のアフィン座標化

楕円曲線の点を射影座標やヤコビ座標で表現しているとき、処理の最後でアフィン座標に変換する場合があります。
アフィン座標に変換するときは、$Z$ 座標の逆数などが必要になるので、複数の楕円曲線の点をアフィン座標化するには上記除算アルゴリズムを使うと高速化されます。

浮動小数点数への適用の可否

今回紹介した手法を有限体ではなく浮動小数点数に適用するのは慎重にした方がよいでしょう。
なぜなら有限体では常に厳密に計算できるので1個ずつ逆数を求めるのと結果は完全に一致します。
しかし、浮動小数点数では $n$ 個の値を全て掛け算したときにオーバーフローやアンダーフローが発生する可能性があり、誤差も発生します。
また、SIMD命令には近似逆数命令を使った逆数計算手法が利用できることも多く、今回の方法を使わなくても高速に計算できる可能性が高いからです。