1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
///
/// Initial setup of sub-elements, starting with last row...
///
template<type T, int row> void det_subelts_setup(const T** elts, const T** subelts)
{
det_subelts_setup<T,row-1>(elts,subelts);
subelts[row] = elts[row]+1;
}
///
/// Specialisation when row index is 0 (no call to previous row)
///
template<type T, 0> void det_subelts_setup(const T** elts, const T** subelts)
{
subelts[0] = elts[0] +1;
}
///
/// Code factorisation to retrieve a sub-matrix * element value
///
template<type T, int n, int row> T det_part(const T** elts, const T** subelts)
{
return (row%2) ?
-elts[row][0] * det<T,n-1>(subelts) :
+elts[row][0] * det<T,n-1>(subelts);
}
///
/// Determinant accumlation for any row (subelts must be already setup
/// for that row).
///
template<type T, int n, int row> T det_acc(const T** elts, const T** subelts)
{
// subelts is already setup to handle row !
T ret = det_part<T,n,row>(elts,subelts);
// replacing previous row
subelts[row-1] = elts[row]+1;
// accumulating (loop)
ret += det_acc<T,n,row-1>(elts,subelts);
return ret;
}
///
/// Specialisation for row 0 (no call to previous row)
///
template<type T, int n, 0> T det_acc(const T** elts, const T** subelts)
{
return det_part<T,n,row>(elts,subelts);
}
///
/// Determinant function:
/// - Prepare sub-rows
/// - Start accumulation
template<type T, int n> det(const T** elts)
{
const T* subelts[n-1];
det_subelts_setup<T,n-2>(elts,subelts);
// making the loop
return det_acc<T,n,n-1>(elts,subelts);
}
/// specialization for n==1
template<type T,1> T det(const T** elts)
{
return elts[0][0];
}
/// specialization for n==2 (optionel)
template<type T,2> T det(const T** elts)
{
return elts[0][0]*elts[1][1] - elts[0][1]*elts[1][0];
} |
