［C++］constexprに双対数（二重数）を実装する - 地面を見下ろす少年の足蹴にされる私

双対数（二重数）

双対数とは、実数に対して $\epsilon^2 = 0$ となるような実数ではない要素を導入し、任意の実数a,bを用いて $a+b\epsilon$ と表わされる数です。見ての通り複素数（ $i^2=-1$ ）と似たもので、性質を考えるときに複素数の考え方が役立ちます。
$\epsilon^2 = 0$ を導入すると何が嬉しいのか？テイラー展開を双対数で見てましょう。
テイラー展開の定義
　 $\displaystyle f(x + \delta) = \sum_{n=0}^{\infty} \frac{f^{(n)}(x)}{n!}\delta^n$
より、微小項 $\delta = \epsilon$ と見て、 $x + \delta$ を双対数 $a+b\epsilon$ に置き換えると
　 $\displaystyle f(a + b\epsilon) = f(a) + f'(a)b\epsilon + \sum_{n=2}^{\infty} \frac{f^{(n)}(a)}{n!}(b\epsilon)^n$
となります。ここで $\epsilon^2 = 0$ を思い出すと、 $n\geqq2$ の項は0になり消えてしまうため
　 $\displaystyle f(a + b\epsilon) = f(a) + f'(a)b\epsilon$
となり、 $\epsilon$ の項の係数が元の関数の微分になっている事が分かります。つまり、任意の関数を双対数で計算すれば、自動的にその関数の微分を得ることが出来ます。また、多項式環 $R[X$ ]を定義多項式 $X^2=0$ で割った剰余環（ $R[X]/X^2$ ）として双対数を定義することが出来、数式・数値的ではなく代数的に微分を計算できているようです。
ちなみに、自動微分界隈では双対数は二重数と呼ばれています（たぶんdual numberの直訳）。

双対数で微分を求めることについての少し詳しい解説、すこし
なぜ双対数（二重数）で微分を求められるのか・・・？ - 地面を見下ろす少年の足蹴にされる私

双対数の性質と実装

双対数の性質を考えながら実装していきましょう。
まずはクラスとしてのガワを整えておきます。ほとんどすべての式はconstexprにできるはずなので、基本的にconstexprを付加しておきます。

template<typename T>
struct dual {
	using this_type= dual<T>;
	using value_type = T;

	//デフォルトコンストラクタ
	constexpr dual()
		: m_a{ 0.0 }
		, m_b{ 0.0 }
	{}

	//値を入れて初期化
	constexpr dual(value_type a, value_type b = value_type{ 0.0 })
		: m_a{ a }
		, m_b{ b }
	{}

	//実部を取得
	constexpr value_type a() const {
		return m_a;
	}

	//虚部？を取得
	constexpr value_type b() const {
		return m_b;
	}
private:
	value_type m_a;
	value_type m_b;
};

ここに各種演算等を入れていく事にします。

四則演算

まずは基本の四則演算。冒頭書いたように複素数での計算を考えることで導出できます。

//加法
constexpr this_type& operator+=(const this_type& rhs) {
	m_a += rhs.m_a;
	m_b += rhs.m_b;

	return *this;
}

//乗法
constexpr this_type& operator*=(const this_type& rhs) {
	m_b *= rhs.m_a;
	m_b += m_a * rhs.m_b;
	m_a *= rhs.m_a;

	return *this;
}

//減法
constexpr this_type& operator-=(const this_type& rhs) {
	m_a -= rhs.m_a;
	m_b -= rhs.m_b;

	return *this;
}

//除法
constexpr this_type& operator/=(const this_type& rhs) {
	m_b *= rhs.m_a;
	m_b -= m_a * rhs.m_b;
	m_b /= rhs.m_a * rhs.m_a;
	m_a /= rhs.m_a;

	return *this;
}

足し算と引き算に関しては悩むところはないでしょう。乗除は $\epsilon^2 = 0$ のために複素数とは少し異なります。
・乗算
$\begin{eqnarray} (a+b\epsilon)(c+d\epsilon) &=&ac+ad\epsilon+bc\epsilon+bd\epsilon^2\\ &=&ac+(ad+bc)\epsilon \end{eqnarray}$
・除算
$\begin{eqnarray} \frac{(a+b\epsilon)}{(c+d\epsilon)} &=&\frac{(a+b\epsilon)}{(c+d\epsilon)} \times \frac{(c-d\epsilon)}{(c-d\epsilon)} \\ &=&\frac{ac-ad\epsilon+bc\epsilon-bd\epsilon^2}{c^2-cd\epsilon+cd\epsilon-(d\epsilon)^2}\\ &=&\frac{ac+(-ad+bc)\epsilon}{c^2}\\ &=&\frac{a}{c}+\frac{-ad+bc}{c^2}\epsilon \end{eqnarray}$

