Zenn
🙄

【NTT(数論変換)入門(1)】DFT(離散フーリエ変換)編

2024/03/29に公開

格子暗号の(高速化)の理解のためにはNTTが避けて通れない。何番煎じかと思われるかもしれないが、フーリエ変換系は多くの人が躓くポイントかと思うので(私も十分理解していないが)、わかりやすくまとめられると良いなと思う。

はじめに(全体像)

まずNTTとは何かを述べる。簡単に言うとNTTはDFTの一般化である。DFTは複素(関)数を複素(関)数へ写す全単射写像であり、要は複素数の表現方法を変換するエンコード処理だとみなしてよいと思う。全単射なので逆写像IDFTが存在し、これはエンコードした複素数を元に戻すデコード処理と考えることができる。

なぜこのような変換が重要なのかというと、変換後のほうが複素数に対する演算を楽に行うことができるからだ。暗号的にはモンゴメリ乗算みたいなものである。工学的には、時間領域より周波数領域が扱いやすいとかという言い方になるが、数論的には値を指数表現することで乗算を加算に変換するような例がわかりやすいかと思う。

下の図の例で、我々は343×49343 \times 49を計算したいとする。しかしながらこの計算は暗算でやるにはちょっとめんどくさい。そこで値を計算しやすい空間に変換することを考える。どうやるかはひとまず置いておいて、343=73,49=72343=7^3, 49=7^2のように表現方法を変換(エンコード)できる。指数表現では乗算は指数の加算に置き換わるため、73×72=757^3 \times 7^2 = 7^5は暗算で楽々計算できるだろう。これまたどうやるかはおいておいて、757^5という表現を16,80716,807に変換(元の空間へデコード)できれば、暗算でやるには少々面倒だった計算がらくらく計算できてうれしい。
気持ち的にはこの例のエンコードに対応するのがDFT (NTT)であり、デコードに対応するのがIDFT (INTT)である。(これはあくまでDFTの気持ちを伝えるための例であり、DFTによって値が指数表現に変換されるわけではない。)
ではDFTとNTTは何が違うんだという話になるが、DFTは複素数含む少数点数の変換が行えるのに対し、NTTは整数を整数へ変換する。こう書くとDFTがNTTの上位互換っぽくなるがそうではない。小数点数というのは厳密には無限の長さの値であり、PCでは値を正確に表現できないため、計算を繰り返していくうち誤差が蓄積して計算誤りが発生する可能性がある。それに対してNTTは整数を扱うため、(ある範囲の値までは)誤差なく計算することができるし、多くの暗号応用上も整数を扱得たほうが便利なことが多い。とはいえDFTとNTTは扱う値が異なるだけでやってることは同じである。そこで本稿ではまずDFTを解説する。

ちなみに、上の例では確かに乗算を加算に変換して易しく計算できるようになったように感じるが、これはエンコード/デコードの計算量を無視していることに注意。343=73343=7^375=16,8077^5=16,807を求めるほうが元の問題より難しいじゃないかといわれると、全くその通りである。これはDFT/IDFTにも同じことが言えて、実はこいつらの計算がかなり難しくて、元の343×49343 \times 49より難しかったりする。すなわち

計算量(計算が難しい世界での乗算)<計算量(エンコード+計算が易しい世界の乗算+デコード)計算量(計算が難しい世界での乗算)< 計算量(エンコード+計算が易しい世界の乗算+デコード)

なのである。意味ないじゃないか。
しかしながらいくつかの工夫を取り入れるとDFT/IDFTの計算量を大きく削減することができ、このアルゴリズムはFFT(高速フーリエ変換)と呼ばれる。暗号の文脈でNTTというとFFTで高速化したNTTを指すことが多く、これによって
計算量(計算が難しい世界での乗算)>計算量(エンコード+計算が易しい世界の乗算+デコード)計算量(計算が難しい世界での乗算)> 計算量(エンコード+計算が易しい世界の乗算+デコード)

とすることができるため、多くの(格子)暗号でNTTが利用されるのである!

