ceph/src/boost/boost/qvm/mat_operations.hpp

   1 #ifndef BOOST_QVM_MAT_OPERATIONS_HPP_INCLUDED
   2 #define BOOST_QVM_MAT_OPERATIONS_HPP_INCLUDED
   3
   4 // Copyright 2008-2022 Emil Dotchevski and Reverge Studios, Inc.
   5 // Copyright 2019 agate-pris
   6
   7 // Distributed under the Boost Software License, Version 1.0. (See accompanying
   8 // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   9
  10 #include <boost/qvm/detail/mat_assign.hpp>
  11 #include <boost/qvm/mat_operations2.hpp>
  12 #include <boost/qvm/mat_operations3.hpp>
  13 #include <boost/qvm/mat_operations4.hpp>
  14 #include <boost/qvm/math.hpp>
  15 #include <boost/qvm/detail/determinant_impl.hpp>
  16 #include <boost/qvm/detail/cofactor_impl.hpp>
  17 #include <boost/qvm/detail/transp_impl.hpp>
  18 #include <boost/qvm/scalar_traits.hpp>
  19 #include <string>
  20
  21 namespace boost { namespace qvm {
  22
  23 namespace
  24 qvm_detail
  25     {
  26     BOOST_QVM_INLINE_CRITICAL
  27     void const *
  28     get_valid_ptr_mat_operations()
  29         {
  30         static int const obj=0;
  31         return &obj;
  32         }
  33     }
  34
  35 ////////////////////////////////////////////////
  36
  37 namespace
  38 qvm_to_string_detail
  39     {
  40     template <class T>
  41     std::string to_string( T const & x );
  42     }
  43
  44 namespace
  45 qvm_detail
  46     {
  47     template <int R,int C>
  48     struct
  49     to_string_m_defined
  50         {
  51         static bool const value=false;
  52         };
  53
  54     template <int I,int SizeMinusOne>
  55     struct
  56     to_string_matrix_elements
  57         {
  58         template <class A>
  59         static
  60         std::string
  61         f( A const & a )
  62             {
  63             using namespace qvm_to_string_detail;
  64             return
  65                 ( (I%mat_traits<A>::cols)==0 ? '(' : ',' ) +
  66                 to_string(mat_traits<A>::template read_element<I/mat_traits<A>::cols,I%mat_traits<A>::cols>(a)) +
  67                 ( (I%mat_traits<A>::cols)==mat_traits<A>::cols-1 ? ")" : "" ) +
  68                 to_string_matrix_elements<I+1,SizeMinusOne>::f(a);
  69             }
  70         };
  71
  72     template <int SizeMinusOne>
  73     struct
  74     to_string_matrix_elements<SizeMinusOne,SizeMinusOne>
  75         {
  76         template <class A>
  77         static
  78         std::string
  79         f( A const & a )
  80             {
  81             using namespace qvm_to_string_detail;
  82             return
  83                 ( (SizeMinusOne%mat_traits<A>::cols)==0 ? '(' : ',' ) +
  84                 to_string(mat_traits<A>::template read_element<SizeMinusOne/mat_traits<A>::cols,SizeMinusOne%mat_traits<A>::cols>(a)) +
  85                 ')';
  86             }
  87         };
  88     }
  89
  90 template <class A>
  91 inline
  92 typename enable_if_c<
  93     is_mat<A>::value  &&
  94     !qvm_detail::to_string_m_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
  95     std::string>::type
  96 to_string( A const & a )
  97     {
  98     return "("+qvm_detail::to_string_matrix_elements<0,mat_traits<A>::rows*mat_traits<A>::cols-1>::f(a)+')';
  99     }
 100
 101 ////////////////////////////////////////////////
 102
 103 template <class A,class B,class Cmp>
 104 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 105 typename enable_if_c<
 106     is_mat<A>::value && is_mat<B>::value &&
 107     mat_traits<A>::rows==mat_traits<B>::rows &&
 108     mat_traits<A>::cols==mat_traits<B>::cols,
 109     bool>::type
 110 cmp( A const & a, B const & b, Cmp f )
 111     {
 112     for( int i=0; i!=mat_traits<A>::rows; ++i )
 113         for( int j=0; j!=mat_traits<A>::cols; ++j )
 114             if( !f(
 115                 mat_traits<A>::read_element_idx(i, j, a),
 116                 mat_traits<B>::read_element_idx(i, j, b)) )
 117                 return false;
 118     return true;
 119     }
 120
 121 ////////////////////////////////////////////////
 122
 123 namespace
 124 qvm_detail
 125     {
 126     template <int M,int N>
 127     struct
 128     convert_to_m_defined
 129         {
 130         static bool const value=false;
 131         };
 132     }
 133
 134 template <class R,class A>
 135 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 136 typename enable_if_c<
 137     is_mat<R>::value && is_mat<A>::value &&
 138     mat_traits<R>::rows==mat_traits<A>::rows &&
 139     mat_traits<R>::cols==mat_traits<A>::cols &&
 140     !qvm_detail::convert_to_m_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 141     R>::type
 142 convert_to( A const & a )
 143     {
 144     R r; assign(r,a);
 145     return r;
 146     }
 147
 148 ////////////////////////////////////////////////
 149
 150 namespace
 151 qvm_detail
 152     {
 153     template <int D>
 154     struct
 155     determinant_defined
 156         {
 157         static bool const value=false;
 158         };
 159     }
 160
 161 template <class A>
 162 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 163 typename enable_if_c<
 164     is_mat<A>::value &&
 165     mat_traits<A>::rows==mat_traits<A>::cols &&
 166     !qvm_detail::determinant_defined<mat_traits<A>::rows>::value,
 167     typename mat_traits<A>::scalar_type>::type
 168 determinant( A const & a )
 169     {
 170     return qvm_detail::determinant_impl(a);
 171     }
 172
 173 ////////////////////////////////////////////////
 174
 175 namespace
 176 qvm_detail
 177     {
 178     template <class T,int Dim>
 179     class
 180     identity_mat_
 181         {
 182         identity_mat_( identity_mat_ const & );
 183         identity_mat_ & operator=( identity_mat_ const & );
 184         ~identity_mat_();
 185
 186         public:
 187
 188         template <class R>
 189         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 190         operator R() const
 191             {
 192             R r;
 193             assign(r,*this);
 194             return r;
 195             }
 196         };
 197     }
 198
 199 template <class T,int Dim>
 200 struct
 201 mat_traits< qvm_detail::identity_mat_<T,Dim> >
 202     {
 203     typedef qvm_detail::identity_mat_<T,Dim> this_matrix;
 204     typedef T scalar_type;
 205     static int const rows=Dim;
 206     static int const cols=Dim;
 207
 208     template <int Row,int Col>
 209     static
 210     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 211     scalar_type
 212     read_element( this_matrix const & /*x*/ )
 213         {
 214         BOOST_QVM_STATIC_ASSERT(Row>=0);
 215         BOOST_QVM_STATIC_ASSERT(Row<Dim);
 216         BOOST_QVM_STATIC_ASSERT(Col>=0);
 217         BOOST_QVM_STATIC_ASSERT(Col<Dim);
 218         return scalar_traits<scalar_type>::value(Row==Col);
 219         }
 220
 221     static
 222     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 223     scalar_type
 224     read_element_idx( int row, int col, this_matrix const & /*x*/ )
 225         {
 226         BOOST_QVM_ASSERT(row>=0);
 227         BOOST_QVM_ASSERT(row<Dim);
 228         BOOST_QVM_ASSERT(col>=0);
 229         BOOST_QVM_ASSERT(col<Dim);
 230         return scalar_traits<scalar_type>::value(row==col);
 231         }
 232     };
 233
 234 template <class T,int Dim>
 235 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 236 qvm_detail::identity_mat_<T,Dim> const &
 237 identity_mat()
 238     {
 239     return *(qvm_detail::identity_mat_<T,Dim> const *)qvm_detail::get_valid_ptr_mat_operations();
 240     }
 241
 242 template <class A>
 243 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 244 typename enable_if_c<
 245     is_mat<A>::value &&
 246     mat_traits<A>::rows==mat_traits<A>::cols,
 247     void>::type
 248 set_identity( A & a )
 249     {
 250     assign(a,identity_mat<typename mat_traits<A>::scalar_type,mat_traits<A>::rows>());
 251     }
 252
 253 ////////////////////////////////////////////////
 254
 255 namespace
 256 qvm_detail
 257     {
 258     template <class T>
 259     struct
 260     projection_
 261         {
 262         T const _00;
 263         T const _11;
 264         T const _22;
 265         T const _23;
 266         T const _32;
 267
 268         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 269         projection_( T _00, T _11, T _22, T _23, T _32 ):
 270             _00(_00),
 271             _11(_11),
 272             _22(_22),
 273             _23(_23),
 274             _32(_32)
 275             {
 276             }
 277
 278         template <class R>
 279         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 280         operator R() const
 281             {
 282             R r;
 283             assign(r,*this);
 284             return r;
 285             }
 286         };
 287
 288     template <int Row,int Col>
 289     struct
 290     projection_get
 291         {
 292         template <class T>
 293         static
 294         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 295         T
 296         get( projection_<T> const & )
 297             {
 298             return scalar_traits<T>::value(0);
 299             }
 300         };
 301
 302     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; } };
 303     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; } };
 304     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; } };
 305     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; } };
 306     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; } };
 307     }
 308
 309 template <class T>
 310 struct
 311 mat_traits< qvm_detail::projection_<T> >
 312     {
 313     typedef qvm_detail::projection_<T> this_matrix;
 314     typedef T scalar_type;
 315     static int const rows=4;
 316     static int const cols=4;
 317
 318     template <int Row,int Col>
 319     static
 320     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 321     scalar_type
 322     read_element( this_matrix const & x )
 323         {
 324         BOOST_QVM_STATIC_ASSERT(Row>=0);
 325         BOOST_QVM_STATIC_ASSERT(Row<rows);
 326         BOOST_QVM_STATIC_ASSERT(Col>=0);
 327         BOOST_QVM_STATIC_ASSERT(Col<cols);
 328         return qvm_detail::projection_get<Row,Col>::get(x);
 329         }
 330     };
 331
 332 template <class T>
 333 qvm_detail::projection_<T>
 334 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 335 perspective_lh( T fov_y, T aspect_ratio, T z_near, T z_far )
 336     {
 337     T const one = scalar_traits<T>::value(1);
 338     T const ys = one/tan(fov_y/scalar_traits<T>::value(2));
 339     T const xs = ys/aspect_ratio;
 340     T const zd = z_far-z_near;
 341     T const z1 = z_far/zd;
 342     T const z2 = -z_near*z1;
 343     return qvm_detail::projection_<T>(xs,ys,z1,z2,one);
 344     }
 345
 346 template <class T>
 347 qvm_detail::projection_<T>
 348 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 349 perspective_rh( T fov_y, T aspect_ratio, T z_near, T z_far )
 350     {
 351     T const one = scalar_traits<T>::value(1);
 352     T const ys = one/tan(fov_y/scalar_traits<T>::value(2));
 353     T const xs = ys/aspect_ratio;
 354     T const zd = z_near-z_far;
 355     T const z1 = z_far/zd;
 356     T const z2 = z_near*z1;
 357     return qvm_detail::projection_<T>(xs,ys,z1,z2,-one);
 358     }
 359
 360 ////////////////////////////////////////////////
 361
 362 namespace
 363 qvm_detail
 364     {
 365     template <class OriginalType,class Scalar>
 366     class
 367     matrix_scalar_cast_
 368         {
 369         matrix_scalar_cast_( matrix_scalar_cast_ const & );
 370         matrix_scalar_cast_ & operator=( matrix_scalar_cast_ const & );
 371         ~matrix_scalar_cast_();
 372
 373         public:
 374
 375         template <class T>
 376         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 377         matrix_scalar_cast_ &
 378         operator=( T const & x )
 379             {
 380             assign(*this,x);
 381             return *this;
 382             }
 383
 384         template <class R>
 385         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 386         operator R() const
 387             {
 388             R r;
 389             assign(r,*this);
 390             return r;
 391             }
 392         };
 393
 394     template <bool> struct scalar_cast_matrix_filter { };
 395     template <> struct scalar_cast_matrix_filter<true> { typedef int type; };
 396     }
 397
 398 template <class OriginalType,class Scalar>
 399 struct
 400 mat_traits< qvm_detail::matrix_scalar_cast_<OriginalType,Scalar> >
 401     {
 402     typedef Scalar scalar_type;
 403     typedef qvm_detail::matrix_scalar_cast_<OriginalType,Scalar> this_matrix;
 404     static int const rows=mat_traits<OriginalType>::rows;
 405     static int const cols=mat_traits<OriginalType>::cols;
 406
 407     template <int Row,int Col>
 408     static
 409     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 410     scalar_type
 411     read_element( this_matrix const & x )
 412         {
 413         BOOST_QVM_STATIC_ASSERT(Row>=0);
 414         BOOST_QVM_STATIC_ASSERT(Row<rows);
 415         BOOST_QVM_STATIC_ASSERT(Col>=0);
 416         BOOST_QVM_STATIC_ASSERT(Col<cols);
 417         return scalar_type(mat_traits<OriginalType>::template read_element<Row,Col>(reinterpret_cast<OriginalType const &>(x)));
 418         }
 419
 420     static
 421     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 422     scalar_type
 423     read_element_idx( int row, int col, this_matrix const & x )
 424         {
 425         BOOST_QVM_ASSERT(row>=0);
 426         BOOST_QVM_ASSERT(row<rows);
 427         BOOST_QVM_ASSERT(col>=0);
 428         BOOST_QVM_ASSERT(col<cols);
 429         return scalar_type(mat_traits<OriginalType>::read_element_idx(col,row,reinterpret_cast<OriginalType const &>(x)));
 430         }
 431     };
 432
 433 template <class OriginalType,class Scalar,int R,int C>
 434 struct
 435 deduce_mat<qvm_detail::matrix_scalar_cast_<OriginalType,Scalar>,R,C>
 436     {
 437     typedef mat<Scalar,R,C> type;
 438     };
 439
 440 template <class Scalar,class T>
 441 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 442 qvm_detail::matrix_scalar_cast_<T,Scalar> const &
 443 scalar_cast( T const & x, typename qvm_detail::scalar_cast_matrix_filter<is_mat<T>::value>::type=0 )
 444     {
 445     return reinterpret_cast<qvm_detail::matrix_scalar_cast_<T,Scalar> const &>(x);
 446     }
 447
 448 ////////////////////////////////////////////////
 449
 450 namespace
 451 qvm_detail
 452     {
 453     template <int M,int N>
 454     struct
 455     div_eq_ms_defined
 456         {
 457         static bool const value=false;
 458         };
 459     }
 460
 461 template <class A,class B>
 462 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 463 typename enable_if_c<
 464     is_mat<A>::value && is_scalar<B>::value &&
 465     !qvm_detail::div_eq_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 466     A &>::type
 467 operator/=( A & a, B b )
 468     {
 469     for( int i=0; i!=mat_traits<A>::rows; ++i )
 470         for( int j=0; j!=mat_traits<A>::cols; ++j )
 471             mat_traits<A>::write_element_idx(i,j,a)/=b;
 472     return a;
 473     }
 474
 475 ////////////////////////////////////////////////
 476
 477 namespace
 478 qvm_detail
 479     {
 480     template <int M,int N>
 481     struct
 482     div_ms_defined
 483         {
 484         static bool const value=false;
 485         };
 486     }
 487
 488 template <class A,class B>
 489 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 490 typename lazy_enable_if_c<
 491     is_mat<A>::value && is_scalar<B>::value &&
 492     !qvm_detail::div_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 493     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
 494 operator/( A const & a, B b )
 495     {
 496     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
 497     R r;
 498     for( int i=0; i!=mat_traits<A>::rows; ++i )
 499         for( int j=0; j!=mat_traits<A>::cols; ++j )
 500             mat_traits<R>::write_element_idx(i,j,r)=mat_traits<A>::read_element_idx(i,j,a)/b;
 501     return r;
 502     }
 503
 504 ////////////////////////////////////////////////
 505
 506 namespace
 507 qvm_detail
 508     {
 509     template <int M,int N>
 510     struct
 511     eq_mm_defined
 512         {
 513         static bool const value=false;
 514         };
 515     }
 516
 517 template <class A,class B>
 518 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 519 typename enable_if_c<
 520     is_mat<A>::value && is_mat<B>::value &&
 521     mat_traits<A>::rows==mat_traits<B>::rows &&
 522     mat_traits<A>::cols==mat_traits<B>::cols &&
 523     !qvm_detail::eq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 524     bool>::type
 525 operator==( A const & a, B const & b )
 526     {
 527     for( int i=0; i!=mat_traits<A>::rows; ++i )
 528         for( int j=0; j!=mat_traits<A>::cols; ++j )
 529             if( mat_traits<A>::read_element_idx(i,j,a)!=mat_traits<B>::read_element_idx(i,j,b) )
 530                 return false;
 531     return true;
 532     }
 533
 534 ////////////////////////////////////////////////
 535
 536 namespace
 537 qvm_detail
 538     {
 539     template <int M,int N>
 540     struct
 541     minus_eq_mm_defined
 542         {
 543         static bool const value=false;
 544         };
 545     }
 546
 547 template <class A,class B>
 548 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 549 typename enable_if_c<
 550     is_mat<A>::value && is_mat<B>::value &&
 551     mat_traits<A>::rows==mat_traits<B>::rows &&
 552     mat_traits<A>::cols==mat_traits<B>::cols &&
 553     !qvm_detail::minus_eq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 554     A &>::type
 555 operator-=( A & a, B const & b )
 556     {
 557     for( int i=0; i!=mat_traits<A>::rows; ++i )
 558         for( int j=0; j!=mat_traits<A>::cols; ++j )
 559             mat_traits<A>::write_element_idx(i,j,a)-=mat_traits<B>::read_element_idx(i,j,b);
 560     return a;
 561     }
 562
 563 ////////////////////////////////////////////////
 564
 565 namespace
 566 qvm_detail
 567     {
 568     template <int M,int N>
 569     struct
 570     minus_m_defined
 571         {
 572         static bool const value=false;
 573         };
 574     }
 575
 576 template <class A>
 577 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 578 typename lazy_enable_if_c<
 579     is_mat<A>::value &&
 580     !qvm_detail::minus_m_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 581     deduce_mat<A> >::type
 582 operator-( A const & a )
 583     {
 584     typedef typename deduce_mat<A>::type R;
 585     R r;
 586     for( int i=0; i!=mat_traits<A>::rows; ++i )
 587         for( int j=0; j!=mat_traits<A>::cols; ++j )
 588             mat_traits<R>::write_element_idx(i,j,r)=-mat_traits<A>::read_element_idx(i,j,a);
 589     return r;
 590     }
 591
 592 ////////////////////////////////////////////////
 593
 594 namespace
 595 qvm_detail
 596     {
 597     template <int M,int N>
 598     struct
 599     minus_mm_defined
 600         {
 601         static bool const value=false;
 602         };
 603     }
 604
 605 template <class A,class B>
 606 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 607 typename lazy_enable_if_c<
 608     is_mat<A>::value && is_mat<B>::value &&
 609     mat_traits<A>::rows==mat_traits<B>::rows &&
 610     mat_traits<A>::cols==mat_traits<B>::cols &&
 611     !qvm_detail::minus_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 612     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
 613 operator-( A const & a, B const & b )
 614     {
 615     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
 616     R r;
 617     for( int i=0; i!=mat_traits<A>::rows; ++i )
 618         for( int j=0; j!=mat_traits<A>::cols; ++j )
 619             mat_traits<R>::write_element_idx(i,j,r)=mat_traits<A>::read_element_idx(i,j,a)-mat_traits<B>::read_element_idx(i,j,b);
 620     return r;
 621     }
 622
 623 ////////////////////////////////////////////////
 624
 625 namespace
 626 qvm_detail
 627     {
 628     template <int D>
 629     struct
 630     mul_eq_mm_defined
 631         {
 632         static bool const value=false;
 633         };
 634     }
 635
 636 template <class A,class B>
 637 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 638 typename enable_if_c<
 639     is_mat<A>::value &&
 640     is_mat<B>::value &&
 641     mat_traits<A>::rows==mat_traits<A>::cols &&
 642     mat_traits<A>::rows==mat_traits<B>::rows &&
 643     mat_traits<A>::cols==mat_traits<B>::cols &&
 644     !qvm_detail::mul_eq_mm_defined<mat_traits<A>::rows>::value,
 645     A &>::type
 646 operator*=( A & r, B const & b )
 647     {
 648     typedef typename mat_traits<A>::scalar_type Ta;
 649     Ta a[mat_traits<A>::rows][mat_traits<A>::cols];
 650     for( int i=0; i<mat_traits<A>::rows; ++i )
 651         for( int j=0; j<mat_traits<B>::cols; ++j )
 652             a[i][j]=mat_traits<A>::read_element_idx(i,j,r);
 653     for( int i=0; i<mat_traits<A>::rows; ++i )
 654         for( int j=0; j<mat_traits<B>::cols; ++j )
 655             {
 656             Ta x(scalar_traits<Ta>::value(0));
 657             for( int k=0; k<mat_traits<A>::cols; ++k )
 658                 x += a[i][k]*mat_traits<B>::read_element_idx(k,j,b);
 659             mat_traits<A>::write_element_idx(i,j,r) = x;
 660             }
 661     return r;
 662     }
 663
 664 ////////////////////////////////////////////////
 665
 666 namespace
 667 qvm_detail
 668     {
 669     template <int M,int N>
 670     struct
 671     mul_eq_ms_defined
 672         {
 673         static bool const value=false;
 674         };
 675     }
 676
 677 template <class A,class B>
 678 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 679 typename enable_if_c<
 680     is_mat<A>::value && is_scalar<B>::value &&
 681     !qvm_detail::mul_eq_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 682     A &>::type
 683 operator*=( A & a, B b )
 684     {
 685     for( int i=0; i!=mat_traits<A>::rows; ++i )
 686         for( int j=0; j!=mat_traits<A>::cols; ++j )
 687             mat_traits<A>::write_element_idx(i,j,a)*=b;
 688     return a;
 689     }
 690
 691 ////////////////////////////////////////////////
 692
 693 namespace
 694 qvm_detail
 695     {
 696     template <int R,int /*CR*/,int C>
 697     struct
 698     mul_mm_defined
 699         {
 700         static bool const value=false;
 701         };
 702     }
 703
 704 template <class A,class B>
 705 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 706 typename lazy_enable_if_c<
 707     is_mat<A>::value && is_mat<B>::value &&
 708     mat_traits<A>::cols==mat_traits<B>::rows &&
 709     !qvm_detail::mul_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols,mat_traits<B>::cols>::value,
 710     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<B>::cols> >::type
 711 operator*( A const & a, B const & b )
 712     {
 713     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<B>::cols>::type R;
 714     R r;
 715     for( int i=0; i<mat_traits<A>::rows; ++i )
 716         for( int j=0; j<mat_traits<B>::cols; ++j )
 717             {
 718             typedef typename mat_traits<A>::scalar_type Ta;
 719             Ta x(scalar_traits<Ta>::value(0));
 720             for( int k=0; k<mat_traits<A>::cols; ++k )
 721                 x += mat_traits<A>::read_element_idx(i,k,a)*mat_traits<B>::read_element_idx(k,j,b);
 722             mat_traits<R>::write_element_idx(i,j,r) = x;
 723             }
 724     return r;
 725     }
 726
 727 ////////////////////////////////////////////////
 728
 729 namespace
 730 qvm_detail
 731     {
 732     template <int M,int N>
 733     struct
 734     mul_ms_defined
 735         {
 736         static bool const value=false;
 737         };
 738     }
 739
 740 template <class A,class B>
 741 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 742 typename lazy_enable_if_c<
 743     is_mat<A>::value && is_scalar<B>::value &&
 744     !qvm_detail::mul_ms_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 745     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
 746 operator*( A const & a, B b )
 747     {
 748     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
 749     R r;
 750     for( int i=0; i!=mat_traits<A>::rows; ++i )
 751         for( int j=0; j!=mat_traits<A>::cols; ++j )
 752             mat_traits<R>::write_element_idx(i,j,r)=mat_traits<A>::read_element_idx(i,j,a)*b;
 753     return r;
 754     }
 755
 756 ////////////////////////////////////////////////
 757
 758 namespace
 759 qvm_detail
 760     {
 761     template <int M,int N>
 762     struct
 763     mul_sm_defined
 764         {
 765         static bool const value=false;
 766         };
 767     }
 768
 769 template <class A,class B>
 770 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 771 typename lazy_enable_if_c<
 772     is_scalar<A>::value && is_mat<B>::value &&
 773     !qvm_detail::mul_sm_defined<mat_traits<B>::rows,mat_traits<B>::cols>::value,
 774     deduce_mat2<A,B,mat_traits<B>::rows,mat_traits<B>::cols> >::type
 775 operator*( A a, B const & b )
 776     {
 777     typedef typename deduce_mat2<A,B,mat_traits<B>::rows,mat_traits<B>::cols>::type R;
 778     R r;
 779     for( int i=0; i!=mat_traits<B>::rows; ++i )
 780         for( int j=0; j!=mat_traits<B>::cols; ++j )
 781             mat_traits<R>::write_element_idx(i,j,r)=a*mat_traits<B>::read_element_idx(i,j,b);
 782     return r;
 783     }
 784
 785 ////////////////////////////////////////////////
 786
 787 namespace
 788 qvm_detail
 789     {
 790     template <int M,int N>
 791     struct
 792     neq_mm_defined
 793         {
 794         static bool const value=false;
 795         };
 796     }
 797
 798 template <class A,class B>
 799 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 800 typename enable_if_c<
 801     is_mat<A>::value && is_mat<B>::value &&
 802     mat_traits<A>::rows==mat_traits<B>::rows &&
 803     mat_traits<A>::cols==mat_traits<B>::cols &&
 804     !qvm_detail::neq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 805     bool>::type
 806 operator!=( A const & a, B const & b )
 807     {
 808     for( int i=0; i!=mat_traits<A>::rows; ++i )
 809         for( int j=0; j!=mat_traits<A>::cols; ++j )
 810             if( mat_traits<A>::read_element_idx(i,j,a)!=mat_traits<B>::read_element_idx(i,j,b) )
 811                 return true;
 812     return false;
 813     }
 814
 815 ////////////////////////////////////////////////
 816
 817 namespace
 818 qvm_detail
 819     {
 820     template <int M,int N>
 821     struct
 822     plus_eq_mm_defined
 823         {
 824         static bool const value=false;
 825         };
 826     }
 827
 828 template <class A,class B>
 829 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 830 typename enable_if_c<
 831     is_mat<A>::value && is_mat<B>::value &&
 832     mat_traits<A>::rows==mat_traits<B>::rows &&
 833     mat_traits<A>::cols==mat_traits<B>::cols &&
 834     !qvm_detail::plus_eq_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 835     A &>::type
 836 operator+=( A & a, B const & b )
 837     {
 838     for( int i=0; i!=mat_traits<A>::rows; ++i )
 839         for( int j=0; j!=mat_traits<A>::cols; ++j )
 840             mat_traits<A>::write_element_idx(i,j,a)+=mat_traits<B>::read_element_idx(i,j,b);
 841     return a;
 842     }
 843
 844 ////////////////////////////////////////////////
 845
 846 namespace
 847 qvm_detail
 848     {
 849     template <int M,int N>
 850     struct
 851     plus_mm_defined
 852         {
 853         static bool const value=false;
 854         };
 855     }
 856
 857 template <class A,class B>
 858 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
 859 typename lazy_enable_if_c<
 860     is_mat<A>::value && is_mat<B>::value &&
 861     mat_traits<A>::rows==mat_traits<B>::rows &&
 862     mat_traits<A>::cols==mat_traits<B>::cols &&
 863     !qvm_detail::plus_mm_defined<mat_traits<A>::rows,mat_traits<A>::cols>::value,
 864     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
 865 operator+( A const & a, B const & b )
 866     {
 867     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type R;
 868     R r;
 869     for( int i=0; i!=mat_traits<A>::rows; ++i )
 870         for( int j=0; j!=mat_traits<A>::cols; ++j )
 871             mat_traits<R>::write_element_idx(i,j,r)=mat_traits<A>::read_element_idx(i,j,a)+mat_traits<B>::read_element_idx(i,j,b);
 872     return r;
 873     }
 874
 875 ////////////////////////////////////////////////
 876
 877 namespace
 878 qvm_detail
 879     {
 880     template <class T>
 881     class
 882     mref_
 883         {
 884         mref_( mref_ const & );
 885         mref_ & operator=( mref_ const & );
 886         ~mref_();
 887
 888         public:
 889
 890         template <class R>
 891         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 892         mref_ &
 893         operator=( R const & x )
 894             {
 895             assign(*this,x);
 896             return *this;
 897             }
 898
 899         template <class R>
 900         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 901         operator R() const
 902             {
 903             R r;
 904             assign(r,*this);
 905             return r;
 906             }
 907         };
 908     }
 909
 910 template <class M>
 911 struct
 912 mat_traits< qvm_detail::mref_<M> >
 913     {
 914     typedef typename mat_traits<M>::scalar_type scalar_type;
 915     typedef qvm_detail::mref_<M> this_matrix;
 916     static int const rows=mat_traits<M>::rows;
 917     static int const cols=mat_traits<M>::cols;
 918
 919     template <int Row,int Col>
 920     static
 921     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 922     scalar_type
 923     read_element( this_matrix const & x )
 924         {
 925         BOOST_QVM_STATIC_ASSERT(Row>=0);
 926         BOOST_QVM_STATIC_ASSERT(Row<rows);
 927         BOOST_QVM_STATIC_ASSERT(Col>=0);
 928         BOOST_QVM_STATIC_ASSERT(Col<cols);
 929         return mat_traits<M>::template read_element<Row,Col>(reinterpret_cast<M const &>(x));
 930         }
 931
 932     template <int Row,int Col>
 933     static
 934     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 935     scalar_type &
 936     write_element( this_matrix & x )
 937         {
 938         BOOST_QVM_STATIC_ASSERT(Row>=0);
 939         BOOST_QVM_STATIC_ASSERT(Row<rows);
 940         BOOST_QVM_STATIC_ASSERT(Col>=0);
 941         BOOST_QVM_STATIC_ASSERT(Col<cols);
 942         return mat_traits<M>::template write_element<Row,Col>(reinterpret_cast<M &>(x));
 943         }
 944
 945     static
 946     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 947     scalar_type
 948     read_element_idx( int row, int col, this_matrix const & x )
 949         {
 950         BOOST_QVM_ASSERT(row>=0);
 951         BOOST_QVM_ASSERT(row<rows);
 952         BOOST_QVM_ASSERT(col>=0);
 953         BOOST_QVM_ASSERT(col<cols);
 954         return mat_traits<M>::read_element_idx(row,col,reinterpret_cast<M const &>(x));
 955         }
 956
 957     static
 958     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
 959     scalar_type &
 960     write_element_idx( int row, int col, this_matrix & x )
 961         {
 962         BOOST_QVM_ASSERT(row>=0);
 963         BOOST_QVM_ASSERT(row<rows);
 964         BOOST_QVM_ASSERT(col>=0);
 965         BOOST_QVM_ASSERT(col<cols);
 966         return mat_traits<M>::write_element_idx(row,col,reinterpret_cast<M &>(x));
 967         }
 968     };
 969
 970 template <class M,int R,int C>
 971 struct
 972 deduce_mat<qvm_detail::mref_<M>,R,C>
 973     {
 974     typedef mat<typename mat_traits<M>::scalar_type,R,C> type;
 975     };
 976
 977 template <class M>
 978 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 979 typename enable_if_c<
 980     is_mat<M>::value,
 981     qvm_detail::mref_<M> const &>::type
 982 mref( M const & a )
 983     {
 984     return reinterpret_cast<qvm_detail::mref_<M> const &>(a);
 985     }
 986
 987 template <class M>
 988 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
 989 typename enable_if_c<
 990     is_mat<M>::value,
 991     qvm_detail::mref_<M> &>::type
 992 mref( M & a )
 993     {
 994     return reinterpret_cast<qvm_detail::mref_<M> &>(a);
 995     }
 996
 997 ////////////////////////////////////////////////
 998
 999 namespace
