Matrix Klasse

[pospiech]

Ich möchte möglichst effizient 2 und 3 dimensionale Matrizen mit double Werten multiplizieren, addieren, Determinante bilden usw.

Dazu suche nach Empfehlungen für Klassen für die Matrizen Berechnung.

[Vevusio]

also hier ist meine matrix klasse, ist nicht vollständig aber naja

//H
#pragma once

#ifdef __DLLBUILD
#define DLL __declspec(dllexport)
#else
#define DLL __declspec(dllimport)
#endif

#include "Vector3.h"
#include "Quaternion.h"
#include <windows.h>

namespace PinkballsOpenGL {
	namespace Math {

	class Quaternion;
	class Vector3;

	class DLL Matrix
	{
#pragma region [Variables]
		/**************************** [Variables] ****************************/
	public:
		static const int DIMENSION = 4;
		static const int FIELDS = DIMENSION * DIMENSION;
		static const int _11 = 0;
		static const int _12 = 1;
		static const int _13 = 2;
		static const int _14 = 3;
		static const int _21 = 4;
		static const int _22 = 5;
		static const int _23 = 6;
		static const int _24 = 7;
		static const int _31 = 8;
		static const int _32 = 9;
		static const int _33 = 10;
		static const int _34 = 11;
		static const int _41 = 12;
		static const int _42 = 13;
		static const int _43 = 14;
		static const int _44 = 15;
	private:
		float matrix[FIELDS];
		/*********************************************************************/
#pragma endregion

#pragma region [Constructors and Destructors]
		/****************** [Methods] ******************/
	public:
		Matrix(void);
		Matrix(const Matrix &aMatrix);

		~Matrix(void);
	private:
		Matrix(float aMatrix[FIELDS]);
		/*********************************************************************/
#pragma endregion

#pragma region [Methods]
		/***************************** [Methods] *****************************/
	public:
		float Determinant();
		static Matrix RotationX(float aAngle);
		static Matrix RotationY(float aAngle);
		static Matrix RotationZ(float aAngle);
		static Matrix RotationQuaternion(float aX, float aY, float aZ, float aW);
		static Matrix RotationQuaternion(const Quaternion &aQuaternion);
		static Matrix Translation(float aX, float aY, float aZ);
		static Matrix Translation(const Vector3 &aTranslation);
		static Matrix Scaling(float aX, float aY, float aZ);
		static Matrix Scaling(const Vector3 &aScaling);
		static Matrix RotationYawPitchRoll(float aYaw, float aPitch, float aRoll);

	private:
		void Initialize(void);
		void Initialize(const Matrix &aMatrix);
		void Initialize(float *aMatrix);
		/*********************************************************************/
#pragma endregion

#pragma region [Properties]
		/*************************** [Properties] ***************************/
	public:
		static Matrix Zero(void) { float newMatrix[FIELDS]; ZeroMemory(newMatrix, sizeof(float) * FIELDS); return Matrix(newMatrix); }
		static Matrix Identity(void) { return Matrix(); }

		const float *GetMatrix(void) const { return matrix; }
		float GetM11(void) { return matrix[_11]; }
		void SetM11(float value) { matrix[_11] = value; }
		float GetM12(void) { return matrix[_12]; }
		void SetM12(float value) { matrix[_12] = value; }
		float GetM13(void) { return matrix[_13]; }
		void SetM13(float value) { matrix[_13] = value; }
		float GetM14(void) { return matrix[_14]; }
		void SetM14(float value) { matrix[_14] = value; }
		float GetM21(void) { return matrix[_21]; }
		void SetM21(float value) { matrix[_21] = value; }
		float GetM22(void) { return matrix[_22]; }
		void SetM22(float value) { matrix[_22] = value; }
		float GetM23(void) { return matrix[_23]; }
		void SetM23(float value) { matrix[_23] = value; }
		float GetM24(void) { return matrix[_24]; }
		void SetM24(float value) { matrix[_24] = value; }
		float GetM31(void) { return matrix[_31]; }
		void SetM31(float value) { matrix[_31] = value; }
		float GetM32(void) { return matrix[_32]; }
		void SetM32(float value) { matrix[_32] = value; }
		float GetM33(void) { return matrix[_33]; }
		void SetM33(float value) { matrix[_33] = value; }
		float GetM34(void) { return matrix[_34]; }
		void SetM34(float value) { matrix[_34] = value; }
		float GetM41(void) { return matrix[_41]; }
		void SetM41(float value) { matrix[_41] = value; }
		float GetM42(void) { return matrix[_42]; }
		void SetM42(float value) { matrix[_42] = value; }
		float GetM43(void) { return matrix[_43]; }
		void SetM43(float value) { matrix[_43] = value; }
		float GetM44(void) { return matrix[_44]; }
		void SetM44(float value) { matrix[_44] = value; }
		/*********************************************************************/
#pragma endregion
#pragma region [Operators]
		/*************************** [Operators] ***************************/
	public:
		Matrix &operator=(const Matrix &aMatrix) 
		{
			Initialize(aMatrix);
			return *this;
		}
		Matrix &operator*=(const Matrix &aMatrix);
		Matrix operator*(const Matrix &aMatrix);
		/*********************************************************************/
#pragma endregion
	};

} }

