数値解析で使う3次元ベクトルクラスの一案

普通の3次元ベクトルクラスはこう書く、と思う。

class qta3dVector
{
public:
    //  成分
    double  x, y, z;

public:
    //  コンストラクタ
    qta3dVector( double rx, double ry, double rz ) :x( rx ), y( ry ), z( rz ){}
    //  長さ
    double  length(){ return ..; }
    //  外積
    qta3dVector crosproduct( const qta3dVector& vec ){ return ..; }
    //  内積
    double  innerproduct( const qta3dVector& vec ){ return ..; }
    //  回転
    bool    rotate( const qtaQuaternion& q ){}
};

（クラス名のqtaはQiitaの略です。）

これは分かりやすいです。でも次のような場合を考えてみます。

可読性や保守に問題がある気もします。

　構造計算などの数値解析において、N次元の連立方程式を解き、N自由度のdouble型の解ベクトルが求まるとします。このベクトルが3次元ベクトルの集まりと考えた時、これを全部ある回転クォータニオンqで回転させるには、最初のクラスを使うと、下記のようになります。

std::vector< double > vecSolution( N ); //  N自由度の解ベクトル
qtaQuaternion   q;  //  回転クォータニオン

for( int i = 0; i < N / 3; ++i ){
    //  値を取り出して3Dベクトルを作成
    const double x = vecSolution[ i * 3 + 0 ];
    const double y = vecSolution[ i * 3 + 1 ];
    const double z = vecSolution[ i * 3 + 2 ];
    qta3dVector vec( x,y,z );   //  ベクトル作成

    vec.rotate( q );    //  回転

    //  回転した値を元に戻す
    vecSolution[ i * 3 + 0 ] = vec.x;
    vecSolution[ i * 3 + 1 ] = vec.y;
    vecSolution[ i * 3 + 2 ] = vec.z;
}

値を元に戻すのでソースも長くなるし、何をやってるのか分かりにくい。xyzがメモリに連続して配置していることを前提にし、非メンバ関数にすれば一行で済みます。

for( int i = 0; i < N / 3; ++i ){
    rotate3dVector( &vecSolution[ i * 3 ]); //  回転
}

しかし、メンバ関数の方が便利なケースもあるため、それを残すと、両方の関数を管理する必要があります。メンバ関数から上記の非メンバ関数を呼ぶと言う方法もあるのですが、以下のようにクラスを２つ作ったほうが汎用的なのではないかと考えております。

クラスを2つに分けたらどうでしょう。

以下のように、3成分の先頭ポインタだけを持つベクトルのクラスと、3成分の実体を持つベクトルの２つに分け、関数を共有するのはどうでしょう。

////////////////////////////////////////////////////////////////////////
//  3成分の先頭ポインタだけを持つベクトル
class qta3dVector_0
{
private:
    double* x;  //  3成分へのポインタ

public:
    double& x(){ return x[ 0 ]; }
    double& y(){ return x[ 1 ]; }
    double& z(){ return x[ 2 ]; }

    //  コンストラクタ
    qta3dVector_0( double* rx ) :x( rx ){}
    //  長さ
    double  length(){ return ..; }
    //  内積
    double  innerproduct( const qta3dVector_0& vec ){ }
    //  回転
    bool    rotate( const qtaQuaternion& q ){}
};

////////////////////////////////////////////////////////////////////////
//  3成分の実体を持つベクトル
class qta3dVector_1 : public qta3dVector_0
{
    //  成分
    double  x[ 3 ];

public:
    //  コンストラクタ
    qta3dVector_1( double rx, double ry, double rz ):
        qta3dVector_0( &x[0] )
    {
        x[ 0 ] = rx;
        x[ 1 ] = ry;
        x[ 2 ] = rz;
    }

    //  外積
    qta3dVector_1 crosproduct( const qta3dVector_0& vec ){  }
};

クラス名末尾の0と1は、そっけないですがおいといて、これにより前述の回転の実装は次のようになります。

for( int i = 0; i < N / 3; ++i ){
    qta3dVector_0   vec( &vecSolution[ i * 3 ]);
    vec.rotate( q );    //  回転
}

これで計算した値を取り出し、計算結果を戻す処理は不要となり、非メンバ関数となることも防ぐことができます。また、3成分の実体を持つベクトルqta3dVector_1もメンバ関数を共有でき、かつqta3dVector_0と_1とで共通の関数を使うことも可能です。データの取り出し、戻す処理が不要となり速度の向上にも繋がると思います。（試してないです）

終わり

以上は、配列や行列に関しても同様のことが適用可能であります。特に行列においては、値のコピーが計算時間に影響がある場合があり、それに対してこの方法は有効ではないか、と考えております。まだ確認はしておりません。またいつか試してみようかと思います。