determinant_impl.hpp 2.5 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879
  1. //Copyright (c) 2008-2016 Emil Dotchevski and Reverge Studios, Inc.
  2. //Distributed under the Boost Software License, Version 1.0. (See accompanying
  3. //file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  4. #ifndef UUID_3DCF6B90AE0E11DE9A315BE555D89593
  5. #define UUID_3DCF6B90AE0E11DE9A315BE555D89593
  6. #include <boost/qvm/inline.hpp>
  7. #include <boost/qvm/mat_traits_array.hpp>
  8. #include <boost/qvm/static_assert.hpp>
  9. namespace
  10. boost
  11. {
  12. namespace
  13. qvm
  14. {
  15. namespace
  16. qvm_detail
  17. {
  18. template <int N>
  19. struct
  20. det_size
  21. {
  22. };
  23. template <class M>
  24. BOOST_QVM_INLINE_TRIVIAL
  25. typename mat_traits<M>::scalar_type
  26. determinant_impl_( M const & a, det_size<2> )
  27. {
  28. return
  29. mat_traits<M>::template read_element<0,0>(a) * mat_traits<M>::template read_element<1,1>(a) -
  30. mat_traits<M>::template read_element<1,0>(a) * mat_traits<M>::template read_element<0,1>(a);
  31. }
  32. template <class M,int N>
  33. BOOST_QVM_INLINE_RECURSION
  34. typename mat_traits<M>::scalar_type
  35. determinant_impl_( M const & a, det_size<N> )
  36. {
  37. typedef typename mat_traits<M>::scalar_type T;
  38. T m[N-1][N-1];
  39. T det=T(0);
  40. for( int j1=0; j1!=N; ++j1 )
  41. {
  42. for( int i=1; i!=N; ++i )
  43. {
  44. int j2 = 0;
  45. for( int j=0; j!=N; ++j )
  46. {
  47. if( j==j1 )
  48. continue;
  49. m[i-1][j2] = mat_traits<M>::read_element_idx(i,j,a);
  50. ++j2;
  51. }
  52. }
  53. T d=determinant_impl_(m,det_size<N-1>());
  54. if( j1&1 )
  55. d=-d;
  56. det += mat_traits<M>::read_element_idx(0,j1,a) * d;
  57. }
  58. return det;
  59. }
  60. template <class M>
  61. BOOST_QVM_INLINE_TRIVIAL
  62. typename mat_traits<M>::scalar_type
  63. determinant_impl( M const & a )
  64. {
  65. BOOST_QVM_STATIC_ASSERT(mat_traits<M>::rows==mat_traits<M>::cols);
  66. return determinant_impl_(a,det_size<mat_traits<M>::rows>());
  67. }
  68. }
  69. }
  70. }
  71. #endif