boost/qvm/mat_operations.hpp

0001 #ifndef BOOST_QVM_MAT_OPERATIONS_HPP_INCLUDED
0002 #define BOOST_QVM_MAT_OPERATIONS_HPP_INCLUDED
0003
0004 // Copyright 2008-2022 Emil Dotchevski and Reverge Studios, Inc.
0005 // Copyright 2019 agate-pris
0006 // Distributed under the Boost Software License, Version 1.0. (See accompanying
0007 // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0008
0009 #include <boost/qvm/detail/mat_assign.hpp>
0010 #include <boost/qvm/mat_operations2.hpp>
0011 #include <boost/qvm/mat_operations3.hpp>
0012 #include <boost/qvm/mat_operations4.hpp>
0013 #include <boost/qvm/math.hpp>
0014 #include <boost/qvm/detail/determinant_impl.hpp>
0015 #include <boost/qvm/detail/cofactor_impl.hpp>
0016 #include <boost/qvm/detail/transp_impl.hpp>
0017 #include <boost/qvm/scalar_traits.hpp>
0018 #include <boost/qvm/to_string.hpp>
0019
0020 namespace boost { namespace qvm {
0021
0022 namespace
0023 qvm_detail
0024     {
0025     BOOST_QVM_INLINE_CRITICAL
0026     void const *
0027     get_valid_ptr_mat_operations()
0028         {
0029         static int const obj=0;
0030         return &obj;
0031         }
0032     }
0033
0034 ////////////////////////////////////////////////
0035
0036
0037 namespace
0038 qvm_detail
0039     {
0040     template <int R,int C>
0041     struct
0042     to_string_m_defined
0043         {
0044         static bool const value=false;
0045         };
0046
0047     template <int I,int SizeMinusOne>
0048     struct
0049     to_string_matrix_elements
0050         {
0051         template <class A>
0052         static
0053         std::string
0054         f( A const & a )
0055             {
0056             using namespace qvm_to_string_detail;
0057             return
0058                 ( (I%mat_traits<A>::cols)==0 ? '(' : ',' ) +
0059                 to_string(mat_traits<A>::template read_element<I/mat_traits<A>::cols,I%mat_traits<A>::cols>(a)) +
0060                 ( (I%mat_traits<A>::cols)==mat_traits<A>::cols-1 ? ")" : "" ) +
0061                 to_string_matrix_elements<I+1,SizeMinusOne>::f(a);
0062             }
0063         };
0064
0065     template <int SizeMinusOne>
0066     struct
0067     to_string_matrix_elements<SizeMinusOne,SizeMinusOne>
0068         {
0069         template <class A>
0070         static
0071         std::string
0072         f( A const & a )
0073             {
0074             using namespace qvm_to_string_detail;
0075             return
0076                 ( (SizeMinusOne%mat_traits<A>::cols)==0 ? '(' : ',' ) +
0077                 to_string(mat_traits<A>::template read_element<SizeMinusOne/mat_traits<A>::cols,SizeMinusOne%mat_traits<A>::cols>(a)) +
0078                 ')';
0079             }
0080         };
0081     }
0082
0083 template <class A>
0084 inline
0085 typename enable_if_c<
0086     is_mat<A>::value  &&
0087     !qvm_detail::to_string_m_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0088     std::string>::type
0089 to_string( A const & a )
0090     {
0091     return "("+qvm_detail::to_string_matrix_elements<0,mat_traits<A>::rows*mat_traits<A>::cols-1>::f(a)+')';
0092     }
0093
0094 ////////////////////////////////////////////////
0095
0096 template <class A,class B,class Cmp>
0097 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0098 typename enable_if_c<
0099     is_mat<A>::value && is_mat<B>::value &&
0100     mat_traits<A>::rows==mat_traits<B>::rows &&
0101     mat_traits<A>::cols==mat_traits<B>::cols,
0102     bool>::type
0103 cmp( A const & a, B const & b, Cmp f )
0104     {
0105     for( int i=0; i!=mat_traits<A>::rows; ++i )
0106         for( int j=0; j!=mat_traits<A>::cols; ++j )
0107             if( !f(
0108                 mat_traits<A>::read_element_idx(i, j, a),
0109                 mat_traits<B>::read_element_idx(i, j, b)) )
0110                 return false;
0111     return true;
0112     }
0113
0114 ////////////////////////////////////////////////
0115
0116 namespace
0117 qvm_detail
0118     {
0119     template <int M,int N>
0120     struct
0121     convert_to_m_defined
0122         {
0123         static bool const value=false;
0124         };
0125     }
0126
0127 template <class R,class A>
0128 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0129 typename enable_if_c<
0130     is_mat<R>::value && is_mat<A>::value &&
0131     mat_traits<R>::rows==mat_traits<A>::rows &&
0132     mat_traits<R>::cols==mat_traits<A>::cols &&
0133     !qvm_detail::convert_to_m_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0134     R>::type
0135 convert_to( A const & a )
0136     {
0137     R r; assign(r,a);
0138     return r;
0139     }
0140
0141 ////////////////////////////////////////////////
0142
0143 namespace
0144 qvm_detail
0145     {
0146     template <int D>
0147     struct
0148     determinant_defined
0149         {
0150         static bool const value=false;
0151         };
0152     }
0153
0154 template <class A>
0155 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0156 typename enable_if_c<
0157     is_mat<A>::value &&
0158     mat_traits<A>::rows==mat_traits<A>::cols &&
0159     !qvm_detail::determinant_defined<mat_traits<A>::rows>::value,
0160     typename mat_traits<A>::scalar_type>::type
0161 determinant( A const & a )
0162     {
0163     return qvm_detail::determinant_impl(a);
0164     }
0165
0166 ////////////////////////////////////////////////
0167
0168 namespace
0169 qvm_detail
0170     {
0171     template <class T,int Dim>
0172     class
0173     identity_mat_
0174         {
0175         identity_mat_( identity_mat_ const & );
0176         identity_mat_ & operator=( identity_mat_ const & );
0177         ~identity_mat_();
0178
0179         public:
0180
0181         template <class R
0182 #if __cplusplus >= 201103L
0183             , class = typename enable_if<is_mat<R> >::type
0184 #endif
0185         >
0186         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0187         operator R() const
0188             {
0189             R r;
0190             assign(r,*this);
0191             return r;
0192             }
0193         };
0194     }
0195
0196 template <class T,int Dim>
0197 struct
0198 mat_traits< qvm_detail::identity_mat_<T,Dim> >
0199     {
0200     typedef qvm_detail::identity_mat_<T,Dim> this_matrix;
0201     typedef T scalar_type;
0202     static int const rows=Dim;
0203     static int const cols=Dim;
0204
0205     template <int Row,int Col>
0206     static
0207     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0208     scalar_type
0209     read_element( this_matrix const & /*x*/ )
0210         {
0211         BOOST_QVM_STATIC_ASSERT(Row>=0);
0212         BOOST_QVM_STATIC_ASSERT(Row<Dim);
0213         BOOST_QVM_STATIC_ASSERT(Col>=0);
0214         BOOST_QVM_STATIC_ASSERT(Col<Dim);
0215         return scalar_traits<scalar_type>::value(Row==Col);
0216         }
0217
0218     static
0219     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0220     scalar_type
0221     read_element_idx( int row, int col, this_matrix const & /*x*/ )
0222         {
0223         BOOST_QVM_ASSERT(row>=0);
0224         BOOST_QVM_ASSERT(row<Dim);
0225         BOOST_QVM_ASSERT(col>=0);
0226         BOOST_QVM_ASSERT(col<Dim);
0227         return scalar_traits<scalar_type>::value(row==col);
0228         }
0229     };
0230
0231 template <class T,int Dim>
0232 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0233 qvm_detail::identity_mat_<T,Dim> const &
0234 identity_mat()
0235     {
0236     return *(qvm_detail::identity_mat_<T,Dim> const *)qvm_detail::get_valid_ptr_mat_operations();
0237     }
0238
0239 template <class A>
0240 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0241 typename enable_if_c<
0242     is_mat<A>::value &&
0243     mat_traits<A>::rows==mat_traits<A>::cols,
0244     void>::type
0245 set_identity( A & a )
0246     {
0247     assign(a,identity_mat<typename mat_traits<A>::scalar_type,mat_traits<A>::rows>());
0248     }
0249
0250 ////////////////////////////////////////////////
0251
0252 namespace
0253 qvm_detail
0254     {
0255     template <class T>
0256     struct
0257     projection_
0258         {
0259         T const _00;
0260         T const _11;
0261         T const _22;
0262         T const _23;
0263         T const _32;
0264
0265         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0266         projection_( T _00, T _11, T _22, T _23, T _32 ):
0267             _00(_00),
0268             _11(_11),
0269             _22(_22),
0270             _23(_23),
0271             _32(_32)
0272             {
0273             }
0274
0275         template <class R
0276 #if __cplusplus >= 201103L
0277             , class = typename enable_if<is_mat<R> >::type
0278 #endif
0279         >
0280         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0281         operator R() const
0282             {
0283             R r;
0284             assign(r,*this);
0285             return r;
0286             }
0287         };
0288
0289     template <int Row,int Col>
0290     struct
0291     projection_get
0292         {
0293         template <class T>
0294         static
0295         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0296         T
0297         get( projection_<T> const & )
0298             {
0299             return scalar_traits<T>::value(0);
0300             }
0301         };
0302
0303     template <> struct projection_get<0,0> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( projection_<T> const & m ) { return m._00; } };
0304     template <> struct projection_get<1,1> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( projection_<T> const & m ) { return m._11; } };
0305     template <> struct projection_get<2,2> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( projection_<T> const & m ) { return m._22; } };
0306     template <> struct projection_get<2,3> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( projection_<T> const & m ) { return m._23; } };
0307     template <> struct projection_get<3,2> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( projection_<T> const & m ) { return m._32; } };
0308     }
0309
0310 template <class T>
0311 struct
0312 mat_traits< qvm_detail::projection_<T> >
0313     {
0314     typedef qvm_detail::projection_<T> this_matrix;
0315     typedef T scalar_type;
0316     static int const rows=4;
0317     static int const cols=4;
0318
0319     template <int Row,int Col>
0320     static
0321     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0322     scalar_type
0323     read_element( this_matrix const & x )
0324         {
0325         BOOST_QVM_STATIC_ASSERT(Row>=0);
0326         BOOST_QVM_STATIC_ASSERT(Row<rows);
0327         BOOST_QVM_STATIC_ASSERT(Col>=0);
0328         BOOST_QVM_STATIC_ASSERT(Col<cols);
0329         return qvm_detail::projection_get<Row,Col>::get(x);
0330         }
0331     };
0332
0333 template <class T>
0334 qvm_detail::projection_<T>
0335 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0336 perspective_lh( T fov_y, T aspect_ratio, T z_near, T z_far )
0337     {
0338     T const one = scalar_traits<T>::value(1);
0339     T const ys = one/tan(fov_y/scalar_traits<T>::value(2));
0340     T const xs = ys/aspect_ratio;
0341     T const zd = z_far-z_near;
0342     T const z1 = z_far/zd;
0343     T const z2 = -z_near*z1;
0344     return qvm_detail::projection_<T>(xs,ys,z1,z2,one);
0345     }
0346
0347 template <class T>
0348 qvm_detail::projection_<T>
0349 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0350 perspective_rh( T fov_y, T aspect_ratio, T z_near, T z_far )
0351     {
0352     T const one = scalar_traits<T>::value(1);
0353     T const ys = one/tan(fov_y/scalar_traits<T>::value(2));
0354     T const xs = ys/aspect_ratio;
0355     T const zd = z_near-z_far;
0356     T const z1 = z_far/zd;
0357     T const z2 = z_near*z1;
0358     return qvm_detail::projection_<T>(xs,ys,z1,z2,-one);
0359     }
0360
0361 ////////////////////////////////////////////////
0362
0363 namespace
0364 qvm_detail
0365     {
0366     template <class OriginalType,class Scalar>
0367     class
0368     matrix_scalar_cast_
0369         {
0370         matrix_scalar_cast_( matrix_scalar_cast_ const & );
0371         matrix_scalar_cast_ & operator=( matrix_scalar_cast_ const & );
0372         ~matrix_scalar_cast_();
0373
0374         public:
0375
0376         template <class T>
0377         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0378         matrix_scalar_cast_ &
0379         operator=( T const & x )
0380             {
0381             assign(*this,x);
0382             return *this;
0383             }
0384
0385         template <class R
0386 #if __cplusplus >= 201103L
0387             , class = typename enable_if<is_mat<R> >::type
0388 #endif
0389         >
0390         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0391         operator R() const
0392             {
0393             R r;
0394             assign(r,*this);
0395             return r;
0396             }
0397         };
0398
0399     template <bool> struct scalar_cast_matrix_filter { };
0400     template <> struct scalar_cast_matrix_filter<true> { typedef int type; };
0401     }
0402
0403 template <class OriginalType,class Scalar>
0404 struct
0405 mat_traits< qvm_detail::matrix_scalar_cast_<OriginalType,Scalar> >
0406     {
0407     typedef Scalar scalar_type;
0408     typedef qvm_detail::matrix_scalar_cast_<OriginalType,Scalar> this_matrix;
0409     static int const rows=mat_traits<OriginalType>::rows;
0410     static int const cols=mat_traits<OriginalType>::cols;
0411
0412     template <int Row,int Col>
0413     static
0414     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0415     scalar_type
0416     read_element( this_matrix const & x )
0417         {
0418         BOOST_QVM_STATIC_ASSERT(Row>=0);
0419         BOOST_QVM_STATIC_ASSERT(Row<rows);
0420         BOOST_QVM_STATIC_ASSERT(Col>=0);
0421         BOOST_QVM_STATIC_ASSERT(Col<cols);
0422         return scalar_type(mat_traits<OriginalType>::template read_element<Row,Col>(reinterpret_cast<OriginalType const &>(x)));
0423         }
0424
0425     static
0426     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0427     scalar_type
0428     read_element_idx( int row, int col, this_matrix const & x )
0429         {
0430         BOOST_QVM_ASSERT(row>=0);
0431         BOOST_QVM_ASSERT(row<rows);
0432         BOOST_QVM_ASSERT(col>=0);
0433         BOOST_QVM_ASSERT(col<cols);
0434         return scalar_type(mat_traits<OriginalType>::read_element_idx(col,row,reinterpret_cast<OriginalType const &>(x)));
0435         }
0436     };
0437
0438 template <class OriginalType,class Scalar,int R,int C>
0439 struct
0440 deduce_mat<qvm_detail::matrix_scalar_cast_<OriginalType,Scalar>,R,C>
0441     {
0442     typedef mat<Scalar,R,C> type;
0443     };
0444
0445 template <class Scalar,class T>
0446 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0447 qvm_detail::matrix_scalar_cast_<T,Scalar> const &
0448 scalar_cast( T const & x, typename qvm_detail::scalar_cast_matrix_filter<is_mat<T>::value>::type=0 )
0449     {
0450     return reinterpret_cast<qvm_detail::matrix_scalar_cast_<T,Scalar> const &>(x);
0451     }
0452
0453 ////////////////////////////////////////////////
0454
0455 namespace
0456 qvm_detail
0457     {
0458     template <int M,int N>
0459     struct
0460     div_eq_ms_defined
0461         {
0462         static bool const value=false;
0463         };
0464     }
0465
0466 template <class A,class B>
0467 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0468 typename enable_if_c<
0469     is_mat<A>::value && is_scalar<B>::value &&
0470     !qvm_detail::div_eq_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0471     A &>::type
0472 operator/=( A & a, B b )
0473     {
0474     for( int i=0; i!=mat_traits<A>::rows; ++i )
0475         for( int j=0; j!=mat_traits<A>::cols; ++j )
0476             write_mat_element_idx(i,j,a,mat_traits<A>::read_element_idx(i,j,a)/b);
0477     return a;
0478     }
0479
0480 ////////////////////////////////////////////////
0481
0482 namespace
0483 qvm_detail
0484     {
0485     template <int M,int N>
0486     struct
0487     div_ms_defined
0488         {
0489         static bool const value=false;
0490         };
0491     }
0492
0493 template <class A,class B>
0494 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0495 typename lazy_enable_if_c<
0496     is_mat<A>::value && is_scalar<B>::value &&
0497     !qvm_detail::div_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0498     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
0499 operator/( A const & a, B b )
0500     {
0501     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
0502     R r;
0503     for( int i=0; i!=mat_traits<A>::rows; ++i )
0504         for( int j=0; j!=mat_traits<A>::cols; ++j )
0505             write_mat_element_idx(i,j,r,mat_traits<A>::read_element_idx(i,j,a)/b);
0506     return r;
0507     }
0508
0509 ////////////////////////////////////////////////
0510
0511 namespace
0512 qvm_detail
0513     {
0514     template <int M,int N>
0515     struct
0516     eq_mm_defined
0517         {
0518         static bool const value=false;
0519         };
0520     }
0521
0522 template <class A,class B>
0523 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0524 typename enable_if_c<
0525     is_mat<A>::value && is_mat<B>::value &&
0526     mat_traits<A>::rows==mat_traits<B>::rows &&
0527     mat_traits<A>::cols==mat_traits<B>::cols &&
0528     !qvm_detail::eq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0529     bool>::type
0530 operator==( A const & a, B const & b )
0531     {
0532     for( int i=0; i!=mat_traits<A>::rows; ++i )
0533         for( int j=0; j!=mat_traits<A>::cols; ++j )
0534             if( mat_traits<A>::read_element_idx(i,j,a)!=mat_traits<B>::read_element_idx(i,j,b) )
0535                 return false;
0536     return true;
0537     }
0538
0539 ////////////////////////////////////////////////
0540
0541 namespace
0542 qvm_detail
0543     {
0544     template <int M,int N>
0545     struct
0546     minus_eq_mm_defined
0547         {
0548         static bool const value=false;
0549         };
0550     }
0551
0552 template <class A,class B>
0553 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0554 typename enable_if_c<
0555     is_mat<A>::value && is_mat<B>::value &&
0556     mat_traits<A>::rows==mat_traits<B>::rows &&
0557     mat_traits<A>::cols==mat_traits<B>::cols &&
0558     !qvm_detail::minus_eq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0559     A &>::type
0560 operator-=( A & a, B const & b )
0561     {
0562     for( int i=0; i!=mat_traits<A>::rows; ++i )
0563         for( int j=0; j!=mat_traits<A>::cols; ++j )
0564             write_mat_element_idx(i,j,a,mat_traits<A>::read_element_idx(i,j,a)-mat_traits<B>::read_element_idx(i,j,b));
0565     return a;
0566     }
0567
0568 ////////////////////////////////////////////////
0569
0570 namespace
0571 qvm_detail
0572     {
0573     template <int M,int N>
0574     struct
0575     minus_m_defined
0576         {
0577         static bool const value=false;
0578         };
0579     }
0580
0581 template <class A>
0582 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0583 typename lazy_enable_if_c<
0584     is_mat<A>::value &&
0585     !qvm_detail::minus_m_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0586     deduce_mat<A> >::type
0587 operator-( A const & a )
0588     {
0589     typedef typename deduce_mat<A>::type R;
0590     R r;
0591     for( int i=0; i!=mat_traits<A>::rows; ++i )
0592         for( int j=0; j!=mat_traits<A>::cols; ++j )
0593             write_mat_element_idx(i,j,r,-mat_traits<A>::read_element_idx(i,j,a));
0594     return r;
0595     }
0596
0597 ////////////////////////////////////////////////
0598
0599 namespace
0600 qvm_detail
0601     {
0602     template <int M,int N>
0603     struct
0604     minus_mm_defined
0605         {
0606         static bool const value=false;
0607         };
0608     }
0609
0610 template <class A,class B>
0611 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0612 typename lazy_enable_if_c<
0613     is_mat<A>::value && is_mat<B>::value &&
0614     mat_traits<A>::rows==mat_traits<B>::rows &&
0615     mat_traits<A>::cols==mat_traits<B>::cols &&
0616     !qvm_detail::minus_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0617     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
0618 operator-( A const & a, B const & b )
0619     {
0620     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
0621     R r;
0622     for( int i=0; i!=mat_traits<A>::rows; ++i )
0623         for( int j=0; j!=mat_traits<A>::cols; ++j )
0624             write_mat_element_idx(i,j,r,mat_traits<A>::read_element_idx(i,j,a)-mat_traits<B>::read_element_idx(i,j,b));
0625     return r;
0626     }
0627
0628 ////////////////////////////////////////////////
0629
0630 namespace
0631 qvm_detail
0632     {
0633     template <int D>
0634     struct
0635     mul_eq_mm_defined
0636         {
0637         static bool const value=false;
0638         };
0639     }
0640
0641 template <class A,class B>
0642 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0643 typename enable_if_c<
0644     is_mat<A>::value &&
0645     is_mat<B>::value &&
0646     mat_traits<A>::rows==mat_traits<A>::cols &&
0647     mat_traits<A>::rows==mat_traits<B>::rows &&
0648     mat_traits<A>::cols==mat_traits<B>::cols &&
0649     !qvm_detail::mul_eq_mm_defined<mat_traits<A>::rows>::value,
0650     A &>::type
0651 operator*=( A & r, B const & b )
0652     {
0653     typedef typename mat_traits<A>::scalar_type Ta;
0654     Ta a[mat_traits<A>::rows][mat_traits<A>::cols];
0655     for( int i=0; i<mat_traits<A>::rows; ++i )
0656         for( int j=0; j<mat_traits<B>::cols; ++j )
0657             a[i][j]=mat_traits<A>::read_element_idx(i,j,r);
0658     for( int i=0; i<mat_traits<A>::rows; ++i )
0659         for( int j=0; j<mat_traits<B>::cols; ++j )
0660             {
0661             Ta x(scalar_traits<Ta>::value(0));
0662             for( int k=0; k<mat_traits<A>::cols; ++k )
0663                 x += a[i][k]*mat_traits<B>::read_element_idx(k,j,b);
0664             write_mat_element_idx(i,j,r,x);
0665             }
0666     return r;
0667     }
0668
0669 ////////////////////////////////////////////////
0670
0671 namespace
0672 qvm_detail
0673     {
0674     template <int M,int N>
0675     struct
0676     mul_eq_ms_defined
0677         {
0678         static bool const value=false;
0679         };
0680     }
0681
0682 template <class A,class B>
0683 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0684 typename enable_if_c<
0685     is_mat<A>::value && is_scalar<B>::value &&
0686     !qvm_detail::mul_eq_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0687     A &>::type
0688 operator*=( A & a, B b )
0689     {
0690     for( int i=0; i!=mat_traits<A>::rows; ++i )
0691         for( int j=0; j!=mat_traits<A>::cols; ++j )
0692             write_mat_element_idx(i,j,a,mat_traits<A>::read_element_idx(i,j,a)*b);
0693     return a;
0694     }
0695
0696 ////////////////////////////////////////////////
0697
0698 namespace
0699 qvm_detail
0700     {
0701     template <int R,int /*CR*/,int C>
0702     struct
0703     mul_mm_defined
0704         {
0705         static bool const value=false;
0706         };
0707     }
0708
0709 template <class A,class B>
0710 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0711 typename lazy_enable_if_c<
0712     is_mat<A>::value && is_mat<B>::value &&
0713     mat_traits<A>::cols==mat_traits<B>::rows &&
0714     !qvm_detail::mul_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols,mat_traits<B>::cols>::value,
0715     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<B>::cols> >::type
0716 operator*( A const & a, B const & b )
0717     {
0718     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<B>::cols>::type R;
0719     R r;
0720     for( int i=0; i<mat_traits<A>::rows; ++i )
0721         for( int j=0; j<mat_traits<B>::cols; ++j )
0722             {
0723             typedef typename mat_traits<A>::scalar_type Ta;
0724             Ta x(scalar_traits<Ta>::value(0));
0725             for( int k=0; k<mat_traits<A>::cols; ++k )
0726                 x += mat_traits<A>::read_element_idx(i,k,a)*mat_traits<B>::read_element_idx(k,j,b);
0727             write_mat_element_idx(i,j,r,x);
0728             }
0729     return r;
0730     }
0731
0732 ////////////////////////////////////////////////
0733
0734 namespace
0735 qvm_detail
0736     {
0737     template <int M,int N>
0738     struct
0739     mul_ms_defined
0740         {
0741         static bool const value=false;
0742         };
0743     }
0744
0745 template <class A,class B>
0746 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0747 typename lazy_enable_if_c<
0748     is_mat<A>::value && is_scalar<B>::value &&
0749     !qvm_detail::mul_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0750     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
0751 operator*( A const & a, B b )
0752     {
0753     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
0754     R r;
0755     for( int i=0; i!=mat_traits<A>::rows; ++i )
0756         for( int j=0; j!=mat_traits<A>::cols; ++j )
0757             write_mat_element_idx(i,j,r,mat_traits<A>::read_element_idx(i,j,a)*b);
0758     return r;
0759     }
0760
0761 ////////////////////////////////////////////////
0762
0763 namespace
0764 qvm_detail
0765     {
0766     template <int M,int N>
0767     struct
0768     mul_sm_defined
0769         {
0770         static bool const value=false;
0771         };
0772     }
0773
0774 template <class A,class B>
0775 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0776 typename lazy_enable_if_c<
0777     is_scalar<A>::value && is_mat<B>::value &&
0778     !qvm_detail::mul_sm_defined<mat_traits<B>::rows,mat_traits<B>::cols>::value,
0779     deduce_mat2<A,B,mat_traits<B>::rows,mat_traits<B>::cols> >::type
0780 operator*( A a, B const & b )
0781     {
0782     typedef typename deduce_mat2<A,B,mat_traits<B>::rows,mat_traits<B>::cols>::type R;
0783     R r;
0784     for( int i=0; i!=mat_traits<B>::rows; ++i )
0785         for( int j=0; j!=mat_traits<B>::cols; ++j )
0786             write_mat_element_idx(i,j,r,a*mat_traits<B>::read_element_idx(i,j,b));
0787     return r;
0788     }
0789
0790 ////////////////////////////////////////////////
0791
0792 namespace
0793 qvm_detail
0794     {
0795     template <int M,int N>
0796     struct
0797     neq_mm_defined
0798         {
0799         static bool const value=false;
0800         };
0801     }
0802
0803 template <class A,class B>
0804 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0805 typename enable_if_c<
0806     is_mat<A>::value && is_mat<B>::value &&
0807     mat_traits<A>::rows==mat_traits<B>::rows &&
0808     mat_traits<A>::cols==mat_traits<B>::cols &&
0809     !qvm_detail::neq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0810     bool>::type
0811 operator!=( A const & a, B const & b )
0812     {
0813     for( int i=0; i!=mat_traits<A>::rows; ++i )
0814         for( int j=0; j!=mat_traits<A>::cols; ++j )
0815             if( mat_traits<A>::read_element_idx(i,j,a)!=mat_traits<B>::read_element_idx(i,j,b) )
0816                 return true;
0817     return false;
0818     }
0819
0820 ////////////////////////////////////////////////
0821
0822 namespace
0823 qvm_detail
0824     {
0825     template <int M,int N>
0826     struct
0827     plus_eq_mm_defined
0828         {
0829         static bool const value=false;
0830         };
0831     }
0832
0833 template <class A,class B>
0834 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0835 typename enable_if_c<
0836     is_mat<A>::value && is_mat<B>::value &&
0837     mat_traits<A>::rows==mat_traits<B>::rows &&
0838     mat_traits<A>::cols==mat_traits<B>::cols &&
0839     !qvm_detail::plus_eq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0840     A &>::type
0841 operator+=( A & a, B const & b )
0842     {
0843     for( int i=0; i!=mat_traits<A>::rows; ++i )
0844         for( int j=0; j!=mat_traits<A>::cols; ++j )
0845             write_mat_element_idx(i,j,a,mat_traits<A>::read_element_idx(i,j,a)+mat_traits<B>::read_element_idx(i,j,b));
0846     return a;
0847     }
0848
0849 ////////////////////////////////////////////////
0850
0851 namespace
0852 qvm_detail
0853     {
0854     template <int M,int N>
0855     struct
0856     plus_mm_defined
0857         {
0858         static bool const value=false;
0859         };
0860     }
0861
0862 template <class A,class B>
0863 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
0864 typename lazy_enable_if_c<
0865     is_mat<A>::value && is_mat<B>::value &&
0866     mat_traits<A>::rows==mat_traits<B>::rows &&
0867     mat_traits<A>::cols==mat_traits<B>::cols &&
0868     !qvm_detail::plus_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
0869     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
0870 operator+( A const & a, B const & b )
0871     {
0872     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
0873     R r;
0874     for( int i=0; i!=mat_traits<A>::rows; ++i )
0875         for( int j=0; j!=mat_traits<A>::cols; ++j )
0876             write_mat_element_idx(i,j,r,mat_traits<A>::read_element_idx(i,j,a)+mat_traits<B>::read_element_idx(i,j,b));
0877     return r;
0878     }
0879
0880 ////////////////////////////////////////////////
0881
0882 namespace
0883 qvm_detail
0884     {
0885     template <class T>
0886     class
0887     mref_
0888         {
0889         mref_( mref_ const & );
0890         mref_ & operator=( mref_ const & );
0891         ~mref_();
0892
0893         public:
0894
0895         template <class R>
0896         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0897         mref_ &
0898         operator=( R const & x )
0899             {
0900             assign(*this,x);
0901             return *this;
0902             }
0903
0904         template <class R
0905 #if __cplusplus >= 201103L
0906             , class = typename enable_if<is_mat<R> >::type
0907 #endif
0908         >
0909         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
0910         operator R() const
0911             {
0912             R r;
0913             assign(r,*this);
0914             return r;
0915             }
0916         };
0917
0918     template <class M,bool WriteElementRef=mat_write_element_ref<M>::value>
0919     struct mref_write_traits;
0920
0921     template <class M>
0922     struct
0923     mref_write_traits<M,true>
0924         {
0925         typedef typename mat_traits<M>::scalar_type scalar_type;
0926         typedef qvm_detail::mref_<M> this_matrix;
0927         static int const rows=mat_traits<M>::rows;
0928         static int const cols=mat_traits<M>::cols;
0929
0930         template <int Row,int Col>
0931         static
0932         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0933         scalar_type &
0934         write_element( this_matrix & x )
0935             {
0936             BOOST_QVM_STATIC_ASSERT(Row>=0);
0937             BOOST_QVM_STATIC_ASSERT(Row<rows);
0938             BOOST_QVM_STATIC_ASSERT(Col>=0);
0939             BOOST_QVM_STATIC_ASSERT(Col<cols);
0940             return mat_traits<M>::template write_element<Row,Col>(reinterpret_cast<M &>(x));
0941             }
0942
0943         static
0944         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0945         scalar_type &
0946         write_element_idx( int row, int col, this_matrix & x )
0947             {
0948             BOOST_QVM_ASSERT(row>=0);
0949             BOOST_QVM_ASSERT(row<rows);
0950             BOOST_QVM_ASSERT(col>=0);
0951             BOOST_QVM_ASSERT(col<cols);
0952             return mat_traits<M>::write_element_idx(row,col,reinterpret_cast<M &>(x));
0953             }
0954         };
0955
0956     template <class M>
0957     struct
0958     mref_write_traits<M,false>
0959         {
0960         typedef typename mat_traits<M>::scalar_type scalar_type;
0961         typedef qvm_detail::mref_<M> this_matrix;
0962         static int const rows=mat_traits<M>::rows;
0963         static int const cols=mat_traits<M>::cols;
0964
0965         template <int Row,int Col>
0966         static
0967         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0968         void
0969         write_element( this_matrix & x, scalar_type s )
0970             {
0971             BOOST_QVM_STATIC_ASSERT(Row>=0);
0972             BOOST_QVM_STATIC_ASSERT(Row<rows);
0973             BOOST_QVM_STATIC_ASSERT(Col>=0);
0974             BOOST_QVM_STATIC_ASSERT(Col<cols);
0975             mat_traits<M>::template write_element<Row,Col>(reinterpret_cast<M &>(x), s);
0976             }
0977
0978         static
0979         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
0980         void
0981         write_element_idx( int row, int col, this_matrix & x, scalar_type s )
0982             {
0983             BOOST_QVM_ASSERT(row>=0);
0984             BOOST_QVM_ASSERT(row<rows);
0985             BOOST_QVM_ASSERT(col>=0);
0986             BOOST_QVM_ASSERT(col<cols);
0987             mat_traits<M>::write_element_idx(row,col,reinterpret_cast<M &>(x), s);
0988             }
0989         };
0990     }
0991
0992 template <class M>
0993 struct
0994 mat_traits< qvm_detail::mref_<M> >:
0995     qvm_detail::mref_write_traits<M>
0996     {
0997     typedef typename mat_traits<M>::scalar_type scalar_type;
0998     typedef qvm_detail::mref_<M> this_matrix;
0999     static int const rows=mat_traits<M>::rows;
1000     static int const cols=mat_traits<M>::cols;
1001
1002     template <int Row,int Col>
1003     static
1004     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1005     scalar_type
1006     read_element( this_matrix const & x )
1007         {
1008         BOOST_QVM_STATIC_ASSERT(Row>=0);
1009         BOOST_QVM_STATIC_ASSERT(Row<rows);
1010         BOOST_QVM_STATIC_ASSERT(Col>=0);
1011         BOOST_QVM_STATIC_ASSERT(Col<cols);
1012         return mat_traits<M>::template read_element<Row,Col>(reinterpret_cast<M const &>(x));
1013         }
1014
1015     static
1016     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1017     scalar_type
1018     read_element_idx( int row, int col, this_matrix const & x )
1019         {
1020         BOOST_QVM_ASSERT(row>=0);
1021         BOOST_QVM_ASSERT(row<rows);
1022         BOOST_QVM_ASSERT(col>=0);
1023         BOOST_QVM_ASSERT(col<cols);
1024         return mat_traits<M>::read_element_idx(row,col,reinterpret_cast<M const &>(x));
1025         }
1026     };
1027
1028 template <class M,int R,int C>
1029 struct
1030 deduce_mat<qvm_detail::mref_<M>,R,C>
1031     {
1032     typedef mat<typename mat_traits<M>::scalar_type,R,C> type;
1033     };
1034
1035 template <class M>
1036 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1037 typename enable_if_c<
1038     is_mat<M>::value,
1039     qvm_detail::mref_<M> const &>::type
1040 mref( M const & a )
1041     {
1042     return reinterpret_cast<qvm_detail::mref_<M> const &>(a);
1043     }
1044
1045 template <class M>
1046 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1047 typename enable_if_c<
1048     is_mat<M>::value,
1049     qvm_detail::mref_<M> &>::type
1050 mref( M & a )
1051     {
1052     return reinterpret_cast<qvm_detail::mref_<M> &>(a);
1053     }
1054
1055 ////////////////////////////////////////////////
1056
1057 namespace
1058 qvm_detail
1059     {
1060     template <class T,int Rows,int Cols>
1061     class
1062     zero_mat_
1063         {
1064         zero_mat_( zero_mat_ const & );
1065         zero_mat_ & operator=( zero_mat_ const & );
1066         ~zero_mat_();
1067
1068         public:
1069
1070         template <class R
1071 #if __cplusplus >= 201103L
1072             , class = typename enable_if<is_mat<R> >::type
1073 #endif
1074         >
1075         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1076         operator R() const
1077             {
1078             R r;
1079             assign(r,*this);
1080             return r;
1081             }
1082         };
1083     }
1084
1085 template <class T,int Rows,int Cols>
1086 struct
1087 mat_traits< qvm_detail::zero_mat_<T,Rows,Cols> >
1088     {
1089     typedef qvm_detail::zero_mat_<T,Rows,Cols> this_matrix;
1090     typedef T scalar_type;
1091     static int const rows=Rows;
1092     static int const cols=Cols;
1093
1094     template <int Row,int Col>
1095     static
1096     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1097     scalar_type
1098     read_element( this_matrix const & )
1099         {
1100         BOOST_QVM_STATIC_ASSERT(Row>=0);
1101         BOOST_QVM_STATIC_ASSERT(Row<rows);
1102         BOOST_QVM_STATIC_ASSERT(Col>=0);
1103         BOOST_QVM_STATIC_ASSERT(Col<cols);
1104         return scalar_traits<scalar_type>::value(0);
1105         }
1106
1107     static
1108     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1109     scalar_type
1110     read_element_idx( int row, int col, this_matrix const & )
1111         {
1112         BOOST_QVM_ASSERT(row>=0); (void)row;
1113         BOOST_QVM_ASSERT(row<rows);
1114         BOOST_QVM_ASSERT(col>=0); (void)col;
1115         BOOST_QVM_ASSERT(col<cols);
1116         return scalar_traits<scalar_type>::value(0);
1117         }
1118     };
1119
1120 template <class T,int Rows,int Cols,int R,int C>
1121 struct
1122 deduce_mat<qvm_detail::zero_mat_<T,Rows,Cols>,R,C>
1123     {
1124     typedef mat<T,R,C> type;
1125     };
1126
1127 template <class T,int Rows,int Cols>
1128 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1129 qvm_detail::zero_mat_<T,Rows,Cols> const &
1130 zero_mat()
1131     {
1132     return *(qvm_detail::zero_mat_<T,Rows,Cols> const *)qvm_detail::get_valid_ptr_mat_operations();
1133     }
1134
1135 template <class T,int Dim>
1136 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1137 qvm_detail::zero_mat_<T,Dim,Dim> const &
1138 zero_mat()
1139     {
1140     return *(qvm_detail::zero_mat_<T,Dim,Dim> const *)qvm_detail::get_valid_ptr_mat_operations();
1141     }
1142
1143 template <class A>
1144 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1145 typename enable_if_c<
1146     is_mat<A>::value,
1147     void>::type
1148 set_zero( A & a )
1149     {
1150     assign(a,zero_mat<typename mat_traits<A>::scalar_type,mat_traits<A>::rows,mat_traits<A>::cols>());
1151     }
1152
1153 ////////////////////////////////////////////////
1154
1155 namespace
1156 qvm_detail
1157     {
1158     template <int D,class S>
1159     struct
1160     rot_mat_
1161         {
1162         typedef S scalar_type;
1163         scalar_type a[3][3];
1164
1165         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1166         rot_mat_(
1167                 scalar_type a00, scalar_type a01, scalar_type a02,
1168                 scalar_type a10, scalar_type a11, scalar_type a12,
1169                 scalar_type a20, scalar_type a21, scalar_type a22 )
1170             {
1171             a[0][0] = a00;
1172             a[0][1] = a01;
1173             a[0][2] = a02;
1174             a[1][0] = a10;
1175             a[1][1] = a11;
1176             a[1][2] = a12;
1177             a[2][0] = a20;
1178             a[2][1] = a21;
1179             a[2][2] = a22;
1180             }
1181
1182         template <class R
1183 #if __cplusplus >= 201103L
1184             , class = typename enable_if<is_mat<R> >::type
1185 #endif
1186         >
1187         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1188         operator R() const
1189             {
1190             R r;
1191             assign(r,*this);
1192             return r;
1193             }
1194         };
1195
1196     template <int Row,int Col>
1197     struct
1198     rot_m_get
1199         {
1200         template <class T>
1201         static
1202         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1203         T
1204         get( T const (&)[3][3] )
1205             {
1206             return scalar_traits<T>::value(Row==Col);
1207             }
1208         };
1209
1210     template <> struct rot_m_get<0,0> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[0][0]; } };
1211     template <> struct rot_m_get<0,1> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[0][1]; } };
1212     template <> struct rot_m_get<0,2> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[0][2]; } };
1213     template <> struct rot_m_get<1,0> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[1][0]; } };
1214     template <> struct rot_m_get<1,1> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[1][1]; } };
1215     template <> struct rot_m_get<1,2> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[1][2]; } };
1216     template <> struct rot_m_get<2,0> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[2][0]; } };
1217     template <> struct rot_m_get<2,1> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[2][1]; } };
1218     template <> struct rot_m_get<2,2> { template <class T> static BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL T get( T const (&a)[3][3] ) { return a[2][2]; } };
1219     }
1220
1221 template <class M>
1222 struct mat_traits;
1223
1224 template <int D,class S>
1225 struct
1226 mat_traits< qvm_detail::rot_mat_<D,S> >
1227     {
1228     typedef qvm_detail::rot_mat_<D,S> this_matrix;
1229     typedef typename this_matrix::scalar_type scalar_type;
1230     static int const rows=D;
1231     static int const cols=D;
1232
1233     template <int Row,int Col>
1234     static
1235     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1236     scalar_type
1237     read_element( this_matrix const & x )
1238         {
1239         BOOST_QVM_STATIC_ASSERT(Row>=0);
1240         BOOST_QVM_STATIC_ASSERT(Row<D);
1241         BOOST_QVM_STATIC_ASSERT(Col>=0);
1242         BOOST_QVM_STATIC_ASSERT(Col<D);
1243         return qvm_detail::rot_m_get<Row,Col>::get(x.a);
1244         }
1245
1246     static
1247     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1248     scalar_type
1249     read_element_idx( int row, int col, this_matrix const & x )
1250         {
1251         BOOST_QVM_ASSERT(row>=0);
1252         BOOST_QVM_ASSERT(row<D);
1253         BOOST_QVM_ASSERT(col>=0);
1254         BOOST_QVM_ASSERT(col<D);
1255         return row<3 && col<3?
1256             x.a[row][col] :
1257             scalar_traits<scalar_type>::value(row==col);
1258         }
1259     };
1260
1261 template <int Dim,class V,class Angle>
1262 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1263 typename enable_if_c<
1264     is_vec<V>::value && vec_traits<V>::dim==3,
1265     qvm_detail::rot_mat_<Dim,Angle> >::type
1266 rot_mat( V const & axis, Angle angle )
1267     {
1268     typedef Angle scalar_type;
1269     scalar_type const x=vec_traits<V>::template read_element<0>(axis);
1270     scalar_type const y=vec_traits<V>::template read_element<1>(axis);
1271     scalar_type const z=vec_traits<V>::template read_element<2>(axis);
1272     scalar_type const m2=x*x+y*y+z*z;
1273     if( m2==scalar_traits<scalar_type>::value(0) )
1274         BOOST_QVM_THROW_EXCEPTION(zero_magnitude_error());
1275     scalar_type const s = sin(angle);
1276     scalar_type const c = cos(angle);
1277     scalar_type const x2 = x*x;
1278     scalar_type const y2 = y*y;
1279     scalar_type const z2 = z*z;
1280     scalar_type const xy = x*y;
1281     scalar_type const xz = x*z;
1282     scalar_type const yz = y*z;
1283     scalar_type const xs = x*s;
1284     scalar_type const ys = y*s;
1285     scalar_type const zs = z*s;
1286     scalar_type const one = scalar_traits<scalar_type>::value(1);
1287     scalar_type const c1 = one-c;
1288     return qvm_detail::rot_mat_<Dim,Angle>(
1289         x2+(one-x2)*c, xy*c1-zs, xz*(one-c)+ys,
1290         xy*c1+zs, y2+(one-y2)*c, yz*c1-xs,
1291         xz*c1-ys, yz*c1+xs, z2+(one-z2)*c );
1292     }
1293
1294 template <class A,class B,class Angle>
1295 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1296 typename enable_if_c<
1297     is_mat<A>::value &&
1298     mat_traits<A>::rows==mat_traits<A>::cols &&
1299     mat_traits<A>::rows>=3 &&
1300     is_vec<B>::value && vec_traits<B>::dim==3,
1301     void>::type
1302 set_rot( A & a, B const & axis, Angle angle )
1303     {
1304     assign(a,rot_mat<mat_traits<A>::rows>(axis,angle));
1305     }
1306
1307 template <class A,class B,class Angle>
1308 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1309 typename enable_if_c<
1310     is_mat<A>::value &&
1311     mat_traits<A>::rows==mat_traits<A>::cols &&
1312     mat_traits<A>::rows>=3 &&
1313     is_vec<B>::value && vec_traits<B>::dim==3,
1314     void>::type
1315 rotate( A & a, B const & axis, Angle angle )
1316     {
1317     a *= rot_mat<mat_traits<A>::rows>(axis,angle);
1318     }
1319
1320 ////////////////////////////////////////////////
1321
1322 template <int Dim,class Angle>
1323 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1324 qvm_detail::rot_mat_<Dim,Angle>
1325 rot_mat_xzy( Angle x1, Angle z2, Angle y3 )
1326     {
1327     typedef Angle scalar_type;
1328     scalar_type const c1 = cos(x1);
1329     scalar_type const s1 = sin(x1);
1330     scalar_type const c2 = cos(z2);
1331     scalar_type const s2 = sin(z2);
1332     scalar_type const c3 = cos(y3);
1333     scalar_type const s3 = sin(y3);
1334     return qvm_detail::rot_mat_<Dim,Angle>(
1335         c2*c3, -s2, c2*s3,
1336         s1*s3 + c1*c3*s2, c1*c2, c1*s2*s3 - c3*s1,
1337         c3*s1*s2 - c1*s3, c2*s1, c1*c3 + s1*s2*s3 );
1338     }
1339
1340 template <class A,class Angle>
1341 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1342 typename enable_if_c<
1343     is_mat<A>::value &&
1344     mat_traits<A>::rows==mat_traits<A>::cols &&
1345     mat_traits<A>::rows>=3,
1346     void>::type
1347 set_rot_xzy( A & a, Angle x1, Angle z2, Angle y3 )
1348     {
1349     assign(a,rot_mat_xzy<mat_traits<A>::rows>(x1,z2,y3));
1350     }
1351
1352 template <class A,class Angle>
1353 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1354 typename enable_if_c<
1355     is_mat<A>::value &&
1356     mat_traits<A>::rows==mat_traits<A>::cols &&
1357     mat_traits<A>::rows>=3,
1358     void>::type
1359 rotate_xzy( A & a, Angle x1, Angle z2, Angle y3 )
1360     {
1361     a *= rot_mat_xzy<mat_traits<A>::rows>(x1,z2,y3);
1362     }
1363
1364 ////////////////////////////////////////////////
1365
1366 template <int Dim,class Angle>
1367 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1368 qvm_detail::rot_mat_<Dim,Angle>
1369 rot_mat_xyz( Angle x1, Angle y2, Angle z3 )
1370     {
1371     typedef Angle scalar_type;
1372     scalar_type const c1 = cos(x1);
1373     scalar_type const s1 = sin(x1);
1374     scalar_type const c2 = cos(y2);
1375     scalar_type const s2 = sin(y2);
1376     scalar_type const c3 = cos(z3);
1377     scalar_type const s3 = sin(z3);
1378     return qvm_detail::rot_mat_<Dim,Angle>(
1379         c2*c3, -c2*s3, s2,
1380         c1*s3 + c3*s1*s2, c1*c3 - s1*s2*s3, -c2*s1,
1381         s1*s3 - c1*c3*s2, c3*s1 + c1*s2*s3, c1*c2 );
1382     }
1383
1384 template <class A,class Angle>
1385 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1386 typename enable_if_c<
1387     is_mat<A>::value &&
1388     mat_traits<A>::rows==mat_traits<A>::cols &&
1389     mat_traits<A>::rows>=3,
1390     void>::type
1391 set_rot_xyz( A & a, Angle x1, Angle y2, Angle z3 )
1392     {
1393     assign(a,rot_mat_xyz<mat_traits<A>::rows>(x1,y2,z3));
1394     }
1395
1396 template <class A,class Angle>
1397 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1398 typename enable_if_c<
1399     is_mat<A>::value &&
1400     mat_traits<A>::rows==mat_traits<A>::cols &&
1401     mat_traits<A>::rows>=3,
1402     void>::type
1403 rotate_xyz( A & a, Angle x1, Angle y2, Angle z3 )
1404     {
1405     a *= rot_mat_xyz<mat_traits<A>::rows>(x1,y2,z3);
1406     }
1407
1408 ////////////////////////////////////////////////
1409
1410 template <int Dim,class Angle>
1411 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1412 qvm_detail::rot_mat_<Dim,Angle>
1413 rot_mat_yxz( Angle y1, Angle x2, Angle z3 )
1414     {
1415     typedef Angle scalar_type;
1416     scalar_type const c1 = cos(y1);
1417     scalar_type const s1 = sin(y1);
1418     scalar_type const c2 = cos(x2);
1419     scalar_type const s2 = sin(x2);
1420     scalar_type const c3 = cos(z3);
1421     scalar_type const s3 = sin(z3);
1422     return qvm_detail::rot_mat_<Dim,Angle>(
1423         c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3, c2*s1,
1424         c2*s3, c2*c3, -s2,
1425         c1*s2*s3 - c3*s1, c1*c3*s2 + s1*s3, c1*c2 );
1426     }
1427
1428 template <class A,class Angle>
1429 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1430 typename enable_if_c<
1431     is_mat<A>::value &&
1432     mat_traits<A>::rows==mat_traits<A>::cols &&
1433     mat_traits<A>::rows>=3,
1434     void>::type
1435 set_rot_yxz( A & a, Angle y1, Angle x2, Angle z3 )
1436     {
1437     assign(a,rot_mat_yxz<mat_traits<A>::rows>(y1,x2,z3));
1438     }
1439
1440 template <class A,class Angle>
1441 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1442 typename enable_if_c<
1443     is_mat<A>::value &&
1444     mat_traits<A>::rows==mat_traits<A>::cols &&
1445     mat_traits<A>::rows>=3,
1446     void>::type
1447 rotate_yxz( A & a, Angle y1, Angle x2, Angle z3 )
1448     {
1449     a *= rot_mat_yxz<mat_traits<A>::rows>(y1,x2,z3);
1450     }
1451
1452 ////////////////////////////////////////////////
1453
1454 template <int Dim,class Angle>
1455 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1456 qvm_detail::rot_mat_<Dim,Angle>
1457 rot_mat_yzx( Angle y1, Angle z2, Angle x3 )
1458     {
1459     typedef Angle scalar_type;
1460     scalar_type const c1 = cos(y1);
1461     scalar_type const s1 = sin(y1);
1462     scalar_type const c2 = cos(z2);
1463     scalar_type const s2 = sin(z2);
1464     scalar_type const c3 = cos(x3);
1465     scalar_type const s3 = sin(x3);
1466     return qvm_detail::rot_mat_<Dim,Angle>(
1467         c1*c2, s1*s3 - c1*c3*s2, c3*s1 + c1*s2*s3,
1468         s2, c2*c3, -c2*s3,
1469         -c2*s1, c1*s3 + c3*s1*s2, c1*c3 - s1*s2*s3 );
1470     }
1471
1472 template <class A,class Angle>
1473 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1474 typename enable_if_c<
1475     is_mat<A>::value &&
1476     mat_traits<A>::rows==mat_traits<A>::cols &&
1477     mat_traits<A>::rows>=3,
1478     void>::type
1479 set_rot_yzx( A & a, Angle y1, Angle z2, Angle x3 )
1480     {
1481     assign(a,rot_mat_yzx<mat_traits<A>::rows>(y1,z2,x3));
1482     }
1483
1484 template <class A,class Angle>
1485 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1486 typename enable_if_c<
1487     is_mat<A>::value &&
1488     mat_traits<A>::rows==mat_traits<A>::cols &&
1489     mat_traits<A>::rows>=3,
1490     void>::type
1491 rotate_yzx( A & a, Angle y1, Angle z2, Angle x3 )
1492     {
1493     a *= rot_mat_yzx<mat_traits<A>::rows>(y1,z2,x3);
1494     }
1495
1496 ////////////////////////////////////////////////
1497
1498 template <int Dim,class Angle>
1499 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1500 qvm_detail::rot_mat_<Dim,Angle>
1501 rot_mat_zyx( Angle z1, Angle y2, Angle x3 )
1502     {
1503     typedef Angle scalar_type;
1504     scalar_type const c1 = cos(z1);
1505     scalar_type const s1 = sin(z1);
1506     scalar_type const c2 = cos(y2);
1507     scalar_type const s2 = sin(y2);
1508     scalar_type const c3 = cos(x3);
1509     scalar_type const s3 = sin(x3);
1510     return qvm_detail::rot_mat_<Dim,Angle>(
1511         c1*c2, c1*s2*s3 - c3*s1, s1*s3 + c1*c3*s2,
1512         c2*s1, c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3,
1513         -s2, c2*s3, c2*c3 );
1514     }
1515
1516 template <class A,class Angle>
1517 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1518 typename enable_if_c<
1519     is_mat<A>::value &&
1520     mat_traits<A>::rows==mat_traits<A>::cols &&
1521     mat_traits<A>::rows>=3,
1522     void>::type
1523 set_rot_zyx( A & a, Angle z1, Angle y2, Angle x3 )
1524     {
1525     assign(a,rot_mat_zyx<mat_traits<A>::rows>(z1,y2,x3));
1526     }
1527
1528 template <class A,class Angle>
1529 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1530 typename enable_if_c<
1531     is_mat<A>::value &&
1532     mat_traits<A>::rows==mat_traits<A>::cols &&
1533     mat_traits<A>::rows>=3,
1534     void>::type
1535 rotate_zyx( A & a, Angle z1, Angle y2, Angle x3 )
1536     {
1537     a *= rot_mat_zyx<mat_traits<A>::rows>(z1,y2,x3);
1538     }
1539
1540 ////////////////////////////////////////////////
1541
1542 template <int Dim,class Angle>
1543 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1544 qvm_detail::rot_mat_<Dim,Angle>
1545 rot_mat_zxy( Angle z1, Angle x2, Angle y3 )
1546     {
1547     typedef Angle scalar_type;
1548     scalar_type const c1 = cos(z1);
1549     scalar_type const s1 = sin(z1);
1550     scalar_type const c2 = cos(x2);
1551     scalar_type const s2 = sin(x2);
1552     scalar_type const c3 = cos(y3);
1553     scalar_type const s3 = sin(y3);
1554     return qvm_detail::rot_mat_<Dim,Angle>(
1555         c1*c3 - s1*s2*s3, -c2*s1, c1*s3 + c3*s1*s2,
1556         c3*s1 + c1*s2*s3, c1*c2, s1*s3 - c1*c3*s2,
1557         -c2*s3, s2, c2*c3 );
1558     }
1559
1560 template <class A,class Angle>
1561 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1562 typename enable_if_c<
1563     is_mat<A>::value &&
1564     mat_traits<A>::rows==mat_traits<A>::cols &&
1565     mat_traits<A>::rows>=3,
1566     void>::type
1567 set_rot_zxy( A & a, Angle z1, Angle x2, Angle y3 )
1568     {
1569     assign(a,rot_mat_zxy<mat_traits<A>::rows>(z1,x2,y3));
1570     }
1571
1572 template <class A,class Angle>
1573 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1574 typename enable_if_c<
1575     is_mat<A>::value &&
1576     mat_traits<A>::rows==mat_traits<A>::cols &&
1577     mat_traits<A>::rows>=3,
1578     void>::type
1579 rotate_zxy( A & a, Angle z1, Angle x2, Angle y3 )
1580     {
1581     a *= rot_mat_zxy<mat_traits<A>::rows>(z1,x2,y3);
1582     }
1583
1584 ////////////////////////////////////////////////
1585
1586 template <int Dim,class Angle>
1587 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1588 qvm_detail::rot_mat_<Dim,Angle>
1589 rot_mat_xzx( Angle x1, Angle z2, Angle x3 )
1590     {
1591     typedef Angle scalar_type;
1592     scalar_type const c1 = cos(x1);
1593     scalar_type const s1 = sin(x1);
1594     scalar_type const c2 = cos(z2);
1595     scalar_type const s2 = sin(z2);
1596     scalar_type const c3 = cos(x3);
1597     scalar_type const s3 = sin(x3);
1598     return qvm_detail::rot_mat_<Dim,Angle>(
1599         c2, -c3*s2, s2*s3,
1600         c1*s2, c1*c2*c3 - s1*s3, -c3*s1 - c1*c2*s3,
1601         s1*s2, c1*s3 + c2*c3*s1, c1*c3 - c2*s1*s3 );
1602     }
1603
1604 template <class A,class Angle>
1605 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1606 typename enable_if_c<
1607     is_mat<A>::value &&
1608     mat_traits<A>::rows==mat_traits<A>::cols &&
1609     mat_traits<A>::rows>=3,
1610     void>::type
1611 set_rot_xzx( A & a, Angle x1, Angle z2, Angle x3 )
1612     {
1613     assign(a,rot_mat_xzx<mat_traits<A>::rows>(x1,z2,x3));
1614     }
1615
1616 template <class A,class Angle>
1617 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1618 typename enable_if_c<
1619     is_mat<A>::value &&
1620     mat_traits<A>::rows==mat_traits<A>::cols &&
1621     mat_traits<A>::rows>=3,
1622     void>::type
1623 rotate_xzx( A & a, Angle x1, Angle z2, Angle x3 )
1624     {
1625     a *= rot_mat_xzx<mat_traits<A>::rows>(x1,z2,x3);
1626     }
1627
1628 ////////////////////////////////////////////////
1629
1630 template <int Dim,class Angle>
1631 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1632 qvm_detail::rot_mat_<Dim,Angle>
1633 rot_mat_xyx( Angle x1, Angle y2, Angle x3 )
1634     {
1635     typedef Angle scalar_type;
1636     scalar_type const c1 = cos(x1);
1637     scalar_type const s1 = sin(x1);
1638     scalar_type const c2 = cos(y2);
1639     scalar_type const s2 = sin(y2);
1640     scalar_type const c3 = cos(x3);
1641     scalar_type const s3 = sin(x3);
1642     return qvm_detail::rot_mat_<Dim,Angle>(
1643         c2, s2*s3, c3*s2,
1644         s1*s2, c1*c3 - c2*s1*s3, -c1*s3 - c2*c3*s1,
1645         -c1*s2, c3*s1 + c1*c2*s3, c1*c2*c3 - s1*s3 );
1646     }
1647
1648 template <class A,class Angle>
1649 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1650 typename enable_if_c<
1651     is_mat<A>::value &&
1652     mat_traits<A>::rows==mat_traits<A>::cols &&
1653     mat_traits<A>::rows>=3,
1654     void>::type
1655 set_rot_xyx( A & a, Angle x1, Angle y2, Angle x3 )
1656     {
1657     assign(a,rot_mat_xyx<mat_traits<A>::rows>(x1,y2,x3));
1658     }
1659
1660 template <class A,class Angle>
1661 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1662 typename enable_if_c<
1663     is_mat<A>::value &&
1664     mat_traits<A>::rows==mat_traits<A>::cols &&
1665     mat_traits<A>::rows>=3,
1666     void>::type
1667 rotate_xyx( A & a, Angle x1, Angle y2, Angle x3 )
1668     {
1669     a *= rot_mat_xyx<mat_traits<A>::rows>(x1,y2,x3);
1670     }
1671
1672 ////////////////////////////////////////////////
1673
1674 template <int Dim,class Angle>
1675 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1676 qvm_detail::rot_mat_<Dim,Angle>
1677 rot_mat_yxy( Angle y1, Angle x2, Angle y3 )
1678     {
1679     typedef Angle scalar_type;
1680     scalar_type const c1 = cos(y1);
1681     scalar_type const s1 = sin(y1);
1682     scalar_type const c2 = cos(x2);
1683     scalar_type const s2 = sin(x2);
1684     scalar_type const c3 = cos(y3);
1685     scalar_type const s3 = sin(y3);
1686     return qvm_detail::rot_mat_<Dim,Angle>(
1687         c1*c3 - c2*s1*s3, s1*s2, c1*s3 + c2*c3*s1,
1688         s2*s3, c2, -c3*s2,
1689         -c3*s1 - c1*c2*s3, c1*s2, c1*c2*c3 - s1*s3 );
1690     }
1691
1692 template <class A,class Angle>
1693 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1694 typename enable_if_c<
1695     is_mat<A>::value &&
1696     mat_traits<A>::rows==mat_traits<A>::cols &&
1697     mat_traits<A>::rows>=3,
1698     void>::type
1699 set_rot_yxy( A & a, Angle y1, Angle x2, Angle y3 )
1700     {
1701     assign(a,rot_mat_yxy<mat_traits<A>::rows>(y1,x2,y3));
1702     }
1703
1704 template <class A,class Angle>
1705 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1706 typename enable_if_c<
1707     is_mat<A>::value &&
1708     mat_traits<A>::rows==mat_traits<A>::cols &&
1709     mat_traits<A>::rows>=3,
1710     void>::type
1711 rotate_yxy( A & a, Angle y1, Angle x2, Angle y3 )
1712     {
1713     a *= rot_mat_yxy<mat_traits<A>::rows>(y1,x2,y3);
1714     }
1715
1716 ////////////////////////////////////////////////
1717
1718 template <int Dim,class Angle>
1719 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1720 qvm_detail::rot_mat_<Dim,Angle>
1721 rot_mat_yzy( Angle y1, Angle z2, Angle y3 )
1722     {
1723     typedef Angle scalar_type;
1724     scalar_type const c1 = cos(y1);
1725     scalar_type const s1 = sin(y1);
1726     scalar_type const c2 = cos(z2);
1727     scalar_type const s2 = sin(z2);
1728     scalar_type const c3 = cos(y3);
1729     scalar_type const s3 = sin(y3);
1730     return qvm_detail::rot_mat_<Dim,Angle>(
1731         c1*c2*c3 - s1*s3, -c1*s2, c3*s1 + c1*c2*s3,
1732         c3*s2, c2, s2*s3,
1733         -c1*s3 - c2*c3*s1, s1*s2, c1*c3 - c2*s1*s3 );
1734     }
1735
1736 template <class A,class Angle>
1737 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1738 typename enable_if_c<
1739     is_mat<A>::value &&
1740     mat_traits<A>::rows==mat_traits<A>::cols &&
1741     mat_traits<A>::rows>=3,
1742     void>::type
1743 set_rot_yzy( A & a, Angle y1, Angle z2, Angle y3 )
1744     {
1745     assign(a,rot_mat_yzy<mat_traits<A>::rows>(y1,z2,y3));
1746     }
1747
1748 template <class A,class Angle>
1749 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1750 typename enable_if_c<
1751     is_mat<A>::value &&
1752     mat_traits<A>::rows==mat_traits<A>::cols &&
1753     mat_traits<A>::rows>=3,
1754     void>::type
1755 rotate_yzy( A & a, Angle y1, Angle z2, Angle y3 )
1756     {
1757     a *= rot_mat_yzy<mat_traits<A>::rows>(y1,z2,y3);
1758     }
1759
1760 ////////////////////////////////////////////////
1761
1762 template <int Dim,class Angle>
1763 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1764 qvm_detail::rot_mat_<Dim,Angle>
1765 rot_mat_zyz( Angle z1, Angle y2, Angle z3 )
1766     {
1767     typedef Angle scalar_type;
1768     scalar_type const c1 = cos(z1);
1769     scalar_type const s1 = sin(z1);
1770     scalar_type const c2 = cos(y2);
1771     scalar_type const s2 = sin(y2);
1772     scalar_type const c3 = cos(z3);
1773     scalar_type const s3 = sin(z3);
1774     return qvm_detail::rot_mat_<Dim,Angle>(
1775         c1*c2*c3 - s1*s3, -c3*s1 - c1*c2*s3, c1*s2,
1776         c1*s3 + c2*c3*s1, c1*c3 - c2*s1*s3, s1*s2,
1777         -c3*s2, s2*s3, c2 );
1778     }
1779
1780 template <class A,class Angle>
1781 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1782 typename enable_if_c<
1783     is_mat<A>::value &&
1784     mat_traits<A>::rows==mat_traits<A>::cols &&
1785     mat_traits<A>::rows>=3,
1786     void>::type
1787 set_rot_zyz( A & a, Angle z1, Angle y2, Angle z3 )
1788     {
1789     assign(a,rot_mat_zyz<mat_traits<A>::rows>(z1,y2,z3));
1790     }
1791
1792 template <class A,class Angle>
1793 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1794 typename enable_if_c<
1795     is_mat<A>::value &&
1796     mat_traits<A>::rows==mat_traits<A>::cols &&
1797     mat_traits<A>::rows>=3,
1798     void>::type
1799 rotate_zyz( A & a, Angle z1, Angle y2, Angle z3 )
1800     {
1801     a *= rot_mat_zyz<mat_traits<A>::rows>(z1,y2,z3);
1802     }
1803
1804 ////////////////////////////////////////////////
1805
1806 template <int Dim,class Angle>
1807 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1808 qvm_detail::rot_mat_<Dim,Angle>
1809 rot_mat_zxz( Angle z1, Angle x2, Angle z3 )
1810     {
1811     typedef Angle scalar_type;
1812     scalar_type const c1 = cos(z1);
1813     scalar_type const s1 = sin(z1);
1814     scalar_type const c2 = cos(x2);
1815     scalar_type const s2 = sin(x2);
1816     scalar_type const c3 = cos(z3);
1817     scalar_type const s3 = sin(z3);
1818     return qvm_detail::rot_mat_<Dim,Angle>(
1819         c1*c3 - c2*s1*s3, -c1*s3 - c2*c3*s1, s1*s2,
1820         c3*s1 + c1*c2*s3, c1*c2*c3 - s1*s3, -c1*s2,
1821         s2*s3, c3*s2, c2 );
1822     }
1823
1824 template <class A,class Angle>
1825 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1826 typename enable_if_c<
1827     is_mat<A>::value &&
1828     mat_traits<A>::rows==mat_traits<A>::cols &&
1829     mat_traits<A>::rows>=3,
1830     void>::type
1831 set_rot_zxz( A & a, Angle z1, Angle x2, Angle z3 )
1832     {
1833     assign(a,rot_mat_zxz<mat_traits<A>::rows>(z1,x2,z3));
1834     }
1835
1836 template <class A,class Angle>
1837 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1838 typename enable_if_c<
1839     is_mat<A>::value &&
1840     mat_traits<A>::rows==mat_traits<A>::cols &&
1841     mat_traits<A>::rows>=3,
1842     void>::type
1843 rotate_zxz( A & a, Angle z1, Angle x2, Angle z3 )
1844     {
1845     a *= rot_mat_zxz<mat_traits<A>::rows>(z1,x2,z3);
1846     }
1847
1848 ////////////////////////////////////////////////
1849
1850 namespace
1851 qvm_detail
1852     {
1853     template <int Dim,class Angle>
1854     struct
1855     rotx_mat_
1856         {
1857         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1858         rotx_mat_()
1859             {
1860             }
1861
1862         template <class R
1863 #if __cplusplus >= 201103L
1864             , class = typename enable_if<is_mat<R> >::type
1865 #endif
1866         >
1867         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1868         operator R() const
1869             {
1870             R r;
1871             assign(r,*this);
1872             return r;
1873             }
1874
1875         private:
1876
1877         rotx_mat_( rotx_mat_ const & );
1878         rotx_mat_ & operator=( rotx_mat_ const & );
1879         ~rotx_mat_();
1880         };
1881
1882     template <int Row,int Col>
1883     struct
1884     rotx_m_get
1885         {
1886         template <class T>
1887         static
1888         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1889         T
1890         get( T const & )
1891             {
1892             return scalar_traits<T>::value(Row==Col);
1893             }
1894         };
1895
1896     template <>
1897     struct
1898     rotx_m_get<1,1>
1899         {
1900         template <class T>
1901         static
1902         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1903         T
1904         get( T const & angle )
1905             {
1906             return cos(angle);
1907             }
1908         };
1909
1910     template <>
1911     struct
1912     rotx_m_get<1,2>
1913         {
1914         template <class T>
1915         static
1916         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1917         T
1918         get( T const & angle )
1919             {
1920             return -sin(angle);
1921             }
1922         };
1923
1924     template <>
1925     struct
1926     rotx_m_get<2,1>
1927         {
1928         template <class T>
1929         static
1930         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1931         T
1932         get( T const & angle )
1933             {
1934             return sin(angle);
1935             }
1936         };
1937
1938     template <>
1939     struct
1940     rotx_m_get<2,2>
1941         {
1942         template <class T>
1943         static
1944         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1945         T
1946         get( T const & angle )
1947             {
1948             return cos(angle);
1949             }
1950         };
1951     }
1952
1953 template <int Dim,class Angle>
1954 struct
1955 mat_traits< qvm_detail::rotx_mat_<Dim,Angle> >
1956     {
1957     typedef qvm_detail::rotx_mat_<Dim,Angle> this_matrix;
1958     typedef Angle scalar_type;
1959     static int const rows=Dim;
1960     static int const cols=Dim;
1961
1962     template <int Row,int Col>
1963     static
1964     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1965     scalar_type
1966     read_element( this_matrix const & x )
1967         {
1968         BOOST_QVM_STATIC_ASSERT(Row>=0);
1969         BOOST_QVM_STATIC_ASSERT(Col>=0);
1970         BOOST_QVM_STATIC_ASSERT(Row<Dim);
1971         BOOST_QVM_STATIC_ASSERT(Col<Dim);
1972         return qvm_detail::rotx_m_get<Row,Col>::get(reinterpret_cast<Angle const &>(x));
1973         }
1974
1975     static
1976     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1977     scalar_type
1978     read_element_idx( int row, int col, this_matrix const & x )
1979         {
1980         BOOST_QVM_ASSERT(row>=0);
1981         BOOST_QVM_ASSERT(col>=0);
1982         BOOST_QVM_ASSERT(row<Dim);
1983         BOOST_QVM_ASSERT(col<Dim);
1984         Angle const & a=reinterpret_cast<Angle const &>(x);
1985         if( row==1 )
1986             {
1987             if( col==1 )
1988                 return cos(a);
1989             if( col==2 )
1990                 return -sin(a);
1991             }
1992         if( row==2 )
1993             {
1994             if( col==1 )
1995                 return sin(a);
1996             if( col==2 )
1997                 return cos(a);
1998             }
1999         return scalar_traits<scalar_type>::value(row==col);
2000         }
2001     };
2002
2003 template <int Dim,class Angle>
2004 struct
2005 deduce_mat<qvm_detail::rotx_mat_<Dim,Angle>,Dim,Dim>
2006     {
2007     typedef mat<Angle,Dim,Dim> type;
2008     };
2009
2010 template <int Dim,class Angle>
2011 struct
2012 deduce_mat2<qvm_detail::rotx_mat_<Dim,Angle>,qvm_detail::rotx_mat_<Dim,Angle>,Dim,Dim>
2013     {
2014     typedef mat<Angle,Dim,Dim> type;
2015     };
2016
2017 template <int Dim,class Angle>
2018 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2019 qvm_detail::rotx_mat_<Dim,Angle> const &
2020 rotx_mat( Angle const & angle )
2021     {
2022     BOOST_QVM_STATIC_ASSERT(Dim>=3);
2023     return reinterpret_cast<qvm_detail::rotx_mat_<Dim,Angle> const &>(angle);
2024     }
2025
2026 template <class A,class Angle>
2027 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2028 typename enable_if_c<
2029     is_mat<A>::value &&
2030     mat_traits<A>::rows>=3 &&
2031     mat_traits<A>::rows==mat_traits<A>::cols,
2032     void>::type
2033 set_rotx( A & a, Angle angle )
2034     {
2035     assign(a,rotx_mat<mat_traits<A>::rows>(angle));
2036     }
2037
2038 template <class A,class Angle>
2039 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2040 typename enable_if_c<
2041     is_mat<A>::value &&
2042     mat_traits<A>::rows>=3 &&
2043     mat_traits<A>::rows==mat_traits<A>::cols,
2044     void>::type
2045 rotate_x( A & a, Angle angle )
2046     {
2047     a *= rotx_mat<mat_traits<A>::rows>(angle);
2048     }
2049
2050 ////////////////////////////////////////////////
2051
2052 namespace
2053 qvm_detail
2054     {
2055     template <int Dim,class Angle>
2056     struct
2057     roty_mat_
2058         {
2059         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2060         roty_mat_()
2061             {
2062             }
2063
2064         template <class R
2065 #if __cplusplus >= 201103L
2066             , class = typename enable_if<is_mat<R> >::type
2067 #endif
2068         >
2069         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2070         operator R() const
2071             {
2072             R r;
2073             assign(r,*this);
2074             return r;
2075             }
2076
2077         private:
2078
2079         roty_mat_( roty_mat_ const & );
2080         roty_mat_ & operator=( roty_mat_ const & );
2081         ~roty_mat_();
2082         };
2083
2084     template <int Row,int Col>
2085     struct
2086     roty_m_get
2087         {
2088         template <class T>
2089         static
2090         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2091         T
2092         get( T const & )
2093             {
2094             return scalar_traits<T>::value(Row==Col);
2095             }
2096         };
2097
2098     template <>
2099     struct
2100     roty_m_get<0,0>
2101         {
2102         template <class T>
2103         static
2104         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2105         T
2106         get( T const & angle )
2107             {
2108             return cos(angle);
2109             }
2110         };
2111
2112     template <>
2113     struct
2114     roty_m_get<0,2>
2115         {
2116         template <class T>
2117         static
2118         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2119         T
2120         get( T const & angle )
2121             {
2122             return sin(angle);
2123             }
2124         };
2125
2126     template <>
2127     struct
2128     roty_m_get<2,0>
2129         {
2130         template <class T>
2131         static
2132         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2133         T
2134         get( T const & angle )
2135             {
2136             return -sin(angle);
2137             }
2138         };
2139
2140     template <>
2141     struct
2142     roty_m_get<2,2>
2143         {
2144         template <class T>
2145         static
2146         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2147         T
2148         get( T const & angle )
2149             {
2150             return cos(angle);
2151             }
2152         };
2153     }
2154
2155 template <int Dim,class Angle>
2156 struct
2157 mat_traits< qvm_detail::roty_mat_<Dim,Angle> >
2158     {
2159     typedef qvm_detail::roty_mat_<Dim,Angle> this_matrix;
2160     typedef Angle scalar_type;
2161     static int const rows=Dim;
2162     static int const cols=Dim;
2163
2164     template <int Row,int Col>
2165     static
2166     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2167     scalar_type
2168     read_element( this_matrix const & x )
2169         {
2170         BOOST_QVM_STATIC_ASSERT(Row>=0);
2171         BOOST_QVM_STATIC_ASSERT(Col>=0);
2172         BOOST_QVM_STATIC_ASSERT(Row<Dim);
2173         BOOST_QVM_STATIC_ASSERT(Col<Dim);
2174         return qvm_detail::roty_m_get<Row,Col>::get(reinterpret_cast<Angle const &>(x));
2175         }
2176
2177     static
2178     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2179     scalar_type
2180     read_element_idx( int row, int col, this_matrix const & x )
2181         {
2182         BOOST_QVM_ASSERT(row>=0);
2183         BOOST_QVM_ASSERT(col>=0);
2184         BOOST_QVM_ASSERT(row<Dim);
2185         BOOST_QVM_ASSERT(col<Dim);
2186         Angle const & a=reinterpret_cast<Angle const &>(x);
2187         if( row==0 )
2188             {
2189             if( col==0 )
2190                 return cos(a);
2191             if( col==2 )
2192                 return sin(a);
2193             }
2194         if( row==2 )
2195             {
2196             if( col==0 )
2197                 return -sin(a);
2198             if( col==2 )
2199                 return cos(a);
2200             }
2201         return scalar_traits<scalar_type>::value(row==col);
2202         }
2203     };
2204
2205 template <int Dim,class Angle>
2206 struct
2207 deduce_mat<qvm_detail::roty_mat_<Dim,Angle>,Dim,Dim>
2208     {
2209     typedef mat<Angle,Dim,Dim> type;
2210     };
2211
2212 template <int Dim,class Angle>
2213 struct
2214 deduce_mat2<qvm_detail::roty_mat_<Dim,Angle>,qvm_detail::roty_mat_<Dim,Angle>,Dim,Dim>
2215     {
2216     typedef mat<Angle,Dim,Dim> type;
2217     };
2218
2219 template <int Dim,class Angle>
2220 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2221 qvm_detail::roty_mat_<Dim,Angle> const &
2222 roty_mat( Angle const & angle )
2223     {
2224     BOOST_QVM_STATIC_ASSERT(Dim>=3);
2225     return reinterpret_cast<qvm_detail::roty_mat_<Dim,Angle> const &>(angle);
2226     }
2227
2228 template <class A,class Angle>
2229 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2230 typename enable_if_c<
2231     is_mat<A>::value &&
2232     mat_traits<A>::rows>=2 &&
2233     mat_traits<A>::rows==mat_traits<A>::cols,
2234     void>::type
2235 set_roty( A & a, Angle angle )
2236     {
2237     assign(a,roty_mat<mat_traits<A>::rows>(angle));
2238     }
2239
2240 template <class A,class Angle>
2241 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2242 typename enable_if_c<
2243     is_mat<A>::value &&
2244     mat_traits<A>::rows>=3 &&
2245     mat_traits<A>::rows==mat_traits<A>::cols,
2246     void>::type
2247 rotate_y( A & a, Angle angle )
2248     {
2249     a *= roty_mat<mat_traits<A>::rows>(angle);
2250     }
2251
2252 ////////////////////////////////////////////////
2253
2254 namespace
2255 qvm_detail
2256     {
2257     template <int Dim,class Angle>
2258     struct
2259     rotz_mat_
2260         {
2261         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2262         rotz_mat_()
2263             {
2264             }
2265
2266         template <class R
2267 #if __cplusplus >= 201103L
2268             , class = typename enable_if<is_mat<R> >::type
2269 #endif
2270         >
2271         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2272         operator R() const
2273             {
2274             R r;
2275             assign(r,*this);
2276             return r;
2277             }
2278
2279         private:
2280
2281         rotz_mat_( rotz_mat_ const & );
2282         rotz_mat_ & operator=( rotz_mat_ const & );
2283         ~rotz_mat_();
2284         };
2285
2286     template <int Row,int Col>
2287     struct
2288     rotz_m_get
2289         {
2290         template <class T>
2291         static
2292         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2293         T
2294         get( T const & )
2295             {
2296             return scalar_traits<T>::value(Row==Col);
2297             }
2298         };
2299
2300     template <>
2301     struct
2302     rotz_m_get<0,0>
2303         {
2304         template <class T>
2305         static
2306         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2307         T
2308         get( T const & angle )
2309             {
2310             return cos(angle);
2311             }
2312         };
2313
2314     template <>
2315     struct
2316     rotz_m_get<0,1>
2317         {
2318         template <class T>
2319         static
2320         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2321         T
2322         get( T const & angle )
2323             {
2324             return -sin(angle);
2325             }
2326         };
2327
2328     template <>
2329     struct
2330     rotz_m_get<1,0>
2331         {
2332         template <class T>
2333         static
2334         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2335         T
2336         get( T const & angle )
2337             {
2338             return sin(angle);
2339             }
2340         };
2341
2342     template <>
2343     struct
2344     rotz_m_get<1,1>
2345         {
2346         template <class T>
2347         static
2348         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2349         T
2350         get( T const & angle )
2351             {
2352             return cos(angle);
2353             }
2354         };
2355     }
2356
2357 template <int Dim,class Angle>
2358 struct
2359 mat_traits< qvm_detail::rotz_mat_<Dim,Angle> >
2360     {
2361     typedef qvm_detail::rotz_mat_<Dim,Angle> this_matrix;
2362     typedef Angle scalar_type;
2363     static int const rows=Dim;
2364     static int const cols=Dim;
2365
2366     template <int Row,int Col>
2367     static
2368     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2369     scalar_type
2370     read_element( this_matrix const & x )
2371         {
2372         BOOST_QVM_STATIC_ASSERT(Row>=0);
2373         BOOST_QVM_STATIC_ASSERT(Col>=0);
2374         BOOST_QVM_STATIC_ASSERT(Row<Dim);
2375         BOOST_QVM_STATIC_ASSERT(Col<Dim);
2376         return qvm_detail::rotz_m_get<Row,Col>::get(reinterpret_cast<Angle const &>(x));
2377         }
2378
2379     static
2380     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2381     scalar_type
2382     read_element_idx( int row, int col, this_matrix const & x )
2383         {
2384         BOOST_QVM_ASSERT(row>=0);
2385         BOOST_QVM_ASSERT(col>=0);
2386         BOOST_QVM_ASSERT(row<Dim);
2387         BOOST_QVM_ASSERT(col<Dim);
2388         Angle const & a=reinterpret_cast<Angle const &>(x);
2389         if( row==0 )
2390             {
2391             if( col==0 )
2392                 return cos(a);
2393             if( col==1 )
2394                 return -sin(a);
2395             }
2396         if( row==1 )
2397             {
2398             if( col==0 )
2399                 return sin(a);
2400             if( col==1 )
2401                 return cos(a);
2402             }
2403         return scalar_traits<scalar_type>::value(row==col);
2404         }
2405     };
2406
2407 template <int Dim,class Angle>
2408 struct
2409 deduce_mat<qvm_detail::rotz_mat_<Dim,Angle>,Dim,Dim>
2410     {
2411     typedef mat<Angle,Dim,Dim> type;
2412     };
2413
2414 template <int Dim,class Angle,int R,int C>
2415 struct
2416 deduce_mat2<qvm_detail::rotz_mat_<Dim,Angle>,qvm_detail::rotz_mat_<Dim,Angle>,R,C>
2417     {
2418     typedef mat<Angle,R,C> type;
2419     };
2420
2421 template <int Dim,class Angle>
2422 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2423 qvm_detail::rotz_mat_<Dim,Angle> const &
2424 rotz_mat( Angle const & angle )
2425     {
2426     BOOST_QVM_STATIC_ASSERT(Dim>=2);
2427     return reinterpret_cast<qvm_detail::rotz_mat_<Dim,Angle> const &>(angle);
2428     }
2429
2430 template <class A,class Angle>
2431 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2432 typename enable_if_c<
2433     is_mat<A>::value &&
2434     mat_traits<A>::rows>=2 &&
2435     mat_traits<A>::rows==mat_traits<A>::cols,
2436     void>::type
2437 set_rotz( A & a, Angle angle )
2438     {
2439     assign(a,rotz_mat<mat_traits<A>::rows>(angle));
2440     }
2441
2442 template <class A,class Angle>
2443 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2444 typename enable_if_c<
2445     is_mat<A>::value &&
2446     mat_traits<A>::rows>=2 &&
2447     mat_traits<A>::rows==mat_traits<A>::cols,
2448     void>::type
2449 rotate_z( A & a, Angle angle )
2450     {
2451     a *= rotz_mat<mat_traits<A>::rows>(angle);
2452     }
2453
2454 ////////////////////////////////////////////////
2455
2456 namespace
2457 qvm_detail
2458     {
2459     template <int D>
2460     struct
2461     inverse_m_defined
2462         {
2463         static bool const value=false;
2464         };
2465     }
2466
2467 template <class A,class B>
2468 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2469 typename lazy_enable_if_c<
2470     is_mat<A>::value && is_scalar<B>::value &&
2471     mat_traits<A>::rows==mat_traits<A>::cols &&
2472     !qvm_detail::inverse_m_defined<mat_traits<A>::rows>::value,
2473     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
2474 inverse( A const & a, B det )
2475     {
2476     typedef typename mat_traits<A>::scalar_type T;
2477     BOOST_QVM_ASSERT(det!=scalar_traits<T>::value(0));
2478     T f=scalar_traits<T>::value(1)/det;
2479     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type cofactor_return_type;
2480     cofactor_return_type c=qvm_detail::cofactor_impl(a);
2481     return reinterpret_cast<qvm_detail::transposed_<cofactor_return_type> const &>(c) * f;
2482     }
2483
2484 template <class A>
2485 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2486 typename lazy_enable_if_c<
2487     is_mat<A>::value &&
2488     mat_traits<A>::rows==mat_traits<A>::cols &&
2489     !qvm_detail::inverse_m_defined<mat_traits<A>::rows>::value,
2490     deduce_mat<A> >::type
2491 inverse( A const & a )
2492     {
2493     typedef typename mat_traits<A>::scalar_type T;
2494     T det=determinant(a);
2495     if( det==scalar_traits<T>::value(0) )
2496         BOOST_QVM_THROW_EXCEPTION(zero_determinant_error());
2497     return inverse(a,det);
2498     }
2499
2500 ////////////////////////////////////////////////
2501
2502 namespace
2503 sfinae
2504     {
2505     using ::boost::qvm::to_string;
2506     using ::boost::qvm::assign;
2507     using ::boost::qvm::determinant;
2508     using ::boost::qvm::cmp;
2509     using ::boost::qvm::convert_to;
2510     using ::boost::qvm::set_identity;
2511     using ::boost::qvm::set_zero;
2512     using ::boost::qvm::scalar_cast;
2513     using ::boost::qvm::operator/=;
2514     using ::boost::qvm::operator/;
2515     using ::boost::qvm::operator==;
2516     using ::boost::qvm::operator-=;
2517     using ::boost::qvm::operator-;
2518     using ::boost::qvm::operator*=;
2519     using ::boost::qvm::operator*;
2520     using ::boost::qvm::operator!=;
2521     using ::boost::qvm::operator+=;
2522     using ::boost::qvm::operator+;
2523     using ::boost::qvm::mref;
2524     using ::boost::qvm::rot_mat;
2525     using ::boost::qvm::set_rot;
2526     using ::boost::qvm::rotate;
2527     using ::boost::qvm::set_rotx;
2528     using ::boost::qvm::rotate_x;
2529     using ::boost::qvm::set_roty;
2530     using ::boost::qvm::rotate_y;
2531     using ::boost::qvm::set_rotz;
2532     using ::boost::qvm::rotate_z;
2533     using ::boost::qvm::inverse;
2534     }
2535
2536 } }
2537
2538 #endif