Re: [boost] [Review][PQS] Review deadline

From: Andy Little (andy_at_[hidden])
Date: 2006-06-08 17:01:10

"Michael Fawcett"  wrote
> What the library would be concerned with is that the result of the
> operation was possible (made sense) given the types that were held
> within the matrix or vector.

Before continuing I better say I'm no expert in linear algebra! I have only
looked at 3D transform matrices.

But yes using Boost.Typeof or  result_of  to run through the expression. IOW the
operations decide the type of the result. Obvious really, but linear algebra
libraries often assume numeric types. However AFAIK a matrix can have matrix or
maybe vector or complex elements so perhaps there are other candidates than
physical quantities for the elements?

> Contrived example using theoretical code:
>
> // Multiply accumulate (Is there a simpler way?)

I think its simpler with Boost.Typeof, You could use it instead of the macc
struct for deducing the result of the expression directly, but gcc wont like it,
so you could use it in the body of macc:

template
<
typename X, typename Y, typename Z,
typename VX, typename VY, typename VZ
>
struct macc
{
    typedef BOOST_TYPEOF_TPL(  X() * VX() + Y() * VY() + Z() * VZ()) type;

};

> template
> <
> typename X, typename Y, typename Z,
> typename VX, typename VY, typename VZ
>>
> typename macc<X, Y, Z, VX, VY, VZ>::type
> dot_product(const vec3<X, Y, Z> &lhs, const vec3<VX, VY, VZ> &rhs)
> {
> return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z;
> }

 template
<
 typename X, typename Y, typename Z,
 typename VX, typename VY, typename VZ

 typename macc<X, Y, Z, VX, VY, VZ>::type
 dot_product(const vec3<X, Y, Z> &lhs, const vec3<VX, VY, VZ> &rhs)
 {
 return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z;
 }

>
> blas::vec3<int, short, float> a(1, 0, 0);
> blas::vec3<double, int, float> b(0, 1, 0);
> double result = dot_product(a, b); // returns double
>
> I left out PQS types there for simplicity, but I think you can see
> where I'm going.

I think so...

 If the result_of struct isn't specialized for the
> types you are trying to combine, then you get a compile-time error.
> So given:
>
> blas::vec3<length::m> len(1, 2, 3);
> blas::vec3<velocity::s> vel(5, 6, 7);
> // Note the result type
> blas::vec3<velocity::m_per_s2> result = len / (vel * vel);
>
> You only have to specialize for the very base operations between
> types.  Everything more complex (matrix multiplies, etc) can be built
> on top of that.

OK. That sounds much better than what I have done, which was to write it out 
long hand. I
implemented a 4 x 4 matrix as a tuple of tuples FWIW. But I'm not much of an 
expert as I said.

The matrix and vect stuff is in  in pqs_3_1_1  in <boost/pqs/three_d/xx.hpp>
directory to see the gory details FWIW

> For the above, operator / would need to have been
> implemented to return boost::result_of_divides<length::m,
> velocity::s2>::type, while operator * would need to have been
> implemented to return boost::result_of_multiplies<velocity::s,
> velocity::s>::type.  I see you were using BOOST_AUTO in some other
> posts which would probably clean some of my code up, but I'm not
> familiar enough with it to use in my examples.

> This would require the math library be more heavily templated than it
> currently is.  Note that a matrix would need to hold mutliple types
> (up to NxN).  For instance, what would be the syntax for a 3D
> transform?  Perhaps something like:
>
> matrix4x4
> <
>   float, float, float, length,
>   float, float, float, length,
>   float, float, float, length,
>   length, length, length, float
>> m;
>
> maybe?

The above shape would  I am reasonably sure be ( or the other way dependent on
if its a R*C or C*R matrix IIRC)

 matrix4x4
 <
   float, float, float, typeof(1/ length()),
   float, float, float, typeof(1/ length()),
   float, float, float,  typeof(1/ length()),
   length, length, length, float
> m;

If you work through the calcs using dimensionally analysis you find that
concatenating two matrices then gives you the same type of matrix (give or take
promotion ), which is kind of nice! Of course you can generalise the float by
replacing it with typeof(length()/length()) too. (Then finally replace length by
T .... Stick  a single template param T on the front, giving matrix4x4<T> and
Bobs your uncle!

That is the basic idea and similarly for complex vect etc. The only one I havent
tried is quaternion, but hopefully I can try it sometime!

regards
Andy Little