ここ、それぞれの虚部をよく見てみると、積の微分法則、商の微分法則、を満たしている事が分かります。これが双対数で微分計算が出来る最大のキモです。

そして、実数との四則演算も定義します。実数は $a+0\epsilon$ という双対数だと考えられるので、これを先ほどの四則演算に当てはめます。

//加法
constexpr this_type& operator+=(const value_type rhs) {
	m_a += rhs;

	return *this;
}

//乗法
constexpr this_type& operator*=(const value_type rhs) {
	m_b *= rhs;
	m_a *= rhs;

	return *this;
}

//減法
constexpr this_type& operator-=(const value_type rhs) {
	m_a -= rhs;

	return *this;
}

//除法
constexpr this_type& operator/=(const ValueType rhs) {
	m_a /= rhs;
	m_b /= rhs;

	return *this;
}

これもやはり加減算に関してはそのままです。乗除は、先ほどの導出時の式のd=0と考えると結局このように直感的な形になります。
クラスの内部で定義するのは自己代入系の演算子のみで、通常の二項演算はクラス外に定義しておきます。その際、先ほど定義した自己代入系の演算子を利用します。

//双対数同士の加算
template<typename T>
constexpr auto operator+(const dual<T>& lhs, const dual<T>& rhs) {
	return dual<T>{lhs} += rhs;
}

//双対数+実数
template<typename T>
constexpr auto operator+(const dual<T>& lhs, const T rhs) {
	return dual<T>{lhs} += rhs;
}

//実数+双対数
template<typename T>
constexpr auto operator+(const T lhs, const dual<T>& rhs) {
	return dual<T>{rhs} += lhs;
}

多いので加算だけ書いておきますが、他の二項演算子も同じように定義できます。

逆元

逆元とはある元との二項演算の結果がその単位元となるような元の事です。二項演算ごとに存在し、存在するなら1つしかありません。実数 $x$ ならば、足し算の $-x$ 、掛け算の $\frac{1}{x}$ がそれにあたります。
双対数の加算の単位元は $0+0\epsilon$ なので、加算の逆元は $a+b\epsilon$ に対して $-a-b\epsilon$ です。
乗算の逆元は少し面倒なので、計算してみましょう。 $d, d^{-1}$ を任意の双対数とその逆元、 $\it{1}$ を双対数の乗法単位元とします。
$\begin{eqnarray} d d^{-1} &=& \it{1}\\ d^{-1}&=&\frac{\it{1}}{d}\\ &=&\frac{\it{1}}{a+b\epsilon}\frac{a-b\epsilon}{a-b\epsilon}\\ &=&\frac{a-b\epsilon}{(a+b\epsilon)(a-b\epsilon)}\\ &=&\frac{a-b\epsilon}{a^2-(b\epsilon)^2}\\ &=&\frac{a-b\epsilon}{a^2}\\ &=&\frac{1}{a}-\frac{b}{a^2}\epsilon \end{eqnarray}$
と求まります。明らかに $a=0$ の時に存在しない事が分かります。

//加法逆元
constexpr this_type operator-() const {
	return type{ -m_a,-m_b };
}

//乗法逆元にする
constexpr void inverse() {
	this->inverted();
}

//乗法逆元を得る
constexpr this_type& inverted() {
	//d^-1 = 1/a - b/a^2
	m_a = T(1.0) / m_a;
	m_b /= -(m_a* m_a);
	return *this;
}

これらの関数は現在の値についてのものなのでメンバで定義。加法逆元に関しては単項-演算子があるのでそれを利用し、乗法逆元は得たい場合としたい場合が想定されるので二つ準備。

関係演算子

まずは取り合えず等値比較。これは何も考えず項毎に比較した結果の論理積を取ればいいでしょう。複素数もそうなっていますし。

//等値比較
constexpr bool operator==(const type& rhs) const {
	return std::tie(m_a, m_b) == std::tie(rhs.m_a, rhs.m_b);
}

//等値比較の否定、メンバ内部の等値比較演算子を利用して実装しクラス外に定義
template<typename T>
constexpr bool operator!=(const dual<T>& lhs, const dual<T>& rhs) {
	return !(*this == rhs);
}

