C++
C++14
テンプレートメタプログラミング
浮動小数点数

非型テンプレートパラメータに浮動小数点数を渡す方法

More than 1 year has passed since last update.

「無理だろそんなん」って言う前にちょっくら話を聞いてくださいなお客さん。

みなさん既に御存知の通り1、非型テンプレートパラメータに浮動小数点数を渡すことはできません。

その事実を知った時、「なんでや!」と思った人も多いのではないでしょうか。私も思いました。

少なくとも、浮動小数点数の演算を定数式で行うことができる以上、制限を緩和することもできそうな気がするのですが、まあ今のところできないものはできないのです。

……と、ここで終わってしまっては話になりません。

「それでも俺は浮動小数点数をテンプレートで扱いたい!」と思った人は今までにも何人もいらっしゃったようで、Qiita内にもいくつかの記事があります。


方法1:浮動小数点定数への参照をとる

非型テンプレート引数は、グローバル変数への参照をとることができます。それを利用して、こんな風にやることができます。

template<const double& f>

class op
{
public:
static void show() { std::cout << f << std::endl; }
};

constexpr double pi = 3.14159;

using op_pi = op<pi>;

クラス名に深い意味はありません。

これなら一応、クラス内からテンプレート引数に渡された浮動小数点数の値を参照することもできます。

いやー良かった良かった、解決だ。

とは言いがたいです。


何がいけないのか

クラスopは、double型定数への参照をとっているように見えます。

ですが実際は、静的なデータ領域にあるdouble型変数であれば、何でも参照をとることができます。

試してみましょう。

#include <iostream>


template<const double& f>
class op
{
public:
static void show() { std::cout << f << std::endl; }
};

int main() {
static double x = 0.0;
using op_x = op<x>;
op_x::show();
x = 2.5;
op_x::show();
return 0;
}


stdout

0

2.5

こんな風に書き換えすらできてしまいます。この方法ではとても、テンプレート引数に浮動小数点数を渡しているとは言えません。

という訳で、この方法は却下です。


方法2:特性クラス内定数を定義して、特性クラスを渡す

直接渡すことができないなら、クラスで包んで渡せばいいじゃない

template<typename T>

class op
{
static constexpr double value = T::value;
public:
static void show() { std::cout << value << std::endl; }
};

struct pi {
static constexpr double value = 3.14159;
};

using op_pi = op<pi>;

なるほど、まあ一応ちゃんと定数としてdouble型の値を渡すことができている。

こういうふうに書けばop<T>::valueがコンパイル時定数であることも保証できる。

少なくとも、参照を渡すよりはだいぶ筋が良さそうに思えます。

でも、やっぱりちょっと……。


何がいけないのか

問題は、一々特性クラスを定義しないといけないという点です。

面倒くさい、ということもありますが、特性クラスがdouble型の値と一対一であるという保証はどこにもないので、例えば、

template<typename T>

class op
{
static constexpr double value = T::value;
public:
static void show() { std::cout << value << std::endl; }
};

struct pi {
static constexpr double value = 3.14159;
};

struct poi {
static constexpr double value = 3.14159;
};

using op_pi = op<pi>;
using op_poi = op<poi>;

全く同じ値を渡しているはずなのに、別の特殊化をされてしまう……なんてこともあるっぽい。ぽーい。

そんなことしないっぽい? 一つのライブラリに依存した別々の機能を別の人が開発しても、そうならないと思えるっぽい?

という訳で、このやり方も推奨しがたいです。


方法3:sprout::stringを使ってfloat⇆文字列の変換をして扱う

力技。でもやってる人がいたので紹介します。詳しいことは解説しません。

template に float 値を渡す

力技を可能にするSproutは凄いとは思いますが、見た目複雑怪奇なコードが出来上がりますし、スマートとは言いがたいと思います。

それと、これは結局、テンプレート引数で型を受け取って中からvalueを取り出す、つまり方法2で言った特性クラスを簡単に定義する方法に過ぎないため、方法2と同じ問題が生じます。