DFT

といわけでまずDFTを解説する。世の中にDFTの解説はごまんとあるが、英語版wikipediaの定義を(適当に訳して)引用する。

DFTはNN個の複素数のシーケンス{xn}:=x0,x1,...,xN1\{\bold{x}_n\}:=x_0,x_1,...,x_{N-1}を、別の複素数シーケンス{Xk}:=X0,X1,...,XN1\{\bold{X}_k\}:=X_0,X_1,...,X_{N-1}に変換するものであり

Xk=n=0N1xnei2πkNnX_k=\sum_{n=0}^{N-1}x_n \cdot e^{-i2\pi\frac{k}{N}n}

で定義される。

これは{xn}\{\bold{x}_n\}{Xk}\{\bold{X}_k\}に変換するときの各要素の計算方法を示している。xnx_nei2πkNne^{-i2\pi\frac{k}{N}n}をそれぞれ長さNNのベクトルと見れば、n=0N1xnei2πkNn\sum_{n=0}^{N-1}x_n \cdot e^{-i2\pi\frac{k}{N}n}はそれら二つのベクトルの内積計算に他ならない。ei2πkNne^{-i2\pi\frac{k}{N}n}が周波数kkの正弦波(のベクトル表示)であることに着目すると、これは周波数kkの正弦波に対する射影すなわち{xn}\{\bold{x}_n\}に含まれている周波数kkの正弦波成分の量となる。

逆変換であるIDFTも適当に引用しておこう

xn=1Nk=0N1Xkei2πkNnx_n=\frac{1}{N}\sum_{k=0}^{N-1}X_k \cdot e^{i2\pi\frac{k}{N}n}

やっていることはDFTと変わらないので特に説明することはない。eeの符号が逆になっているところと、係数に1N\frac{1}{N}がついているところがDFTとは異なる。1N\frac{1}{N}は正規化係数なので本質的にはDFT/IDFTどちらにつけても構わないし、両方に1N\frac{1}{\sqrt{N}}をつけてもよい。
シーケンス{xn}\{\bold{x}_n\}が時間波形、{Xk}\{\bold{X}_k\}が周波数波形と考えるとこれは確かに時間領域から周波数領域の波形である。オイラーの公式的に周波数が何となく出てくるのは想像できる。
ただし時間とか周波数とかいうのは、工学応用上それが重要なだけであって、フーリエ変換の本質をとらえていない。実際NTTには時間とか周波数とか言った概念は出てこない。要はei2πkNne^{-i2\pi\frac{k}{N}n}が周波数を持つ三角関数であることはそんなに重要でないのである。DFT(およびNTT)の本質は、ei2πkNne^{-i2\pi\frac{k}{N}n}回転という性質を持っていることにある。

回転子

ei2πkNn      (i)e^{-i2\pi\frac{k}{N}n} \;\;\; (i)

の部分は、オイラーの公式

eix=cosx+isinxe^{ix}=\text{cos}x+i\text{sin}x

を使うと、複素平面上の単位円として幾何的に表すことができる。
k,n,Nk, n, Nは全てインデックスなので自然数であり、従って(i)式は、単位円をNN等分したときの(マイナスがついているので逆方向に)knkn番目の点であるとみなせる。
これはkn=Nkn=Nで1周する回転を表すものであるから、knknはmod NNの世界で考える。

このように(i)式は単位円NN等分点のどれかをknknによって示すものだった。そこでi2πN\frac{-i2\pi}{N}を定数とみなしei2πN=ωNe^{\frac{-i2\pi}{N}}=\omega_Nと置くと(i)式=ωNkn=\omega_N^{kn}と書くことができ、ずいぶん見通しがよくなる。ωN\omega_Nを回転子と呼んだりする。

