Low dimensional matrices, vectors, quaternions and a cubic equation solver

Gareth Charnock gareth.tpc at gmail.com
Sun Apr 18 16:11:35 PDT 2010


Andrei Alexandrescu wrote:
> On 04/16/2010 04:25 PM, Gareth Charnock wrote:
>> Okay, here goes. I've collected together the basic functionality that
>> would probably make a good starting point. As I've mentioned my code is
>> very messy and has bits missing (e.g. I never had a use for the cross
>> product but it's pretty important in general). I guess a good way to
>> begin would be to write the pubic interfaces then start on the
>> implementation.
>>
>> Cubic Solvers:
>> General complex cubic solver with two algorithms (one requiring a
>> complex cosine and and arcosine one using only +-*/ and roots). A
>> special case cubic solver for the reduced cubic x^^3 + px - q = 0.
> 
> This sounds a candidate for std.numeric. How popular/needed are cubic 
> solvers?

They're probably don't come up that frequently but they do they're 
rather fiddly. Lots of operations to get right. But the solution doesn't 
depend on anything but basic math operators so once it's written, it's 
written. I guess the question is whether Phobos is meant to be a small 
library or a kitchen sink library.

>> Quaternions:
>> opAdd, opSub, opMult(quaternion), opMult(vector), opDiv, Normalise,
>> Normalized, conjugate, conjugated, toEulerAngles*, fromEulerAngles,
>> this(real,i,j,k), this(angle,axis), getAngle(), getAxis()
> 
> Sounds good, but you'd need to convert the code to the new overloaded 
> operators approach.

Fair enough, and this will be a good opportunity to show off why the new 
overloading scheme is more powerful (e.g. opAdd and opSub can be combined).

>> *Currently I've only got a quaternion->euler angles routine that works
>> in the ZYZ convention but I have read a paper that generalises my method
>> to all valid axis conventions. Will probably impliment as something like:
>> toEulerAngles(string convention="XYZ")()
>>
>> Vectors:
>> opAdd, opSub, opMult(scalar), opMult(vector)*, cross**, Normalise,
>> Normalized, Length
> 
> What is the representation of vectors? I'm afraid the design above would 
> be too limited for what we need.
A fixed sized array where V[0] ~ x, V[1] ~ y and V[2] ~ z. The field the 
vector is defined over is templated.

What other operators are needed? I'd defiantly want to add swizzling. 
http://www.ogre3d.org/docs/api/html/classOgre_1_1Vector3.html looks like 
it could be a good source of ideas.

>> * dot product. Would this be better named as dot()?

> We already have dot product and normalization routines that work with 
> general ranges.
> 
> http://www.digitalmars.com/d/2.0/phobos/std_numeric.html#dotProduct
> http://www.digitalmars.com/d/2.0/phobos/std_numeric.html#normalize
> 
> Generally I'd strongly suggest making operations free generic functions 
> instead of members.

I've not really thought about operators as members vs operators as free 
functions. I just tend to put them as members because it feels more 
organised. But looking at other implementations, I seem to be in the 
minority.

>> ** 3D vectors only. Perhaps defining a cross product on the
>>
>> Matrices:
>> opAdd, opSub, opMult(scalar), opMult(vector), opMult(matrix)**, Invert,
>> Inverted, Orthogonalize, Orthogonalized, Reorthogonalize***,
>> Reorthogonalized***, Det, Transpose, Transposed, Dagger*, Daggered*,
>> Eigenvalues****, Eigenvectors****
> 
> What is the representation of matrices?
private:
   F[n*n] mat;
where F is the type of the field and n is the dimension of the matrix.

>> *The hermitian conjugate/conjugate transpose. Reduces to the transpose
>> for a real matrix
> 
> Transposition should also be handled in a static manner, e.g. define a 
> transposed view of a matrix that doesn't actually move elements.

Do you mean that a bool should be stored to count the number of 
transpositions or that there is a type that behaves like a matrix but 
actually just presents a view of another matrix e.g.

Matrix A;
...
Matrix B = A.Transposed();
//Changes to B now affect A

>> ** Matrix-matrix multiplication doesn't commute. Could this be a problem
>> when using operator notation?
> 
> Should be fine.
> 
>> *** Othogonalize assuming the matrix is nearly orthogonal already
>> (possibly using some quick, approximate method such as a Taylor series)
>> **** I have a eigenvalue/vector solver for 3x3 matrices which seems
>> reasonably stable but needs more testing.
>>
>> Free functions:
>> MatrixToQuaternion
>> QuaternionToMatrix
>> + code to allow easy printing to stdout/streams and such.
> 
> Sounds encouraging.
> 
> I think a good next step is to go through a community scrutiny process 
> by dropping the code somewhere on the Web so people can review it.
Couldn't agree more because I'm sure I'll miss tricks and conventions. I 
would have never thought of that funky swizzling idea.

I've also got another question: should matrices, vectors and quaternions 
be classes or structs? My gut reaction is that they should be structs 
and thus act like value types. But matrices might be too big and should 
be passed by reference which would imply they should be a class. Anyone 
know any rules of thumb that might apply?

Gareth Charnock



More information about the Digitalmars-d mailing list