このように、std::tieを使うとこういう複数の値の比較を簡単に書くことが出来ます（2つ程度ならあまり旨みはないけれども・・・）。

全順序となるような大小関係については複素数同様に定義することが出来ませんが、辞書式順序（実部を比較の後虚部を比較、するような比較）を導入すれば狭義の弱順序を満たし、C++の要求的には並べることができます。おそらく実用的にはそれで充分ですので辞書式順序で実装しておきます。

//operator<、クラス内部で定義
constexpr bool operator<(const type& rhs) const {
	return std::tie(m_a, m_b) < std::tie(rhs.m_a, rhs.m_b);
}


//残りの関係演算子、クラス外部で定義
template<typename T>
constexpr bool operator<=(const dual<T>& lhs, const dual<T>& rhs) {
	return (lhs == rhs) || (lhs < rhs);
}

template<typename T>
constexpr bool operator>(const dual<T>& lhs, const dual<T>& rhs) {
	return rhs < lhs;
}

template<typename T>
constexpr bool operator>=(const dual<T>& lhs, const dual<T>& rhs) {
	return (lhs == rhs) || (lhs > rhs);
}

ここでもstd::tieを使います。こっちの場合は複雑な入れ子比較を書かなくてよくなるのでより旨みが増します。
そして、クラス内に定義した二つの比較演算子を用いて残りの演算子をクラス外に定義しておきます。

cmathの各種数学関数

三角関数に代表される各関数も双対数に対して定義しておくと便利です（というかないと使い物にならない）。
実部はそのままですが、虚部はどう定義すればよいのでしょうか？
ある関数に双対数を通すと、実部にその値が得られ、虚部にはその微分が得られます。つまりは、関数の導関数を求め、入力の実部をそこに通したものをもとの虚部の値にかけてやれば無事に微分を得ることが出来ます。

template<typename T>
auto sin(const dual<T>& d) {
	using std::sin;
	using std::cos;

	return dual<T>{sin(d.a()), d.b()*cos(d.a())};
}

template<typename T>
auto cos(const dual<T>& d) {
	using std::sin;
	using std::cos;

	return dual<T>{cos(d.a()), -d.b()*sin(d.a())};
}

//二変数関数の導関数はその全微分を求める
template<typename T>
auto atan2(const dual<T>& y, const dual<T>& x) {
	using std::atan2;

	auto sumsq_inv = T(1.0) / (x.a() * x.a() + y.a() * y.a());
	return dual<T>{atan2(y.a(), x.a()), sumsq_inv*(-y.a()*x.b() + x.a()*y.b())};
}

sinとcosを示しておきます。SLTのsin,cosをusingしているのは、Tにプリミティブでない型が与えられた時を考慮しての事です。各関数の導関数はググれば出てくるので、他の関数の定義も同じようにすればいいでしょう。
二変数関数となる一部の関数（atanやpow）はその全微分を求めてやります。
がんばってベッセル関数も定義したけど、clangでは使えない様子。残念・・・

C++的事情による残りの演算子

//実部だけを取り出す暗黙変換、explicit付けた方がいいのかしら
constexpr operator value_type() const {
	return m_a;
}

//単項-の対になる物、必要性・・・？
constexpr this_type operator+() const {
	return *this;
}

//インクリメント演算子4つ、後置はクラス外で定義
constexpr this_type& operator++() {
	++m_a;
	return *this;
}

constexpr this_type& operator--() {
	--m_a;
	return *this;
}

template<typename T>
constexpr auto operator++(dual<T>& dual, int) {
	auto copy = dual;
	++dual;
	return copy;
}

template<typename T>
constexpr auto operator--(dual<T>& dual, int) {
	auto copy = dual;
	--dual;
	return copy;
}

//ストリーム出力演算子、これはクラス外で定義、出力形式はお好み
template<typename T>
std::ostream& operator<<(std::ostream& ostream, const dual<T>& rhs) {
	ostream << rhs.a() << " + " << rhs.b() << "e";
	return ostream;
}

//ユーザー定義リテラル、サフィックスはお好みで
constexpr dual<double> operator"" _d(long double real) {
	return {static_cast<double>(real), 0.0};
}

constexpr dual<double> operator"" _eps(long double eps) {
	return {0.0, static_cast<double>(eps)};
}

