ceph/src/boost/boost/qvm/detail/determinant_impl.hpp

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