2024年04月01日(月) [長年日記]
■ [c++][dev][math] 高速フーリエ変換
高速フーリエ変換についてメモ。
高速フーリエ変換とは、離散フーリエ変換を高速に計算するアルゴリズムのこと。
離散フーリエ変換
離散フーリエ変換とは、項数 N の複素数列 \((a_n)_{n=0}^{N-1}\) から、次で定める項数 N の複素数列 \((b_n)_{n=0}^{N-1}\) を求める操作のこととする。
\[ \ b_n = \sum_{x=0}^{N-1}{a_x \exp(-i \frac{2 \pi nx}{N})} \]
\(W_N=\exp\left(-i\frac{2 \pi}{N}\right)\) とおくと、この式は次のように書ける。
\[ \ b_n = \sum_{x=0}^{N-1}{W_N^{nx} a_x} \]
この式の掛け算の回数は N 回で、すべての 0≦n≦N に対して計算するには N×N 回の掛け算が必要になる。
高速フーリエ変換の解説は N を2のべき乗に制限するものが多いが、ここでは非素数の N に対して（素朴に計算するのに比べて）高速に計算できるアルゴリズムを考える。なお、以下の話は N が素数であっても適用できる（R=N、S=1 とすればよい）が、その場合は計算量が減らない。
計算量を減らすための式変形
R を N の約数、S を N/R とする。求めたい複素数列 \((b_n)_{n=0}^{N-1}\) を項数 S の R 個の複素数列 \((b_{Rm})_{m=0}^{S-1}, (b_{Rm+1})_{m=0}^{S-1}, ..., (b_{Rm+R-1})_{m=0}^{S-1} \) に分割する。数列 \((b_{Rm+k})_{m=0}^{S-1}\) の一般項は次のように表せる（この式変形の導出は本稿の最後に示す）。
\[ \ b_{Rm+k} = \sum_{y=0}^{S-1}{W_S^{my} \left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)} \]
この式の掛け算の回数を考えると、
- 括弧内の計算に R+1 回の掛け算が必要で、これをすべての 0≦y≦S-1 と 0≦k≦R-1 に対して計算するので全体では (R+1)×S×R 回
- \(W_S^{my}\)との積が S 回でこれをすべての 0≦m≦S-1 と 0≦k≦R-1 に対して計算するので全体では S×S×R 回
になる（と思う）。合わせて (R+1)×S×R + S×S×R = (R+S+1)×N 回の掛け算が必要で、これはほとんどのケースで元の式の N×N 回より少ない回数になる。括弧内の計算がmに依存しておらず計算結果を再利用できることが効いている。ただし、元の式には無かった \(W_S^{my}\) と \(W_R^{kx}\) の計算が新たに発生していることには注意。また、\((b_n)_{n=0}^{N-1}\) を分割する際に並び順が変わっているので、並び替えのコストも発生する。
ここで、\(b_{Rm+k}\) の式が項数 S の複素数列 \(\left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)_{y=0}^{S-1}\) に対する離散フーリエ変換の形になっていることに注目してほしい。これにより、R個の数列 \((b_{Rm+k})_{m=0}^{S-1}\) それぞれの計算に対しても上記と同じ式変形で計算量を減らすことができ、この操作は数列の要素数が1になるまで再帰的に行えることになる。
高速フーリエ変換の実装
上のようにして計算量を減らしたコードと減らさず素朴に計算するコードの比較を行うプログラムを書いた。200行くらいあるので全体はGistに置きここにコードを抜粋する。コンパイルにはC++20が必要だが、もっと古いC++でコンパイルできるように直すのは難しくない。
高速フーリエ変換の計算をする関数 fft() は次の通り。
// inputの値を離散フーリエ変換した結果を返す
template<typename T>
std::vector<std::complex<T>> fft(
std::vector<std::complex<T>> input,
const twiddle_factor_map<T>& tf_map)
{
std::vector<std::complex<double>> work(input.size());
detail::fft(
input.begin(),
input.end(),
work.begin(),
detail::find_divisor(input.size()),
tf_map);
return input;
}
- 入力をvectorで与え結果をvectorで受け取る。
- \( W_N \) のべき乗については事前に計算したものを tf_map というマップに設定して渡すようにした。
- 計算の作業領域用に入力と同じサイズのvectorを使用している。
上で呼び出している実際の処理を行う関数 detail::fft() は次の通り。
// 範囲[first, last)の値を離散フーリエ変換した結果を範囲[first, last)に上書きする。
// work_firstは[first, last)と同じサイズの作業用範囲の先頭。
// radixには範囲[first, last)の要素数の約数を指定する。
template<typename T, typename Iter>
void fft(
Iter first,
Iter last,
Iter work_first,
std::size_t radix,
const twiddle_factor_map<T>& tf_map)
{
const auto n = static_cast<std::size_t>(std::distance(first, last));
if (n <= 1) return;
const auto& tfs_r = tf_map.at(radix);
const auto& tfs_n = tf_map.at(n);
const auto s = n / radix;
for (std::size_t i = 0; i < s; ++i) {
auto it = work_first + i;
for (std::size_t j = 0; j < radix; ++j) {
*it = 0;
for (std::size_t k = 0; k < radix; ++k) {
*it += *(first + i + s * k) * tfs_r[j * k % radix];
}
*it *= tfs_n[i * j];
it += s;
}
}
for (std::size_t i = 0; i < radix; ++i) {
const auto it = work_first + s * i;
fft(it, it + s, first, find_divisor(s), tf_map);
}
for (std::size_t i = 0; i < s; ++i) {
auto it = first + radix * i;
for (std::size_t j = 0; j < radix; ++j) {
*it = *(work_first + i + s * j);
++it;
}
}
}
- detail::fft() は自分自身を呼び出している。再帰呼び出しの中でvectorを作り直さなくていいように、イテレータ範囲でデータを渡すようにしている。
- 作業領域内にR個の複素数列 \(\left( W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}} \right)_{y=0}^{S-1}\) を設定している。
- 再帰呼び出しでは作業領域を入力データ、入力データ範囲を作業領域として渡している。
- 再帰呼び出しが終わったら作業領域内のデータを並び替えたものを計算結果として入力データ範囲に設定している。
実行結果
Gistに置いたファイルの実行結果は以下のような感じ。素朴に計算した場合の実行時間（dft）を高速に実行した場合の実行時間（fft）で割った値を ratio として出力している。delta は計算結果の差の絶対値の最大値。
なお、NumPyのFFTルーチンと比べると遅い（項数 129140163 = 3の17乗 のときに上のコードは23秒、NumPyは10秒で計算）ので、上のコードもしくは本アルゴリズムは改善の余地が大いにある。
入力の項数が少ないうちは素朴に計算する方が速い。
size 2 ratio 0.379 dft 0.000 fft 0.000 delta 0.000e+00 size 3 ratio 0.858 dft 0.000 fft 0.000 delta 0.000e+00 size 4 ratio 0.897 dft 0.000 fft 0.000 delta 0.000e+00 size 5 ratio 0.950 dft 0.000 fft 0.000 delta 0.000e+00 size 6 ratio 0.400 dft 0.000 fft 0.000 delta 7.589e-15 size 7 ratio 0.907 dft 0.000 fft 0.000 delta 0.000e+00 size 8 ratio 0.640 dft 0.000 fft 0.000 delta 3.794e-15 size 9 ratio 2.072 dft 0.000 fft 0.000 delta 0.000e+00 size 10 ratio 0.621 dft 0.000 fft 0.000 delta 7.944e-15 size 11 ratio 0.966 dft 0.000 fft 0.000 delta 0.000e+00 size 12 ratio 0.600 dft 0.000 fft 0.000 delta 2.197e-14 size 13 ratio 1.021 dft 0.000 fft 0.000 delta 0.000e+00 size 14 ratio 0.704 dft 0.000 fft 0.000 delta 1.507e-14 size 15 ratio 0.524 dft 0.000 fft 0.000 delta 1.888e-14
項数2000くらいで、大きくて40倍くらいの差が出てくる。項数が素数の場合は速くならない。
size 1987 ratio 0.998 dft 0.010 fft 0.010 delta 0.000e+00 size 1988 ratio 22.035 dft 0.010 fft 0.000 delta 3.122e-12 size 1989 ratio 46.875 dft 0.010 fft 0.000 delta 2.890e-12 size 1990 ratio 8.610 dft 0.010 fft 0.001 delta 2.621e-12 size 1991 ratio 10.168 dft 0.009 fft 0.001 delta 2.390e-12 size 1992 ratio 19.410 dft 0.010 fft 0.000 delta 2.912e-12 size 1993 ratio 1.000 dft 0.009 fft 0.009 delta 0.000e+00 size 1994 ratio 1.997 dft 0.009 fft 0.005 delta 3.429e-12 size 1995 ratio 46.528 dft 0.010 fft 0.000 delta 2.430e-12 size 1996 ratio 3.932 dft 0.010 fft 0.002 delta 4.590e-12 size 1997 ratio 1.039 dft 0.010 fft 0.009 delta 0.000e+00 size 1998 ratio 33.008 dft 0.010 fft 0.000 delta 2.993e-12 size 1999 ratio 0.892 dft 0.010 fft 0.011 delta 0.000e+00
項数が2のべき乗の場合。数が大きくなるとどんどん差がつく。
--- 2のべき乗 --- size 2 ratio 1.911 dft 0.000 fft 0.000 delta 0.000e+00 size 4 ratio 0.982 dft 0.000 fft 0.000 delta 0.000e+00 size 8 ratio 0.741 dft 0.000 fft 0.000 delta 7.105e-15 size 16 ratio 0.908 dft 0.000 fft 0.000 delta 1.137e-14 size 32 ratio 1.188 dft 0.000 fft 0.000 delta 6.071e-14 size 64 ratio 2.473 dft 0.000 fft 0.000 delta 2.050e-13 size 128 ratio 4.226 dft 0.000 fft 0.000 delta 3.422e-13 size 256 ratio 8.059 dft 0.000 fft 0.000 delta 1.833e-12 size 512 ratio 15.466 dft 0.001 fft 0.000 delta 2.854e-11 size 1024 ratio 28.213 dft 0.003 fft 0.000 delta 8.084e-11 size 2048 ratio 56.931 dft 0.012 fft 0.000 delta 3.748e-10 size 4096 ratio 104.531 dft 0.044 fft 0.000 delta 7.867e-10 size 8192 ratio 204.091 dft 0.172 fft 0.001 delta 1.531e-09 size 16384 ratio 361.594 dft 0.644 fft 0.002 delta 2.899e-08 size 32768 ratio 715.049 dft 2.648 fft 0.004 delta 6.994e-08 size 65536 ratio 1349.311 dft 11.721 fft 0.009 delta 3.444e-07
項数が3のべき乗の場合。
--- 3のべき乗 --- size 3 ratio 1.306 dft 0.000 fft 0.000 delta 0.000e+00 size 9 ratio 1.336 dft 0.000 fft 0.000 delta 0.000e+00 size 27 ratio 1.078 dft 0.000 fft 0.000 delta 4.550e-14 size 81 ratio 2.865 dft 0.000 fft 0.000 delta 2.558e-13 size 243 ratio 8.203 dft 0.000 fft 0.000 delta 2.752e-12 size 729 ratio 23.335 dft 0.001 fft 0.000 delta 1.347e-11 size 2187 ratio 57.941 dft 0.013 fft 0.000 delta 1.113e-10 size 6561 ratio 169.265 dft 0.116 fft 0.001 delta 2.806e-09 size 19683 ratio 466.513 dft 1.041 fft 0.002 delta 4.331e-08 size 59049 ratio 1307.727 dft 9.614 fft 0.007 delta 8.864e-08
項数の素因数が異なる数になる場合。
--- 素数の積 --- size 2 ratio 1.211 dft 0.000 fft 0.000 delta 0.000e+00 size 6 ratio 0.452 dft 0.000 fft 0.000 delta 7.119e-15 size 30 ratio 1.100 dft 0.000 fft 0.000 delta 3.197e-14 size 210 ratio 6.740 dft 0.000 fft 0.000 delta 2.975e-12 size 2310 ratio 61.381 dft 0.014 fft 0.000 delta 6.370e-11 size 30030 ratio 576.535 dft 2.444 fft 0.004 delta 1.515e-07
項数が素数の場合は同じくらいの速さ。
--- 素数 --- size 2 ratio 0.943 dft 0.000 fft 0.000 delta 0.000e+00 size 3 ratio 0.777 dft 0.000 fft 0.000 delta 0.000e+00 size 5 ratio 0.982 dft 0.000 fft 0.000 delta 0.000e+00 size 7 ratio 1.054 dft 0.000 fft 0.000 delta 0.000e+00 size 11 ratio 1.069 dft 0.000 fft 0.000 delta 0.000e+00 size 13 ratio 1.436 dft 0.000 fft 0.000 delta 0.000e+00 size 1009 ratio 0.999 dft 0.003 fft 0.003 delta 0.000e+00 size 10009 ratio 0.967 dft 0.269 fft 0.278 delta 0.000e+00 size 100003 ratio 0.988 dft 28.695 fft 29.054 delta 0.000e+00
参考
次のページが参考になった。
- FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法 （京都大学数理解析研究所）
上で行った変換は混合基数アルゴリズムに該当する。
式変形の詳細
最後に、先に示した \((b_{Rm+k})_{m=0}^{S-1}\) の一般項の導出をメモしておく。
計算の過程で次を使う。xは整数、Nは自然数、RはNの約数、Sは N/R とする。
\[ \ W_N^{Nx} = 1, \quad \ W_N^{Rx} = W_S^x, \quad \ \sum_{x=0}^{N-1}{f(x)} = \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{f(Sx+y)}} \]
これらを使うことで次のように変形できる。
\begin{align} \ b_{Rm+k} &= \sum_{x=0}^{N-1}{W_N^{(Rm+k)x} a_x} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_N^{(Rm+k)(Sx+y)} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_N^{Nmx} W_N^{Rmy} W_N^{kSx} W_N^{ky} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{\sum_{x=0}^{R-1}{W_S^{my} W_R^{kx} W_N^{ky} a_{Sx+y}}} \\ \ &= \sum_{y=0}^{S-1}{W_S^{my} W_N^{ky} \sum_{x=0}^{R-1}{W_R^{kx} a_{Sx+y}}} \\ \end{align}
2024年04月02日(火) [長年日記]
■ [windows][net][howto] Hyper-V上のVMを内部仮想スイッチ経由で繋ぎポートフォワードで外部からアクセス可能にする
Hyper-V上のVMに外部からアクセスできるようにしたいが、ホストの物理ネットワークに直接接続する（「外部仮想スイッチ」経由で繋ぐ）のではなく、仮想ネットワークに接続して必要なサービスだけをポートフォワードにより外部からアクセスできるようにしたい。
そのためには「内部仮想スイッチ」と「NATネットワーク」を作成して両者を関連付け、VMを内部仮想スイッチに接続してIPアドレスを設定し、更にポートマッピングの設定をすればいいようだ。
内部仮想スイッチの作成
Hyper-Vにはデフォルトの仮想スイッチとして「Default Switch」が用意されているが、このスイッチはホストの再起動でIPアドレスが変わりポートフォワードの設定がうまくできないようだった。なので、それとは別の仮想スイッチを作成する。
仮想スイッチはGUIのHyper-Vマネージャから次のようにして作成できる。
- 対象の仮想化サーバーに接続する。
- 仮想化サーバーに対するメニューから［操作］→［仮想スイッチマネージャー］を選び仮想スイッチマネージャーの画面を開く。
- 左ペインの「新しい仮想ネットワークスイッチ」を選択し、右ペインで「内部」を選択して「仮想スイッチの作成」ボタンをクリックする。
- 新しい仮想スイッチが追加されるので、名前を付けて「OK」ボタンをクリックする。
なお、GUIを使わずNew-VMSwitchコマンドを使うことでも作れそう。
次に仮想スイッチへIPアドレスを設定する。New-NetIPAddressコマンドで設定できる。
New-NetIPAddress -InterfaceIndex 「仮想スイッチインデックス」 `
-IPAddress 「IPアドレス」 -PrefixLength 「プレフィックス長」
ここで、「仮想スイッチインデックス」はGet-NetAdapterコマンドやGet-NetIPAddressコマンドにより取得できる。
NATネットワークの作成
NATネットワークはNew-NetNatコマンドにより作成できる。
New-NetNat -Name 「NATネットワーク名」 -InternalIPInterfaceAddressPrefix 「IPアドレス/プレフィックス長」
ここで、「IPアドレス/プレフィックス長」には仮想スイッチへ付与したIPアドレスが属するネットワークを指定する。
VMがNATを使えるようにする
仮想化サーバー上にVMを作成し、Hyper-Vマネージャから次を行う。
- 仮想マシンに対するメニューから［操作］→［設定］を選び仮想マシンの設定画面を開く。
- 左ペインの「ネットワークアダプター」を選択し、右ペインで「仮想スイッチ」を上で作成した内部仮想スイッチに変更して「OK」ボタンをクリックする。
これにより、NATを使ってVMからインターネットなどホスト側のネットワークへ接続できるようになる。
はじめ、既存のネットワークアダプタの設定はそのまま（Default Switchを使用）にしてネットワークアダプターを新規追加してそちらに新しい内部仮想スイッチを設定するようにしていたのだけれど、そうするとポートフォワードが機能せず時間を無駄にした。Default SwitchもNATを使っており、NATを使うネットワークアダプタが複数あると機能しないのだと思う。
ポートマッピングを設定する
まず、VMの当該ネットワークアダプターに固定のIPアドレスを付与する。これはVM上で動くOSの機能を使って行う。
その上で、Add-NetNatStaticMappingコマンドによりポートマッピングを設定する。
Add-NetNatStaticMapping 「NATネットワーク名」 -Protocol TCP `
-ExternalIPAddress 「VMのホストのIPアドレス」 -ExternalPort 「VMのホストが公開するポート」 `
-InternalIPAddress 「VMのIPアドレス」 -InternalPort 「VM上のサービスのポート」
これにより、ホスト側のネットワークからVM上の当該サービスへアクセスできるようになる。
参考
- Hyper-V を使用して仮想スイッチを作成して構成する （Microsoft Learn）
- NAT ネットワークの設定 （Microsoft Learn）
2024年04月15日(月) [長年日記]
■ [howto][unix][security] SSHで未知のホストの確認プロンプトを表示しないようにする
SSHで未知のホストに接続しようとすると次のようなプロンプトが表示されるが、これを表示したくないことがある（スクリプト内でSSHを実行するときなど）。
The authenticity of host 'localhost (::1)' can't be established. ED25519 key fingerprint is SHA256:V1NLamwjGYVi4hjxTg1Ko96YVh4xhsoYI2sc1rPGdno. This key is not known by any other names Are you sure you want to continue connecting (yes/no/[fingerprint])?
検索するとSSHのStrictHostKeyCheckingオプションを no にする方法がヒットするが、それよりも、対象ホストの公開鍵を事前に ~/.ssh/known_hosts に登録しておく方が自然だと思う。
公開鍵の取得は ssh-keyscanコマンドでできる。
$ ssh-keyscan -t ed25519 localhost # localhost:22 SSH-2.0-OpenSSH_8.7 localhost ssh-ed25519 AAAAC3NzaC1lZDI1NTE5AAAAIHO8iGPxukxiQPrGGJqtTDdfaRYxNvf+CegVOvK2JyBy
ssh-keyscanの出力はそのまま ~/.ssh/known_hosts ファイルに追加できる。~/.ssh/known_hosts に行を追加することで件のプロンプトは表示されなくなる。
ただし、ssh-keyscan の接続先が正しいホストでない可能性があるので、別途公開鍵またはそのフィンガープリントの確認は必要。
接続先ホストの公開鍵はホストの /etc/ssh/ あたりにあるはず。フィンガープリントはssh-keygen -lで表示できる。
$ cat /etc/ssh/ssh_host_ed25519_key.pub ssh-ed25519 AAAAC3NzaC1lZDI1NTE5AAAAIHO8iGPxukxiQPrGGJqtTDdfaRYxNvf+CegVOvK2JyBy $ ssh-keygen -lf /etc/ssh/ssh_host_ed25519_key.pub 256 SHA256:V1NLamwjGYVi4hjxTg1Ko96YVh4xhsoYI2sc1rPGdno no comment (ED25519) $ ssh-keygen -lf <(ssh-keyscan -t ed25519 localhost) # localhost:22 SSH-2.0-OpenSSH_8.7 256 SHA256:V1NLamwjGYVi4hjxTg1Ko96YVh4xhsoYI2sc1rPGdno localhost (ED25519)
2024年04月22日(月) [長年日記]
■ [howto][unix] dnfコマンドでDVDからパッケージをインストールする
Rocky Linux 9を使っている。Rocky Linux 9のパッケージマネージャはDNF。
DNFでDVDからパッケージをインストールするにはどうするのか。
Rocky Linux 9のDVDイメージには media.repo というファイルが含まれているのでこれを使って簡単にできるのかと思ったが、どういう使い方が想定されているのかよくわからなかった。このファイルを /etc/yum.repo.d/ にコピーして baseurl 指定を追加すれば使えそうだけど、余計なファイルを作るのは避けたい。
dnfコマンドのリファレンスを見ると--repofrompathオプションと--setoptオプションがあるので、これらを使って次のようにすれば .repo ファイルを用意せずにインストールができた。
dnf --disablerepo="*" \ --repofrompath dvd1,/media/dvd/BaseOS \ --repofrompath dvd2,/media/dvd/AppStream \ --setopt=dvd1.gpgkey=/etc/pki/rpm-gpg/RPM-GPG-KEY-Rocky-9 \ --setopt=dvd2.gpgkey=/etc/pki/rpm-gpg/RPM-GPG-KEY-Rocky-9 \ install boost tar
上記は /media/dvd/ にDVDをマウントしている前提。