2018-12-23

なぜ双対数（二重数）で微分を求められるのか・・・？

双対数二重数数学自動微分

前回記事の冒頭でサラッと流した部分の補完。
［C++］constexprに双対数（二重数）を実装する - 地面を見下ろす少年の足蹴にされる私

$\epsilon^2 = 0$ を導入するだけで、なぜに2項目に微分が計算されるのでしょうか。純に数学的には微分環との関連が云々・・・、と日本語でおｋな感じになるので、数学弱めな自分なりに考えてみました。

双対数の四則演算

後で使うので、双対数の四則演算をここにまとめておきます。二つの双対数 $d1 = a+b\epsilon$ と $d2 = c+d\epsilon$ について

$\begin{eqnarray} d1+d2&=& (a+b\epsilon) + (c+d\epsilon)\\ &=&(a+c) + (b+d)\epsilon \end{eqnarray}$

$\begin{eqnarray} d1 - d2 &=& (a+b\epsilon) - (c+d\epsilon)\\ &=&(a-c) + (b-d)\epsilon \end{eqnarray}$

$\begin{eqnarray} d1 \times d2 &=& (a+b\epsilon) \times (c+d\epsilon)\\ &=&ac+ad\epsilon+bc\epsilon+bd\epsilon^2\\ &=&ac+(ad+bc)\epsilon \end{eqnarray}$

$\begin{eqnarray} d1 \div d2 &=&\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}$

双対数の虚部

虚部と呼ぶのも変ですが、双対数の2項目（ $a+b\epsilon$ の $b\epsilon$ ）をそう呼ぶことにします。
双対数で微分を求めていくためには、この2項目の係数bをまず1にしておく必要があります。なぜでしょうか？計算していけば分かるのですが、1にしておかないと微分が正しく求まりません。
では、これは何を示しているのでしょうか？双対数の一般形（ $d = a+b\epsilon$ ）を少し書き換えてみます、
$\begin{eqnarray} d = x+1\epsilon \end{eqnarray}$
aには任意の実数を入れられるので、そこを変数xと見ます。
つまり、実部を任意の変数xと置いたときの微分値が虚部です。実際に使う時はaに具体的な値を入れるためにイメージしにくいですが、 $\frac{d}{dx}x = 1$ なのでb=1を入れておくわけです。
この様に見れば、任意の実数aは双対数で
$a = a+0\epsilon$
と表わされることも違和感がないでしょう（ $\frac{d}{dx}a = 0$ なので）。

加減算

次に足し算をやってみましょう。
$\begin{eqnarray} (x + 1\epsilon) + (x + 1\epsilon)&=&(x+x) + (1+1)\epsilon\\ &=&2x + 2\epsilon \end{eqnarray}$
当たり前のように、虚部は2になります。 $\frac{d}{dx}2x = 2$ であるので、足し算しても虚部は実部に対して微分を表し続けます。引き算をしても
$\begin{eqnarray} (2x + 2\epsilon) - (x + 1\epsilon)&=&(2x-x) + (2-1)\epsilon\\ &=&x + 1\epsilon \end{eqnarray}$
のように、引き算に対しても虚部は実部に対して微分となります。

また実数aとの計算でも
$\begin{eqnarray} (a + 0\epsilon) - (x + 1\epsilon)&=&(a-x) + (0-1)\epsilon\\ &=&(a-x) - 1\epsilon \end{eqnarray}$
同じように虚部は実部の関数に対して微分を示し続けます（ $\frac{d}{dx}(a-x) = -1$ ）。

乗除算

次は掛け算をやってみます。
$\begin{eqnarray} (x + 1\epsilon) \times (x + 1\epsilon)&=&(x\times x) + (x\times 1+1\times x)\epsilon\\ &=&x^2 + (x+x)\epsilon\\ &=&x^2 + 2x\epsilon \end{eqnarray}$

$\frac{d}{dx}(x2) = 2x$ なので、積に関しても虚部は微分を示しています。
ここで、今までは双対数の実部を変数と置いてきましたが、一歩進んで関数と見てそれを $f,g$ とおき、その微分（虚部）を $f',g'$ とおいて双対数の積を考えると

$\begin{eqnarray} (f + f'\epsilon) \times (g + g'\epsilon)&=&(f\times g) + (f\times g'+g\times f')\epsilon\\ &=&fg + (fg'+gf')\epsilon \end{eqnarray}$

実部はそのまま関数の積になります。
そして虚部は、積の微分法則（ライプニッツ則）となっています。

次に割り算です。
$\begin{eqnarray} (x^2 + 2x\epsilon) \div (x + 1\epsilon)&=&\frac{x^2}{x} + \left( \frac{-x^2\times1 + 2x\times x }{x\times x} \right)\epsilon\\ &=&x + \left( \frac{-x^2 + 2x^2 }{x^2} \right)\epsilon\\ &=&x + \left( \frac{x^2}{x^2} \right)\epsilon\\ &=&x + 1\epsilon \end{eqnarray}$

割り算に関しても、虚部は実部の微分になる関係を維持しています。
もうお気づきかもしれませんが、先ほどと同じように関数に置き換えて計算してみると

