Boost.uBLAS の補完ライブラリの使い方 - その4

boost C++ uBLAS

今回はベクトルの外積と，前回紹介したuBLASの補完ライブラリの外積の使い方について説明します．
今回説明するのは次の機能です．

  • cross_prod

uBLAS には外積を計算する関数がありません．outer_prod というそれらしい名前のものがありますが，これは外積ではなく直積(Outer product)と呼ばれるものを計算します．3次元ベクトルの外積は本来クロス積やベクトル積(Cross product)と呼ばれるものです．
補完ライブラリではこの外積（クロス積）を式テンプレートによる遅延評価で実装してあります．
しかし，もしこの補完ライブラリを導入するのが面倒で（numeric.hpp を include するだけですが），かつ先行評価(Eager evaluation)で十分だというのであれば，以下の関数を使うことができます．

// "Eager" evaluation cross product
template <class V1, class V2>
boost::numeric::ublas::vector<typename boost::numeric::ublas::promote_traits<typename V1::value_type,
                                                                             typename V2::value_type>::promote_type>
cross_prod(const V1& lhs, const V2& rhs)
{
    BOOST_UBLAS_CHECK(lhs.size() == 3, boost::numeric::ublas::external_logic());
    BOOST_UBLAS_CHECK(rhs.size() == 3, boost::numeric::ublas::external_logic());

    typedef typename boost::numeric::ublas::promote_traits<typename V1::value_type,
                                                           typename V2::value_type>::promote_type promote_type;

    boost::numeric::ublas::vector<promote_type> temporary(3);

    temporary(0) = lhs(1) * rhs(2) - lhs(2) * rhs(1);
    temporary(1) = lhs(2) * rhs(0) - lhs(0) * rhs(2);
    temporary(2) = lhs(0) * rhs(1) - lhs(1) * rhs(0);

    return temporary;
}

実際，外積を遅延評価で実装する意味はあまりないかもしれません．外積の結果は3次元ベクトルでしかないですし，結果の一部の成分しか必要としないなんてことも稀であるだろうからです．しかしだからといって積極的に先行評価を用いる理由もあまりありませんが．


以下のテストコードで補完ライブラリの cross_prod の動作を確認できます．

test.cpp

#include "numeric.hpp" // for "lazy" evaluation cross product
#include <boost/numeric/ublas/io.hpp>


int main()
{
    boost::numeric::ublas::vector<double> v1(3);
    boost::numeric::ublas::vector<double> v2(3);

    v1(0) = 1; v1(1) = 2; v1(2) = 3;
    v2(0) = 4; v2(1) = 5; v2(2) = 6;

    // 計算結果を http://en.wikipedia.org/wiki/Cross_product#Examples と比較してみましょう
    std::cout << numeric::ublas::cross_prod(v1, v2) << std::endl;

    return 0;
}
  • VC9
  • g++ (GCC) 3.4.4 (cygming special)
  • g++ (GCC) 4.1.3

でのコンパイルを確認しています．

実行してみると，正しく外積を計算できていることが分かります．


一つだけ注意が必要なのは，prod(A, prod(B, C)) が禁止されているのと同じ理由で，cross_prod(a, cross_prod(b, c)) が禁止されているということです．つまり，cross_prod は複雑(complexity)な式を受け取ることができず，かつ cross_prod は複雑な式であると定義しています．そのため上記のベクトル三重積のような計算をしたい場合には，prod のときと同様に次のように一旦 vector 等の複雑でない式を表すコンテナに代入し，すべての要素を評価してから cross_prod に渡す必要があります．

ublas::vector<double> v1(3), v2(3), v3(3);

v1(0) = 1; v1(1) = 2; v1(2) = 3;
v2(0) = 4; v2(1) = 5; v2(2) = 6;
v3(0) = 7; v3(1) = 8; v3(2) = 9;

// numeric::ublas::cross_prod(v1, numeric::ublas::cross_prod(v2, v3));                 これはコンパイルエラー
numeric::ublas::cross_prod(v1, vector<double>(numeric::ublas::cross_prod(v2, v3))); // 正しくはこう呼び出す


次回は Boost による正規乱数について説明する予定です．
それでは．