//CPP
#include "StdAfx.h"
#include "Matrix.h"
#include "Geometry.h"

PinkballsOpenGL::Math::Matrix::Matrix(void)
{
	Initialize();
}

PinkballsOpenGL::Math::Matrix::Matrix(const Matrix &aMatrix)
{
	*this = aMatrix;
}

PinkballsOpenGL::Math::Matrix::Matrix(float aMatrix[FIELDS])
{
	Initialize(aMatrix);
}

PinkballsOpenGL::Math::Matrix::~Matrix(void)
{
}

void PinkballsOpenGL::Math::Matrix::Initialize(void)
{
	ZeroMemory(matrix, sizeof(float) * FIELDS);
	matrix[_11] = matrix[_22] = matrix[_33] = matrix[_44] = 1.0f;
}

void PinkballsOpenGL::Math::Matrix::Initialize(const Matrix &aMatrix)
{
	for(int i = 0; i < FIELDS; i++)
	{
		matrix[i] = aMatrix.matrix[i];
	}
}

void PinkballsOpenGL::Math::Matrix::Initialize(float aMatrix[FIELDS])
{
	for(int i = 0; i < FIELDS; i++)
	{
		matrix[i] = aMatrix[i];
	}
}

float PinkballsOpenGL::Math::Matrix::Determinant()
{
	return
		matrix[_41] * matrix[_32] * matrix[_23] * matrix[_14] -
		matrix[_31] * matrix[_42] * matrix[_23] * matrix[_14] -
		matrix[_41] * matrix[_22] * matrix[_33] * matrix[_14] +
		matrix[_21] * matrix[_42] * matrix[_33] * matrix[_14] +
		matrix[_31] * matrix[_22] * matrix[_43] * matrix[_14] -
		matrix[_21] * matrix[_32] * matrix[_43] * matrix[_14] -
		matrix[_41] * matrix[_32] * matrix[_13] * matrix[_24] +
		matrix[_31] * matrix[_42] * matrix[_13] * matrix[_24] +
		matrix[_41] * matrix[_12] * matrix[_33] * matrix[_24] -
		matrix[_11] * matrix[_42] * matrix[_33] * matrix[_24] -
		matrix[_31] * matrix[_12] * matrix[_43] * matrix[_24] +
		matrix[_11] * matrix[_32] * matrix[_43] * matrix[_24] +
		matrix[_41] * matrix[_22] * matrix[_13] * matrix[_34] -
		matrix[_21] * matrix[_42] * matrix[_13] * matrix[_34] -
		matrix[_41] * matrix[_12] * matrix[_23] * matrix[_34] +
		matrix[_11] * matrix[_42] * matrix[_23] * matrix[_34] +
		matrix[_21] * matrix[_12] * matrix[_43] * matrix[_34] -
		matrix[_11] * matrix[_22] * matrix[_43] * matrix[_34] -
		matrix[_31] * matrix[_22] * matrix[_13] * matrix[_44] +
		matrix[_21] * matrix[_32] * matrix[_13] * matrix[_44] +
		matrix[_31] * matrix[_12] * matrix[_23] * matrix[_44] -
		matrix[_11] * matrix[_32] * matrix[_23] * matrix[_44] -
		matrix[_21] * matrix[_12] * matrix[_33] * matrix[_44] +
		matrix[_11] * matrix[_22] * matrix[_33] * matrix[_44] ;
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::operator*(const Matrix &aMatrix) 
{
	float newMatrix[FIELDS];
	float sum;
	int x, y, i;
	for (x = 0; x < DIMENSION; x++)
	{
		for (y = 0; y < DIMENSION; y++)
		{
			sum = 0;
			for (i = 0; i < 4; i++)
			{
				sum += matrix[y + 4*i] * aMatrix.matrix[4*x + i];
			}
			newMatrix[y + 4*x] = sum;
		}
	}
	return Matrix(newMatrix);
}

PinkballsOpenGL::Math::Matrix &PinkballsOpenGL::Math::Matrix::operator*=(const Matrix &aMatrix) 
{
	float newMatrix[FIELDS];
	float sum;
	int x, y, i;
	for (x = 0; x < DIMENSION; x++)
	{
		for (y = 0; y < DIMENSION; y++)
		{
			sum = 0;
			for (i = 0; i < 4; i++)
			{
				sum += matrix[y + 4*i] * aMatrix.matrix[4*x + i];
			}
			newMatrix[y + 4*x] = sum;
		}
	}
	Initialize(newMatrix);
	return *this;
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::RotationX(float aAngle)
{
	float matrix[FIELDS];

	//  |  1  0       0       0 |
	//M=|  0  cos(A) -sin(A)  0 |
	//  |  0  sin(A)  cos(A)  0 |
	//  |  0  0       0       1 |
	matrix[_11] = 1; matrix[_21] = 0;          matrix[_31] = 0;           matrix[_41] = 0;
	matrix[_12] = 0; matrix[_22] = Geometry::Cos(aAngle); matrix[_32] = -Geometry::Sin(aAngle); matrix[_42] = 0;
	matrix[_13] = 0; matrix[_23] = Geometry::Sin(aAngle); matrix[_33] = Geometry::Cos(aAngle); matrix[_43] = 0;
	matrix[_14] = 0; matrix[_24] = 0;          matrix[_34] = 0;          matrix[_44] = 1;

	return Matrix(matrix);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::RotationY(float aAngle)
{
	float matrix[FIELDS];

	//  |  cos(A)  0  sin(A)  0 |
	//M=|  0       1  0       0 |
	//  | -sin(A)  0  cos(A)  0 |
	//  |  0       0  0       1 |
	matrix[_11] = Geometry::Cos(aAngle);  matrix[_21] = 0; matrix[_31] = Geometry::Sin(aAngle);  matrix[_41] = 0;
	matrix[_12] = 0;           matrix[_22] = 1; matrix[_32] = 0;           matrix[_42] = 0;
	matrix[_13] = -Geometry::Sin(aAngle); matrix[_23] = 0; matrix[_33] = Geometry::Cos(aAngle); matrix[_43] = 0;
	matrix[_14] = 0;           matrix[_24] = 0; matrix[_34] = 0;          matrix[_44] = 1;

	return Matrix(matrix);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::RotationZ(float aAngle)
{
	float matrix[FIELDS];

	//  |  cos(A)  -sin(A) 0  0 |
	//M=|  sin(A)   cos(A) 0  0 |
	//  |  0        0      1  0 |
	//  |  0        0      0  1 |
	matrix[_11] = Geometry::Cos(aAngle); matrix[_21] = -Geometry::Sin(aAngle); matrix[_31] = 0;  matrix[_41] = 0;
	matrix[_12] = Geometry::Sin(aAngle); matrix[_22] = Geometry::Cos(aAngle);  matrix[_32] = 0;  matrix[_42] = 0;
	matrix[_13] = 0;          matrix[_23] = 0;           matrix[_33] = 1; matrix[_43] = 0;
	matrix[_14] = 0;          matrix[_24] = 0;           matrix[_34] = 0; matrix[_44] = 1;

	return Matrix(matrix);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::RotationQuaternion(float aX, float aY, float aZ, float aW)
{
	float matrix[FIELDS];

	double xx = aX * aX;
	double xy = aX * aY;
	double xz = aX * aZ;
	double xw = aX * aW;
	double yy = aY * aY;
	double yz = aY * aZ;
	double yw = aY * aW;
	double zz = aZ * aZ;
	double zw = aZ * aW;

	matrix[_11] = (float)(1.0 - 2.0 * ( yy + zz ));
	matrix[_12] = (float)(		2.0 * ( xy + zw ));
	matrix[_13] = (float)(		2.0 * ( xz - yw ));
	matrix[_21] = (float)(		2.0 * ( xy - zw ));
	matrix[_22] = (float)(1.0 - 2.0 * ( xx + zz ));
	matrix[_23] = (float)(		2.0 * ( yz + xw ));
	matrix[_31] = (float)(		2.0 * ( xz + yw ));
	matrix[_32] = (float)(		2.0 * ( yz - xw ));
	matrix[_33] = (float)(1.0 - 2.0 * ( xx + yy ));
	matrix[_14] = matrix[_24] = matrix[_34] = matrix[_41] = matrix[_42] = matrix[_43] = 0.0f;
	matrix[_44] = 1.0f;

	return Matrix(matrix);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::RotationQuaternion(const Quaternion &aQuaternion)
{
	return RotationQuaternion(aQuaternion.X, aQuaternion.Y, aQuaternion.Z, aQuaternion.W);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::Translation(float aX, float aY, float aZ)
{
	float matrix[FIELDS];

	//  | 1  0  0  X |
	//M=| 0  1  0  Y |
	//  | 0  0  1  Z |
	//  | 0  0  0  1 |
	matrix[_11] = 1; matrix[_21] = 0; matrix[_31] = 0; matrix[_41] = aX;
	matrix[_12] = 0; matrix[_22] = 1; matrix[_32] = 0; matrix[_42] = aY;
	matrix[_13] = 0; matrix[_23] = 0; matrix[_33] = 1; matrix[_43] = aZ;
	matrix[_14] = 0; matrix[_24] = 0; matrix[_34] = 0; matrix[_44] = 1;

	return Matrix(matrix);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::Translation(const Vector3 &aTranslation)
{
	return Translation(aTranslation.X, aTranslation.Y, aTranslation.Z);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::Scaling(float aX, float aY, float aZ)
{
	float matrix[FIELDS];

	//  | X  0  0  0 |
	//M=| 0  Y  0  0 |
	//  | 0  0  Z  0 |
	//  | 0  0  0  1 |
	matrix[_11] = aX; matrix[_21] = 0;    matrix[_31] = 0;     matrix[_41] = 0;
	matrix[_12] = 0;    matrix[_22] = aY; matrix[_32] = 0;     matrix[_42] = 0;
	matrix[_13] = 0;    matrix[_23] = 0;    matrix[_33] = aZ; matrix[_43] = 0;
	matrix[_14] = 0;    matrix[_24] = 0;    matrix[_34] = 0;    matrix[_44] = 1;

	return Matrix(matrix);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::Scaling(const Vector3 &aScaling)
{
	return Scaling(aScaling.X, aScaling.Y, aScaling.Z);
}

PinkballsOpenGL::Math::Matrix PinkballsOpenGL::Math::Matrix::RotationYawPitchRoll(float aYaw, float aPitch, float aRoll)
{
	//TODO : This should be faster (recalculate it though cause its late and you screwed up somewhere probably):
	//
	// cosY * cosZ							| cosy * -sinz							| siny			| 0
	// -sinx * -siny * cosz + cosx * sinz	| -sinx * -siny * -sinz + cosx * cosz	| -sinx * cosy	| 0
	// cosx * -siny * cosz + sinx * sinz	| cosx * -siny * -sinz + sinx * cosz	| cosx * cosy	| 0
	// 0									| 0										| 0				| 1
	return ( Matrix::RotationX(aYaw) * Matrix::RotationY(aPitch) * Matrix::RotationZ(aRoll) );
}

das ganze ist hauptsächlich ne objektorientierte version von einer c library für opengl die du dir hier runterladen kannst, ich implementiere immer nur nach und nach die funktionen die ich gerade brauche also wird nicht alles drin sein

die c library gibts hier:
http://members.tripod.com/markus_ilmola/glm.htm

natürlich lässt sich das alles halt vermutlich beschleunigen wenn man es in assembler umschreibt teilweise aber darum kümmere ich mich erstmal nicht

//EDIT
achja noch etwas am rande, openGL verwendet COL mayor matrices und directx ROW mayor matrices, darauf musst du halt achten, also dieser code ist col-mayor das heißt die indices der matrix werden so durchnummeriert

col mayor    row mayor
 0  4  8 12    0  1  2  3
 1  5  9 13    4  5  6  7
 2  6 10 14    8  9 10 11
 3  7 11 15   12 13 14 15

[Tyrdal]

Vielleicht ist boost::uBlas was für dich.

[unskilled]

"Ich möchte möglichst effizient 2 und 3 dimensionale Matrizen mit double Werten multiplizieren, addieren, Determinante bilden usw."

Das kannste eigtl auch relativ leicht selbst schreiben...

Die Klasse von Vevusio bringt dich glaube nicht sooooo sehr viel weiter ^^ (kann eigtl nix von dem, was du willst - weil sie einfach mal nicht fürs rechnen im mathematischen Sinne gedacht ist ^^)
Ich find die Klasse auch nicht wirklich hübsch : P

Wenn du wirklich so extrem viel Wert auf Performance legst und nur 2d und 3d matrizen betrachtest, dann würde ich einzelne Klassen dafür schreiben und auf std::vector verzichten

class Matrix_2D
{
public:
double var[2][2]; //würde eher ins private gehören und dann noch Getter und Setter-Methoden schreiben...
Matrix_2D (bool benutzereingabe = false)
{
if (benutzereingabe)
{
//benutzer auffordern, werte einzugeben und einlesen...
}
else
{
//var[0][0] ... var[1][1] = 0;
}
}

Matrix_2D (const double (& arr)[2][2])
{
var[0][0] = arr[0][0];
//... usw ^^
}

Matrix_2D operator + (const Matrix_2D &a)
{
Matrix_2D temp;
temp[0][0] = var[0][0] + a.var[0][0];
temp[0][1] = var[0][1] + a.var[0][1];
temp[1][0] = var[1][0] + a.var[1][0];
temp[1][1] = var[1][1] + a.var[1][1];
return tmp;
}
Matrix_2D operator += (const Matrix_2D &a);
bool operator == (const Matrix_2D &a);
bool operator != (const Matrix_2D &a);
//todo...

Matrix_2D operator * (const Matrix_2D &a);
Matrix_2D operator *= (const Matrix_2D &a);
/*so einfach geht das nur, wenn du wirklich immer quadratische Matrizen hast (2 * 2 bzw 3 * 3) - aber das ist hier bei mehreren Bsp. so
Das liegt daran, dass die Multiplikation von Determinaten nicht kommutativ ist (falls du nicht weist, was das heißt:
Matrix c = Matrix A * Matrix B
Matrix d = Matrix B * Matrix A

Matrix c und d sind in der Regel NICHT gleich - es könnte passieren, aber das passiert nur bei ausreichend trivialen Matrizen A und B...)

und die erste nur gleich viele Zeilen haben muss, wie die zweite Spalten hat und dementsprechend bei Mat (2*3) * Mat (3*2) eine Mat (3*3) rauskommt (vll ist das auch genau falsch rum-ich hab jz ne nachgedacht ^^) - sollte dich aber nicht weiter interessieren, wenn du nur quadratische betrachtest
*/

double GetDeterminate (void)
{
/*
 A    B
 C    D
Determinante = A*D - B*C
*/
return (var[0][0]*var[1][1] - var[0][1] * var[1][0]);
}
}

Die 3d-Variante geht fast genau so - nur, dass es mehr getippe ist und die Determinante sich so berechnen lässt:

Matrix:
 A B C
 D E F
 G H I

Determinate: (kannst ja ma nach 'Regel von Sarrus' (glaube ich ^^) suchen)

det = A*E*I+F*F*G+C*D*H-C*E*G-B*D*I-A*F*H

sry, dass ich nirgendwo eingerückt habe ^^

bye

[Vevusio]

unskilled schrieb:

(..)weil sie einfach mal nicht fürs rechnen im mathematischen Sinne gedacht ist(..)

naja wie gesagt man müsste halt noch mehr operatoren implementieren und so aber ich seh nicht wieso sie nicht auch für "mathematische sachen" einsetzbar sein sollte

btw, was findest du daran schlecht? bin für jegliche kritik dankbar die mir hilft die klasse zu verbessern (dass die matrix selbst als 1d und nicht 2d array abgespeichert ist liegt daran dass opengl ein 1d float array erwartet also hab ich es gleich direkt so implementiert falls du mit unter das gemeint hast)

[*D*Evil]

naja generell gilt das pointer schneller als index-operatoren sind. Dann kannst du vllt. noch was rauskitzeln, wenn du dir GPG 5 anguckst ... ("Foating Points Tricks")

[pospiech]

Kennt einer 'Eigen' bzw was haltet ihr davon:
http://eigen.tuxfamily.org/

Enthält auch 'fixed size' Matrizen.

Matthias

[unskilled]

geht alles @vevusio:
(hab alles noch ma editiert, weil selbst ich ne mehr durchsah ^^)

Ja, das meinte ich auch - aber dann ist das ja ok ^^ ansonsten find ich sie zwar nicht schlecht aber eben unpassend für ihn ^^

trotzdem noch 2 pkt:

1. als bsp mal dein *= operator:

PinkballsOpenGL::Math::Matrix &PinkballsOpenGL::Math::Matrix::operator*=(const Matrix &aMatrix) 
{
    float newMatrix[FIELDS];
    for (size_t x = 0; x < DIMENSION; ++x) //unsigned char x würde auch reichen - könnte man mit nem typedef am anfang schreiben - da wär es zwar statisch aber man könnte es schnell ändern...
    {
        for (size_t y = 0; y < DIMENSION; ++y)
        {
            float sum = matrix[y] * aMatrix.matrix[4*x]; //nicht erst als 0 deklarieren
            for (size_t i = 1; i < 4; ++i)
            {
                sum += matrix[y + 4*i] * aMatrix.matrix[4*x + i];
            }
            newMatrix[y + 4*x] = sum;
        }
    }
    Initialize (newMatrix);
    return *this;
}

2. zur Determiante:
Ich weiß nicht, was für ne Formel du genommen hast, aber:

Erst mal deine Fkt:

Du hast 24 Zeilen á 
+ 3 Multiplikationen 
+ 1/2 Addition 
+ 1/2 Subtraktion 
====================== 
+ 72 Multiplikationen //ich betrachte jz nur noch die Multiplikationen - Additionen und Subtraktionen sollten wir immer fast gleich viel haben...
+ 12 Subtraktionen
+ 12 Additionen

Also wenn ich den Entwicklungssatz von Laplace nehm und dann die Regel von Sarrus:
(ich habs auch noch 2 Mal mit dem Laplaceschen Entwicklungssatz gemacht und dann ein wenig zusammengefasst - kommt das gleiche raus...)

A (F (KP - LO) + G (LN - JP) + H (JO - KN))

ist die erste (von 4) (3x3) Matrizen
=>
4*10 * Multiplikation

==>
40 Multiplikationen : 72 Multiplikationen
                (ME : YOU)

;o)

Bei dir würde das also so anfangen:

matrix[_11] * ( matrix[_22] * (matrix[_33] * matrix[_44] - matrix[_34] * matrix[_43])+...

ist mir zu viel getippe ^^

Da mir langweilig war hab ich mal ne Funktion draus gemacht:

Der Laplacesche-Entwicklungssatz in ne Fkt verpackt würde bei mir in etwa so aussehen (ok - nen reservere wäre vll noch gut ^^). So kannst du damit auf jeden Fall dann auch Determinanten von jeder beliebigen quadratischen Matrix ausrechnen (natürlich nur so lang, wie keine parameter drin sind ^^ aber auch das sollte nicht unmöglich sein ^^)...

double Determinante (const std::vector <std::vector <double> > &matrix)
			{ 
				const size_t dim = matrix.size ();
				if (dim >= 3)
					{
						const size_t nextdim = dim-1; 
						std::vector <std::vector <double> > nextmatrix;
						double sum (0);
						for (size_t m (0); m < nextdim; ++m) 
							{
								for (size_t n (0); n < nextdim; ++n) 
									{ 
										std::vector <double> ceil;
										for (size_t i (0); i < n; ++i) 
											{
												ceil.push_back (matrix[n+1][i]);
											} 
										for (size_t i (n); i < nextdim; ++i) 
											{ 
												ceil.push_back (matrix[n+1][i+1]); 
											}
										nextmatrix.push_back (ceil);
									}
								sum += matrix[0][m] * Determinante (nextmatrix);
							}
						return sum;
					}
				else if (dim == 2)
					{
						return matrix[0][0]*matrix[1][1]-matrix[0][1]*matrix[1][0];
					}
				else if (dim == 1)
					{
						return matrix[0][0];
					}
				else return 0;
			}

Ich bin mir auch relativ sicher, dass es geht - hab es aber nur einmal ausprobiert... (ok - ich hab jz noch 2 Fehler gefunden gehabt - also so sicher bin ich mir dann doch nicht mehr ^^)
Schade, dass das nicht so recht in den Thread hier reinpasst, aber ich dachte, ich poste es einfach mal mit um bissl rumzuprahlen ;o)

Also bye

[Vevusio]

ah ok, danke für die vorschläge