$\begin{eqnarray} (f + f'\epsilon) \div (g + g'\epsilon)&=&\frac{f}{g} + \left( \frac{-f\times g' + g\times f' }{g\times g} \right)\epsilon\\ &=&\frac{f}{g} + \left( \frac{gf'-fg'}{g^2} \right)\epsilon \end{eqnarray}$

この様に、虚部は商の微分法則を構成しています。

関数係数双対数の四則演算

任意の関数を係数とした双対数の加減算についても求めて、まとめておきます。二つの双対数 $d1 = f+f'\epsilon$ と $d2 = g+g'\epsilon$ について

$\begin{eqnarray} d1+d2&=& (f+f'\epsilon) + (g+g'\epsilon)\\ &=&(f+g) + (f'+g')\epsilon \end{eqnarray}$

$\begin{eqnarray} d1 - d2 &=& (f+f'\epsilon) - (g+g'\epsilon)\\ &=&(f-g) + (f'-g')\epsilon \end{eqnarray}$

$\begin{eqnarray} d1 \times d2 &=& (f + f'\epsilon) \times (g + g'\epsilon)\\ &=&fg + (fg'+gf')\epsilon \end{eqnarray}$

$\begin{eqnarray} d1 \div d2 &=&(f + f'\epsilon) \div (g + g'\epsilon)\\ &=&\frac{f}{g} + \left( \frac{gf'-fg'}{g^2} \right)\epsilon \end{eqnarray}$

$\epsilon^2 = 0$

ここまで見てきたように、双対数の虚部は常に実部の関数に対しての導関数となっています。少なくとも四則演算の範囲ではその関係は崩れません。
双対数の加減算は複素数と同じ形になっているので、加減算が虚部と実部の微分についての関係を保存する性質は複素数でも成り立ちます。しかし、乗除算についてはどうでしょうか？
複素数では $i^2 = -1$ となるため $i^2$ の項は再び実数に戻ります。そのため、積に関しては実部が、商に関しては両方が、虚部と実部の微分についての関係を保存しません。
実数に付け加える任意の数を $\chi$ と置いて乗除を計算してみると（ $\chi$ を付け加えて実数を拡大した数をd1とd2とします）
$\begin{eqnarray} d1 \times d2 &=& (a+b\chi) \times (c+d\chi)\\ &=&ac+ad\chi+bc\chi+bd\chi^2\\ &=&ac+(ad+bc)\chi+bd\chi^2 \end{eqnarray}$

$\begin{eqnarray} d1 \div d2 &=&\frac{(a+b\chi)}{(c+d\chi)}\\ &=&\frac{(a+b\chi)}{(c+d\chi)} \times \frac{(c-d\chi)}{(c-d\chi)} \\ &=&\frac{ac-ad\chi+bc\chi-bd\chi^2}{c^2-cd\chi+cd\chi-(d\chi)^2}\\ &=&\frac{ac-(ad+bc)\chi-bd\chi^2}{c^2-(d\chi)^2} \end{eqnarray}$

ここで、 $\chi^2 = -1$ とすれば複素数となり、 $\chi^2 = 0$ とすると双対数になります。積と商の微分法則を虚部で成り立たせつつ実部では素のままの計算となるようにするためには $\chi^2 = \epsilon^2 = 0$ が必要で、そうでなければいけない事が分かるでしょう。

複素数が $i^2 = -1$ の導入によって二次元ベクトル空間の元と同値になるように、双対数は $\epsilon^2 = 0$ の導入によって実部と虚部の微分についての関係を四則演算で保存するようになります。その結果、実部と虚部が関数と導関数の関係を保持する双対数同士の計算において実部に計算結果が、虚部にその微分値が、自動で計算されます。
双対数で微分を計算するためには $\epsilon^2 = 0$ だけでは不十分で、入力となる双対数が実部と虚部で微分関係を維持している（例：実部 $a$ に対して虚部 $b=1$ ）ことも必要になります。虚部係数bに1か0以外の値を入れる場合、実部がその値に対応した値でないと微分は正しく求まらず、実部も虚部もでたらめな値が計算されます。

つまりは、双対数は微分が自動計算できる魔法の数というわけではなく、双対数の性質によって入力双対数の実部と虚部の微分関係が保持されたまま四則演算を実行可能で、結果として微分を得ることが出来るわけです。

参考文献

積の微分法則 - Wikipedia
商の微分法則 - Wikipedia
複素数の四則演算

2018-12-22

［C++］constexprに双対数（二重数）を実装する

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

2018-12-14

［C++］配列の参照、その書き方

C++ 配列参照

C++ではおよそ全ての型に対して、その実体へのエイリアス（別名）として参照を宣言できます。少しだけ安全なポインタとして扱うことが出来、コンパイル後はポインタと同じコードになります（ならない処理系もあるかもしれません）。そして、配列もまたその参照を書くことが出来ます。

変数宣言

int array[]{ 1,2,3,4,5 };

//配列の参照宣言
int(&r2)[5] = array;

//decltype(r1) = int(&)[5]
auto& r1 = array;

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

