determinant_impl.hpp 2.1 KB

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