リンク先では部分特殊化を使っているのである程度はコンパイルエラーで弾けますが、部分特殊化を一々書かないといけないのもマイナスポイントです。

あとは、邪悪なマクロを使っているのも、個人的には好きません。


方法4:IEEE754の内部表現を分解してテンプレートクラスに渡す

こちらの記事で紹介されているやり方。

浮動小数点数をテンプレートメタプログラミングで扱う

個人的には、今までの中では一番好みです。IEEE 754の内部表現を分解しよう、という考え方もいい。

ただ、このやり方も、方法2における特性クラスを簡単に定義しようという方向です。

ですので、方法2と同じ問題が生じる点は方法3とも変わりません。

また、やはり邪悪なマクロが使われています。


方法5:浮動小数点数を列挙型に変換してしまう

という訳で本題。

今まで見たやり方は全て、浮動小数点数を直接渡すことを諦めて、間接的に渡すパターンでした。

それでもまあ、それなりに問題なく動くとは思います。そもそも浮動小数点数をテンプレート引数にしたいなんて要求はあまり多くはないでしょうし。

ですが、もうちょっと見た目をなんとかしたい。

template<double f>

class op;

こういう書き方ができれば最高。でも流石にそれは無理。無理でも、そこに少しは近づけようじゃないか。

一応言っておきますが、ここで解説するやり方も決してスマートとは言いがたいです。内部的にはかなり泥臭い。

ですが、泥臭い部分を独立したライブラリにまとめてしまえば、外部からの扱いやすさの面では一番ではないかなあ、と思っております。

どうやるのかは上の見出しでも書きましたが、それだけでは何をしたいのか分からない方もいるでしょう。

キャストするの? 小数点以下は切り捨てられるし、桁が大きい場合もオーバーフローだよ? なんて思われる方もいるかもしれません。

もちろん、単純にキャストするわけではありません。

IEEE 754の浮動小数点数の内部表現を使います、と言えば察しがつく方もいるでしょう。

今更説明するまでもないことかもしれませんが――、IEEE 754は浮動小数点数の表現方法や計算方法の標準規格です。

昔はIBMのSystem/360なんかでは違う形式の浮動小数点数が使われていたこともあるようですが、現在ではほぼ全てのアーキテクチャがIEEE 754を採用しています。ハードには詳しくないので、そうじゃないアーキテクチャは何だと問われても答えることはできませんが。

C++の標準規格において、浮動小数点数の形式は処理系定義です。とはいえ、IEEE 754の単精度および倍精度浮動小数点数がほぼデファクトスタンダードになっていますが。

とはいえ互換性を考えるならなるべく処理系依存にはしないほうがいい。

という訳で、floatdoubleの内部表現が何であっても互換性があるように作ろうと思います。


浮動小数点数型の代わりの型

まず用意いたしますのは、今回の操作の核となる列挙型。

列挙型を使うのは、strong typedef代わりになるのと、演算子を定義できるためです。

C++11から入ったenum classで、任意の整数型と同一サイズの列挙型を作ることができます。

namespace alt_float {

enum class f32_t : std::uint32_t {};
enum class f64_t : std::uint64_t {};

}

浮動小数点数の代わりなので、名前空間はalt_floatとしました。

f32_tはIEEE 754の単精度浮動小数点数を、f64_tはIEEE 754の倍精度浮動小数点数を表す型です。

とは言え、現時点では単なる列挙型に過ぎません。

これからこの型に対して、


  • 浮動小数点数からのキャスト演算

  • 加減乗除

  • 浮動小数点数への精度を失わない形での変換

といった操作を行えるようにしようと思います。

具体的には、以下のようなクラスや関数を定義します。