//略記のためのエイリアス、これもクラス外
using dual_f = dual<float>;
using dual_d = dual<double>;
using dual_ld = dual<long double>;

なるべくプリミティブな浮動小数点型と同じようにあつかえることを目指しました。
最後に省略した部分も含めた実装を置いておきます。ヘッダオンリーかつSTLがあれば使えます。C++14以降に対応したコンパイラが必要になるでしょう。
github.com

本当に動くの？？

双対数の定義を見てもなぜあれで微分が出来るのかさっぱりな上、天下り的に実装されても本当に動くのかは疑問なところ。
実際に動かして実験してみます。

以下のように、簡単な多項式とその微分、sinとその微分を返す関数を用意して、双対数と普通の数を渡して結果を見てみます。

//f(x) = 4x^3 + 3x^2 + 2x + 1
template<typename T>
auto f(const T& x) {
    return 4.0*x*x*x + 3.0*x*x + 2.0*x + 1.0;
}

//df(x)/dx = 12x^2 + 6x + 2
template<typename T>
auto df(const T& x) {
    return 12.0*x*x + 6.0*x + 2.0;
}

//sin(x)
template<typename T>
auto fsin(const T& x) {
    using std::sin;
    return sin(x);
}

//dsin(x)/dx = cos(x)
template<typename T>
auto dfsin(const T& x) {
    using std::cos;
    return cos(x);
}

実行結果
[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

f(x) = 10 + 20e
df(x)/dx = 20
sin(π) = 1.22465e-16 + -1e
dsin(π)/dx = -1

双対数の出力は2項目を見てやりますと、確かに結果が一致している事が分かります。どうやら微分を求められているようですね。

ニュートン法

では早速、実装した双対数を使ってニュートン法をやってみましょう。分かりやすいので2乗根を求めるのにニュートン法を使ってみます。
2乗根を求めるには次のような二次関数の零点となる $x$ を求めてやればいいです。
$f(x)=x^2-n$
ちなみに、2乗根を求めるだけならもっといいアルゴリズムがあるのでそちらをお使いください。これはあくまでサンプルです。

template<typename T>
constexpr auto fabs_static(T x) {
	return (x < 0.0) ? (-x) : (x);
}

template<typename Func>
constexpr double newtonMethod(double x0, Func&& f) {
	using namespace DualNumbers;

	double xn = x0;
	double diff{};

	do {
		dual<double> d = f({ xn, 1.0 });
		diff = d.a() / d.b();
		xn -= diff;
	} while (1.0E-15 < fabs_static(diff));

	return xn;
}

まずは一般的なニュートン法の実装部分。初期値と関数（f:dual<double> → dual<double>）を受け取りオーソドックスにニュートン法を実行します。std::fabsがconstexprではないので、constexprなfabsを定義しておきます。

struct squareroot{
	DualNumbers::dual<double> N;

	constexpr squareroot(double n)
		: N{ n, 1.0 }
	{}

	constexpr DualNumbers::dual<double> operator()(const DualNumbers::dual<double>& x) const {
		return x * x - N;
	}
};

constexpr double sqrt_newton(double n) {
	using namespace DualNumbers;

	//return newtonMethod(n, [N = dual<double>{ n, 1.0 }](const dual<double>& x) constexpr {return x * x - N; });
	return newtonMethod(n, squareroot{ n });
}

次にn乗根を求めるための関数を計算するファンクタを定義し、それと受け取った初期値からニュートン法を実行するインターフェースとなる関数を定義します。なぜラムダで渡さないのかというと、VS2017(15.8.9)ではエラーが出てconstexpr実行できないからです（少なくともclangはコメントアウトしてある方をコンパイル可能）。

int main() {
	constexpr auto sqrt_n = sqrt_newton(10.0);

	std::cout << std::setprecision(16) << sqrt_n << std::endl;   //3.162277660168379
}

最後に呼び出すところです。constexprで受けてもエラーが出ないことからconstexpr実行できている事が分かるでしょう。
wandboxでのコードと実行結果を置いておきます。
[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

この様に、双対数を用いれば数値微分に頼らずとも任意の関数の微分を計算することが出来ます。しかも数値微分のような刻みの概念を気にしなくてもよくなります、凄い。
（ろくにテストしてないのでバグを見つけたらご報告ください・・・）

参考文献

二重数 - Wikipedia
微分環と双対数 - 記号の世界ゟ
 自動微分を実装して理解する（前編） - Qiita
atan2 - Wikipedia