この記事はC++ Advent Calendar 2016の19日目の記事です．
初めに
普段、計算科学系の研究をしていると行列ライブラリおよび複素数系のライブラリは必須になるかと思います。
多くの方はEigenを使用していると考えられますが、今回は自分向けにc++の勉強がてらに少し便利に記述できる行列ライブラリと超複素数系ライブラリを作ってみました。
以下にgithubのリポジトリをのけっておきます。
READMEにも少し記述しましたが、今のところc++14及びgcc5上のみでテストしています。
また依存ライブラリはboostのpreprocessorのみです。
サンプルのプログラムにはテストをそのまま乗っけていますが、かなりガバガバなテストケースだと思います・・・。今後自分で使っていく上で追加できたらと思います！
行列ライブラリ
行列ライブラリは*1の書籍を参考に作りました。今のところ最低限の機能しか揃っていませんが、以下にサンプルコードを記述します。
matrix_test.cpp
#include <gtest/gtest.h>
#include <cstddef>
#include "../../src/Matrix/Matrix.h"
using namespace omniclid;
TEST(MatrixTest, Matrix)
{
//initializerlistを使って挿入
Matrix<int,1> mi{1,2,3,4};
ASSERT_EQ(1, mi(0));
ASSERT_EQ(4, mi(3));
//普通に挿入。8要素文確保
Matrix<int,3> mi2(2,2,2);
mi2(0,0,0) = 2;
ASSERT_EQ(2, mi2(0,0,0));
ASSERT_EQ(0, mi2(0,0,1));
//initalizer_listの木構造。
Matrix<int,3> m3 { { {{1},{2}},{{3},{4}},{{5},{6}} }, { {{7},{8}},{{9},{10}},{{11},{12}} } , { {{11},{12}},{{13},{14}},{{15},{16}} } , { {{17},{18}},{{19},{110}},{{111},{112}} } };
Matrix<int,2> m3_copy(m3[2]);
ASSERT_EQ(14,m3_copy(1,1));
Matrix<int,4> m4 {
{ { {{1},{2}},{{3},{4}},{{5},{6}} }, { {{7},{8}},{{9},{10}},{{11},{12}} } , { {{1},{2}},{{3},{4}},{{5},{6}} } , { {{7},{8}},{{9},{10}},{{11},{12}} } } ,
{ { {{1},{2}},{{3},{4}},{{5},{6}} }, { {{7},{8}},{{9},{10}},{{11},{12}} } , { {{1},{2}},{{3},{4}},{{5},{6}} } , { {{7},{8}},{{9},{10}},{{11},{12}} } } ,
{ { {{1},{2}},{{3},{4}},{{5},{6}} }, { {{7},{8}},{{9},{10}},{{11},{12}} } , { {{1},{2}},{{3},{4}},{{5},{6}} } , { {{7},{8}},{{9},{10}},{{11},{12}} } } ,
{ { {{1},{2}},{{3},{4}},{{5},{6}} }, { {{7},{8}},{{9},{10}},{{11},{12}} } , { {{1},{2}},{{3},{4}},{{5},{6}} } , { {{7},{8}},{{9},{10}},{{11},{12}} } } ,
{ { {{1},{2}},{{3},{4}},{{5},{6}} }, { {{7},{8}},{{9},{10}},{{11},{12}} } , { {{1},{2}},{{3},{4}},{{5},{6}} } , { {{7},{8}},{{9},{10}},{{11},{12}} } }
};
Matrix<int,3> m4_copy(m4[1]);
ASSERT_EQ(1,m4_copy(0,0,0));
ASSERT_EQ(7,m4_copy(1,0,0));
Matrix<int,2> mult1 { {1,2},{3,4},{5,6} };
Matrix<int,2> mult2 { {1,2,3},{4,5,6} };
auto mult3 = mult1 * mult2;
ASSERT_EQ(9,mult3(0,0));
ASSERT_EQ(51,mult3(2,2));
mult1 += 2;
mult1 += mult1;
auto add3 = mult1 + mult1;
ASSERT_EQ(12,add3(0,0));
ASSERT_EQ(32,add3(2,1));
Matrix<int, 1> single {1,2,3,4,5};
ASSERT_EQ(3,single[2]());
single[2]() = 10;
Matrix<int,0> single_copy(single[2]);
ASSERT_EQ(single_copy(),10);
}
インターフェースはシンプルなのでサンプルに書いてあるのが大体全部です。
便利な点はより汎用性を高めるため行列（2階テンソル）だけではなく階数を自由に指定できる点です。これにより、ベクトルや実数の表現も可能になります。
それに加えて、階数を指定すればinitalizer_listを使った感覚的な要素な格納ができるというところも魅力の一つです！
内部の実装はsliceとreferenceを用いたc++の標準ライブラリっぽい設計になっています。
超複素数
計算科学、特に3DCG系をやっていると複素数などだけではなく、四元数や双体四元数といったトリッキーな代数系が結構出てきます。
外部のライブラリは結構ありますが、勉強がてらに自分で作ってみました！
以下にサンプルコードを記載します。
quaternion_map.cpp
#include "hypercomplex_map.h"
#include <boost/preprocessor.hpp>
namespace omniclid{
//四元数の例
struct quaternion_map : public hypercomplex_map<4>{
public:
//ここに演算表を描く。0はdummy. 実数が1から始まる点に注意（掛けあわせた時の負を表現するため）
using opmap = meta_list<0,
1, 2, 3, 4,
2,-1, 4,-3,
3,-4,-1, 2,
4, 3,-2,-1
>;
//割り算の表を定義
CREATE(16,quaternion_map)
};
}
imagine_map.h
namespace omniclid{
//虚数の例
struct imagine_map : public hypercomplex_map<2> {
public:
/*元々はkeyvalueで演算表を表現していた。
using opmap = meta_map< 0,
kv< x<1,1>, 1> , kv< x<1,2> , 2>,
kv< x<2,1>, 2> , kv< x<2,2> , -1>
>;
*/
using opmap = meta_list<0,
1, 2,
2,-1
>;
CREATE(4,imagine_map)
};
}
hypercomplex_test.cpp
#include <gtest/gtest.h>
#include <cstddef>
#include "../../src/HyperComplex/HyperComplex.h"
#include "../../src/HyperComplex/imagine_map.h"
#include "../../src/HyperComplex/quaternion_map.h"
using namespace omniclid;
TEST(QUATERNION, HyperComplex)
{
//r,i,j,kの型を単体で取り出すこともできる
using QR = typename decltype(make_hypercomplex<0,quaternion_map>())::type;
using QI = typename decltype(make_hypercomplex<1,quaternion_map>())::type;
using QJ = typename decltype(make_hypercomplex<2,quaternion_map>())::type;
using QK = typename decltype(make_hypercomplex<3,quaternion_map>())::type;
QR q(2.0);
QR q2(3);
q += 2;
q += q2;
ASSERT_EQ(q.get_dim(),0);
ASSERT_DOUBLE_EQ(q.get_value(),7.0);
HyperComplex<quaternion_map> complex{1,2,3,4};
HyperComplex<quaternion_map> complex2{1,2,3,4};
complex2 += complex;
QK q3(3);
q3 *= 10.5;
complex2 += q3;
ASSERT_DOUBLE_EQ(complex2.get<3>().get_value(),39.5);
QI i(2);
QJ j(2);
//コンパイル時にqx, qx2の型を推論できる
auto qx = i * j;
auto qx2 = j * i;
ASSERT_EQ(i.get_dim(), 1);
ASSERT_EQ(j.get_dim(), 2);
ASSERT_EQ(qx.get_dim(),3);
ASSERT_EQ(qx2.get_dim(),3);
ASSERT_EQ(i.get_value(), 2);
ASSERT_EQ(j.get_value(), 2);
ASSERT_EQ(qx.get_value(),4);
ASSERT_EQ(qx2.get_value(),-4);
HyperComplex<quaternion_map> q_copy(qx2);
ASSERT_EQ(q_copy.get<1>().get_value(), 0);
ASSERT_EQ(q_copy.get<3>().get_value(),-4);
HyperComplex<quaternion_map> complex3{1,2,3,4};
HyperComplex<quaternion_map> complex4{1,2,3,4};
complex3 *= complex4;
ASSERT_EQ(complex3.get<2>().get_value(),6);
ASSERT_EQ(complex4.get<3>().get_value(),4);
complex4 += qx2;
ASSERT_EQ(complex4.get<3>().get_value(),0);
auto complex5 = qx - qx2;
ASSERT_EQ(complex5.get<3>().get_value(),8);
}
TEST(COMPLEX, HyperComplex)
{
double delta = 0.001;
using IR = typename decltype(make_hypercomplex<0,imagine_map>())::type;
using II = typename decltype(make_hypercomplex<1,imagine_map>())::type;
HyperComplex<imagine_map> imagine1{2,3};
HyperComplex<imagine_map> imagine2{4,5};
auto imagine3 = imagine1/imagine2;
ASSERT_NEAR(imagine3.get<0>().get_value(), 0.560976,delta);
ASSERT_NEAR(imagine3.get<1>().get_value(), 0.0487805,delta);
}
こんな感じに演算表を書けば、勝手に演算が定義できている面白いライブラリに出来たと思います！演算表から自動的に割り算用の演算表も作りたかったのでpreprocessorをどうしても使用してしまいました。
最後に
まだまだ自己満足ライブラリから抜けきれてないので今後はちゃんとテストケースを増やすこと、及びSSE命令などを用いた高速化に挑戦できたらと思います。プルリクやご意見お待ちしております。
明日は I(@wx257osn2) 様の エラーハンドリングを綺麗にこなすためのライブラリ・expectedの紹介と応用 です。
参考文系
*1 プログラミング言語c++（第4版）（ビャーネ・ストラウストラップ著 柴田望洋訳 28章）
https://www.amazon.co.jp/プログラミング言語C-第4版-ビャーネ・ストラウストラップ/dp/4797375957