この様に、配列の参照を変数として書くときは
T(&name)[N]のように書き、Tに型名、nameに変数名、Nに要素数、となります。
ただし、配列の型Tと要素数Nは参照先の配列と一致している必要があります（Nはコンパイル時に確定する必要があります）。
しかし、これは書き辛いし見辛いので、auto&を使ってコンパイラさんに推論してもらうのをお勧めします。

引数宣言

//int型の5要素の配列を受け取る
void f(int(&array)[5]) {
    for(auto n : array) {
        std::cout << n << std::endl;
    }
}

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

関数の引数として配列を取るときは、変数宣言時と同じ書き方で書いてあげます。見辛いのでテンプレート使いたくなりますが、その場合はポインタに推論されるため、配列への参照にはなりません（ほぼ同じ意味にはなりますが）。

配列への参照を返す関数

一体いつこれを使うことになるのかは定かではありませんが、もちろん書くことが出来ます。しかし、奇妙さが倍増しになります・・・

//int型5要素の配列への参照を返すretarray関数
int (&retarray())[5] {
    static int array[]{ 1,2,3,4,5 };
    return array;
}

//cv、参照修飾、noexcept指定などは引数宣言の閉じかっこ後
int (&retarray2(int n) noexcept)[5] {
    static int array[]{ n,n+1,n+2,n+3,n+4 };
    return array;
}

//auto&を使って簡便に
auto& retarray3() {
    static int array[]{ 1,2,3,4,5 };
    return array;
}

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

この様に、配列の参照を返す関数を書くときは
T(&name(args...) const noexcept...)[N] {...}のように書き、Tに型名、nameに関数名、args...に引数列、Nに要素数、となり、変数宣言の時と同じように、TとNは参照先のものと一致していなければなりません。
CV、参照修飾や例外指定は普通の関数と同じ位置、つまりは引数宣言の閉じかっこの後に置きます。
こちらもauto&を使えるのでこちらで書くと少し見やすくなります（が、型が分かり辛くなります）。

しかし、戻り値型の後置き構文を使う事でかなり見やすく書くことが出来ます。

auto retarray() -> int(&)[5] {
    static int array[]{ 1,2,3,4,5 };
    return array;
}

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

後置きの書式は変数宣言時と同じ書き方です（名前はいりません）。これが一番見やすいのでこの書き方で書くといいでしょう。
関係ない注意点ですが、関数ローカルの配列への参照を返してしまわないように気を付けてください。

配列の参照のメリット

一見すると構文の奇妙さが目立ち、使いどころも分からない配列への参照ですが、その要素数をコンパイル時に確定できるというメリットがあります。
どういう事でしょうか？

void func(int array[5]);
//decltype(func) = void (int *)

int array[2]{};
func(array);

