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

[Lars T. Kyllingstad]

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?

IMO, general vectors/matrices should be structs wrapping a pointer to 
the data:

   struct Vector(T)
   {
      T* ptr;
      size_t length;
      size_t stride;
   }

Low-dimensional fixed-size vectors should probably be value types.

-Lars