メモの日々


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

参考

次のページが参考になった。

上で行った変換は混合基数アルゴリズムに該当する。

式変形の詳細

最後に、先に示した \((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マネージャから次のようにして作成できる。

  1. 対象の仮想化サーバーに接続する。
  2. 仮想化サーバーに対するメニューから[操作]→[仮想スイッチマネージャー]を選び仮想スイッチマネージャーの画面を開く。
  3. 左ペインの「新しい仮想ネットワークスイッチ」を選択し、右ペインで「内部」を選択して「仮想スイッチの作成」ボタンをクリックする。
  4. 新しい仮想スイッチが追加されるので、名前を付けて「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マネージャから次を行う。

  1. 仮想マシンに対するメニューから[操作]→[設定]を選び仮想マシンの設定画面を開く。
  2. 左ペインの「ネットワークアダプター」を選択し、右ペインで「仮想スイッチ」を上で作成した内部仮想スイッチに変更して「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上の当該サービスへアクセスできるようになる。

参考


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をマウントしている前提。