関数の引数に配列を取るときはこのように書き、呼び出し側では配列の名前を渡します。引数宣言では要素数を書かなくてもいいのにあえて書くことで配列の要素数に制限を加えている気分になります。
なりますが、実際はこれは単なるポインタ渡しになっており、配列でなくてもint*となる物なら何でも渡すことが出来てしまいます。
[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

しかし、配列の参照を受け取るようにすると

void func(int(&array)[5]);
//decltype(func) = void (int (&)[5])

int array[2]{};
func(array);   //コンパイルエラー

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ
引数の型名に要素数がしっかりと入り、その要素数以外の配列を渡そうとするとコンパイルエラーとなります。もちろん、int*を渡すこともできません。しかも、関数の中でその要素数を超えてアクセスしている場合に警告を出してくれます（clangの場合）。
ある意味Conceptの様な制約を加えることが出来ているわけです。しかも、参照は実質ポインタなので実行時にはポインタ渡しと同じコードになっているはずです。

また、これを利用して静的配列の要素数を得るテクニックがあります。

template<typename T, size_t N>
constexpr auto size(const T(&)[N]) noexcept -> size_t {
   return N;
}

int array[5] = {5, 4, 3, 2, 1};
auto N = size(array); // N = 5

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ
 C言語由来のマクロで書く方法もありますが、C++ならばこちらを使うべきでしょう。ちなみにC++17でこれはstd::size関数として標準ライブラリ入りしています（iteratorヘッダ）。

多次元配列の場合

配列は別に1次元だけでなく、多次元配列を宣言できます。そして多次元配列の参照も当然あります。さらに複雑な構文を要求されることはなく、普通に次元を書き足してやるだけです。

//int型3行4列の配列への参照を返す、retarray関数
int(&retarray() noexcept)[3][4] {
    static int array[3][4] = {{4, 3, 2, 1}, {8, 7, 6, 5}, {12, 11, 10, 9}};
    return array;
}

//int型3行4列の配列への参照を宣言
int(&array2)[3][4] = retarray();

//二次元配列の要素数を取得するsize関数
template<typename T, size_t N, size_t M>
constexpr auto size(const T(&)[N][M]) noexcept -> std::tuple<size_t, size_t> {
    return {N, M};
}

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

参考文献

「要素数 4 の配列を渡してください」ではなく「要素数 4 の配列を渡さなければならない」にする - Qiita
size - cpprefjp C++日本語リファレンス
この記事を書くきっかけ
https://twitter.com/non_archimedean/status/1069039533688115200

2018-11-23

［C++］final指定子と最適化

C++ final 最適化仮想関数テーブル

final specifier

final指定子はC++11にてoverride指定子とともに導入されたもので、指定した関数がオーバーライドされないことを明示し、オーバーライドされた場合にコンパイルエラーを起こすものです。
また、クラスに対しても指定でき、継承できないことを表します。継承した場合はコンパイルエラーとなります。
クラスの名前の後ろ、関数の引数定義の後ろ（overrideと同じところ）に書くことが出来ます。
overrideと違ってあまり利用されていないような気がする子です・・・

仮想関数テーブル参照のスキップ

このfinal指定子ですが、指定しておくとある条件の下で最適化に利用される可能性があります。
以下のようなコードで実験してみます。

struct Base {
    virtual int Num() = 0;
};

struct Derived1 : Base {
    int Num() override {
        return 10;
    }
};

struct Derived2 : Base {
    int Num() override {
        return 100;
    }
};

auto output(Base* ptr) {
    return ptr->Num();
}

auto output(Derived2* ptr) {
    return ptr->Num();
}

このコードをコンパイルしてアセンブリ出力を見てみます。一番見やすかったのでclang7.0.0のものを見てみます（両output関数の部分のみコピペします）。

output(Base*): # @output(Base*)
  mov rax, qword ptr [rdi]
  jmp qword ptr [rax] # TAILCALL
output(Derived2*): # @output(Derived2*)
  mov rax, qword ptr [rdi]
  jmp qword ptr [rax] # TAILCALL

アセンブラが読めなくても、なんとなくポインタでジャンプしてんなあというのが分かると思います。
私も詳しくはないですが簡単に説明すると
mov rax, QWORD PTR [rdi]→rdiの値を64bitポインタとしてraxにコピー
jmp qword ptr [rax]→raxの値を64bitポインタとしてジャンプ
両方ともこんな感じで仮想関数テーブルを参照しているらしい事が分かります。

次にDerived2を修正します。

struct Derived2 : Base {
    int Num() override final {
        return 100;
    }
};

このようにDerived2のNum()にfinalをつけておきます、すると

output(Base*): # @output(Base*)
  mov rax, qword ptr [rdi]
  jmp qword ptr [rax] # TAILCALL
output(Derived2*): # @output(Derived2*)
  mov eax, 100
  ret

お分かりいただけたでしょうか？
Derived2のポインタからNum()を呼び出す方は、仮想関数テーブルの参照がスキップされ、関数がインライン展開されています。
何が起きたのかというと
final指定を追加したことによって、コンパイラはDerived2::Num()がこれ以上オーバーライドされないことを知っているため、Derived2のポインタからの関数呼び出しに関しては仮想関数テーブルを参照して関数呼び出しをする必要がない、と判断することが出来ます。
結果、Derived2*からのNum()の呼び出し先は静的に確定し、インライン展開することが出来ます。
ちなみに、Derived2クラス自体にfinal指定をして、オーバーライドした関数にfinalを付けなくても結果は同じです。これ以上継承されないという事は関数がオーバーライドされることはありません。

その他のコンパイラ

clangだけではあれなのでgccやmsvcでも結果を見てみましょう。
Compiler Explorer
上記のclangの出力も含めてこちらのサイトを利用しました。

GCC 8.1 x64

finalなし（10行目付近）

output(Derived2*):
  mov rax, QWORD PTR [rdi]
  mov rax, QWORD PTR [rax]
  cmp rax, OFFSET FLAT:Derived2::Num()
  jne .L7
  mov eax, 100
  ret

finalあり（7行目付近）

output(Derived2*):
  mov eax, 100
  ret

MSVC 19.14 x64

finalなし（95行目付近）

int output(Derived2 * __ptr64) PROC ; output, COMDAT
  mov rax, QWORD PTR [rcx]
  rex_jmp QWORD PTR [rax]
int output(Derived2 * __ptr64) ENDP ; output

finalあり（95行目付近）

int output(Derived2 * __ptr64) PROC ; output, COMDAT
  mov eax, 100 ; 00000064H
  ret 0
int output(Derived2 * __ptr64) ENDP ; output

icc 19.0.0 x64

finalなし（43行目付近）

output(Derived2*):
  mov rax, QWORD PTR [rdi] #22.12
  mov rdx, QWORD PTR [rax] #22.12
  jmp rdx #22.12

finalあり（36行目付近）

output(Derived2*):
  mov eax, 100 #22.17
  ret #22.17

どのコンパイラでも同じようにfinalを指定すると仮想関数テーブル参照がスキップされている事が分かるでしょう。
特筆すべきはfinalなしの時のGCCの出力で、渡されたアドレスがDerived2のものであるかを調べて、そうであれば仮想関数テーブル参照をスキップするコードを出力しています。なんという貪欲さでしょう・・・

まとめ

以上のように、final指定子を追加しておくだけでC++使いの悩みの種？である仮想関数テーブル参照をスキップする最適化を促進できる場合があります。
その条件とは以下の二つ

あるクラスが仮想関数テーブルを持っていて（virtual指定された関数を持つ別クラスを継承している）、かつ
そのクラスのポインタからメンバ関数を呼び出す時（動的ポリモーフィズムしないとき）

派生する予定がない、してほしくない、するつもりがない、場合には積極的にfinalを指定してやるとよいかと思われます。その場合、関数毎につけるよりはクラスに直接つけた方が楽でしょう。

参考文献

overrideとfinal - cpprefjp C++日本語リファレンス
 c++11 - Is final used for optimization in C++? - Stack Overflow

2018-11-14

［C#］Visual Studio 2015でC#7.0以降を利用する

C# nuget

Visual Studio 2019のプレビュー版も出そうな時にまだ2015なんて使ってるんですか？？奇遇ですね。

Visual Studio 2015を強いられている間に、C#はもう8.0を伺うところまで進歩しています。タプルやらspanやら便利なものや構文が追加されています。しかしVS2015に同梱のC#コンパイラはC#6.0までのコードしかコンパイルすることが出来ません・・・。
とはいえ、.NetFramework4.7以降同梱のコンパイラ（csc.exe）ならばC#7.0のコードもコンパイルできるはずです。なんとかVS2015のコンパイラを差し替えられないものでしょうか・・・

Microsoft.Net.Compilers

www.nuget.org
こちらのパッケージをインストールすると、その悩みを解決することが出来ます。
これはMicrosoftのRoslynチームがリリースしているもので、このパッケージをインストールしたプロジェクトをビルドするときに使用するコンパイラをパッケージ内臓のものを使うようにするものです。
まさにVS2015のC#コンパイラを差し替えてくれます。

ちなみに、このパッケージC#だけでなくVBにも同じ効果をもたらします。

しかし、このパッケージはコンパイラを差し替えるだけのため、ターゲットフレームワークによっては実装が無く使えない機能がいくつかあります。

ValueTuple

タプルを扱うにはValueTuple構造体等が必要ですが、これは.Net4.7以降にしか含まれていません。しかし、同名同機能のクラスを自前で用意すれば同じように使うことが出来ます。そして、そのような実装がMicrosoftから提供されています。
www.nuget.org

System.Memory

同じ理由でSpan構造体関連の機能を使うにはSpanの実装が必要です。こちらも提供があります。
www.nuget.org

パッケージバージョンとC#バージョンの対応

パッケージバージョン	C#バージョン
1.x	6.0
2.0	7.0
2.1	7.0
2.2	7.0
2.3	7.1
2.4	7.2
2.6	7.2
2.7	7.2
2.8	7.3
2.9	7.3

バージョンが被っているものはバグ修正です。新しいバージョンが使いたい場合は最新をインストールすればいいと思われます。CIサーバ等にVSインストールせずにコンパイラを最新にする、あえて古いバージョンのコンパイラを使う、コンパイラのバージョンを環境問わず一致させる、等の使い方ではバージョンを見極める必要があるでしょう（どうやらこちらが本来の使い方らしいですが）。

インテリセンス

いいことずくめに思えるMicrosoft.Net.Compilersパッケージですが、インテリセンスには適用されません。もともとインテリセンスに利用するコンパイラは別に用意されているらしいので、単純に置き換えできないのかもしれません。
なので、赤線で怒られまくるけどビルドは通る、という状況で使用することになります・・・

参考文献

NuGet packages · dotnet/roslyn Wiki · GitHub

2018-10-27

［C++］surrogate call function(代理呼び出し関数)なるもの

C++ 関数ポインタオーバーロード

surrogate call functionなる物をたまたま見かけたけども、特に日本語の記事とかなかったので調べてみました。

surrogate call function？？

以下のような動作をする関数オブジェクトのような何かのことです。

template<typename Func1, typename Func2>
struct Surrogate {

    constexpr Surrogate(Func1* f1, Func2* f2)
        : m_f1{f1}
        , m_f2{f2}
    {}
    
    constexpr operator Func1*(){
        return m_f1;
    }
    
    constexpr operator Func2*(){
        return m_f2;
    }
    
private:
    Func1* m_f1;
    Func2* m_f2;
};

void func1(int n) {
    std::cout << "n = " << n << std::endl;
}

void func2(double f) {
    std::cout << std::setprecision(8) << "f = " << f << std::endl;
}

int main()
{
    Surrogate callable{func1, func2};
        
    callable(128);
    callable(1024.2048);
}

//出力
//n = 128
//f = 1024.2048

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

Surrogateクラスは関数呼び出し演算子を一切実装していないはずなのに関数オブジェクトのように扱えているうえに、オーバーロードが正しく処理されています。
何が起きてるのか分かり辛いですが、よく見ると実装されている二つの演算子オーバーロードは暗黙の型変換演算子です。
つまりは、Surrogate型のオブジェクトに対し関数呼び出しがなされた時、そのような演算子が無いので暗黙の型変換を行い、そこから得られた関数ポインタに対して関数呼び出しを試みてくれます。しかもその際に、引数型に応じてオーバーロード解決までしてくれています。コンパイラさん素敵・・・。
分かってしまえば簡単なのですが、絶対に気付かないしやろうとも思わないでしょう・・・。でもこの挙動は少し悪用できそうな気がしないでもないですね。

ちなみに、関数ポインタでなくてラムダ式を渡すこともできます。ただし、キャプチャをしてはいけません。

//省略

int main()
{
    Surrogate blender{+[](int n){func1(n + n);}, +[](double f){func2(f + f);}};
        
    blender(128);
    blender(1024.2048);
}

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

キャプチャをしていないラムダ式は、同じシグネチャ（引数、戻り値の型が同じ）の関数ポインタへ暗黙変換することが出来ます。しかし、この場合テンプレートパラメータの推論時に関数オブジェクトとしての型が推論されてしまいエラーになるので、単項+演算子(decay operator)を頭につけて、明示的に関数ポインタに変換してやります。
では、Surrogateクラスは関数オブジェクトをメンバに持つようにして、暗黙の型変換はそれの参照を返すようにしても上手くいくのでは？と思いますが、こちらはどうやら許されないご様子、関数呼び出し演算子は実装されていないとエラーになります。

合法？

とても怪しいこの挙動ですが、実は規格にちゃんと記載があります。
N4659 §16.3.1.1.2より

16.3.1.1.2 Call to object of class type
1 If the primary-expression E in the function call syntax evaluates to a class object of type “cv T”, then the set of candidate functions includes at least the function call operators of T. The function call operators of T are obtained by ordinary lookup of the name operator() in the context of (E).operator().
2 In addition, for each non-explicit conversion function declared in T of the form
operator conversion-type-id () cv-qualiﬁer ref-qualiﬁeropt noexcept-speciﬁeropt attribute-speciﬁer-seqopt ;
where cv-qualiﬁer is the same cv-qualiﬁcation as, or a greater cv-qualiﬁcation than, cv, and where conversiontype-id denotes the type “pointer to function of (P1,...,Pn) returning R”, or the type “reference to pointer to function of (P1,...,Pn) returning R”, or the type “reference to function of (P1,...,Pn) returning R”, a surrogate call function with the unique name call-function and having the form
R call-function ( conversion-type-id F, P1 a1, ..., Pn an) { return F (a1, ..., an); }
is also considered as a candidate function. Similarly, surrogate call functions are added to the set of candidate functions for each non-explicit conversion function declared in a base class of T provided the function is not hidden within T by another intervening declaration.
3 If such a surrogate call function is selected by overload resolution, the corresponding conversion function will be called to convert E to the appropriate function pointer or reference, and the function will then be invoked with the arguments of the call. If the conversion function cannot be called (e.g., because of an ambiguity), the program is ill-formed.
4 The argument list submitted to overload resolution consists of the argument expressions present in the function call syntax preceded by the implied object argument (E).

以下google翻訳片手の超意訳

16.3.1.1.2 クラスのオブジェクトに対する関数呼び出し
1. obj(P1...Pn)と書かれたときにobjがクラスT(そのCV修飾含む)のオブジェクトであれば、Tの関数呼び出し演算子を（オーバーロードの）候補関数に加える。
2. さらに、Tに暗黙の型変換関数が（1つ以上）宣言されているとき、その変換先の型が
Rを返し引数(P1...Pn)を受け取る関数ポインタ
Rを返し引数(P1...Pn)を受け取る関数ポインタへの参照
Rを返し引数(P1...Pn)を受け取る関数への参照
であれば（CV修飾は同じかより厳しいものである必要がある）
次のような形式の一意の名前(call-function)を持つ代理呼び出し関数として、それらも候補関数に加える。
R call-function ( conversion-type-id F, P1 a1, ..., Pn an) { return F (a1, ..., an); }
※conversion-type-idは型変換関数の変換先の型
そのような暗黙の型変換演算子がTの基底クラスで宣言されているときも同様に扱われる。ただし、そのような型変換演算子が別の宣言によってT内に隠蔽されていない場合に限る（同じシグネチャでTとその基底で型変換演算子が宣言されている場合はTの宣言が優先される）。
3.そのような代理呼び出し関数がオーバーロード解決の結果選択されたとき、対応する暗黙の型変換関数が呼び出され、適切な関数ポインタまたは参照に変換し、与えられた引数で呼び出される。型変換関数を呼び出すことが出来ない場合（例えば、曖昧さのため）はそのプログラムは不適格（ill-formed）となる。
4. オーバーロード解決に用いられる引数リストは、関数呼び出し構文の引数リストの先頭に暗黙のオブジェクト引数を置いた形で構成される（メンバ関数ならthis、代理呼び出し関数ならばconversion-type-id）

つまりは、クラス内の関数ポインタへの暗黙変換演算子は一旦関数の形に置き換えられオーバーロード解決を行うようです。
そしてオーバーロード解決後に型変換等を行う様子。そして、その一旦置き換えられる関数をsurrogate call functionと呼ぶようです。

ちなみにこの挙動、テンプレート完全ガイドという本には記載があったり、stackoverflowの質問が2012年だったりするので、C++03の時代から存在している様子？

なんだかC++の闇を覗いてしまったような気がします。ラムダや関数ポインタを一つのクラスにごちゃ混ぜにして、統一的に呼び出すのは確かにどこかで使うことが出来そうではありますが・・・

おまけ、任意のラムダで同じことをする

やはり時代はラムダなので、ラムダでも同じようなことがしたいです。しかし、surrogate call functionとなれるのは関数ポインタへの変換だけなので、同じ手は使えません。
2つのラムダの関数呼び出し演算子を一所に集め、オーバーロード解決をしたうえで呼び出す、それが出来るのは継承しかありません（多分）。

template<typename Func1, typename Func2>
struct Blender : Func1, Func2 {

    constexpr Blender(Func1& f1, Func2& f2)
        : Func1{f1}
        , Func2{f2}
    {}

    constexpr Blender(Func1&& f1, Func2&& f2)
        : Func1{std::move(f1)}
        , Func2{std::move(f2)}
    {}

    using Func1::operator();
    using Func2::operator();
};

int main()
{
    int bi = 10;
    double bd = 11.0;
    
    Blender blender{[bi](int n){func1(bi + n);}, [bd](double f){func2(bd + f);}};
        
    blender(128);
    blender(1024.2048);
}

[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

C++におけるラムダの正体は関数オブジェクトなので、当然継承できます。公開継承しておけば関数呼び出し演算子も継承先で公開されるため、複数のラムダを継承することで一つのクラスに複数の関数呼び出し演算子を実装できます。
そのクラスのオブジェクトで関数呼び出しを行えば、あとは普通のオーバーロード解決と同じプロセスで適切なオーバーロードが選択されます。
ちなみに、C++17のコンストラクタ引数からのテンプレート引数推定が無いと、いったんラムダをautoで受けてからdecltypeするという手順を踏まなければならず、あまり実用的ではありません。
（もう少し素敵な実装→C++ - C++における関数の受け渡しについて｜teratail）

参考文献

c++ - How does surrogate call function work? - Stack Overflow
N4659

2018-10-23

［C++］メンバに参照型を持つクラス（構造体）の取り扱い

C++ std::launder placment new

C++のクラスはそのメンバとして任意の参照型を持つことができます。その初期化はコンストラクタ初期化子のみで行えますが、それ以降参照そのものを変更することはできません（その変数に対する代入等の操作は全て参照先に対する操作になるため）。そのようなクラスの取り扱いについてのまとめです。
以下のような簡単な実装で考えていきます。

struct Test {
    Test(int& value) : m_ref{value}
    {}
    
    operator int() const {
        return m_ref;
    }
        
private:
    int& m_ref;
};

特性

先ほど定義したTest型はコンストラクタのみを明示的に定義していますが、その他の特殊関数を宣言も定義もしていないため、コンパイラによって暗黙的に定義されていることが期待されます。
しかし、その場合に構築や代入は何が許されて何が許されないのでしょうか？type_traitsヘッダのメタ関数を用いて先ほどのTest型に対して以下のような特性を見てみます。

    std::cout << std::is_trivially_copyable_v<Test> << std::endl;
    std::cout << std::is_trivial_v<Test> << std::endl;
    std::cout << std::is_pod_v<Test> << std::endl;
    std::cout << std::is_copy_constructible_v<Test> << std::endl;
    std::cout << std::is_move_constructible_v<Test> << std::endl;
    std::cout << std::is_copy_assignable_v<Test> << std::endl;
    std::cout << std::is_move_assignable_v<Test> << std::endl;

コードと実行結果
[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

結果の一覧

特性	type_traits	結果
トリビアルコピー可能	is_trivially_copyable	true
トリビアル型	is_trivial	false
POD型	is_pod	false
コピー構築可	is_copy_constructible	true
ムーブ構築可	is_move_constructible	true
コピー代入可	is_copy_assignable	false
ムーブ代入可	is_move_assignable	false

構築に関してはコピー/ムーブともに可能ですが、代入は一切できません。自分で定義したとしてもコピー/ムーブコンストラクタではコンストラクタ初期化子を使えるので、言われてみれば納得の結果です。
また、デフォルトコンストラクタが無いのでトリビアル型（コンストラクタとコピーコンストラクタを自分で定義していない型、default指定はok）にはなりえず、そのためPOD型（C言語のstruct/unionと互換性のある型）にもなれないのも分かります。
トリビアルコピー可能（memcopyでコピーしても問題ない）がtrueになっているのはその条件を満たしているためですが、どう考えてもmemcopyしない方が良いでしょう。
※トリビアル～、というのはコンパイラが暗黙的に生成してくれる各特殊関数について、自分で定義していない状態の事です。ただし、default指定をした場合はトリビアルであるとして扱われます。

コードにすれば以下のようになります。

int n = 1024;
Test t1{n};
Test t2{t1};   //コピー構築、ok
Test t3{std::move(t1)};   //ムーブ構築、ok

int m = 2048;
Test t4{m};
t4 = t3;   //コピー代入、ng
t4 = std::move(t3);   //ムーブ代入、ng

集成体(Aggregate)の場合

参照型メンバの初期化はコンストラクタ初期化子のみで行えると言いましたが、実は集成体にすることもできます。

struct Test2 {
    int& ref;
    
    operator int() const {
        return ref;
    }
};

//集成体初期化が可能
int v = 1;
Test2 t2 = {v};

こうしたとしても性質はほとんど変わりません。

特性	type_traits	結果
トリビアルコピー可能	is_trivially_copyable	true
トリビアル型	is_trivial	true
POD型	is_pod	false
コピー構築可	is_copy_constructible	true
ムーブ構築可	is_move_constructible	true
コピー代入可	is_copy_assignable	false
ムーブ代入可	is_move_assignable	false
集成体	is_aggregate	true
デフォルト構築可	is_default_constructible	false

コードと実行結果
[Wandbox]三へ( へ՞ਊ ՞)へﾊｯﾊｯ

構築や代入に関してはそのままです。集成体なので、is_aggregateがtrueになります。そして、コンストラクタを一切定義していないためトリビアル型であると判定されるようになりますが、参照型が初期化必須のためデフォルト構築不可、という不思議な状態になります。そして、トリビアル型ではありますがPOD型にはなれません、C言語のライブラリ等には渡さないようにしましょう。

注意点

以上のように、メンバとして参照を持つ場合、何も考えずに定義しても少し変わった扱いになる事が分かります。
そして、そこに罠が潜んでいたりするものです・・・

コンストラクタで参照を初期化する場合の引数の型

これは罠というほどではない軽度な物ですが、割とミスりやすいものです。
先ほどのTest型で言うなら以下のようにしてしまった場合です。

struct Test {
    Test(int value) : m_ref{value}
    {}
    //省略
};

コンストラクタ引数を参照で受け取らずに初期化を試みています。多分コンパイルは通るでしょう。しかし実行されたとき、コンストラクタのvalueは渡された変数がコピーされて構築されます、そしてその参照でm_refを初期化します。ところが、valueはこのコンストラクタ内のローカル変数ですので、コンストラクタの処理が終わればその寿命も尽きます。つまり、m_refは不正な参照になってしまいます。これは未定義動作ですので何が起こるかはわかりません。

そのストレージ領域を再利用する場合（placement new）

こちらは完全に罠です、知らなければ気付きようもありません。しかし、これを行おうと思う事はめったにないのではないかとも思います。
placement newを使うとすでに存在している変数の領域を再初期化することが出来ます。その際にはコンストラクタを呼び出して初期化します。
これはどのような型の領域に対しても可能ですが、例外があります。
非静的メンバにconst修飾されたメンバ変数、もしくは参照型を持つ型に対するplacement new は未定義動作となります。すなわち、何が起こるかわかりません・・・
なぜでしょうか？
それは、この二つはコンストラクタで初期化されて以降は内容が変化しないとみなすことが出来るからです。それは最適化に利用されます。つまり、constもしくは参照型のメンバは値が変わらないため、定数として埋め込むことが出来ます。その結果、placement newで値を変更したつもりでも、埋め込まれた定数が帰り続けることが起こりえます。

int n = 256;
Test* t1 = new Test{n};   //通常のnew

int m = 512;
new(t1) Test{m};   //placement new、戻り値を捨てる

n = 1024;
m = 2048;

//ここで、t1はnとmのどちらを指すことになるかは未定義。
std::cout << t1 << std::endl;

これを解決する手段は主に二つあります。
一つは、placement newの戻り値を正しく利用することです。

int n = 256;
Test* t1 = new Test{n};   //通常のnew

int m = 512;
t1 = new(t1) Test{m};   //placement new、戻り値で更新（値としては変わらない）

n = 1024;
m = 2048;

//t1はmを指す。
std::cout << t1 << std::endl;

これはC++のバージョンに関係なく利用可能です、可能ならこうするべきです。

もう一つはC++17で追加されたstd::launderを利用することです。

int n = 256;
Test* t1 = new Test{n};   //通常のnew

int m = 512;
new(t1) Test{m};   //placement new、戻り値は捨てる。

n = 1024;
m = 2048;

//t1はmを指す。
std::cout << *std::launder(t1) << std::endl;

std::launderはポインタを受け取りそのポインタをそのまま返します。役割はコンパイラへ最適化抑制を指示することです。constメンバを含む型に対しても同様に利用できます。

自分の手でこんなことをすることはないでしょうが、もしかしたら知らぬ間にハマる危険があります。
例えば、デフォルトでないアロケータを指定された標準コンテナを利用する場合や、中身のよく分かっていないライブラリにこのような型を渡す場合等、自分のあずかり知らない所で引っかかるかもしれません。
その場合上記の解決策ではどうしようもないわけですが・・・

参考文献

type_traits - cpprefjp C++日本語リファレンス
 http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2017/p0532r0.pdf

パッケージバージョン	C#バージョン
1.x	6.0
2.0	7.0
2.1	7.0
2.2	7.0
2.3	7.1
2.4	7.2
2.6	7.2
2.7	7.2
2.8	7.3
2.9	7.3

パッケージバージョン	C#バージョン
1.x	6.0
2.0	7.0
2.1	7.0
2.2	7.0
2.3	7.1
2.4	7.2
2.6	7.2
2.7	7.2
2.8	7.3
2.9	7.3

パッケージバージョン	C#バージョン
1.x	6.0
2.0	7.0
2.1	7.0
2.2	7.0
2.3	7.1
2.4	7.2
2.6	7.2
2.7	7.2
2.8	7.3
2.9	7.3