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