namespace alt_float {

template<typename T>
struct is_alt_float;

template <typename T>
constexpr auto is_alt_float_v = is_alt_float<T>::value;

template<typename T, typename U>
constexpr T cast(U);

template<typename T, typename U>
constexpr T operator +(T, U);
template<typename T, typename U>
constexpr T operator -(T, U);
template<typename T, typename U>
constexpr T operator *(T, U);
template<typename T, typename U>
constexpr T operator /(T, U);
template<typename T>
constexpr T operator -(T);
template<typename T>
constexpr auto operator *(T) -> long double; // 実際は十分な大きさの浮動小数点数型

}

定義の大半が演算子ですね。


特性クラスを定義する

上で散々特性クラスを使うやり方を批判しておいてなんですが、特性クラスを定義します。

どうか怒らないで聞いてほしい。特性クラスを使うやり方が駄目という訳ではないんです。TMPにおいて、特性クラスはとても便利な考え方です。ただ今回の場合、特性クラスを直接渡すだけでは、私が理想とする書き方に辿りつけないのです。

そう言えば、特性クラス特性クラスと言っていますが、用語的に正しいのか、というかそもそもそういう言い方はするのか不安になってきました。

ここらへんのテクニックって、ネット漁ってるだけだと情報が少ないんですよね……。

毎回偉そうに記事書いてますが、C++に関する書籍(例えばEffective modern C++とか)は私ほとんど読んだことがないので、私の常識は世間の非常識になっている可能性が……

まあ、それはともかく。私がここで使っている「特性クラス」というのは、「TMPで必要な定数や型の集合として定義されたクラス」のことです。何かもうちょっと良い言い方があれば教えて下さい。

サイズの違う型をTMPで扱う場合、操作内容はほぼ同じでも、シフト演算の回数だけ異なる、なんてことはよくあります。

え? そんな場面に出くわしたことがない? もっとTMPしようぜ!

そんな時、まあstd::numeric_limitsなんかで解決できる場合もありますが、特性クラスを定義しておくと、それを介して異なる部分を解決することができます。

ごちゃごちゃと書いてもイメージが掴みづらい気がするので、実際のコードを見てみたほうがいいかもしれません。

namespace alt_float {

template<int n>
using float_least_t = first_enabled_t<
std::enable_if<std::numeric_limits<float>::digits >= n, float>,
std::enable_if<std::numeric_limits<double>::digits >= n, double>,
long double
>;

template<typename T>
struct f_traits_base;

template<>
struct f_traits_base<f32_t>
{
using u_type = std::uint32_t;
static constexpr auto fraction_bits = 23;
static constexpr auto exponent_bits = 8;
};
template<>
struct f_traits_base<f64_t>
{
using u_type = std::uint64_t;
static constexpr auto fraction_bits = 52;
static constexpr auto exponent_bits = 11;
};

// 定数式でx*2^nを得る関数
template<typename T>
constexpr T ldexp(T x, int n) noexcept {
constexpr auto lshift64 = (uint64_t{} -1) + T{1.};
constexpr auto rshift64 = T{1.} / lshift64;
if(n>=0) {
while(n>=64) {
x *= lshift64;
n-=64;
}
if(n) x *= uint64_t{1} << n;
}
else {
while(n<=-64) {
x *= rshift64;
n+=64;
}
if(n) x *= rshift64 * 2 * (uint64_t{1} << (63+n));
}
return x;
}

template<typename T>
constexpr T pow2(int n) noexcept { return ldexp(T{1.}, n); }

template<typename T>
struct f_traits
{
using base = f_traits_base<T>;
using u_type = typename base::u_type; // 同一サイズの符号なし整数型
static constexpr auto fraction_bits = base::fraction_bits; // IEEE 754の内部表現における基数部分のビット数
static constexpr auto exponent_bits = base::exponent_bits; // IEEE 754の内部表現における指数部分のビット数
using f_type = float_least_t<fraction_bits+1>; // 少なくとも(fraction+1)bitの精度を持つ浮動小数点数型
static constexpr auto fraction_mask = (u_type{1} << fraction_bits) -1; // and演算で基数部のみを取り出すためのマスク
static constexpr auto exponent_mask = (u_type{1} << exponent_bits) -1; // and演算で指数部のみを取り出すためのマスク
static constexpr auto sign_mask = u_type{1} << (fraction_bits + exponent_bits); // and演算で符号部のみを取り出すためのマスク
static constexpr auto bias = (1 << (exponent_bits-1))-1; // 指数部のバイアス
static constexpr auto denorm_min = pow2<f_type>(-fraction_bits-bias+1); // 最小の非正規化数の値
static constexpr auto norm_min = pow2<f_type>(-bias+1);
};

}