さて、ωN\omega_Nがフーリエ変換の本質といってよいほど重要なものであるが、図を見れば明らかなようにωNkn\omega_N^{kn}は複素平面の相違なるNN個の値をとり、かつωNN\omega_N^{N}の時に初めて11となる数である。このようなNN乗して初めて1になる数は実は三角関数以外にも存在し、1の原始NN乗根と呼ぶ。フーリエ変換の性質の多くは1の原始NN乗根の性質に依存しており、1の原始NN乗根が存在する空間であればフーリエ変換を定義することができる。具体的には、複素数体から一般の環上へフーリエ変換を一般化できる。NTTは一般化DFTの一種であり、有限体上の1の原始NN乗根を使うものであるが、解説のためには複素平面が使えるとありがたいので、まずこのまま複素数体上のDFTで解説を進める。

1の原始NN乗根の性質

(複素数体上の)1の原始nn乗根=ωN=\omega_Nには、次のような性質がある。

  • ωNN=1\omega_N^{N}=1
  • ωN0,ωN1,ωN2,...ωNN1\omega_N^{0},\omega_N^{1},\omega_N^{2},...\omega_N^{N-1}は全て相違なる。
  • n=0N1ωNkn=0(1k<N)(ii)\sum_{n=0}^{N-1} \omega_N^{kn} = 0 \quad (1 \leq k < N) \qquad (ii)

一番上は定義より自明、二つ目も上の複素平面の図を見れば、NN等分点がすべて異なる複素数となるのは明らか、問題は3つ目(ii)式である。
(ii)式はk,nk,nの2変数あるため見通しが悪いが、行列的に書くと次のような意味である。

[000]=[ωN0ωN1ωN2ωNN1ωN0ωN2ωN4ωN2(N1)ωN0ωNn1ωN2(n1)ωN(N1)(N1)][111]\begin{bmatrix}0\\0\\\vdots\\0\end{bmatrix} = \begin{bmatrix} \omega_N^0&\omega_N^1&\omega_N^2&\cdots&\omega_N^{N-1} \\ \omega_N^0&\omega_N^2&\omega_N^4&\cdots&\omega_N^{2(N-1)}\\ \vdots&\vdots&\vdots&&\vdots\\ \omega_N^0&\omega_N^{n-1}&\omega_N^{2(n-1)}&\cdots&\omega_N^{(N-1)(N-1)}\\ \end{bmatrix} \begin{bmatrix}1\\1\\\vdots\\1\end{bmatrix}

この行列のjj行目は、ωN0=1\omega_N^0=1から始めて、ωNj\omega_N^j回転を繰り返したものの和が00である、といっている。これに関しては、等比数列の和の公式を使った証明や代数的な証明がよく見られて、確かにそうなるのだがいまいち釈然としない。なので正確さをある程度犠牲にして幾何的に考えてみよう。

複素平面で見ると下図のようになる。この図はN=6N=6なのでωN\omega_Nは-60度回転を表す。よって行列の1行目(n=1n=1)は-60度単位の回転をすべて足したものの和である。これは赤で示した位置ベクトルの和に等しく、あるベクトルに対してその逆向きのベクトルが存在することから和が0になることは容易にわかる。
2行目(n=2n=2)は-120度単位の回転ωN2\omega_N^2の和であり、3回で1周するが、和としては2周分とっている。しかしながらωN2+ωN4=ωN3\omega_N^2+\omega_N^4=\omega_N^3であり、これとωN0\omega_N^0が逆向きであるから、総和は0となる。

また、次のように考えることもできる。行列の1行目(n=1n=1)の和は、原点からスタートして距離1進み-60度回転を6回繰り返すものであるから、その軌跡(和)は正6角形となり原点(0ベクトル)に戻ってくるのである。2行目(n=2n=2)では-120度回転になるから正三角形になり、同じように原点へ戻ってくる。このように考えると(ii)式が成り立つことがわかるかと思う。

直交性

(複素数体上の)1の原始NN乗根の最も重要な性質が、直交性である。
(ii)式は実は直交性を表している。
(ii)式をk=abk=a-bとして変形すると

