ルギア君の戯言

雑多な記事。

ベクトルの計算

このあいだ提示した部分を実装しました。
間違っていたら御指摘願います。


著作権さえ消さなければ、何をしてくれてもかまいません。
(あくまで、何かあった時のために僕に連絡できるようにするためであり、著作権は主張する気はあまりないです)
煮るなり焼くなりなんでもどうぞ。

//
// C++ Implementation: mvector.h
//
// Description: Mathematical Vector Calculation
//
//
// Author: Lugia Kun <lugia.kun@gmail.com>, (C) 2008
//
// Copyright: See COPYING file that comes with this distribution
//
//

#ifndef _MVECTOR_H_
#define _MVECTOR_H_

#ifndef __cplusplus
#error mvector.h is only for C++.
#endif // __cplusplus

#include <vector>
#include <list>
#include <string>

namespace std {

// Error Message
#define MISS_SIZE 1	// Two operands have different sizes.
#define DIV_ZERO  2	// Calculating division by 0.
#define EUNKNOWN  3	// Unknown error.

// Error Location
#define LOC_ADD  1	// Error when Adding two vectors.
#define LOC_SUB  2	// Error when Substracting two vectors.
#define LOC_SMUL 3	// Error when Multipling vector and scalor.
#define LOC_IPRD 4	// Error when calculating inner product of vectors.
#define LOC_OPRD 5	// Error when calculating outer product of vectors.
#define LOC_UNKNOWN 0	// Some Error happened. But couldn't specify the place.
typedef class _mvec_ex {
	int  data;
	int  loc;
    public:
	_mvec_ex(int d, int l) { data = d; loc = l; }
	int  getErrorNum() { return data; }
} mvector_exception;

template <class _Tp, class _A = allocator<_Tp> >
class mvector : public vector<_Tp, _A> {
	typedef vector<_Tp, _A> vector_type;
	typedef mvector<_Tp, _A> mvector_type;

    public:	
	mvector() : vector_type() {}
	mvector(int n) : vector_type(n) {}
	mvector(int n, _Tp x) : vector_type(n, x) {}
	
	template<class _InputIterator>
	mvector(_InputIterator first, _InputIterator last) : vector_type(first, last) {}
	mvector(const mvector& x) : vector_type((vector_type) x) {}

	// addition
	mvector_type  operator+(mvector_type b) {
		if(this->size() != b.size()) throw _mvec_ex(MISS_SIZE, LOC_ADD);
		mvector_type tmp;
		for(int i = 0; i < this->size(); i++) {
			tmp.push_back((*this)[i] + b[i]);
		}
		return tmp;
	}

	// substraction
	mvector_type  operator-(mvector_type b){
		if(this->size() != b.size()) throw _mvec_ex(MISS_SIZE, LOC_SUB);
		mvector_type tmp;
		for(int i = 0; i < this->size(); i++) {
			tmp.push_back((*this)[i] - b[i]);
		}
		return tmp;
	}

	// Multiply -1.
	mvector_type  operator-() {
		return operator*(-1);
	}
	
	// Scalor Multiplication
	mvector_type  operator*(_Tp b){
		mvector_type tmp;
		for(int i = 0; i < this->size(); i++) {
			tmp.push_back((*this)[i] * b);
		}
		return tmp;
	}

	// Scalor Division
	mvector_type  operator/(_Tp b){
		mvector_type tmp;
		for(int i = 0; i < this->size(); i++) {
			tmp.push_back((*this)[i] / b);
		}
		return tmp;
	}
	
	// Inner Product.
	_Tp  operator*(mvector_type b){
		if(this->size() != b.size()) throw _mvec_ex(MISS_SIZE, LOC_IPRD);
		_Tp sum;
		for(int i = 0; i < this->size(); i++) {
			sum += (*this)[i] * b[i];
		}
		return sum;
	}

	// Outer Product. a ** b will calc outer product.
	mvector_type operator*(mvector_type *b) {
		if(this->size() != 3 || b->size() != 3) throw _mvec_ex(MISS_SIZE, LOC_IPRD);
		// this outer product can use only for 3rd dimension.
		// if you want to use in 1st dimension, use InnerProduct.
		// if you want to use in 7th dimension, please implement new function. I'm sorry.

		mvector_type tmp;
		mvector_type& l = *this;
		mvector_type& m = *b;
		int y, z;
		for(int i = 0; i < 3; i++) {
			y = (i + 1) % 3;
			z = (i + 2) % 3;
			tmp.push_back(l[y] * m[z] - l[z] * m[y]);
		}
		return tmp;
	}

	mvector_type* operator*(){ return this; }

	mvector_type  operator+=(mvector_type b) {
		(*this) = operator+(b);
		return *this;
	}
};

}

#endif // _MVECTOR_H_

前回に紹介した通りですが、一部変わっています。


a, b, c はベクトル、s はスカラーです。

c = a + b;  // 加算
c = a - b;  // 減算
c = a * s;  // スカラー
c = a / s;  // スカラー倍 (1/s)
c = -a;     // 逆方向ベクトル
s = a * b;  // 内積
c = a ** b; // 外積 (3次元のみ)
c = a * &b; // 外積 (同)
c += b;     // 加算代入 (間違えて実装)

7次元の外積を計算したい場合は、別途実装のこと(工学系的には不要だと思われる)。
1次元の外積は内積の数値をベクトル化することで可能。
それ以外の次元は外積を計算できない*1
(1, 3, 7次元以外は、ウェッジ積*2を用いるが、これも成分計算に関して問題があるので、未実装。必要な方は、各自で実装してみてください。あるいは僕がそれを学習した上で実装することとします。)


また、サイズが間違っているなどの場合には、mvector_exception 例外を送出します。


+= は継承前の vector の動作(後ろに追加)をさせるようにします (そのほうが使い易いでしょう)


これ以外に必要な機能

  • 一次独立判定 (行列側で実装する可能性あり)
  • 演算子の関数版
  • ベクトルの大きさ

注意点

mvector の配列 a と b があった場合、外積の計算は、C++ 本来の計算順序の関係上、

c = a[0] ** b[0];

c = a[0] *(* b)[0];

とみなされます。
したがって、外積でなく、スカラー倍になります。
ですから、この場合で、外積を使用したい場合は、

c = (a[0]) ** (b[0]);

とかっこで括ってやる必要があります。
それがいやな場合には、

c = a[0].OuterProduct(b[0]);

としてください。(実装予定)

*1:http://www12.plala.or.jp/ksp/vectoranalysis/SevenDCrossProd/ を参照のこと。

*2:http://www12.plala.or.jp/ksp/differentialforms/ExteriorAlgebra/ を参照のこと。同ページにも書いてあるが、多くの教科書で「ウェッジ積」のことを「外積」と書いているとのこと。