2002年04月01日(月) 入社式の日
■ 会社
- 11:00 久しぶりに入社式に出席。会社紹介のビデオが流され高幣さん村地さん児玉くんなどが登場。新人は一言スピーチ。服部さんの話は退屈。
- hnsに再び脆弱性。クロスサイトスクリプティングだそうです。またバージョンアップしなきゃ。
- 16:55 kondaraをダウンロードしてみたがなかなか終わらない。うーんまずかったかなあ。今2枚目のCDROMイメージをダウンロードしているところ。
- kondaraを落とす上で参考になったサイト。.isoという拡張子のファイルをダウンロードしているのだ。本屋に行っても売ってなかったから。
- そういえばDatulaの新しいバージョンもリリースされているのだった。Windows Updateはもう実行したのだ。そのうちNetscapeもバージョンアップしなくちゃいけないんだろうなあ。
- 19:40 kondaraのダウンロードは終わっている。が、SFAの見積もりと実行計画、WineとGIGAの完了報告をやらねばならないのでインストールをやる時間がとれなそう。残念。
■ テレビ／本
- パパパパパフィースペシャル。終わってしまった。
- ガキの使いやあらへんで。トークのみ。仮面ライダーエロスの話とか。
- 日経Windowsプロ4月号読み終わり。どうということはない。
2003年04月01日(火) 名古屋で自転車女による通り魔続発中
■ [life] 告別式
長い闘病の末亡くなってしまった同期の告別式へ行った。時間があったので、少し遠い駅で降りて会場まで歩いた。暖かく晴れた日の下、広い桜並木の道をゆっくり歩く。
会場に着いて席に座ろうとしたら、小坂さんに一番前の席へと導かれた。会場にはキーボードを弾く人がいて音楽が流れている。程なく牧師さん（？）の進行で式が始まる。キリスト教式なので、黙祷の後賛美歌を歌ったり聖書を読んだり。故人の紹介などはなかった。
そして献花。一番前に座っていたことで、一般会葬者の中では一番最初に献花することになってしまった。なにやら説明があったがよく聞き取れない。ので、しどろもどろになりあまり落ち着いて献花できず残念。
献花後そのまま退席し、会場後部で待つ。全員の献花が終わると椅子が片付けられ、棺が中央へ運ばれた。棺が開かれ、そこへ献花された花を添え最後のお別れ。とても綺麗な顔だった。泣きそうになったが涙はこらえた。どうぞ安らかに。
喪主の挨拶の後、棺が運び出され終了・解散。皆はそのまま会社へ行くがおれには行く所が無い。同期と一緒に食事をした後家へ帰る。帰りも少し遠めの駅から歩いた。電車賃をケチるため。そのためか、またもやヘトヘトになり家に着くなり寝てしまった。
■ やること
- 退職後の手続きをまとめる
2022年04月01日(金)
■ [windows][howto] PowerShellでJSONを手軽に整形
Windowsで1行のJSONを整形して表示したかった。
フォーマットにこだわらなければ、PowerShell上でConvertFrom-JsonとConvertTo-Jsonを組合せて次のようにするのが一番手軽に思えたがどうだろうか。
convertfrom-json '{"a": 1, "b1": {"b2": {"b3": {"b4": 2}}}, "c": ["hello", "world"]}' `
| convertto-json
{
"a": 1,
"b1": {
"b2": {
"b3": "@{b4=2}"
}
},
"c": [
"hello",
"world"
]
}
ただし、よく見ると "b3" の所がおかしい。
これは ConvertTo-Json に -Depth オプションを指定すれば解決する。
convertfrom-json '{"a": 1, "b1": {"b2": {"b3": {"b4": 2}}}, "c": ["hello", "world"]}' `
| convertto-json -depth 3
{
"a": 1,
"b1": {
"b2": {
"b3": {
"b4": 2
}
}
},
"c": [
"hello",
"world"
]
}
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}
● 私の名前は高田といいます [小坂さんの教会、一度行って見てみたいです。埼玉ですか？都会はあまりしらないのですが初めての人でも行かれますか？]