namespace GA3D
{

struct mvec
{
	float e[8];
};

static const mvec e1 =   { { 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f } };
static const mvec e2 =   { { 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f } };
static const mvec e3 =   { { 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f } };
static const mvec e23 =  { { 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f } };
static const mvec e31 =  { { 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f } };
static const mvec e12 =  { { 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f } };
static const mvec e123 = { { 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f } };

inline mvec real(float r) {
	mvec v;
	v.e[0] = r;
	v.e[1] = 0.0f;
	v.e[2] = 0.0f;
	v.e[3] = 0.0f;
	v.e[4] = 0.0f;
	v.e[5] = 0.0f;
	v.e[6] = 0.0f;
	v.e[7] = 0.0f;
	return v;
}
inline mvec vec(float x, float y, float z) {
	mvec v;
	v.e[0] = 0.0f;
	v.e[1] = x;
	v.e[2] = y;
	v.e[3] = z;
	v.e[4] = 0.0f;
	v.e[5] = 0.0f;
	v.e[6] = 0.0f;
	v.e[7] = 0.0f;
	return v;
}
inline mvec bivec(float x, float y, float z) {
	mvec v;
	v.e[0] = 0.0f;
	v.e[1] = 0.0f;
	v.e[2] = 0.0f;
	v.e[3] = 0.0f;
	v.e[4] = z;
	v.e[5] = y;
	v.e[6] = x;
	v.e[7] = 0.0f;
	return v;
}
inline mvec rotor(float x, float y, float z, float theta) {
	mvec v;
	float s = sinf(theta/2.0f);
	v.e[0] = cosf(theta/2.0f);
	v.e[1] = 0.0f;
	v.e[2] = 0.0f;
	v.e[3] = 0.0f;
	v.e[4] = s*z;
	v.e[5] = s*y;
	v.e[6] = s*x;
	v.e[7] = 0.0f;
	return v;
}



inline mvec operator+(const mvec &l, const mvec &r) {
	mvec u;
	u.e[0] = l.e[0] + r.e[0];
	u.e[1] = l.e[1] + r.e[1];
	u.e[2] = l.e[2] + r.e[2];
	u.e[3] = l.e[3] + r.e[3];
	u.e[4] = l.e[4] + r.e[4];
	u.e[5] = l.e[5] + r.e[5];
	u.e[6] = l.e[6] + r.e[6];
	u.e[7] = l.e[7] + r.e[7];
	return u;
}
inline mvec operator-(const mvec &l, const mvec &r) {
	mvec u;
	u.e[0] = l.e[0] - r.e[0];
	u.e[1] = l.e[1] - r.e[1];
	u.e[2] = l.e[2] - r.e[2];
	u.e[3] = l.e[3] - r.e[3];
	u.e[4] = l.e[4] - r.e[4];
	u.e[5] = l.e[5] - r.e[5];
	u.e[6] = l.e[6] - r.e[6];
	u.e[7] = l.e[7] - r.e[7];
	return u;
}
inline mvec operator*(const mvec &v, float r) {
	mvec u;
	u.e[0] = v.e[0] * r;
	u.e[1] = v.e[1] * r;
	u.e[2] = v.e[2] * r;
	u.e[3] = v.e[3] * r;
	u.e[4] = v.e[4] * r;
	u.e[5] = v.e[5] * r;
	u.e[6] = v.e[6] * r;
	u.e[7] = v.e[7] * r;
	return u;
}
inline mvec operator*(float r, const mvec &v) { return v * r; }
inline mvec operator/(const mvec &v, float r) { return v * (1.0f/r); }
inline mvec operator/(float r, const mvec &v) { return v * (1.0f/r); }

inline float normsq(const mvec &v){
	return v.e[0]*v.e[0] + v.e[1]*v.e[1] + v.e[2]*v.e[2] + v.e[3]*v.e[3]
		+ v.e[4]*v.e[4] + v.e[5]*v.e[5] + v.e[6]*v.e[6] + v.e[7]*v.e[7];
}
inline float norm(const mvec &v){ return sqrtf(normsq(v)); }
inline mvec normalized(const mvec &v){ return v / norm(v); }
inline mvec operator*(const mvec &l, const mvec &r) {
	mvec v;
v.e[0] = + l.e[0]*r.e[0] + l.e[1]*r.e[1] + l.e[2]*r.e[2] + l.e[3]*r.e[3] - l.e[4]*r.e[4] - l.e[5]*r.e[5] - l.e[6]*r.e[6] - l.e[7]*r.e[7];
v.e[1] = + l.e[0]*r.e[1] + l.e[1]*r.e[0] - l.e[2]*r.e[4] + l.e[3]*r.e[5] + l.e[4]*r.e[2] - l.e[5]*r.e[3] - l.e[6]*r.e[7] - l.e[7]*r.e[6];
v.e[2] = + l.e[0]*r.e[2] + l.e[1]*r.e[4] + l.e[2]*r.e[0] - l.e[3]*r.e[6] - l.e[4]*r.e[1] - l.e[5]*r.e[7] + l.e[6]*r.e[3] - l.e[7]*r.e[5];
v.e[3] = + l.e[0]*r.e[3] - l.e[1]*r.e[5] + l.e[2]*r.e[6] + l.e[3]*r.e[0] - l.e[4]*r.e[7] + l.e[5]*r.e[1] - l.e[6]*r.e[2] - l.e[7]*r.e[4];
v.e[4] = + l.e[0]*r.e[4] + l.e[1]*r.e[2] - l.e[2]*r.e[1] + l.e[3]*r.e[7] + l.e[4]*r.e[0] + l.e[5]*r.e[6] - l.e[6]*r.e[5] + l.e[7]*r.e[3];
v.e[5] = + l.e[0]*r.e[5] - l.e[1]*r.e[3] + l.e[2]*r.e[7] + l.e[3]*r.e[1] - l.e[4]*r.e[6] + l.e[5]*r.e[0] + l.e[6]*r.e[4] + l.e[7]*r.e[2];
v.e[6] = + l.e[0]*r.e[6] + l.e[1]*r.e[7] + l.e[2]*r.e[3] - l.e[3]*r.e[2] + l.e[4]*r.e[5] - l.e[5]*r.e[4] + l.e[6]*r.e[0] + l.e[7]*r.e[1];
v.e[7] = + l.e[0]*r.e[7] + l.e[1]*r.e[6] + l.e[2]*r.e[5] + l.e[3]*r.e[4] + l.e[4]*r.e[3] + l.e[5]*r.e[2] + l.e[6]*r.e[1] + l.e[7]*r.e[0];
	return v;
}
inline mvec operator^(const mvec &l, const mvec &r) {
	mvec v;
v.e[0] = + l.e[0]*r.e[0];
v.e[1] = + l.e[0]*r.e[1] + l.e[1]*r.e[0];
v.e[2] = + l.e[0]*r.e[2] + l.e[2]*r.e[0];
v.e[3] = + l.e[0]*r.e[3] + l.e[3]*r.e[0];
v.e[4] = + l.e[0]*r.e[4] + l.e[1]*r.e[2] - l.e[2]*r.e[1] + l.e[4]*r.e[0];
v.e[5] = + l.e[0]*r.e[5] - l.e[1]*r.e[3] + l.e[3]*r.e[1] + l.e[5]*r.e[0];
v.e[6] = + l.e[0]*r.e[6] + l.e[2]*r.e[3] - l.e[3]*r.e[2] + l.e[6]*r.e[0];
v.e[7] = + l.e[0]*r.e[7] + l.e[1]*r.e[6] + l.e[2]*r.e[5] + l.e[3]*r.e[4] + l.e[4]*r.e[3] + l.e[5]*r.e[2] + l.e[6]*r.e[1] + l.e[7]*r.e[0];
	return v;
}
inline mvec operator|(const mvec &l, const mvec &r) {
	mvec v;
v.e[0] = + l.e[0]*r.e[0] + l.e[1]*r.e[1] + l.e[2]*r.e[2] + l.e[3]*r.e[3] - l.e[4]*r.e[4] - l.e[5]*r.e[5] - l.e[6]*r.e[6] - l.e[7]*r.e[7];
v.e[1] = + l.e[0]*r.e[1] + l.e[1]*r.e[0] - l.e[2]*r.e[4] + l.e[3]*r.e[5] + l.e[4]*r.e[2] - l.e[5]*r.e[3] - l.e[6]*r.e[7] - l.e[7]*r.e[6];
v.e[2] = + l.e[0]*r.e[2] + l.e[1]*r.e[4] + l.e[2]*r.e[0] - l.e[3]*r.e[6] - l.e[4]*r.e[1] - l.e[5]*r.e[7] + l.e[6]*r.e[3] - l.e[7]*r.e[5];
v.e[3] = + l.e[0]*r.e[3] - l.e[1]*r.e[5] + l.e[2]*r.e[6] + l.e[3]*r.e[0] - l.e[4]*r.e[7] + l.e[5]*r.e[1] - l.e[6]*r.e[2] - l.e[7]*r.e[4];
v.e[4] = + l.e[0]*r.e[4] + l.e[3]*r.e[7] + l.e[4]*r.e[0] + l.e[7]*r.e[3];
v.e[5] = + l.e[0]*r.e[5] + l.e[2]*r.e[7] + l.e[5]*r.e[0] + l.e[7]*r.e[2];
v.e[6] = + l.e[0]*r.e[6] + l.e[1]*r.e[7] + l.e[6]*r.e[0] + l.e[7]*r.e[1];
v.e[7] = + l.e[0]*r.e[7] + l.e[7]*r.e[0];
	return v;
}
inline mvec operator~(const mvec &v) {
	mvec u;
	u.e[0] = v.e[0];
	u.e[1] = v.e[1];
	u.e[2] = v.e[2];
	u.e[3] = v.e[3];
	u.e[4] = - v.e[4];
	u.e[5] = - v.e[5];
	u.e[6] = - v.e[6];
	u.e[7] = - v.e[7];
	return u;
}
inline mvec operator!(const mvec &v) {
	mvec u;
	u.e[0] = v.e[7];
	u.e[1] = v.e[6];
	u.e[2] = v.e[5];
	u.e[3] = v.e[4];
	u.e[4] = v.e[3];
	u.e[5] = v.e[2];
	u.e[6] = v.e[1];
	u.e[7] = v.e[0];
	return u;
}
inline mvec operator&(const mvec &l, const mvec &r) { return !(!l ^ !r); }

inline mvec grade(int k, const mvec &v) {
	mvec u;
	switch(k) {
	default:
	case 0:
		u.e[0] = v.e[0];
		u.e[1] = 0.0f;
		u.e[2] = 0.0f;
		u.e[3] = 0.0f;
		u.e[4] = 0.0f;
		u.e[5] = 0.0f;
		u.e[6] = 0.0f;
		u.e[7] = 0.0f;
		break;
	case 1:
		u.e[0] = 0.0f;
		u.e[1] = v.e[1];
		u.e[2] = v.e[2];
		u.e[3] = v.e[3];
		u.e[4] = 0.0f;
		u.e[5] = 0.0f;
		u.e[6] = 0.0f;
		u.e[7] = 0.0f;
		break;
	case 2:
		u.e[0] = 0.0f;
		u.e[1] = 0.0f;
		u.e[2] = 0.0f;
		u.e[3] = 0.0f;
		u.e[4] = v.e[4];
		u.e[5] = v.e[5];
		u.e[6] = v.e[6];
		u.e[7] = 0.0f;
		break;
	case 3:
		u.e[0] = 0.0f;
		u.e[1] = 0.0f;
		u.e[2] = 0.0f;
		u.e[3] = 0.0f;
		u.e[4] = 0.0f;
		u.e[5] = 0.0f;
		u.e[6] = 0.0f;
		u.e[7] = v.e[7];
		break;
	}
	return u;
}

}