n=0N1ωN(ab)n=n=0N1ωNanωNbn=0(k0)\sum_{n=0}^{N-1} \omega_N^{(a-b)n} = \sum_{n=0}^{N-1} \omega_N^{an}\omega_N^{-bn} = 0 \quad (k \neq 0)

これは複素ベクトル(ωN0,ωNa,ωN2a,...\omega_N^{0}, \omega_N^{a}, \omega_N^{2a},...)と(ωN0,ωNb,ωN2b,...\omega_N^{0}, \omega_N^{b}, \omega_N^{2b},...)の内積に他ならない。
つまり、上に示した行列の相違なる2行の内積は0になる(すなわち直交しているのである)。

また、k=0k=0すなわちa=ba=bの場合は明らかに

n=0N1ωN(ab)n=n=0N11=N\sum_{n=0}^{N-1} \omega_N^{(a-b)n} = \sum_{n=0}^{N-1} 1 = N

であるから、これらをまとめて

n=0N1ωN(ab)n={0(ab)N(a=b) \sum_{n=0}^{N-1} \omega_N^{(a-b)n} = \begin{cases} 0 &(a \neq b) \\ N &(a = b) \end{cases}

と表すことが多い。両辺を1/N1/Nで正規化することで右辺をクロネッカーのデルタとし、規格直交系としてもいい。

DFTとIDFT

直交性がわかれば、DFTとIDFTがなぜ逆変換の関係にあるかを理解することができる。
DFTとIDFTの式をもう一度示す。

Xk=n=0N1xnωNkn(DFT)X_k=\sum_{n=0}^{N-1}x_n \cdot \omega_N^{-kn} \qquad \text{(DFT)}
xn=1Nk=0N1XkωNkn(IDFT)x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_k \cdot \omega_N^{kn} \qquad \text{(IDFT)}

IDFTの式へDFTの式を代入すると

xn=1Nk=0N1(n=0N1xnωNkn)ωNkn=1Nk=0N1(n=0N1xnωNk(nn))=1Nn=0N1(k=0N1xnωNk(nn))=1Nn=0N1xnk=0N1ωNk(nn) \begin{align} x_{n'} &= \frac{1}{N}\sum_{k=0}^{N-1} (\sum_{n=0}^{N-1}x_n \cdot \omega_N^{-kn}) \cdot \omega_N^{kn'} \\ &= \frac{1}{N}\sum_{k=0}^{N-1}(\sum_{n=0}^{N-1} x_n \cdot \omega_N^{k(n'-n)}) \\ &= \frac{1}{N}\sum_{n=0}^{N-1}(\sum_{k=0}^{N-1} x_n \cdot \omega_N^{k(n'-n)}) \\ &= \frac{1}{N}\sum_{n=0}^{N-1} x_n \sum_{k=0}^{N-1} \omega_N^{k(n'-n)} \end{align}

直交性を利用するとn=nn=n'のときにxn=xnx_{n'}=x_{n}となり元に戻ることがわかる。

このようにDFT/IDFTによって、あるシーケンスを同じ長さの別のシーケンスに表現変更できることがわかった。ただこの段階ではDFT/IDFTともにO(N2)O(N^2)の計算量がかかる。多項式乗算の計算量も次数NNに対してO(N2)O(N^2)であることを考えると、全然早くならないのであるが、FFTという素晴らしいアルゴリズムによってO(NlogN)O(N\text{log}N)に抑えられる。それはまたいつか・・・

次回

https://zenn.dev/ankoma/articles/72212c58d3f6ef

参考

ある程度厳密なDFT
https://mathlog.info/articles/2579
http://bygzam.seesaa.net/article/150579016.html
https://qiita.com/AngrySadEight/items/0dfde26060daaf6a2fda
https://qiita.com/TumoiYorozu/items/5855d75a47ef2c7e62c8
https://www.creativ.xyz/fast-fourier-transform/
https://drive.google.com/file/d/1gMdA1Krok9qpxeEPAyFSnU8ekFi5tYis/view
http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/fs_comp.html

Discussion

ログインするとコメントできます