EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CbmLitMatrixMath.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CbmLitMatrixMath.cxx
1 #include "utils/CbmLitMatrixMath.h"
2 
3 #include <iostream>
4 
5 /*
6  * Matrix operations
7  */
8 
9 // SMij are indices for a 5x5 matrix.
10 
11 #define SM00 0
12 #define SM01 1
13 #define SM02 2
14 #define SM03 3
15 #define SM04 4
16 
17 #define SM10 1
18 #define SM11 5
19 #define SM12 6
20 #define SM13 7
21 #define SM14 8
22 
23 #define SM20 2
24 #define SM21 6
25 #define SM22 9
26 #define SM23 10
27 #define SM24 11
28 
29 #define SM30 3
30 #define SM31 7
31 #define SM32 10
32 #define SM33 12
33 #define SM34 13
34 
35 #define SM40 4
36 #define SM41 8
37 #define SM42 11
38 #define SM43 13
39 #define SM44 14
40 
41 bool InvSym15(
42  std::vector<litfloat>& a)
43 {
44  if (a.size() != 15) {
45  std::cout << "-E- InvSym15: size is not correct" << std::endl;
46  return false;
47  }
48 
49  litfloat* pM = &a[0];
50 
51  // Find all NECESSARY 2x2 dets: (25 of them)
52 
53  const litfloat mDet2_23_01 = pM[SM20]*pM[SM31] - pM[SM21]*pM[SM30];
54  const litfloat mDet2_23_02 = pM[SM20]*pM[SM32] - pM[SM22]*pM[SM30];
55  const litfloat mDet2_23_03 = pM[SM20]*pM[SM33] - pM[SM23]*pM[SM30];
56  const litfloat mDet2_23_12 = pM[SM21]*pM[SM32] - pM[SM22]*pM[SM31];
57  const litfloat mDet2_23_13 = pM[SM21]*pM[SM33] - pM[SM23]*pM[SM31];
58  const litfloat mDet2_23_23 = pM[SM22]*pM[SM33] - pM[SM23]*pM[SM32];
59  const litfloat mDet2_24_01 = pM[SM20]*pM[SM41] - pM[SM21]*pM[SM40];
60  const litfloat mDet2_24_02 = pM[SM20]*pM[SM42] - pM[SM22]*pM[SM40];
61  const litfloat mDet2_24_03 = pM[SM20]*pM[SM43] - pM[SM23]*pM[SM40];
62  const litfloat mDet2_24_04 = pM[SM20]*pM[SM44] - pM[SM24]*pM[SM40];
63  const litfloat mDet2_24_12 = pM[SM21]*pM[SM42] - pM[SM22]*pM[SM41];
64  const litfloat mDet2_24_13 = pM[SM21]*pM[SM43] - pM[SM23]*pM[SM41];
65  const litfloat mDet2_24_14 = pM[SM21]*pM[SM44] - pM[SM24]*pM[SM41];
66  const litfloat mDet2_24_23 = pM[SM22]*pM[SM43] - pM[SM23]*pM[SM42];
67  const litfloat mDet2_24_24 = pM[SM22]*pM[SM44] - pM[SM24]*pM[SM42];
68  const litfloat mDet2_34_01 = pM[SM30]*pM[SM41] - pM[SM31]*pM[SM40];
69  const litfloat mDet2_34_02 = pM[SM30]*pM[SM42] - pM[SM32]*pM[SM40];
70  const litfloat mDet2_34_03 = pM[SM30]*pM[SM43] - pM[SM33]*pM[SM40];
71  const litfloat mDet2_34_04 = pM[SM30]*pM[SM44] - pM[SM34]*pM[SM40];
72  const litfloat mDet2_34_12 = pM[SM31]*pM[SM42] - pM[SM32]*pM[SM41];
73  const litfloat mDet2_34_13 = pM[SM31]*pM[SM43] - pM[SM33]*pM[SM41];
74  const litfloat mDet2_34_14 = pM[SM31]*pM[SM44] - pM[SM34]*pM[SM41];
75  const litfloat mDet2_34_23 = pM[SM32]*pM[SM43] - pM[SM33]*pM[SM42];
76  const litfloat mDet2_34_24 = pM[SM32]*pM[SM44] - pM[SM34]*pM[SM42];
77  const litfloat mDet2_34_34 = pM[SM33]*pM[SM44] - pM[SM34]*pM[SM43];
78 
79  // Find all NECESSARY 3x3 dets: (30 of them)
80 
81  const litfloat mDet3_123_012 = pM[SM10]*mDet2_23_12 - pM[SM11]*mDet2_23_02 + pM[SM12]*mDet2_23_01;
82  const litfloat mDet3_123_013 = pM[SM10]*mDet2_23_13 - pM[SM11]*mDet2_23_03 + pM[SM13]*mDet2_23_01;
83  const litfloat mDet3_123_023 = pM[SM10]*mDet2_23_23 - pM[SM12]*mDet2_23_03 + pM[SM13]*mDet2_23_02;
84  const litfloat mDet3_123_123 = pM[SM11]*mDet2_23_23 - pM[SM12]*mDet2_23_13 + pM[SM13]*mDet2_23_12;
85  const litfloat mDet3_124_012 = pM[SM10]*mDet2_24_12 - pM[SM11]*mDet2_24_02 + pM[SM12]*mDet2_24_01;
86  const litfloat mDet3_124_013 = pM[SM10]*mDet2_24_13 - pM[SM11]*mDet2_24_03 + pM[SM13]*mDet2_24_01;
87  const litfloat mDet3_124_014 = pM[SM10]*mDet2_24_14 - pM[SM11]*mDet2_24_04 + pM[SM14]*mDet2_24_01;
88  const litfloat mDet3_124_023 = pM[SM10]*mDet2_24_23 - pM[SM12]*mDet2_24_03 + pM[SM13]*mDet2_24_02;
89  const litfloat mDet3_124_024 = pM[SM10]*mDet2_24_24 - pM[SM12]*mDet2_24_04 + pM[SM14]*mDet2_24_02;
90  const litfloat mDet3_124_123 = pM[SM11]*mDet2_24_23 - pM[SM12]*mDet2_24_13 + pM[SM13]*mDet2_24_12;
91  const litfloat mDet3_124_124 = pM[SM11]*mDet2_24_24 - pM[SM12]*mDet2_24_14 + pM[SM14]*mDet2_24_12;
92  const litfloat mDet3_134_012 = pM[SM10]*mDet2_34_12 - pM[SM11]*mDet2_34_02 + pM[SM12]*mDet2_34_01;
93  const litfloat mDet3_134_013 = pM[SM10]*mDet2_34_13 - pM[SM11]*mDet2_34_03 + pM[SM13]*mDet2_34_01;
94  const litfloat mDet3_134_014 = pM[SM10]*mDet2_34_14 - pM[SM11]*mDet2_34_04 + pM[SM14]*mDet2_34_01;
95  const litfloat mDet3_134_023 = pM[SM10]*mDet2_34_23 - pM[SM12]*mDet2_34_03 + pM[SM13]*mDet2_34_02;
96  const litfloat mDet3_134_024 = pM[SM10]*mDet2_34_24 - pM[SM12]*mDet2_34_04 + pM[SM14]*mDet2_34_02;
97  const litfloat mDet3_134_034 = pM[SM10]*mDet2_34_34 - pM[SM13]*mDet2_34_04 + pM[SM14]*mDet2_34_03;
98  const litfloat mDet3_134_123 = pM[SM11]*mDet2_34_23 - pM[SM12]*mDet2_34_13 + pM[SM13]*mDet2_34_12;
99  const litfloat mDet3_134_124 = pM[SM11]*mDet2_34_24 - pM[SM12]*mDet2_34_14 + pM[SM14]*mDet2_34_12;
100  const litfloat mDet3_134_134 = pM[SM11]*mDet2_34_34 - pM[SM13]*mDet2_34_14 + pM[SM14]*mDet2_34_13;
101  const litfloat mDet3_234_012 = pM[SM20]*mDet2_34_12 - pM[SM21]*mDet2_34_02 + pM[SM22]*mDet2_34_01;
102  const litfloat mDet3_234_013 = pM[SM20]*mDet2_34_13 - pM[SM21]*mDet2_34_03 + pM[SM23]*mDet2_34_01;
103  const litfloat mDet3_234_014 = pM[SM20]*mDet2_34_14 - pM[SM21]*mDet2_34_04 + pM[SM24]*mDet2_34_01;
104  const litfloat mDet3_234_023 = pM[SM20]*mDet2_34_23 - pM[SM22]*mDet2_34_03 + pM[SM23]*mDet2_34_02;
105  const litfloat mDet3_234_024 = pM[SM20]*mDet2_34_24 - pM[SM22]*mDet2_34_04 + pM[SM24]*mDet2_34_02;
106  const litfloat mDet3_234_034 = pM[SM20]*mDet2_34_34 - pM[SM23]*mDet2_34_04 + pM[SM24]*mDet2_34_03;
107  const litfloat mDet3_234_123 = pM[SM21]*mDet2_34_23 - pM[SM22]*mDet2_34_13 + pM[SM23]*mDet2_34_12;
108  const litfloat mDet3_234_124 = pM[SM21]*mDet2_34_24 - pM[SM22]*mDet2_34_14 + pM[SM24]*mDet2_34_12;
109  const litfloat mDet3_234_134 = pM[SM21]*mDet2_34_34 - pM[SM23]*mDet2_34_14 + pM[SM24]*mDet2_34_13;
110  const litfloat mDet3_234_234 = pM[SM22]*mDet2_34_34 - pM[SM23]*mDet2_34_24 + pM[SM24]*mDet2_34_23;
111 
112  // Find all NECESSARY 4x4 dets: (15 of them)
113 
114  const litfloat mDet4_0123_0123 = pM[SM00]*mDet3_123_123 - pM[SM01]*mDet3_123_023
115  + pM[SM02]*mDet3_123_013 - pM[SM03]*mDet3_123_012;
116  const litfloat mDet4_0124_0123 = pM[SM00]*mDet3_124_123 - pM[SM01]*mDet3_124_023
117  + pM[SM02]*mDet3_124_013 - pM[SM03]*mDet3_124_012;
118  const litfloat mDet4_0124_0124 = pM[SM00]*mDet3_124_124 - pM[SM01]*mDet3_124_024
119  + pM[SM02]*mDet3_124_014 - pM[SM04]*mDet3_124_012;
120  const litfloat mDet4_0134_0123 = pM[SM00]*mDet3_134_123 - pM[SM01]*mDet3_134_023
121  + pM[SM02]*mDet3_134_013 - pM[SM03]*mDet3_134_012;
122  const litfloat mDet4_0134_0124 = pM[SM00]*mDet3_134_124 - pM[SM01]*mDet3_134_024
123  + pM[SM02]*mDet3_134_014 - pM[SM04]*mDet3_134_012;
124  const litfloat mDet4_0134_0134 = pM[SM00]*mDet3_134_134 - pM[SM01]*mDet3_134_034
125  + pM[SM03]*mDet3_134_014 - pM[SM04]*mDet3_134_013;
126  const litfloat mDet4_0234_0123 = pM[SM00]*mDet3_234_123 - pM[SM01]*mDet3_234_023
127  + pM[SM02]*mDet3_234_013 - pM[SM03]*mDet3_234_012;
128  const litfloat mDet4_0234_0124 = pM[SM00]*mDet3_234_124 - pM[SM01]*mDet3_234_024
129  + pM[SM02]*mDet3_234_014 - pM[SM04]*mDet3_234_012;
130  const litfloat mDet4_0234_0134 = pM[SM00]*mDet3_234_134 - pM[SM01]*mDet3_234_034
131  + pM[SM03]*mDet3_234_014 - pM[SM04]*mDet3_234_013;
132  const litfloat mDet4_0234_0234 = pM[SM00]*mDet3_234_234 - pM[SM02]*mDet3_234_034
133  + pM[SM03]*mDet3_234_024 - pM[SM04]*mDet3_234_023;
134  const litfloat mDet4_1234_0123 = pM[SM10]*mDet3_234_123 - pM[SM11]*mDet3_234_023
135  + pM[SM12]*mDet3_234_013 - pM[SM13]*mDet3_234_012;
136  const litfloat mDet4_1234_0124 = pM[SM10]*mDet3_234_124 - pM[SM11]*mDet3_234_024
137  + pM[SM12]*mDet3_234_014 - pM[SM14]*mDet3_234_012;
138  const litfloat mDet4_1234_0134 = pM[SM10]*mDet3_234_134 - pM[SM11]*mDet3_234_034
139  + pM[SM13]*mDet3_234_014 - pM[SM14]*mDet3_234_013;
140  const litfloat mDet4_1234_0234 = pM[SM10]*mDet3_234_234 - pM[SM12]*mDet3_234_034
141  + pM[SM13]*mDet3_234_024 - pM[SM14]*mDet3_234_023;
142  const litfloat mDet4_1234_1234 = pM[SM11]*mDet3_234_234 - pM[SM12]*mDet3_234_134
143  + pM[SM13]*mDet3_234_124 - pM[SM14]*mDet3_234_123;
144 
145  // Find the 5x5 det:
146 
147  const litfloat det = pM[SM00]*mDet4_1234_1234 - pM[SM01]*mDet4_1234_0234 + pM[SM02]*mDet4_1234_0134
148  - pM[SM03]*mDet4_1234_0124 + pM[SM04]*mDet4_1234_0123;
149 
150  if (det == 0) {
151  std::cout << "-E- InvSym15: zero determinant" << std::endl;
152  return false;
153  }
154 
155  const litfloat oneOverDet = 1.0/det;
156  const litfloat mn1OverDet = - oneOverDet;
157 
158  pM[SM00] = mDet4_1234_1234 * oneOverDet;
159  pM[SM01] = mDet4_1234_0234 * mn1OverDet;
160  pM[SM02] = mDet4_1234_0134 * oneOverDet;
161  pM[SM03] = mDet4_1234_0124 * mn1OverDet;
162  pM[SM04] = mDet4_1234_0123 * oneOverDet;
163 
164  pM[SM11] = mDet4_0234_0234 * oneOverDet;
165  pM[SM12] = mDet4_0234_0134 * mn1OverDet;
166  pM[SM13] = mDet4_0234_0124 * oneOverDet;
167  pM[SM14] = mDet4_0234_0123 * mn1OverDet;
168 
169  pM[SM22] = mDet4_0134_0134 * oneOverDet;
170  pM[SM23] = mDet4_0134_0124 * mn1OverDet;
171  pM[SM24] = mDet4_0134_0123 * oneOverDet;
172 
173  pM[SM33] = mDet4_0124_0124 * oneOverDet;
174  pM[SM34] = mDet4_0124_0123 * mn1OverDet;
175 
176  pM[SM44] = mDet4_0123_0123 * oneOverDet;
177 
178  return true;
179 }
180 
181 
182 bool Mult25(
183  const std::vector<litfloat>& a,
184  const std::vector<litfloat>& b,
185  std::vector<litfloat>& c)
186 {
187  if (a.size() != 25 || b.size() != 25 || c.size() != 25) {
188  std::cout << "-E- Mult25: size is not correct" << std::endl;
189  return false;
190  }
191 
192  c[0] = a[0] * b[0] + a[1] * b[5] + a[2] * b[10] + a[3] * b[15] + a[4] * b[20];
193  c[1] = a[0] * b[1] + a[1] * b[6] + a[2] * b[11] + a[3] * b[16] + a[4] * b[21];
194  c[2] = a[0] * b[2] + a[1] * b[7] + a[2] * b[12] + a[3] * b[17] + a[4] * b[22];
195  c[3] = a[0] * b[3] + a[1] * b[8] + a[2] * b[13] + a[3] * b[18] + a[4] * b[23];
196  c[4] = a[0] * b[4] + a[1] * b[9] + a[2] * b[14] + a[3] * b[19] + a[4] * b[24];
197  c[5] = a[5] * b[0] + a[6] * b[5] + a[7] * b[10] + a[8] * b[15] + a[9] * b[20];
198  c[6] = a[5] * b[1] + a[6] * b[6] + a[7] * b[11] + a[8] * b[16] + a[9] * b[21];
199  c[7] = a[5] * b[2] + a[6] * b[7] + a[7] * b[12] + a[8] * b[17] + a[9] * b[22];
200  c[8] = a[5] * b[3] + a[6] * b[8] + a[7] * b[13] + a[8] * b[18] + a[9] * b[23];
201  c[9] = a[5] * b[4] + a[6] * b[9] + a[7] * b[14] + a[8] * b[19] + a[9] * b[24];
202  c[10] = a[10] * b[0] + a[11] * b[5] + a[12] * b[10] + a[13] * b[15] + a[14] * b[20];
203  c[11] = a[10] * b[1] + a[11] * b[6] + a[12] * b[11] + a[13] * b[16] + a[14] * b[21];
204  c[12] = a[10] * b[2] + a[11] * b[7] + a[12] * b[12] + a[13] * b[17] + a[14] * b[22];
205  c[13] = a[10] * b[3] + a[11] * b[8] + a[12] * b[13] + a[13] * b[18] + a[14] * b[23];
206  c[14] = a[10] * b[4] + a[11] * b[9] + a[12] * b[14] + a[13] * b[19] + a[14] * b[24];
207  c[15] = a[15] * b[0] + a[16] * b[5] + a[17] * b[10] + a[18] * b[15] + a[19] * b[20];
208  c[16] = a[15] * b[1] + a[16] * b[6] + a[17] * b[11] + a[18] * b[16] + a[19] * b[21];
209  c[17] = a[15] * b[2] + a[16] * b[7] + a[17] * b[12] + a[18] * b[17] + a[19] * b[22];
210  c[18] = a[15] * b[3] + a[16] * b[8] + a[17] * b[13] + a[18] * b[18] + a[19] * b[23];
211  c[19] = a[15] * b[4] + a[16] * b[9] + a[17] * b[14] + a[18] * b[19] + a[19] * b[24];
212  c[20] = a[20] * b[0] + a[21] * b[5] + a[22] * b[10] + a[23] * b[15] + a[24] * b[20];
213  c[21] = a[20] * b[1] + a[21] * b[6] + a[22] * b[11] + a[23] * b[16] + a[24] * b[21];
214  c[22] = a[20] * b[2] + a[21] * b[7] + a[22] * b[12] + a[23] * b[17] + a[24] * b[22];
215  c[23] = a[20] * b[3] + a[21] * b[8] + a[22] * b[13] + a[23] * b[18] + a[24] * b[23];
216  c[24] = a[20] * b[4] + a[21] * b[9] + a[22] * b[14] + a[23] * b[19] + a[24] * b[24];
217 
218  return true;
219 }
220 
221 
222 bool Transpose25(
223  std::vector<litfloat>& a)
224 {
225  if (a.size() != 25) {
226  std::cout << "-E- Transpose25: size is not correct" << std::endl;
227  return true;
228  }
229  std::vector<litfloat> b(a);
230  a[0] = b[0];
231  a[1] = b[5];
232  a[2] = b[10];
233  a[3] = b[15];
234  a[4] = b[20];
235  a[5] = b[1];
236  a[6] = b[6];
237  a[7] = b[11];
238  a[8] = b[16];
239  a[9] = b[21];
240  a[10] = b[2];
241  a[11] = b[7];
242  a[12] = b[12];
243  a[13] = b[17];
244  a[14] = b[22];
245  a[15] = b[3];
246  a[16] = b[8];
247  a[17] = b[13];
248  a[18] = b[18];
249  a[19] = b[23];
250  a[20] = b[4];
251  a[21] = b[9];
252  a[22] = b[14];
253  a[23] = b[19];
254  a[24] = b[24];
255  return true;
256 }
257 
258 
259 bool Mult25On5(
260  const std::vector<litfloat>& a,
261  const std::vector<litfloat>& b,
262  std::vector<litfloat>& c)
263 {
264  if (a.size() != 25 || b.size() != 5 || c.size() != 5) {
265  std::cout << "-E- Mult25On5: size is not correct" << std::endl;
266  return false;
267  }
268  c[0] = a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3] + a[4] * b[4];
269  c[1] = a[5] * b[0] + a[6] * b[1] + a[7] * b[2] + a[8] * b[3] + a[9] * b[4];
270  c[2] = a[10] * b[0] + a[11] * b[1] + a[12] * b[2] + a[13] * b[3] + a[14] * b[4];
271  c[3] = a[15] * b[0] + a[16] * b[1] + a[17] * b[2] + a[18] * b[3] + a[19] * b[4];
272  c[4] = a[20] * b[0] + a[21] * b[1] + a[22] * b[2] + a[23] * b[3] + a[24] * b[4];
273  return true;
274 }
275 
276 bool Mult15On5(
277  const std::vector<litfloat>& a,
278  const std::vector<litfloat>& b,
279  std::vector<litfloat>& c)
280 {
281  if (a.size() != 15 || b.size() != 5 || c.size() != 5) {
282  std::cout << "-E- Mult15On5: size is not correct" << std::endl;
283  return false;
284  }
285  c[0] = a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3] + a[4] * b[4];
286  c[1] = a[1] * b[0] + a[5] * b[1] + a[6] * b[2] + a[7] * b[3] + a[8] * b[4];
287  c[2] = a[2] * b[0] + a[6] * b[1] + a[9] * b[2] + a[10] * b[3] + a[11] * b[4];
288  c[3] = a[3] * b[0] + a[7] * b[1] + a[10] * b[2] + a[12] * b[3] + a[13] * b[4];
289  c[4] = a[4] * b[0] + a[8] * b[1] + a[11] * b[2] + a[13] * b[3] + a[14] * b[4];
290  return true;
291 }
292 
293 bool Subtract(
294  const std::vector<litfloat>& a,
295  const std::vector<litfloat>& b,
296  std::vector<litfloat>& c)
297 {
298  if (a.size() != b.size() || a.size() != c.size()) {
299  std::cout << "-E- Subtract: size is not correct" << std::endl;
300  return false;
301  }
302  for (unsigned int i = 0; i < a.size(); ++i) { c[i] = a[i] - b[i]; }
303  return true;
304 }
305 
306 
307 
308 bool Add(
309  const std::vector<litfloat>& a,
310  const std::vector<litfloat>& b,
311  std::vector<litfloat>& c)
312 {
313  if (a.size() != b.size() || a.size() != c.size()) {
314  std::cout << "-E- Add: size is not correct" << std::endl;
315  return false;
316  }
317  for (unsigned int i = 0; i < a.size(); ++i) { c[i] = a[i] + b[i]; }
318  return true;
319 }
320 
321 
322 /* a*b*a^T */
323 bool Similarity(
324  const std::vector<litfloat>& a,
325  const std::vector<litfloat>& b,
326  std::vector<litfloat>& c)
327 {
328  if (a.size() != 25 || b.size() != 15 || c.size() != 15) {
329  std::cout << "-E- Similarity: size is not correct" << std::endl;
330  return false;
331  }
332 
333  litfloat A = a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3] + a[4] * b[4];
334  litfloat B = a[0] * b[1] + a[1] * b[5] + a[2] * b[6] + a[3] * b[7] + a[4] * b[8];
335  litfloat C = a[0] * b[2] + a[1] * b[6] + a[2] * b[9] + a[3] * b[10] + a[4] * b[11];
336  litfloat D = a[0] * b[3] + a[1] * b[7] + a[2] * b[10] + a[3] * b[12] + a[4] * b[13];
337  litfloat E = a[0] * b[4] + a[1] * b[8] + a[2] * b[11] + a[3] * b[13] + a[4] * b[14];
338 
339  litfloat F = a[5] * b[0] + a[6] * b[1] + a[7] * b[2] + a[8] * b[3] + a[9] * b[4];
340  litfloat G = a[5] * b[1] + a[6] * b[5] + a[7] * b[6] + a[8] * b[7] + a[9] * b[8];
341  litfloat H = a[5] * b[2] + a[6] * b[6] + a[7] * b[9] + a[8] * b[10] + a[9] * b[11];
342  litfloat I = a[5] * b[3] + a[6] * b[7] + a[7] * b[10] + a[8] * b[12] + a[9] * b[13];
343  litfloat J = a[5] * b[4] + a[6] * b[8] + a[7] * b[11] + a[8] * b[13] + a[9] * b[14];
344 
345  litfloat K = a[10] * b[0] + a[11] * b[1] + a[12] * b[2] + a[13] * b[3] + a[14] * b[4];
346  litfloat L = a[10] * b[1] + a[11] * b[5] + a[12] * b[6] + a[13] * b[7] + a[14] * b[8];
347  litfloat M = a[10] * b[2] + a[11] * b[6] + a[12] * b[9] + a[13] * b[10] + a[14] * b[11];
348  litfloat N = a[10] * b[3] + a[11] * b[7] + a[12] * b[10] + a[13] * b[12] + a[14] * b[13];
349  litfloat O = a[10] * b[4] + a[11] * b[8] + a[12] * b[11] + a[13] * b[13] + a[14] * b[14];
350 
351  litfloat P = a[15] * b[0] + a[16] * b[1] + a[17] * b[2] + a[18] * b[3] + a[19] * b[4];
352  litfloat Q = a[15] * b[1] + a[16] * b[5] + a[17] * b[6] + a[18] * b[7] + a[19] * b[8];
353  litfloat R = a[15] * b[2] + a[16] * b[6] + a[17] * b[9] + a[18] * b[10] + a[19] * b[11];
354  litfloat S = a[15] * b[3] + a[16] * b[7] + a[17] * b[10] + a[18] * b[12] + a[19] * b[13];
355  litfloat T = a[15] * b[4] + a[16] * b[8] + a[17] * b[11] + a[18] * b[13] + a[19] * b[14];
356 
357  c[0] = A * a[0] + B * a[1] + C * a[2] + D * a[3] + E * a[4];
358  c[1] = A * a[5] + B * a[6] + C * a[7] + D * a[8] + E * a[9];
359  c[2] = A * a[10] + B * a[11] + C * a[12] + D * a[13] + E * a[14];
360  c[3] = A * a[15] + B * a[16] + C * a[17] + D * a[18] + E * a[19];
361  c[4] = A * a[20] + B * a[21] + C * a[22] + D * a[23] + E * a[24];
362 
363  c[5] = F * a[5] + G * a[6] + H * a[7] + I * a[8] + J * a[9];
364  c[6] = F * a[10] + G * a[11] + H * a[12] + I * a[13] + J * a[14];
365  c[7] = F * a[15] + G * a[16] + H * a[17] + I * a[18] + J * a[19];
366  c[8] = F * a[20] + G * a[21] + H * a[22] + I * a[23] + J * a[24];
367 
368  c[9] = K * a[10] + L * a[11] + M * a[12] + N * a[13] + O * a[14];
369  c[10] = K * a[15] + L * a[16] + M * a[17] + N * a[18] + O * a[19];
370  c[11] = K * a[20] + L * a[21] + M * a[22] + N * a[23] + O * a[24];
371 
372  c[12] = P * a[15] + Q * a[16] + R * a[17] + S * a[18] + T * a[19];
373  c[13] = P * a[20] + Q * a[21] + R * a[22] + S * a[23] + T * a[24];
374 
375  c[14] = (a[20] * b[0] + a[21] * b[1] + a[22] * b[2] + a[23] * b[3] + a[24] * b[4]) * a[20] +
376  (a[20] * b[1] + a[21] * b[5] + a[22] * b[6] + a[23] * b[7] + a[24] * b[8]) * a[21] +
377  (a[20] * b[2] + a[21] * b[6] + a[22] * b[9] + a[23] * b[10] + a[24] * b[11]) * a[22] +
378  (a[20] * b[3] + a[21] * b[7] + a[22] * b[10] + a[23] * b[12] + a[24] * b[13]) * a[23] +
379  (a[20] * b[4] + a[21] * b[8] + a[22] * b[11] + a[23] * b[13] + a[24] * b[14]) * a[24];
380  return true;
381 }
382 
383 
384 
385 bool Mult15On25(
386  const std::vector<litfloat>& a,
387  const std::vector<litfloat>& b,
388  std::vector<litfloat>& c)
389 {
390  if (a.size() != 15 || b.size() != 25 || c.size() != 25) {
391  std::cout << "-E- Mult15On25: size is not correct" << std::endl;
392  return false;
393  }
394  c[0] = a[0] * b[0] + a[1] * b[5] + a[2] * b[10] + a[3] * b[15] + a[4] * b[20];
395  c[1] = a[0] * b[1] + a[1] * b[6] + a[2] * b[11] + a[3] * b[16] + a[4] * b[21];
396  c[2] = a[0] * b[2] + a[1] * b[7] + a[2] * b[12] + a[3] * b[17] + a[4] * b[22];
397  c[3] = a[0] * b[3] + a[1] * b[8] + a[2] * b[13] + a[3] * b[18] + a[4] * b[23];
398  c[4] = a[0] * b[4] + a[1] * b[9] + a[2] * b[14] + a[3] * b[19] + a[4] * b[24];
399  c[5] = a[1] * b[0] + a[5] * b[5] + a[6] * b[10] + a[7] * b[15] + a[8] * b[20];
400  c[6] = a[1] * b[1] + a[5] * b[6] + a[6] * b[11] + a[7] * b[16] + a[8] * b[21];
401  c[7] = a[1] * b[2] + a[5] * b[7] + a[6] * b[12] + a[7] * b[17] + a[8] * b[22];
402  c[8] = a[1] * b[3] + a[5] * b[8] + a[6] * b[13] + a[7] * b[18] + a[8] * b[23];
403  c[9] = a[1] * b[4] + a[5] * b[9] + a[6] * b[14] + a[7] * b[19] + a[8] * b[24];
404  c[10] = a[2] * b[0] + a[6] * b[5] + a[9] * b[10] + a[10] * b[15] + a[11] * b[20];
405  c[11] = a[2] * b[1] + a[6] * b[6] + a[9] * b[11] + a[10] * b[16] + a[11] * b[21];
406  c[12] = a[2] * b[2] + a[6] * b[7] + a[9] * b[12] + a[10] * b[17] + a[11] * b[22];
407  c[13] = a[2] * b[3] + a[6] * b[8] + a[9] * b[13] + a[10] * b[18] + a[11] * b[23];
408  c[14] = a[2] * b[4] + a[6] * b[9] + a[9] * b[14] + a[10] * b[19] + a[11] * b[24];
409  c[15] = a[3] * b[0] + a[7] * b[5] + a[10] * b[10] + a[12] * b[15] + a[13] * b[20];
410  c[16] = a[3] * b[1] + a[7] * b[6] + a[10] * b[11] + a[12] * b[16] + a[13] * b[21];
411  c[17] = a[3] * b[2] + a[7] * b[7] + a[10] * b[12] + a[12] * b[17] + a[13] * b[22];
412  c[18] = a[3] * b[3] + a[7] * b[8] + a[10] * b[13] + a[12] * b[18] + a[13] * b[23];
413  c[19] = a[3] * b[4] + a[7] * b[9] + a[10] * b[14] + a[12] * b[19] + a[13] * b[24];
414  c[20] = a[4] * b[0] + a[8] * b[5] + a[11] * b[10] + a[13] * b[15] + a[14] * b[20];
415  c[21] = a[4] * b[1] + a[8] * b[6] + a[11] * b[11] + a[13] * b[16] + a[14] * b[21];
416  c[22] = a[4] * b[2] + a[8] * b[7] + a[11] * b[12] + a[13] * b[17] + a[14] * b[22];
417  c[23] = a[4] * b[3] + a[8] * b[8] + a[11] * b[13] + a[13] * b[18] + a[14] * b[23];
418  c[24] = a[4] * b[4] + a[8] * b[9] + a[11] * b[14] + a[13] * b[19] + a[14] * b[24];
419 
420  return true;
421 }
422 
423 bool Mult25On15(
424  const std::vector<litfloat>& a,
425  const std::vector<litfloat>& b,
426  std::vector<litfloat>& c)
427 {
428  if (a.size() != 25 || b.size() != 15 || c.size() != 25) {
429  std::cout << "-E- Mult15On25: size is not correct" << std::endl;
430  return false;
431  }
432  c[0] = a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3] + a[4] * b[4];
433  c[1] = a[0] * b[1] + a[1] * b[5] + a[2] * b[6] + a[3] * b[7] + a[4] * b[8];
434  c[2] = a[0] * b[2] + a[1] * b[6] + a[2] * b[9] + a[3] * b[10] + a[4] * b[11];
435  c[3] = a[0] * b[3] + a[1] * b[7] + a[2] * b[10] + a[3] * b[12] + a[4] * b[13];
436  c[4] = a[0] * b[4] + a[1] * b[8] + a[2] * b[11] + a[3] * b[13] + a[4] * b[14];
437  c[5] = a[5] * b[0] + a[6] * b[1] + a[7] * b[2] + a[8] * b[3] + a[9] * b[4];
438  c[6] = a[5] * b[1] + a[6] * b[5] + a[7] * b[6] + a[8] * b[7] + a[9] * b[8];
439  c[7] = a[5] * b[2] + a[6] * b[6] + a[7] * b[9] + a[8] * b[10] + a[9] * b[11];
440  c[8] = a[5] * b[3] + a[6] * b[7] + a[7] * b[10] + a[8] * b[12] + a[9] * b[13];
441  c[9] = a[5] * b[4] + a[6] * b[8] + a[7] * b[11] + a[8] * b[13] + a[9] * b[14];
442  c[10] = a[10] * b[0] + a[11] * b[1] + a[12] * b[2] + a[13] * b[3] + a[14] * b[4];
443  c[11] = a[10] * b[1] + a[11] * b[5] + a[12] * b[6] + a[13] * b[7] + a[14] * b[8];
444  c[12] = a[10] * b[2] + a[11] * b[6] + a[12] * b[9] + a[13] * b[10] + a[14] * b[11];
445  c[13] = a[10] * b[3] + a[11] * b[7] + a[12] * b[10] + a[13] * b[12] + a[14] * b[13];
446  c[14] = a[10] * b[4] + a[11] * b[8] + a[12] * b[11] + a[13] * b[13] + a[14] * b[14];
447  c[15] = a[15] * b[0] + a[16] * b[1] + a[17] * b[2] + a[18] * b[3] + a[19] * b[4];
448  c[16] = a[15] * b[1] + a[16] * b[5] + a[17] * b[6] + a[18] * b[7] + a[19] * b[8];
449  c[17] = a[15] * b[2] + a[16] * b[6] + a[17] * b[9] + a[18] * b[10] + a[19] * b[11];
450  c[18] = a[15] * b[3] + a[16] * b[7] + a[17] * b[10] + a[18] * b[12] + a[19] * b[13];
451  c[19] = a[15] * b[4] + a[16] * b[8] + a[17] * b[11] + a[18] * b[13] + a[19] * b[14];
452  c[20] = a[20] * b[0] + a[21] * b[1] + a[22] * b[2] + a[23] * b[3] + a[24] * b[4];
453  c[21] = a[20] * b[1] + a[21] * b[5] + a[22] * b[6] + a[23] * b[7] + a[24] * b[8];
454  c[22] = a[20] * b[2] + a[21] * b[6] + a[22] * b[9] + a[23] * b[10] + a[24] * b[11];
455  c[23] = a[20] * b[3] + a[21] * b[7] + a[22] * b[10] + a[23] * b[12] + a[24] * b[13];
456  c[24] = a[20] * b[4] + a[21] * b[8] + a[22] * b[11] + a[23] * b[13] + a[24] * b[14];
457 
458  return true;
459 }
460