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_3DCF6B90AE0E11DE9A315BE555D89593
7 #define UUID_3DCF6B90AE0E11DE9A315BE555D89593
8 
9 #include <boost/qvm/inline.hpp>
10 #include <boost/qvm/mat_traits_array.hpp>
11 #include <boost/qvm/static_assert.hpp>
12 
13 namespace
14 boost
15     {
16     namespace
17     qvm
18         {
19         namespace
20         qvm_detail
21             {
22             template <int N>
23             struct
24             det_size
25                 {
26                 };
27 
28             template <class M>
29             BOOST_QVM_INLINE_TRIVIAL
30             typename mat_traits<M>::scalar_type
determinant_impl_(M const & a,det_size<2>)31             determinant_impl_( M const & a, det_size<2> )
32                 {
33                 return
34                     mat_traits<M>::template read_element<0,0>(a) * mat_traits<M>::template read_element<1,1>(a) -
35                     mat_traits<M>::template read_element<1,0>(a) * mat_traits<M>::template read_element<0,1>(a);
36                 }
37 
38             template <class M,int N>
39             BOOST_QVM_INLINE_RECURSION
40             typename mat_traits<M>::scalar_type
determinant_impl_(M const & a,det_size<N>)41             determinant_impl_( M const & a, det_size<N> )
42                 {
43                 typedef typename mat_traits<M>::scalar_type T;
44                 T m[N-1][N-1];
45                 T det=T(0);
46                 for( int j1=0; j1!=N; ++j1 )
47                     {
48                     for( int i=1; i!=N; ++i )
49                         {
50                         int j2 = 0;
51                         for( int j=0; j!=N; ++j )
52                             {
53                             if( j==j1 )
54                                 continue;
55                             m[i-1][j2] = mat_traits<M>::read_element_idx(i,j,a);
56                             ++j2;
57                             }
58                         }
59                     T d=determinant_impl_(m,det_size<N-1>());
60                     if( j1&1 )
61                         d=-d;
62                     det += mat_traits<M>::read_element_idx(0,j1,a) * d;
63                     }
64                 return det;
65                 }
66 
67             template <class M>
68             BOOST_QVM_INLINE_TRIVIAL
69             typename mat_traits<M>::scalar_type
determinant_impl(M const & a)70             determinant_impl( M const & a )
71                 {
72                 BOOST_QVM_STATIC_ASSERT(mat_traits<M>::rows==mat_traits<M>::cols);
73                 return determinant_impl_(a,det_size<mat_traits<M>::rows>());
74                 }
75             }
76         }
77     }
78 
79 #endif
80