first_enabled_tに関してはここあたり参照。ところで私の拙い英語力では果たしてこの名前でいいのだろうかと悩んでいます。

f_typeに関しては、浮動小数点型の内部表現は処理系依存なので、要求する精度を満たす最小の型を使おうとした結果こうなりました。とは言え、大抵の環境では、f_traits<f32_t>::f_typefloatだし、f_traits<f64_t>::f_typedoubleだと思います。

個人的には、どうしようもないけど気に入らないのがldexpです。どうしたってループが入る。

浮動小数点数の演算では2のn乗は指数部を増減させるだけでいいのに!

苦肉の策として、シフト演算と合わせることでループ回数を削減しています。でももしかしたらコンパイラの最適化に任せたほうがいいのかも……(未計測 & 未逆アセンブル)。

まあ、非型テンプレート引数以外の場所では素直に浮動小数点数型使えばいいと思いますけどね。


castの実装

必要となるのは以下の操作です。


  • プリミティブな型からalt_floatの型へのキャスト


  • alt_floatの型からプリミティブな型へのキャスト

あとついでに、



  • alt_floatの型からalt_floatの型へのキャスト

も必要かな。

static_castを使うと、ビット表現そのままキャストされてしまうので、先ほども書いたようにcastを定義します。

namespace alt_float {

template<bool cond>
using enable_when = typename std::enable_if<cond, std::nullptr_t>::type;

template<typename T>
struct is_alt_float : std::conditional_t<std::is_same<T, f32_t>{} || std::is_same<T, f64_t>{}, std::true_type, std::false_type> {};
template<typename T>
constexpr auto is_alt_float_v = is_alt_float<T>::value;

template<typename T, typename U>
struct is_compatible : std::false_type {};
template<typename T>
struct is_compatible<T, f32_t> : std::conditional_t<std::is_same<T, typename f_traits<f32_t>::f_type>{}, std::true_type, std::false_type> {};
template<typename T>
struct is_compatible<T, f64_t> : std::conditional_t<std::is_same<T, typename f_traits<f64_t>::f_type>{}, std::true_type, std::false_type> {};
template<typename T, typename U>
const auto is_compatible_v = is_compatible<T, U>::value;

template<typename T, typename U, enable_when<is_compatible_v<U, T>> = nullptr>
constexpr T cast(U value) noexcept {
using t = f_traits<T>;
using f_type = typename t::f_type;
using u_type = typename t::u_type;

constexpr auto fb = t::fraction_bits;
constexpr auto fm = t::fraction_mask;
constexpr auto em = t::exponent_mask;
constexpr auto sm = t::sign_mask;
constexpr auto nm = t::norm_min;
constexpr auto bi = t::bias;

if(value != value) return static_cast<T>((em << fb) | fm); // NaNの場合。全部+qNaNにされます
auto sign = value < 0;
if(sign) value = -value;
u_type result{};
if(value == std::numeric_limits<f_type>::infinity()) result = em << fb; // infinityの場合
else if(value < nm) result = static_cast<u_type>(value * pow2<U>(fb + bi - 1)); // 非正規化数の場合
else {
auto e = bi + fb;
while(value >= pow2<U>(fb + 1)) { value *= 0.5; e++; }
while(value < pow2<U>(fb)) { value *= 2; e--; }
result = (static_cast<u_type>(e) << fb) | (static_cast<u_type>(value) & fm);
}
return static_cast<T>(sign ? sm | result : result);
}

template<typename T, typename U, enable_when<is_alt_float_v<T> && !is_alt_float_v<U> && !is_compatible_v<U, T>> = nullptr>
constexpr T cast(U value) noexcept(noexcept(static_cast<typename f_traits<T>::f_type>(value))) {
return cast<T>(static_cast<typename f_traits<T>::f_type>(value));
}

template<typename T, typename U, enable_when<is_compatible_v<T, U>> = nullptr>
constexpr T cast(U value) noexcept {
using t = f_traits<U>;
using u_type = typename t::u_type;
using limits = std::numeric_limits<T>;

constexpr auto fb = t::fraction_bits;
constexpr auto fm = t::fraction_mask;
constexpr auto em = t::exponent_mask;
constexpr auto sm = t::sign_mask;
constexpr auto bi = t::bias;
constexpr auto dm = t::denorm_min;

auto s = !!(static_cast<u_type>(value) & sm);
auto e = (static_cast<u_type>(value) >> fb) & em;
auto f = static_cast<u_type>(value) & fm;
T result{};
if(e == em) { // infinityもしくはNaNの場合
if(!f) result = limits::infinity();
else if(f & (u_type{1} << (fb -1))) result = limits::quiet_NaN();
else result = limits::signaling_NaN();
}
else if(!e) { //非正規化数の場合
result = f *dm;
}
else { //正規化数の場合
result = ldexp<T>(f | (u_type{1} << fb), e-fb-bi);
}
return s ? -result: result;
}
template<typename T, typename U, enable_when<is_alt_float_v<U> && !is_alt_float_v<T> && !is_compatible_v<T, U>> = nullptr>
constexpr T cast(U value) noexcept(noexcept(static_cast<T>(typename f_traits<U>::f_type{}))) {
return static_cast<T>(cast<typename f_traits<U>::f_type>(value));
}

template<typename T, typename U, enable_when<is_alt_float_v<U> && is_alt_float_v<T>> = nullptr>
constexpr T cast(U value) noexcept {
return cast<T>(cast<typename f_traits<U>::f_type>(value));
}

enable_whenはお馴染みのアレ。お馴染みのアレと言って一体どれほどの人に通じるのか分かりませんが。

SFINAEを多用しているので、複雑なコードになっているのは否定できません。


演算子の定義

さて、まあしかし、これで一番核となる部分はできました。後はおまけのようなものです。

template<typename T, typename U> struct common_type : std::common_type<T, U> {};

template<> struct common_type<f32_t, f64_t> { using type = f64_t; };
template<> struct common_type<f64_t, f32_t> { using type = f64_t; };
template<typename T, typename U>
using common_type_t = typename common_type<T, U>::type;

template<typename T, enable_when<is_alt_float_v<T>> = nullptr>
constexpr T operator -(T value) noexcept {
return value ^ f_traits<T>::sign_mask;
}
template<typename T, typename U, enable_when<is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr auto operator +(T l, U r) noexcept {
using f1 = typename f_traits<T>::f_type;
using f2 = typename f_traits<U>::f_type;
return cast<common_type_t<T, U>>(cast<f1>(l) + cast<f2>(r));
}
template<typename T, typename U, enable_when<is_alt_float_v<T> && !is_alt_float_v<U>> = nullptr>
constexpr T operator +(T l, U r) noexcept(noexcept(cast<T>(r))) {
using f_type = typename f_traits<T>::f_type;
return cast<T>(cast<f_type>(l) + static_cast<f_type>(r));
}
template<typename T, typename U, enable_when<!is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr U operator +(T l, U r) noexcept(noexcept(r + l)) { return r + l; }
template<typename T, typename U, enable_when<is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr auto operator -(T l, U r) noexcept {
using f1 = typename f_traits<T>::f_type;
using f2 = typename f_traits<U>::f_type;
return cast<common_type_t<T, U>>(cast<f1>(l) - cast<f2>(r));
}
template<typename T, typename U, enable_when<is_alt_float_v<T> && !is_alt_float_v<U>> = nullptr>
constexpr T operator -(T l, U r) noexcept(noexcept(cast<T>(r))) {
using f_type = typename f_traits<T>::f_type;
return cast<T>(cast<f_type>(l) - static_cast<f_type>(r));
}
template<typename T, typename U, enable_when<!is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr U operator -(T l, U r) noexcept(noexcept(cast<U>(l))) {
using f_type = typename f_traits<U>::f_type;
return cast<U>(static_cast<f_type>(l) - cast<f_type>(r));
}
template<typename T, typename U, enable_when<is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr auto operator *(T l, U r) noexcept {
using f1 = typename f_traits<T>::f_type;
using f2 = typename f_traits<U>::f_type;
return cast<common_type_t<T, U>>(cast<f1>(l) * cast<f2>(r));
}
template<typename T, typename U, enable_when<is_alt_float_v<T> && !is_alt_float_v<U>> = nullptr>
constexpr T operator *(T l, U r) noexcept(noexcept(cast<T>(r))) {
using f_type = typename f_traits<T>::f_type;
return cast<T>(cast<f_type>(l) * static_cast<f_type>(r));
}
template<typename T, typename U, enable_when<!is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr U operator *(T l, U r) noexcept(noexcept(r * l)) { return r * l; }
template<typename T, typename U, enable_when<is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr auto operator /(T l, U r) noexcept {
using f1 = typename f_traits<T>::f_type;
using f2 = typename f_traits<U>::f_type;
return cast<common_type_t<T, U>>(cast<f1>(l) / cast<f2>(r));
}
template<typename T, typename U, enable_when<is_alt_float_v<T> && !is_alt_float_v<U>> = nullptr>
constexpr T operator /(T l, U r) noexcept(noexcept(cast<T>(r))) {
using f_type = typename f_traits<T>::f_type;
return cast<T>(cast<f_type>(l) / static_cast<f_type>(r));
}
template<typename T, typename U, enable_when<!is_alt_float_v<T> && is_alt_float_v<U>> = nullptr>
constexpr U operator /(T l, U r) noexcept(noexcept(cast<U>(l))) {
using f_type = typename f_traits<U>::f_type;
return cast<U>(static_cast<f_type>(l) / cast<f_type>(r));
}
template<typename T, enable_when<is_alt_float_v<T>> = nullptr>
constexpr auto operator *(T value) noexcept {
return cast<typename f_traits<T>::f_type>(value);
}

加減乗除が定義できました。あと、単項*演算子で変換ができるようになっています。

さて、ここまででやりたかったことの大半はできました。

template<alt_float::f32_t f>

struct op {
static void show() {
std::cout << *f << std::endl;
}
};

constexpr auto pi = alt_float::cast<f32_t>(3.14159);
using op_pi = op<pi>;

int main() {
op_pi::show();
}

どうでしょう。実装方法を気にしなければ見た感じ綺麗じゃないでしょうか。

ここで終わってしまってもいいのですが、せっかくなのでユーザー定義リテラルを作ってもっといい感じにしましょう。

namespace alt_float {

namespace literal {
constexpr f32_t operator ""_f32(long double f) noexcept { return cast<f32_t>(f); }
constexpr f64_t operator ""_f64(long double f) noexcept { return cast<f64_t>(f); }
}
}

こうしておけば、


using namespace alt_float::literal;

template<alt_float::f32_t f>
struct op {
static void show() {
std::cout << *f << std::endl;
}
};

constexpr auto pi = 3.14159_f32;
using op_pi = op<pi>;

int main() {
op_pi::show();
}

なんて書けます。すっきりしましたね。


おしまい

という訳で特に使う予定はありません。

誰か使ってくれ。

いかがだったでしょーか。

例によって説明がおざなりな記事だなあと思いますが、多少なりとも何かの役に立てば幸いです。

以前の記事で書きなぐった奴も合わせて、近いうちにライブラリにまとめて公開しようかななんて考えています。

それでは今回はこの辺で。





  1. どの程度までご存知としてしまって良いものかどうか毎回悩むのですが