Low dimensional matrices, vectors, quaternions and a cubic equation

[#ponce]

Gareth Charnock Wrote:

> 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?

I definately think low-dimensional vectors/matrices must be structs.

The cost of writing/reading some float/double/real values sitting next in memory is nowhere near the cost of allocating such a new area for each opAdd. Structs can be pooled too, and referred by pointers.