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

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