1000 qvm_detail
1001     {
1002     template <class T,int Rows,int Cols>
1003     class
1004     zero_mat_
1005         {
1006         zero_mat_( zero_mat_ const & );
1007         zero_mat_ & operator=( zero_mat_ const & );
1008         ~zero_mat_();
1009
1010         public:
1011
1012         template <class R>
1013         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1014         operator R() const
1015             {
1016             R r;
1017             assign(r,*this);
1018             return r;
1019             }
1020         };
1021     }
1022
1023 template <class T,int Rows,int Cols>
1024 struct
1025 mat_traits< qvm_detail::zero_mat_<T,Rows,Cols> >
1026     {
1027     typedef qvm_detail::zero_mat_<T,Rows,Cols> this_matrix;
1028     typedef T scalar_type;
1029     static int const rows=Rows;
1030     static int const cols=Cols;
1031
1032     template <int Row,int Col>
1033     static
1034     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1035     scalar_type
1036     read_element( this_matrix const & )
1037         {
1038         BOOST_QVM_STATIC_ASSERT(Row>=0);
1039         BOOST_QVM_STATIC_ASSERT(Row<Rows);
1040         BOOST_QVM_STATIC_ASSERT(Col>=0);
1041         BOOST_QVM_STATIC_ASSERT(Col<Cols);
1042         return scalar_traits<scalar_type>::value(0);
1043         }
1044
1045     static
1046     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1047     scalar_type
1048     read_element_idx( int row, int col, this_matrix const & )
1049         {
1050         BOOST_QVM_ASSERT(row>=0);
1051         BOOST_QVM_ASSERT(row<rows);
1052         BOOST_QVM_ASSERT(col>=0);
1053         BOOST_QVM_ASSERT(col<cols);
1054         return scalar_traits<scalar_type>::value(0);
1055         }
1056     };
1057
1058 template <class T,int Rows,int Cols,int R,int C>
1059 struct
1060 deduce_mat<qvm_detail::zero_mat_<T,Rows,Cols>,R,C>
1061     {
1062     typedef mat<T,R,C> type;
1063     };
1064
1065 template <class T,int Rows,int Cols>
1066 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1067 qvm_detail::zero_mat_<T,Rows,Cols> const &
1068 zero_mat()
1069     {
1070     return *(qvm_detail::zero_mat_<T,Rows,Cols> const *)qvm_detail::get_valid_ptr_mat_operations();
1071     }
1072
1073 template <class T,int Dim>
1074 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1075 qvm_detail::zero_mat_<T,Dim,Dim> const &
1076 zero_mat()
1077     {
1078     return *(qvm_detail::zero_mat_<T,Dim,Dim> const *)qvm_detail::get_valid_ptr_mat_operations();
1079     }
1080
1081 template <class A>
1082 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1083 typename enable_if_c<
1084     is_mat<A>::value,
1085     void>::type
1086 set_zero( A & a )
1087     {
1088     assign(a,zero_mat<typename mat_traits<A>::scalar_type,mat_traits<A>::rows,mat_traits<A>::cols>());
1089     }
1090
1091 ////////////////////////////////////////////////
1092
1093 namespace
1094 qvm_detail
1095     {
1096     template <int D,class S>
1097     struct
1098     rot_mat_
1099         {
1100         typedef S scalar_type;
1101         scalar_type a[3][3];
1102
1103         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1104         rot_mat_(
1105                 scalar_type a00, scalar_type a01, scalar_type a02,
1106                 scalar_type a10, scalar_type a11, scalar_type a12,
1107                 scalar_type a20, scalar_type a21, scalar_type a22 )
1108             {
1109             a[0][0] = a00;
1110             a[0][1] = a01;
1111             a[0][2] = a02;
1112             a[1][0] = a10;
1113             a[1][1] = a11;
1114             a[1][2] = a12;
1115             a[2][0] = a20;
1116             a[2][1] = a21;
1117             a[2][2] = a22;
1118             }
1119
1120         template <class R>
1121         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1122         operator R() const
1123             {
1124             R r;
1125             assign(r,*this);
1126             return r;
1127             }
1128         };
1129
1130     template <int Row,int Col>
1131     struct
1132     rot_m_get
1133         {
1134         template <class T>
1135         static
1136         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1137         T
1138         get( T const (&)[3][3] )
1139             {
1140             return scalar_traits<T>::value(Row==Col);
1141             }
1142         };
1143
1144     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]; } };
1145     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]; } };
1146     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]; } };
1147     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]; } };
1148     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]; } };
1149     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]; } };
1150     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]; } };
1151     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]; } };
1152     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]; } };
1153     }
1154
1155 template <class M>
1156 struct mat_traits;
1157
1158 template <int D,class S>
1159 struct
1160 mat_traits< qvm_detail::rot_mat_<D,S> >
1161     {
1162     typedef qvm_detail::rot_mat_<D,S> this_matrix;
1163     typedef typename this_matrix::scalar_type scalar_type;
1164     static int const rows=D;
1165     static int const cols=D;
1166
1167     template <int Row,int Col>
1168     static
1169     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1170     scalar_type
1171     read_element( this_matrix const & x )
1172         {
1173         BOOST_QVM_STATIC_ASSERT(Row>=0);
1174         BOOST_QVM_STATIC_ASSERT(Row<D);
1175         BOOST_QVM_STATIC_ASSERT(Col>=0);
1176         BOOST_QVM_STATIC_ASSERT(Col<D);
1177         return qvm_detail::rot_m_get<Row,Col>::get(x.a);
1178         }
1179
1180     static
1181     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1182     scalar_type
1183     read_element_idx( int row, int col, this_matrix const & x )
1184         {
1185         BOOST_QVM_ASSERT(row>=0);
1186         BOOST_QVM_ASSERT(row<D);
1187         BOOST_QVM_ASSERT(col>=0);
1188         BOOST_QVM_ASSERT(col<D);
1189         return row<3 && col<3?
1190             x.a[row][col] :
1191             scalar_traits<scalar_type>::value(row==col);
1192         }
1193     };
1194
1195 template <int Dim,class V,class Angle>
1196 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1197 typename enable_if_c<
1198     is_vec<V>::value && vec_traits<V>::dim==3,
1199     qvm_detail::rot_mat_<Dim,Angle> >::type
1200 rot_mat( V const & axis, Angle angle )
1201     {
1202     typedef Angle scalar_type;
1203     scalar_type const x=vec_traits<V>::template read_element<0>(axis);
1204     scalar_type const y=vec_traits<V>::template read_element<1>(axis);
1205     scalar_type const z=vec_traits<V>::template read_element<2>(axis);
1206     scalar_type const m2=x*x+y*y+z*z;
1207     if( m2==scalar_traits<scalar_type>::value(0) )
1208         BOOST_QVM_THROW_EXCEPTION(zero_magnitude_error());
1209     scalar_type const s = sin(angle);
1210     scalar_type const c = cos(angle);
1211     scalar_type const x2 = x*x;
1212     scalar_type const y2 = y*y;
1213     scalar_type const z2 = z*z;
1214     scalar_type const xy = x*y;
1215     scalar_type const xz = x*z;
1216     scalar_type const yz = y*z;
1217     scalar_type const xs = x*s;
1218     scalar_type const ys = y*s;
1219     scalar_type const zs = z*s;
1220     scalar_type const one = scalar_traits<scalar_type>::value(1);
1221     scalar_type const c1 = one-c;
1222     return qvm_detail::rot_mat_<Dim,Angle>(
1223         x2+(one-x2)*c, xy*c1-zs, xz*(one-c)+ys,
1224         xy*c1+zs, y2+(one-y2)*c, yz*c1-xs,
1225         xz*c1-ys, yz*c1+xs, z2+(one-z2)*c );
1226     }
1227
1228 template <class A,class B,class Angle>
1229 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1230 typename enable_if_c<
1231     is_mat<A>::value &&
1232     mat_traits<A>::rows==mat_traits<A>::cols &&
1233     mat_traits<A>::rows>=3 &&
1234     is_vec<B>::value && vec_traits<B>::dim==3,
1235     void>::type
1236 set_rot( A & a, B const & axis, Angle angle )
1237     {
1238     assign(a,rot_mat<mat_traits<A>::rows>(axis,angle));
1239     }
1240
1241 template <class A,class B,class Angle>
1242 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1243 typename enable_if_c<
1244     is_mat<A>::value &&
1245     mat_traits<A>::rows==mat_traits<A>::cols &&
1246     mat_traits<A>::rows>=3 &&
1247     is_vec<B>::value && vec_traits<B>::dim==3,
1248     void>::type
1249 rotate( A & a, B const & axis, Angle angle )
1250     {
1251     a *= rot_mat<mat_traits<A>::rows>(axis,angle);
1252     }
1253
1254 ////////////////////////////////////////////////
1255
1256 template <int Dim,class Angle>
1257 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1258 qvm_detail::rot_mat_<Dim,Angle>
1259 rot_mat_xzy( Angle x1, Angle z2, Angle y3 )
1260     {
1261     typedef Angle scalar_type;
1262     scalar_type const c1 = cos(x1);
1263     scalar_type const s1 = sin(x1);
1264     scalar_type const c2 = cos(z2);
1265     scalar_type const s2 = sin(z2);
1266     scalar_type const c3 = cos(y3);
1267     scalar_type const s3 = sin(y3);
1268     return qvm_detail::rot_mat_<Dim,Angle>(
1269         c2*c3, -s2, c2*s3,
1270         s1*s3 + c1*c3*s2, c1*c2, c1*s2*s3 - c3*s1,
1271         c3*s1*s2 - c1*s3, c2*s1, c1*c3 + s1*s2*s3 );
1272     }
1273
1274 template <class A,class Angle>
1275 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1276 typename enable_if_c<
1277     is_mat<A>::value &&
1278     mat_traits<A>::rows==mat_traits<A>::cols &&
1279     mat_traits<A>::rows>=3,
1280     void>::type
1281 set_rot_xzy( A & a, Angle x1, Angle z2, Angle y3 )
1282     {
1283     assign(a,rot_mat_xzy<mat_traits<A>::rows>(x1,z2,y3));
1284     }
1285
1286 template <class A,class Angle>
1287 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1288 typename enable_if_c<
1289     is_mat<A>::value &&
1290     mat_traits<A>::rows==mat_traits<A>::cols &&
1291     mat_traits<A>::rows>=3,
1292     void>::type
1293 rotate_xzy( A & a, Angle x1, Angle z2, Angle y3 )
1294     {
1295     a *= rot_mat_xzy<mat_traits<A>::rows>(x1,z2,y3);
1296     }
1297
1298 ////////////////////////////////////////////////
1299
1300 template <int Dim,class Angle>
1301 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1302 qvm_detail::rot_mat_<Dim,Angle>
1303 rot_mat_xyz( Angle x1, Angle y2, Angle z3 )
1304     {
1305     typedef Angle scalar_type;
1306     scalar_type const c1 = cos(x1);
1307     scalar_type const s1 = sin(x1);
1308     scalar_type const c2 = cos(y2);
1309     scalar_type const s2 = sin(y2);
1310     scalar_type const c3 = cos(z3);
1311     scalar_type const s3 = sin(z3);
1312     return qvm_detail::rot_mat_<Dim,Angle>(
1313         c2*c3, -c2*s3, s2,
1314         c1*s3 + c3*s1*s2, c1*c3 - s1*s2*s3, -c2*s1,
1315         s1*s3 - c1*c3*s2, c3*s1 + c1*s2*s3, c1*c2 );
1316     }
1317
1318 template <class A,class Angle>
1319 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1320 typename enable_if_c<
1321     is_mat<A>::value &&
1322     mat_traits<A>::rows==mat_traits<A>::cols &&
1323     mat_traits<A>::rows>=3,
1324     void>::type
1325 set_rot_xyz( A & a, Angle x1, Angle y2, Angle z3 )
1326     {
1327     assign(a,rot_mat_xyz<mat_traits<A>::rows>(x1,y2,z3));
1328     }
1329
1330 template <class A,class Angle>
1331 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1332 typename enable_if_c<
1333     is_mat<A>::value &&
1334     mat_traits<A>::rows==mat_traits<A>::cols &&
1335     mat_traits<A>::rows>=3,
1336     void>::type
1337 rotate_xyz( A & a, Angle x1, Angle y2, Angle z3 )
1338     {
1339     a *= rot_mat_xyz<mat_traits<A>::rows>(x1,y2,z3);
1340     }
1341
1342 ////////////////////////////////////////////////
1343
1344 template <int Dim,class Angle>
1345 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1346 qvm_detail::rot_mat_<Dim,Angle>
1347 rot_mat_yxz( Angle y1, Angle x2, Angle z3 )
1348     {
1349     typedef Angle scalar_type;
1350     scalar_type const c1 = cos(y1);
1351     scalar_type const s1 = sin(y1);
1352     scalar_type const c2 = cos(x2);
1353     scalar_type const s2 = sin(x2);
1354     scalar_type const c3 = cos(z3);
1355     scalar_type const s3 = sin(z3);
1356     return qvm_detail::rot_mat_<Dim,Angle>(
1357         c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3, c2*s1,
1358         c2*s3, c2*c3, -s2,
1359         c1*s2*s3 - c3*s1, c1*c3*s2 + s1*s3, c1*c2 );
1360     }
1361
1362 template <class A,class Angle>
1363 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1364 typename enable_if_c<
1365     is_mat<A>::value &&
1366     mat_traits<A>::rows==mat_traits<A>::cols &&
1367     mat_traits<A>::rows>=3,
1368     void>::type
1369 set_rot_yxz( A & a, Angle y1, Angle x2, Angle z3 )
1370     {
1371     assign(a,rot_mat_yxz<mat_traits<A>::rows>(y1,x2,z3));
1372     }
1373
1374 template <class A,class Angle>
1375 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1376 typename enable_if_c<
1377     is_mat<A>::value &&
1378     mat_traits<A>::rows==mat_traits<A>::cols &&
1379     mat_traits<A>::rows>=3,
1380     void>::type
1381 rotate_yxz( A & a, Angle y1, Angle x2, Angle z3 )
1382     {
1383     a *= rot_mat_yxz<mat_traits<A>::rows>(y1,x2,z3);
1384     }
1385
1386 ////////////////////////////////////////////////
1387
1388 template <int Dim,class Angle>
1389 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1390 qvm_detail::rot_mat_<Dim,Angle>
1391 rot_mat_yzx( Angle y1, Angle z2, Angle x3 )
1392     {
1393     typedef Angle scalar_type;
1394     scalar_type const c1 = cos(y1);
1395     scalar_type const s1 = sin(y1);
1396     scalar_type const c2 = cos(z2);
1397     scalar_type const s2 = sin(z2);
1398     scalar_type const c3 = cos(x3);
1399     scalar_type const s3 = sin(x3);
1400     return qvm_detail::rot_mat_<Dim,Angle>(
1401         c1*c2, s1*s3 - c1*c3*s2, c3*s1 + c1*s2*s3,
1402         s2, c2*c3, -c2*s3,
1403         -c2*s1, c1*s3 + c3*s1*s2, c1*c3 - s1*s2*s3 );
1404     }
1405
1406 template <class A,class Angle>
1407 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1408 typename enable_if_c<
1409     is_mat<A>::value &&
1410     mat_traits<A>::rows==mat_traits<A>::cols &&
1411     mat_traits<A>::rows>=3,
1412     void>::type
1413 set_rot_yzx( A & a, Angle y1, Angle z2, Angle x3 )
1414     {
1415     assign(a,rot_mat_yzx<mat_traits<A>::rows>(y1,z2,x3));
1416     }
1417
1418 template <class A,class Angle>
1419 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1420 typename enable_if_c<
1421     is_mat<A>::value &&
1422     mat_traits<A>::rows==mat_traits<A>::cols &&
1423     mat_traits<A>::rows>=3,
1424     void>::type
1425 rotate_yzx( A & a, Angle y1, Angle z2, Angle x3 )
1426     {
1427     a *= rot_mat_yzx<mat_traits<A>::rows>(y1,z2,x3);
1428     }
1429
1430 ////////////////////////////////////////////////
1431
1432 template <int Dim,class Angle>
1433 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1434 qvm_detail::rot_mat_<Dim,Angle>
1435 rot_mat_zyx( Angle z1, Angle y2, Angle x3 )
1436     {
1437     typedef Angle scalar_type;
1438     scalar_type const c1 = cos(z1);
1439     scalar_type const s1 = sin(z1);
1440     scalar_type const c2 = cos(y2);
1441     scalar_type const s2 = sin(y2);
1442     scalar_type const c3 = cos(x3);
1443     scalar_type const s3 = sin(x3);
1444     return qvm_detail::rot_mat_<Dim,Angle>(
1445         c1*c2, c1*s2*s3 - c3*s1, s1*s3 + c1*c3*s2,
1446         c2*s1, c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3,
1447         -s2, c2*s3, c2*c3 );
1448     }
1449
1450 template <class A,class Angle>
1451 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1452 typename enable_if_c<
1453     is_mat<A>::value &&
1454     mat_traits<A>::rows==mat_traits<A>::cols &&
1455     mat_traits<A>::rows>=3,
1456     void>::type
1457 set_rot_zyx( A & a, Angle z1, Angle y2, Angle x3 )
1458     {
1459     assign(a,rot_mat_zyx<mat_traits<A>::rows>(z1,y2,x3));
1460     }
1461
1462 template <class A,class Angle>
1463 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1464 typename enable_if_c<
1465     is_mat<A>::value &&
1466     mat_traits<A>::rows==mat_traits<A>::cols &&
1467     mat_traits<A>::rows>=3,
1468     void>::type
1469 rotate_zyx( A & a, Angle z1, Angle y2, Angle x3 )
1470     {
1471     a *= rot_mat_zyx<mat_traits<A>::rows>(z1,y2,x3);
1472     }
1473
1474 ////////////////////////////////////////////////
1475
1476 template <int Dim,class Angle>
1477 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1478 qvm_detail::rot_mat_<Dim,Angle>
1479 rot_mat_zxy( Angle z1, Angle x2, Angle y3 )
1480     {
1481     typedef Angle scalar_type;
1482     scalar_type const c1 = cos(z1);
1483     scalar_type const s1 = sin(z1);
1484     scalar_type const c2 = cos(x2);
1485     scalar_type const s2 = sin(x2);
1486     scalar_type const c3 = cos(y3);
1487     scalar_type const s3 = sin(y3);
1488     return qvm_detail::rot_mat_<Dim,Angle>(
1489         c1*c3 - s1*s2*s3, -c2*s1, c1*s3 + c3*s1*s2,
1490         c3*s1 + c1*s2*s3, c1*c2, s1*s3 - c1*c3*s2,
1491         -c2*s3, s2, c2*c3 );
1492     }
1493
1494 template <class A,class Angle>
1495 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1496 typename enable_if_c<
1497     is_mat<A>::value &&
1498     mat_traits<A>::rows==mat_traits<A>::cols &&
1499     mat_traits<A>::rows>=3,
1500     void>::type
1501 set_rot_zxy( A & a, Angle z1, Angle x2, Angle y3 )
1502     {
1503     assign(a,rot_mat_zxy<mat_traits<A>::rows>(z1,x2,y3));
1504     }
1505
1506 template <class A,class Angle>
1507 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1508 typename enable_if_c<
1509     is_mat<A>::value &&
1510     mat_traits<A>::rows==mat_traits<A>::cols &&
1511     mat_traits<A>::rows>=3,
1512     void>::type
1513 rotate_zxy( A & a, Angle z1, Angle x2, Angle y3 )
1514     {
1515     a *= rot_mat_zxy<mat_traits<A>::rows>(z1,x2,y3);
1516     }
1517
1518 ////////////////////////////////////////////////
1519
1520 template <int Dim,class Angle>
1521 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1522 qvm_detail::rot_mat_<Dim,Angle>
1523 rot_mat_xzx( Angle x1, Angle z2, Angle x3 )
1524     {
1525     typedef Angle scalar_type;
1526     scalar_type const c1 = cos(x1);
1527     scalar_type const s1 = sin(x1);
1528     scalar_type const c2 = cos(z2);
1529     scalar_type const s2 = sin(z2);
1530     scalar_type const c3 = cos(x3);
1531     scalar_type const s3 = sin(x3);
1532     return qvm_detail::rot_mat_<Dim,Angle>(
1533         c2, -c3*s2, s2*s3,
1534         c1*s2, c1*c2*c3 - s1*s3, -c3*s1 - c1*c2*s3,
1535         s1*s2, c1*s3 + c2*c3*s1, c1*c3 - c2*s1*s3 );
1536     }
1537
1538 template <class A,class Angle>
1539 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1540 typename enable_if_c<
1541     is_mat<A>::value &&
1542     mat_traits<A>::rows==mat_traits<A>::cols &&
1543     mat_traits<A>::rows>=3,
1544     void>::type
1545 set_rot_xzx( A & a, Angle x1, Angle z2, Angle x3 )
1546     {
1547     assign(a,rot_mat_xzx<mat_traits<A>::rows>(x1,z2,x3));
1548     }
1549
1550 template <class A,class Angle>
1551 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1552 typename enable_if_c<
1553     is_mat<A>::value &&
1554     mat_traits<A>::rows==mat_traits<A>::cols &&
1555     mat_traits<A>::rows>=3,
1556     void>::type
1557 rotate_xzx( A & a, Angle x1, Angle z2, Angle x3 )
1558     {
1559     a *= rot_mat_xzx<mat_traits<A>::rows>(x1,z2,x3);
1560     }
1561
1562 ////////////////////////////////////////////////
1563
1564 template <int Dim,class Angle>
1565 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1566 qvm_detail::rot_mat_<Dim,Angle>
1567 rot_mat_xyx( Angle x1, Angle y2, Angle x3 )
1568     {
1569     typedef Angle scalar_type;
1570     scalar_type const c1 = cos(x1);
1571     scalar_type const s1 = sin(x1);
1572     scalar_type const c2 = cos(y2);
1573     scalar_type const s2 = sin(y2);
1574     scalar_type const c3 = cos(x3);
1575     scalar_type const s3 = sin(x3);
1576     return qvm_detail::rot_mat_<Dim,Angle>(
1577         c2, s2*s3, c3*s2,
1578         s1*s2, c1*c3 - c2*s1*s3, -c1*s3 - c2*c3*s1,
1579         -c1*s2, c3*s1 + c1*c2*s3, c1*c2*c3 - s1*s3 );
1580     }
1581
1582 template <class A,class Angle>
1583 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1584 typename enable_if_c<
1585     is_mat<A>::value &&
1586     mat_traits<A>::rows==mat_traits<A>::cols &&
1587     mat_traits<A>::rows>=3,
1588     void>::type
1589 set_rot_xyx( A & a, Angle x1, Angle y2, Angle x3 )
1590     {
1591     assign(a,rot_mat_xyx<mat_traits<A>::rows>(x1,y2,x3));
1592     }
1593
1594 template <class A,class Angle>
1595 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1596 typename enable_if_c<
1597     is_mat<A>::value &&
1598     mat_traits<A>::rows==mat_traits<A>::cols &&
1599     mat_traits<A>::rows>=3,
1600     void>::type
1601 rotate_xyx( A & a, Angle x1, Angle y2, Angle x3 )
1602     {
1603     a *= rot_mat_xyx<mat_traits<A>::rows>(x1,y2,x3);
1604     }
1605
1606 ////////////////////////////////////////////////
1607
1608 template <int Dim,class Angle>
1609 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1610 qvm_detail::rot_mat_<Dim,Angle>
1611 rot_mat_yxy( Angle y1, Angle x2, Angle y3 )
1612     {
1613     typedef Angle scalar_type;
1614     scalar_type const c1 = cos(y1);
1615     scalar_type const s1 = sin(y1);
1616     scalar_type const c2 = cos(x2);
1617     scalar_type const s2 = sin(x2);
1618     scalar_type const c3 = cos(y3);
1619     scalar_type const s3 = sin(y3);
1620     return qvm_detail::rot_mat_<Dim,Angle>(
1621         c1*c3 - c2*s1*s3, s1*s2, c1*s3 + c2*c3*s1,
1622         s2*s3, c2, -c3*s2,
1623         -c3*s1 - c1*c2*s3, c1*s2, c1*c2*c3 - s1*s3 );
1624     }
1625
1626 template <class A,class Angle>
1627 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1628 typename enable_if_c<
1629     is_mat<A>::value &&
1630     mat_traits<A>::rows==mat_traits<A>::cols &&
1631     mat_traits<A>::rows>=3,
1632     void>::type
1633 set_rot_yxy( A & a, Angle y1, Angle x2, Angle y3 )
1634     {
1635     assign(a,rot_mat_yxy<mat_traits<A>::rows>(y1,x2,y3));
1636     }
1637
1638 template <class A,class Angle>
1639 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1640 typename enable_if_c<
1641     is_mat<A>::value &&
1642     mat_traits<A>::rows==mat_traits<A>::cols &&
1643     mat_traits<A>::rows>=3,
1644     void>::type
1645 rotate_yxy( A & a, Angle y1, Angle x2, Angle y3 )
1646     {
1647     a *= rot_mat_yxy<mat_traits<A>::rows>(y1,x2,y3);
1648     }
1649
1650 ////////////////////////////////////////////////
1651
1652 template <int Dim,class Angle>
1653 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1654 qvm_detail::rot_mat_<Dim,Angle>
1655 rot_mat_yzy( Angle y1, Angle z2, Angle y3 )
1656     {
1657     typedef Angle scalar_type;
1658     scalar_type const c1 = cos(y1);
1659     scalar_type const s1 = sin(y1);
1660     scalar_type const c2 = cos(z2);
1661     scalar_type const s2 = sin(z2);
1662     scalar_type const c3 = cos(y3);
1663     scalar_type const s3 = sin(y3);
1664     return qvm_detail::rot_mat_<Dim,Angle>(
1665         c1*c2*c3 - s1*s3, -c1*s2, c3*s1 + c1*c2*s3,
1666         c3*s2, c2, s2*s3,
1667         -c1*s3 - c2*c3*s1, s1*s2, c1*c3 - c2*s1*s3 );
1668     }
1669
1670 template <class A,class Angle>
1671 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1672 typename enable_if_c<
1673     is_mat<A>::value &&
1674     mat_traits<A>::rows==mat_traits<A>::cols &&
1675     mat_traits<A>::rows>=3,
1676     void>::type
1677 set_rot_yzy( A & a, Angle y1, Angle z2, Angle y3 )
1678     {
1679     assign(a,rot_mat_yzy<mat_traits<A>::rows>(y1,z2,y3));
1680     }
1681
1682 template <class A,class Angle>
1683 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1684 typename enable_if_c<
1685     is_mat<A>::value &&
1686     mat_traits<A>::rows==mat_traits<A>::cols &&
1687     mat_traits<A>::rows>=3,
1688     void>::type
1689 rotate_yzy( A & a, Angle y1, Angle z2, Angle y3 )
1690     {
1691     a *= rot_mat_yzy<mat_traits<A>::rows>(y1,z2,y3);
1692     }
1693
1694 ////////////////////////////////////////////////
1695
1696 template <int Dim,class Angle>
1697 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1698 qvm_detail::rot_mat_<Dim,Angle>
1699 rot_mat_zyz( Angle z1, Angle y2, Angle z3 )
1700     {
1701     typedef Angle scalar_type;
1702     scalar_type const c1 = cos(z1);
1703     scalar_type const s1 = sin(z1);
1704     scalar_type const c2 = cos(y2);
1705     scalar_type const s2 = sin(y2);
1706     scalar_type const c3 = cos(z3);
1707     scalar_type const s3 = sin(z3);
1708     return qvm_detail::rot_mat_<Dim,Angle>(
1709         c1*c2*c3 - s1*s3, -c3*s1 - c1*c2*s3, c1*s2,
1710         c1*s3 + c2*c3*s1, c1*c3 - c2*s1*s3, s1*s2,
1711         -c3*s2, s2*s3, c2 );
1712     }
1713
1714 template <class A,class Angle>
1715 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1716 typename enable_if_c<
1717     is_mat<A>::value &&
1718     mat_traits<A>::rows==mat_traits<A>::cols &&
1719     mat_traits<A>::rows>=3,
1720     void>::type
1721 set_rot_zyz( A & a, Angle z1, Angle y2, Angle z3 )
1722     {
1723     assign(a,rot_mat_zyz<mat_traits<A>::rows>(z1,y2,z3));
1724     }
1725
1726 template <class A,class Angle>
1727 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1728 typename enable_if_c<
1729     is_mat<A>::value &&
1730     mat_traits<A>::rows==mat_traits<A>::cols &&
1731     mat_traits<A>::rows>=3,
1732     void>::type
1733 rotate_zyz( A & a, Angle z1, Angle y2, Angle z3 )
1734     {
1735     a *= rot_mat_zyz<mat_traits<A>::rows>(z1,y2,z3);
1736     }
1737
1738 ////////////////////////////////////////////////
1739
1740 template <int Dim,class Angle>
1741 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE
1742 qvm_detail::rot_mat_<Dim,Angle>
1743 rot_mat_zxz( Angle z1, Angle x2, Angle z3 )
1744     {
1745     typedef Angle scalar_type;
1746     scalar_type const c1 = cos(z1);
1747     scalar_type const s1 = sin(z1);
1748     scalar_type const c2 = cos(x2);
1749     scalar_type const s2 = sin(x2);
1750     scalar_type const c3 = cos(z3);
1751     scalar_type const s3 = sin(z3);
1752     return qvm_detail::rot_mat_<Dim,Angle>(
1753         c1*c3 - c2*s1*s3, -c1*s3 - c2*c3*s1, s1*s2,
1754         c3*s1 + c1*c2*s3, c1*c2*c3 - s1*s3, -c1*s2,
1755         s2*s3, c3*s2, c2 );
1756     }
1757
1758 template <class A,class Angle>
1759 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1760 typename enable_if_c<
1761     is_mat<A>::value &&
1762     mat_traits<A>::rows==mat_traits<A>::cols &&
1763     mat_traits<A>::rows>=3,
1764     void>::type
1765 set_rot_zxz( A & a, Angle z1, Angle x2, Angle z3 )
1766     {
1767     assign(a,rot_mat_zxz<mat_traits<A>::rows>(z1,x2,z3));
1768     }
1769
1770 template <class A,class Angle>
1771 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1772 typename enable_if_c<
1773     is_mat<A>::value &&
1774     mat_traits<A>::rows==mat_traits<A>::cols &&
1775     mat_traits<A>::rows>=3,
1776     void>::type
1777 rotate_zxz( A & a, Angle z1, Angle x2, Angle z3 )
1778     {
1779     a *= rot_mat_zxz<mat_traits<A>::rows>(z1,x2,z3);
1780     }
1781
1782 ////////////////////////////////////////////////
1783
1784 namespace
1785 qvm_detail
1786     {
1787     template <int Dim,class Angle>
1788     struct
1789     rotx_mat_
1790         {
1791         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1792         rotx_mat_()
1793             {
1794             }
1795
1796         template <class R>
1797         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1798         operator R() const
1799             {
1800             R r;
1801             assign(r,*this);
1802             return r;
1803             }
1804
1805         private:
1806
1807         rotx_mat_( rotx_mat_ const & );
1808         rotx_mat_ & operator=( rotx_mat_ const & );
1809         ~rotx_mat_();
1810         };
1811
1812     template <int Row,int Col>
1813     struct
1814     rotx_m_get
1815         {
1816         template <class T>
1817         static
1818         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1819         T
1820         get( T const & )
1821             {
1822             return scalar_traits<T>::value(Row==Col);
1823             }
1824         };
1825
1826     template <>
1827     struct
1828     rotx_m_get<1,1>
1829         {
1830         template <class T>
1831         static
1832         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1833         T
1834         get( T const & angle )
1835             {
1836             return cos(angle);
1837             }
1838         };
1839
1840     template <>
1841     struct
1842     rotx_m_get<1,2>
1843         {
1844         template <class T>
1845         static
1846         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1847         T
1848         get( T const & angle )
1849             {
1850             return -sin(angle);
1851             }
1852         };
1853
1854     template <>
1855     struct
1856     rotx_m_get<2,1>
1857         {
1858         template <class T>
1859         static
1860         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1861         T
1862         get( T const & angle )
1863             {
1864             return sin(angle);
1865             }
1866         };
1867
1868     template <>
1869     struct
1870     rotx_m_get<2,2>
1871         {
1872         template <class T>
1873         static
1874         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1875         T
1876         get( T const & angle )
1877             {
1878             return cos(angle);
1879             }
1880         };
1881     }
1882
1883 template <int Dim,class Angle>
1884 struct
1885 mat_traits< qvm_detail::rotx_mat_<Dim,Angle> >
1886     {
1887     typedef qvm_detail::rotx_mat_<Dim,Angle> this_matrix;
1888     typedef Angle scalar_type;
1889     static int const rows=Dim;
1890     static int const cols=Dim;
1891
1892     template <int Row,int Col>
1893     static
1894     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1895     scalar_type
1896     read_element( this_matrix const & x )
1897         {
1898         BOOST_QVM_STATIC_ASSERT(Row>=0);
1899         BOOST_QVM_STATIC_ASSERT(Col>=0);
1900         BOOST_QVM_STATIC_ASSERT(Row<Dim);
1901         BOOST_QVM_STATIC_ASSERT(Col<Dim);
1902         return qvm_detail::rotx_m_get<Row,Col>::get(reinterpret_cast<Angle const &>(x));
1903         }
1904
1905     static
1906     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
1907     scalar_type
1908     read_element_idx( int row, int col, this_matrix const & x )
1909         {
1910         BOOST_QVM_ASSERT(row>=0);
1911         BOOST_QVM_ASSERT(col>=0);
1912         BOOST_QVM_ASSERT(row<Dim);
1913         BOOST_QVM_ASSERT(col<Dim);
1914         Angle const & a=reinterpret_cast<Angle const &>(x);
1915         if( row==1 )
1916             {
1917             if( col==1 )
1918                 return cos(a);
1919             if( col==2 )
1920                 return -sin(a);
1921             }
1922         if( row==2 )
1923             {
1924             if( col==1 )
1925                 return sin(a);
1926             if( col==2 )
1927                 return cos(a);
1928             }
1929         return scalar_traits<scalar_type>::value(row==col);
1930         }
1931     };
1932
1933 template <int Dim,class Angle>
1934 struct
1935 deduce_mat<qvm_detail::rotx_mat_<Dim,Angle>,Dim,Dim>
1936     {
1937     typedef mat<Angle,Dim,Dim> type;
1938     };
1939
1940 template <int Dim,class Angle>
1941 struct
1942 deduce_mat2<qvm_detail::rotx_mat_<Dim,Angle>,qvm_detail::rotx_mat_<Dim,Angle>,Dim,Dim>
1943     {
1944     typedef mat<Angle,Dim,Dim> type;
1945     };
1946
1947 template <int Dim,class Angle>
1948 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1949 qvm_detail::rotx_mat_<Dim,Angle> const &
1950 rotx_mat( Angle const & angle )
1951     {
1952     BOOST_QVM_STATIC_ASSERT(Dim>=3);
1953     return reinterpret_cast<qvm_detail::rotx_mat_<Dim,Angle> const &>(angle);
1954     }
1955
1956 template <class A,class Angle>
1957 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1958 typename enable_if_c<
1959     is_mat<A>::value &&
1960     mat_traits<A>::rows>=3 &&
1961     mat_traits<A>::rows==mat_traits<A>::cols,
1962     void>::type
1963 set_rotx( A & a, Angle angle )
1964     {
1965     assign(a,rotx_mat<mat_traits<A>::rows>(angle));
1966     }
1967
1968 template <class A,class Angle>
1969 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
1970 typename enable_if_c<
1971     is_mat<A>::value &&
1972     mat_traits<A>::rows>=3 &&
1973     mat_traits<A>::rows==mat_traits<A>::cols,
1974     void>::type
1975 rotate_x( A & a, Angle angle )
1976     {
1977     a *= rotx_mat<mat_traits<A>::rows>(angle);
1978     }
1979
1980 ////////////////////////////////////////////////
1981
1982 namespace
1983 qvm_detail
1984     {
1985     template <int Dim,class Angle>
1986     struct
1987     roty_mat_
1988         {
1989         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1990         roty_mat_()
1991             {
1992             }
1993
1994         template <class R>
1995         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
1996         operator R() const
1997             {
1998             R r;
1999             assign(r,*this);
2000             return r;
2001             }
2002
2003         private:
2004
2005         roty_mat_( roty_mat_ const & );
2006         roty_mat_ & operator=( roty_mat_ const & );
2007         ~roty_mat_();
2008         };
2009
2010     template <int Row,int Col>
2011     struct
2012     roty_m_get
2013         {
2014         template <class T>
2015         static
2016         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2017         T
2018         get( T const & )
2019             {
2020             return scalar_traits<T>::value(Row==Col);
2021             }
2022         };
2023
2024     template <>
2025     struct
2026     roty_m_get<0,0>
2027         {
2028         template <class T>
2029         static
2030         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2031         T
2032         get( T const & angle )
2033             {
2034             return cos(angle);
2035             }
2036         };
2037
2038     template <>
2039     struct
2040     roty_m_get<0,2>
2041         {
2042         template <class T>
2043         static
2044         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2045         T
2046         get( T const & angle )
2047             {
2048             return sin(angle);
2049             }
2050         };
2051
2052     template <>
2053     struct
2054     roty_m_get<2,0>
2055         {
2056         template <class T>
2057         static
2058         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2059         T
2060         get( T const & angle )
2061             {
2062             return -sin(angle);
2063             }
2064         };
2065
2066     template <>
2067     struct
2068     roty_m_get<2,2>
2069         {
2070         template <class T>
2071         static
2072         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2073         T
2074         get( T const & angle )
2075             {
2076             return cos(angle);
2077             }
2078         };
2079     }
2080
2081 template <int Dim,class Angle>
2082 struct
2083 mat_traits< qvm_detail::roty_mat_<Dim,Angle> >
2084     {
2085     typedef qvm_detail::roty_mat_<Dim,Angle> this_matrix;
2086     typedef Angle scalar_type;
2087     static int const rows=Dim;
2088     static int const cols=Dim;
2089
2090     template <int Row,int Col>
2091     static
2092     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2093     scalar_type
2094     read_element( this_matrix const & x )
2095         {
2096         BOOST_QVM_STATIC_ASSERT(Row>=0);
2097         BOOST_QVM_STATIC_ASSERT(Col>=0);
2098         BOOST_QVM_STATIC_ASSERT(Row<Dim);
2099         BOOST_QVM_STATIC_ASSERT(Col<Dim);
2100         return qvm_detail::roty_m_get<Row,Col>::get(reinterpret_cast<Angle const &>(x));
2101         }
2102
2103     static
2104     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2105     scalar_type
2106     read_element_idx( int row, int col, this_matrix const & x )
2107         {
2108         BOOST_QVM_ASSERT(row>=0);
2109         BOOST_QVM_ASSERT(col>=0);
2110         BOOST_QVM_ASSERT(row<Dim);
2111         BOOST_QVM_ASSERT(col<Dim);
2112         Angle const & a=reinterpret_cast<Angle const &>(x);
2113         if( row==0 )
2114             {
2115             if( col==0 )
2116                 return cos(a);
2117             if( col==2 )
2118                 return sin(a);
2119             }
2120         if( row==2 )
2121             {
2122             if( col==0 )
2123                 return -sin(a);
2124             if( col==2 )
2125                 return cos(a);
2126             }
2127         return scalar_traits<scalar_type>::value(row==col);
2128         }
2129     };
2130
2131 template <int Dim,class Angle>
2132 struct
2133 deduce_mat<qvm_detail::roty_mat_<Dim,Angle>,Dim,Dim>
2134     {
2135     typedef mat<Angle,Dim,Dim> type;
2136     };
2137
2138 template <int Dim,class Angle>
2139 struct
2140 deduce_mat2<qvm_detail::roty_mat_<Dim,Angle>,qvm_detail::roty_mat_<Dim,Angle>,Dim,Dim>
2141     {
2142     typedef mat<Angle,Dim,Dim> type;
2143     };
2144
2145 template <int Dim,class Angle>
2146 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2147 qvm_detail::roty_mat_<Dim,Angle> const &
2148 roty_mat( Angle const & angle )
2149     {
2150     BOOST_QVM_STATIC_ASSERT(Dim>=3);
2151     return reinterpret_cast<qvm_detail::roty_mat_<Dim,Angle> const &>(angle);
2152     }
2153
2154 template <class A,class Angle>
2155 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2156 typename enable_if_c<
2157     is_mat<A>::value &&
2158     mat_traits<A>::rows>=2 &&
2159     mat_traits<A>::rows==mat_traits<A>::cols,
2160     void>::type
2161 set_roty( A & a, Angle angle )
2162     {
2163     assign(a,roty_mat<mat_traits<A>::rows>(angle));
2164     }
2165
2166 template <class A,class Angle>
2167 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2168 typename enable_if_c<
2169     is_mat<A>::value &&
2170     mat_traits<A>::rows>=3 &&
2171     mat_traits<A>::rows==mat_traits<A>::cols,
2172     void>::type
2173 rotate_y( A & a, Angle angle )
2174     {
2175     a *= roty_mat<mat_traits<A>::rows>(angle);
2176     }
2177
2178 ////////////////////////////////////////////////
2179
2180 namespace
2181 qvm_detail
2182     {
2183     template <int Dim,class Angle>
2184     struct
2185     rotz_mat_
2186         {
2187         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2188         rotz_mat_()
2189             {
2190             }
2191
2192         template <class R>
2193         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2194         operator R() const
2195             {
2196             R r;
2197             assign(r,*this);
2198             return r;
2199             }
2200
2201         private:
2202
2203         rotz_mat_( rotz_mat_ const & );
2204         rotz_mat_ & operator=( rotz_mat_ const & );
2205         ~rotz_mat_();
2206         };
2207
2208     template <int Row,int Col>
2209     struct
2210     rotz_m_get
2211         {
2212         template <class T>
2213         static
2214         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2215         T
2216         get( T const & )
2217             {
2218             return scalar_traits<T>::value(Row==Col);
2219             }
2220         };
2221
2222     template <>
2223     struct
2224     rotz_m_get<0,0>
2225         {
2226         template <class T>
2227         static
2228         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2229         T
2230         get( T const & angle )
2231             {
2232             return cos(angle);
2233             }
2234         };
2235
2236     template <>
2237     struct
2238     rotz_m_get<0,1>
2239         {
2240         template <class T>
2241         static
2242         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2243         T
2244         get( T const & angle )
2245             {
2246             return -sin(angle);
2247             }
2248         };
2249
2250     template <>
2251     struct
2252     rotz_m_get<1,0>
2253         {
2254         template <class T>
2255         static
2256         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2257         T
2258         get( T const & angle )
2259             {
2260             return sin(angle);
2261             }
2262         };
2263
2264     template <>
2265     struct
2266     rotz_m_get<1,1>
2267         {
2268         template <class T>
2269         static
2270         BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2271         T
2272         get( T const & angle )
2273             {
2274             return cos(angle);
2275             }
2276         };
2277     }
2278
2279 template <int Dim,class Angle>
2280 struct
2281 mat_traits< qvm_detail::rotz_mat_<Dim,Angle> >
2282     {
2283     typedef qvm_detail::rotz_mat_<Dim,Angle> this_matrix;
2284     typedef Angle scalar_type;
2285     static int const rows=Dim;
2286     static int const cols=Dim;
2287
2288     template <int Row,int Col>
2289     static
2290     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2291     scalar_type
2292     read_element( this_matrix const & x )
2293         {
2294         BOOST_QVM_STATIC_ASSERT(Row>=0);
2295         BOOST_QVM_STATIC_ASSERT(Col>=0);
2296         BOOST_QVM_STATIC_ASSERT(Row<Dim);
2297         BOOST_QVM_STATIC_ASSERT(Col<Dim);
2298         return qvm_detail::rotz_m_get<Row,Col>::get(reinterpret_cast<Angle const &>(x));
2299         }
2300
2301     static
2302     BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_CRITICAL
2303     scalar_type
2304     read_element_idx( int row, int col, this_matrix const & x )
2305         {
2306         BOOST_QVM_ASSERT(row>=0);
2307         BOOST_QVM_ASSERT(col>=0);
2308         BOOST_QVM_ASSERT(row<Dim);
2309         BOOST_QVM_ASSERT(col<Dim);
2310         Angle const & a=reinterpret_cast<Angle const &>(x);
2311         if( row==0 )
2312             {
2313             if( col==0 )
2314                 return cos(a);
2315             if( col==1 )
2316                 return -sin(a);
2317             }
2318         if( row==1 )
2319             {
2320             if( col==0 )
2321                 return sin(a);
2322             if( col==1 )
2323                 return cos(a);
2324             }
2325         return scalar_traits<scalar_type>::value(row==col);
2326         }
2327     };
2328
2329 template <int Dim,class Angle>
2330 struct
2331 deduce_mat<qvm_detail::rotz_mat_<Dim,Angle>,Dim,Dim>
2332     {
2333     typedef mat<Angle,Dim,Dim> type;
2334     };
2335
2336 template <int Dim,class Angle,int R,int C>
2337 struct
2338 deduce_mat2<qvm_detail::rotz_mat_<Dim,Angle>,qvm_detail::rotz_mat_<Dim,Angle>,R,C>
2339     {
2340     typedef mat<Angle,R,C> type;
2341     };
2342
2343 template <int Dim,class Angle>
2344 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2345 qvm_detail::rotz_mat_<Dim,Angle> const &
2346 rotz_mat( Angle const & angle )
2347     {
2348     BOOST_QVM_STATIC_ASSERT(Dim>=2);
2349     return reinterpret_cast<qvm_detail::rotz_mat_<Dim,Angle> const &>(angle);
2350     }
2351
2352 template <class A,class Angle>
2353 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2354 typename enable_if_c<
2355     is_mat<A>::value &&
2356     mat_traits<A>::rows>=2 &&
2357     mat_traits<A>::rows==mat_traits<A>::cols,
2358     void>::type
2359 set_rotz( A & a, Angle angle )
2360     {
2361     assign(a,rotz_mat<mat_traits<A>::rows>(angle));
2362     }
2363
2364 template <class A,class Angle>
2365 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_OPERATIONS
2366 typename enable_if_c<
2367     is_mat<A>::value &&
2368     mat_traits<A>::rows>=2 &&
2369     mat_traits<A>::rows==mat_traits<A>::cols,
2370     void>::type
2371 rotate_z( A & a, Angle angle )
2372     {
2373     a *= rotz_mat<mat_traits<A>::rows>(angle);
2374     }
2375
2376 ////////////////////////////////////////////////
2377
2378 namespace
2379 qvm_detail
2380     {
2381     template <int D>
2382     struct
2383     inverse_m_defined
2384         {
2385         static bool const value=false;
2386         };
2387     }
2388
2389 template <class A,class B>
2390 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2391 typename lazy_enable_if_c<
2392     is_mat<A>::value && is_scalar<B>::value &&
2393     mat_traits<A>::rows==mat_traits<A>::cols &&
2394     !qvm_detail::inverse_m_defined<mat_traits<A>::rows>::value,
2395     deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols> >::type
2396 inverse( A const & a, B det )
2397     {
2398     typedef typename mat_traits<A>::scalar_type T;
2399     BOOST_QVM_ASSERT(det!=scalar_traits<T>::value(0));
2400     T f=scalar_traits<T>::value(1)/det;
2401     typedef typename deduce_mat2<A,B,mat_traits<A>::rows,mat_traits<A>::cols>::type cofactor_return_type;
2402     cofactor_return_type c=qvm_detail::cofactor_impl(a);
2403     return reinterpret_cast<qvm_detail::transposed_<cofactor_return_type> const &>(c) * f;
2404     }
2405
2406 template <class A>
2407 BOOST_QVM_CONSTEXPR BOOST_QVM_INLINE_TRIVIAL
2408 typename lazy_enable_if_c<
2409     is_mat<A>::value &&
2410     mat_traits<A>::rows==mat_traits<A>::cols &&
2411     !qvm_detail::inverse_m_defined<mat_traits<A>::rows>::value,
2412     deduce_mat<A> >::type
2413 inverse( A const & a )
2414     {
2415     typedef typename mat_traits<A>::scalar_type T;
2416     T det=determinant(a);
2417     if( det==scalar_traits<T>::value(0) )
2418         BOOST_QVM_THROW_EXCEPTION(zero_determinant_error());
2419     return inverse(a,det);
2420     }
2421
2422 ////////////////////////////////////////////////
2423
2424 namespace
2425 sfinae
2426     {
2427     using ::boost::qvm::to_string;
2428     using ::boost::qvm::assign;
2429     using ::boost::qvm::determinant;
2430     using ::boost::qvm::cmp;
2431     using ::boost::qvm::convert_to;
2432     using ::boost::qvm::set_identity;
2433     using ::boost::qvm::set_zero;
2434     using ::boost::qvm::scalar_cast;
2435     using ::boost::qvm::operator/=;
2436     using ::boost::qvm::operator/;
2437     using ::boost::qvm::operator==;
2438     using ::boost::qvm::operator-=;
2439     using ::boost::qvm::operator-;
2440     using ::boost::qvm::operator*=;
2441     using ::boost::qvm::operator*;
2442     using ::boost::qvm::operator!=;
2443     using ::boost::qvm::operator+=;
2444     using ::boost::qvm::operator+;
2445     using ::boost::qvm::mref;
2446     using ::boost::qvm::rot_mat;
2447     using ::boost::qvm::set_rot;
2448     using ::boost::qvm::rotate;
2449     using ::boost::qvm::set_rotx;
2450     using ::boost::qvm::rotate_x;
2451     using ::boost::qvm::set_roty;
2452     using ::boost::qvm::rotate_y;
2453     using ::boost::qvm::set_rotz;
2454     using ::boost::qvm::rotate_z;
2455     using ::boost::qvm::inverse;
2456     }
2457
2458 } }
2459
2460 #endif