この記事はD言語 Advent Calendarの18日目の記事です.
この記事はTUT Advent Calendar 2018の17日目の記事です.
この記事はD言語くんAdvent Calendarの16日目の記事です.
D言語くんの生態に学術・産業界から注目が集まっている.
本報告では,D言語くんが支配する都市の存在について,またその都市における人々やD言語くんの生態について報告する.
世の中には専門家が書いた精度保証付きの高機能・高性能な数値積分ライブラリが多数あります.
しかし,実際はそこまで機能・性能を求めていなかったり,関数の形がキレイで積分が簡単だったりで,そのようなライブラリがオーバースペックな場合が結構あります.
しかも,そのようなライブラリは高度な技法を使っているため,専門的な知識がないとコードの理解ができないという問題があります.
でも,自分で台形公式やシンプソン則を書いてもあまり精度が出ないような場合もあります.
そのようなときに役に立つ二重指数関数型積分公式(Double Exponential: DE公式)を用いた数値積分ライブラリdeint
を少し前に公開しましたので,利用方法を紹介します.
私が現在勉強しながら研究している非線形性に関する理論を簡単にまとめておきます.
もちろん,まだ一回も発表もしていない研究の内容をそのまま書くのはさすがにいけないので,ここではすでに世間ではよく知られているような内容のみをまとめておきます.
ここで,\(S(f)\)を\(s(t)\)のスペクトル,\(\left\lvert S(f) \right\rvert^2\)を電力スペクトル密度とそれぞれ呼びます.
\(\delta(x)\):ディラックのデルタ関数
\(\text{sinc}(x)=\frac{\sin \pi x}{\pi x}\):sinc関数
まず,フーリエ変換について復習しておきます. 正弦波は周期関数なので,正弦波のフーリエ変換(スペクトル)はフーリエ級数に展開したときの係数で表せます.
ここで,周期\(T_0 = \frac{1}{f_0}\)の任意の信号\(s(t)\)の複素フーリエ係数が\(c_n\)であれば,\(s(t)\)のスペクトルは次のように表せます.
\[S(f) = \sum_{n=-\infty}^{\infty} c_n \delta(f + n f_0)\]そのため,まずは複素フーリエ係数\(c_n\)を求めてみます. ここで,\(c_n\)とは複素正弦波\(e^{j 2 \pi n f_0 t }\)と信号\(s(t)\)の内積で計算されます.
\[c_n = f_0 \int_{-\frac{1}{2 f_0}}^{\frac{1}{2 f_0}} s(t) e^{-j 2\pi nf_0 t} dt\]\(s(t)=e^{j2\pi f_0 t}\)のとき,\(c_1=1\)で\(c_{n\not=1}=0\)になりますね?
では,話を戻して正弦波,つまり\(s(t)=\cos(2\pi f_0 t)\)のときにはどうなるでしょうか? ここでは,オイラーの公式より
\[\cos(2\pi f_0 t) = \frac{1}{2} \left(e^{j2\pi f_0 t} + e^{-j2\pi f_0 t} \right)\]になることから直ちに
\[c_n = \begin{cases} \frac{1}{2}, & n = 1 \\ \frac{1}{2},& n = -1 \\ 0, & \text{otherwise} \end{cases}\]がわかります. つまり,正弦波のフーリエ変換はは次のようになります.
\[\mathcal{F}[\cos(2 \pi f_0 t)] = \frac{1}{2} \delta(f-f_0) + \frac{1}{2} \delta(f+f_0)\]グラフで表すと以下の図のようになります. ここでは,ディラックのデルタ関数を表現するために矢印を用いて,その高さを係数\(c_n\)にしています.
次に,よくある矩形波のスペクトルを考えてみます. ここでは周期\(T_0=\frac{1}{f_0}\)の次のような矩形波\(r(t)\)を考えます.
\[r(-\frac{1}{2f_0}<t<\frac{1}{2f_0}) = \begin{cases} 1, & \lvert 2 f_0 t \rvert < 0.5 \\ -1, & |2 f_0 t| > 0.5 \\ 0, & 2 f_0 t = 0, 1, -1 \\ \end{cases}\]ただし,上の式の定義外については\(r(t)\)が周期\(T_0=\frac{1}{f_0}\)の周期関数となるように定義します.
この矩形波は\(0\)を中心に半分の区間では\(1\)に,残り半分の区間では\(-1\)になっています.
この矩形波のスペクトルを調べるために,まずは複素フーリエ係数を調べます. \(r(t)\)が偶関数であることに注目すると,
\[\begin{align} c_n &= f_0 \int_{-\frac{1}{2f_0}}^{\frac{1}{2f_0}} r(t) e^{-j2\pi nf_0 t} dt \\ &= 2 f_0 \int_{0}^{\frac{1}{2f_0}} r(t) \cos(2\pi n f_0 t) dt \\ &= 2 f_0 \left( \int_{0}^{\frac{1}{4f_0}} \cos(2\pi n f_0 t) dt - \int_{\frac{1}{4f_0}}^{\frac{1}{2f_0}} \cos(2\pi n f_0 t) dt \right) \end{align}\]\(n\not=0\)においては
\[\begin{align} c_n &= 2f_0 \left(\left[\frac{\sin(2\pi n f_0 t)}{2\pi n f_0}\right]_{0}^{\frac{1}{4f_0}} - \left[\frac{\sin(2\pi n f_0 t)}{2\pi n f_0}\right]_{\frac{1}{4f_0}}^{\frac{1}{2f_0}} \right) \\ &= \frac{1}{\pi n} (2\sin(0.5 \pi n) - \sin(\pi n)) \\ &= \frac{\sin(\frac{\pi}{2} n)}{\frac{\pi}{2} n} = \text{sinc}(n/2) \end{align}\]を得ます. また,\(n=0\)では\(c_0=0\)です. そのため,矩形波のスペクトルは
\[S(f) = \mathcal{F}[r(t)] = \sum_{n=-\infty,n\not=0}^{\infty} \text{sinc}(n/2) \delta(f - n f_0)\]となります.
ここで,次のような非線形関数を考えます.
\[\begin{align} g(x) = \begin{cases} 1, & x > 0\\ -1, & x < 0\\ 0, & x= 0 \end{cases} \end{align}\]この非線形関数\(g(x)\)は\(x\)が少しでも正に値を取れば\(1\)に,少しでも負に値を取れば\(-1\)に,\(0\)であれば\(0\)になるような関数です.
この関数は矩形波の定義に似ているかもしれませんが,周期関数ではありません.
この非線形関数を用いて矩形関数は次のように定義できます.
\[r(t) = g(\cos(2\pi f_0 t))\]これを周波数領域で考えてみると,正弦波の2本の線スペクトルは\(g(x)\)な非線形特性を有するシステムを通過することで無限個の線スペクトルに化ける,と考えられます.
では,任意の非線形関数が与えられたとき,その非線形関数のどのような特徴が出力のスペクトルを決定づけるのでしょうか?
つまり,\(g(x)\)の特性だけがわかっているときに,正弦波を実際に適用することなくその出力のスペクトルを予想することはできるのでしょうか?
この問に関する答えは,意外かもしれませんが,チェビシェフ多項式という直交多項式にあります.
実は,任意の関数\(g(x)\)をチェビシェフ多項式\(T_k(x)\)で展開したときの係数と,\(g(\cos(2\pi f_0 t))\)のスペクトルにはある関係があります.
ここで,チェビシェフ多項式\(T_k(x)\)とは\(m\not=n\)において次の関係を満たすものです.
\[I_{m,n} = \int_{-1}^{1} T_m(x) T_n(x) \frac{1}{\sqrt{1-x^2}} dx = 0\]具体的には
\[\begin{align} T_0(x) &= 1 \\ T_1(x) &= x\\ T_2(x) &= 2x^2 - 1\\ T_3(x) &= 4x^3 - 3x\\ T_4(x) &= 8x^4 - 8x^2 + 1 \end{align}\]となっています. また,面白い特性として
\[T_n(\cos(x)) = \cos(nx)\]となります.
このチェビシェフ多項式を用いて\(g(\cos(2\pi f_0 t))\)のスペクトルは次のように記述できます.
\[\mathcal{F}[g(\cos(2\pi f_0 t))] = \frac{a_0}{2} \delta(f) + \sum_{n=0}^{\infty} \frac{a_n}{2} \left(\delta(f-nf_0) + \delta(f+nf_0) \right)\] \[a_n = \frac{1}{\pi} \int_{-1}^{1} g(x) T_n(x) \frac{1}{\sqrt{1-x^2}} dx\]なぜこんなことをするのでしょうか?
それは,\(a_n\)の積分で\(x=\cos(2\pi f_0 t)\)と置換積分することでわかります.
\(\frac{dx}{dt} = -2\pi f_0 \sin(2\pi f_0 t)\) で積分区間は\(\frac{1}{2f_0} \rightarrow 0\)です. また,\(dx/\sqrt{1-x^2} = -2\pi f_0 dt\)になるので,
\[\begin{align} a_n &= \frac{1}{\pi} \int_{\frac{1}{2f_0}}^{0} g(\cos(2\pi f_0 t)) T_n(\cos(2\pi f_0 t)) (-2 \pi f_0) dt \\ &= 2f_0 \int_{0}^{\frac{1}{2f_0}} g(\cos(2\pi f_0 t)) \cos(2\pi n f_0 t) dt \end{align}\]となります. よく観察してみれば,これは複素フーリエ係数\(c_n\)そのものだということがわかります.
ん?じゃあこれは単に複素フーリエ係数を求める積分を変数変換しただけでは?と思うかもしれません.
そのとおりです!
ですが,もう少し考えていきましょう.
チェビシェフ多項式\(T_n(x)\)と係数\(a_n\)を用いることで,任意の非線形関数\(g(x)\)を表現することができます. つまり,
\[g(x) = \sum_{n=0}^{\infty} a_n T_n(x)\]となる\(a_n\)があります. そしてこの\(a_n\)はフーリエ係数なのです.
何が嬉しいかというと, ある非線形な入出力特性をもつ素子や電気回路の特性をチェビシェフ多項式で近似したいとします.
この素子や回路の入力に正弦波を入れ,素子の出力をフーリエ変換すると,その係数が素子をチェビシェフ多項式で近似したときの係数になる,というのです!
さらに,逆も成り立ちます. チェビシェフ多項式を用いて特性が近似されたある素子に正弦波を入れたときの出力信号のフーリエ係数はチェビシェフ多項式の係数なのです.
以上までは,入力信号を正弦波にしたときの結果を述べました. ここではもう少し拡張してみて,入力信号が逆正弦分布という確率分布に従う場合を考えてみます. この逆正弦分布の確率密度関数は
\[p(x) = \frac{1}{\pi} \frac{1}{\sqrt{1-x^2}}\]となります. \(\theta \sim U(0,2\pi)\)のとき\(\sin \theta\),や\(\cos \theta\)は逆正弦分布に従います.
\(a_n\)を求めるときに用いた積分は,実は逆正弦分布に従う入力信号\(x\)について\(\mathbb{E}[g(x)T_n(x)]\)を求めていただなのです. そして,逆正弦分布に従う入力信号についてはチェビシェフ多項式が直交多項式となり,その係数はスペクトルと面白い関係を持ちます.
Rust入門3日目のけーさんです.
Rustを使ってスタックを実装したいとき, 再帰的データ構造を意識すれば,次のようにスタックを実装できると思います.
struct Stack<T> {
head: Option<Box<Node<T>>>
}
struct Node<T> {
elem: T,
next: Option<Box<Node<T>>>
}
impl<T> Stack<T> {
pub fn new() -> Self {
Stack { head: None }
}
pub fn pop(&mut self) -> Option<T> {
// ...省略
}
pub fn push(&mut self, e: T) {
// ...省略
}
}
このStack<T>
はもちろんちゃんと次のテストが通ります.
(そのようにpush, popを実装します)
let mut stack = Stack::new();
assert_eq!(stack.pop(), None);
stack.push(123);
assert_eq!(stack.pop(), Some(123));
stack.push(1);
stack.push(2);
stack.push(3);
assert_eq!(stack.pop(), Some(3));
assert_eq!(stack.pop(), Some(2));
assert_eq!(stack.pop(), Some(1));
assert_eq!(stack.pop(), None);
しかし,次のコードはfatal runtime error: stack overflow
でプログラムが異常終了してしまいます.
let mut stack = Stack::new();
for i in 0..20_000 {
stack.push(i)
}
たったこれだけのコードで,このStack
の何がいけないのでしょうか?
結論からいえば,タイトルにもある通りデストラクタが原因です.
この場合はstack
のデストラクタの呼び出しを皮切りに,再帰的にNode
のデストラクタが呼び出されることでスタックオーバーフローになります.
もう少し噛み砕くと,まずstack
のデストラクタが行われますが,stack
のデストラクタの内部ではstack.head
のデストラクタが行われます.
同様にstack.head
のデストラクタではstack.head.elem
とstack.head.next
がデストラクタが実行されますが,stack.head.next.next
も...というように再帰的にデストラクタが実行されるのです.
この問題はStack
に数個要素を追加するだけでは発生しませんが,しかし2万個程度要素があるだけでデストラクタはスタック領域を食いつぶしてしまうのです.
僕がこの問題に直面したとき,なぜスタックオーバーフローとなっているのか全く理解が出来ませんでした.
この問題に気づきにくい原因は2つあると思います.
特に2つ目が主な原因だと思います. C++やD言語でもそうですが,構造体のデストラクタはコンパイラが自動的に生成してくれるので,プログラマがわざわざ自分で書く必要がない場合が多いと思います.
しかし,再帰的なデータ構造ではこのデフォルトのデストラクタが悪さをしてプログラムをスタックオーバーフローへと導きます.
この問題はRustのみの問題ではなく,C++でも同様に起こります.
template <typename T>
struct StackNode {
StackNode(T e, std::unique_ptr<StackNode<T>> && n)
: elem(e), next(std::move(n)) {}
T elem;
std::unique_ptr<StackNode<T>> next;
};
template <typename T>
struct Stack {
std::unique_ptr<StackNode<T>> head;
void push(T elem) { /* ...*/ }
};
int main() {
Stack<int> stack{};
for(int i = 0; i < 1000000; ++i) {
stack.push(i);
}
}
自分でちゃんとデストラクタ,つまりRustではDrop
トレイトのdrop
を実装すれば問題は起こりません.
impl<T> Drop for List<T> {
fn drop(&mut self) {
let mut node = mem::replace(&mut self.head, None);
while let Some(x) = node {
node = x.next;
}
}
}
C++でも同じです.
template <typename T>
Stack<T>::~Stack()
{
while(head) {
auto h = std::move(head);
head = std::move(h->next);
}
}