1#![allow(non_snake_case)]
2use ndarray::prelude::*;
3use ndarray::{ScalarOperand, Zip};
4use num_traits::{Float, One, Zero};
5use std::fmt::Debug;
6use std::iter::Sum;
7use std::ops::{Mul, Neg, Sub, AddAssign};
8
9
10pub trait Inverse<T: Float> {
11 fn det(&self) -> T;
12
13 fn inv(&self) -> Option<Self>
14 where
15 Self: Sized;
16
17 fn inv_diag(&self) -> Option<Self>
18 where
19 Self: Sized;
20
21 fn lu_inv(&self) -> Option<Self>
22 where
23 Self: Sized;
24
25 fn inv_lt(&self) -> Option<Self>
26 where
27 Self: Sized;
28
29 fn inv_ut(&self) -> Option<Self>
30 where
31 Self: Sized;
32
33 fn cholesky(&self) -> Self
34 where
35 Self: Sized;
36
37 fn lu(&self) -> Option<(Self, Self, Self)>
38 where
39 Self: Sized;
40}
41
42impl<T> Inverse<T> for Array2<T>
43where
44 T: Float + Debug + ScalarOperand + Sum<T> + AddAssign
45{
46 fn det(&self) -> T {
47 let s = self.raw_dim();
48 assert!(s[0] == s[1]);
49 let vm: Vec<T> = self.iter().copied().collect();
51 determinant(&vm, s[0])
52 }
53
54 fn inv(&self) -> Option<Self> {
55 fn submat<T: Float>(res: &mut [T], m: &Array2<T>, sr: usize, sc: usize) {
56 let s = m.raw_dim();
57 let mut i: usize = 0;
58 (0..s[0]).for_each(|r| (0..s[1]).for_each(|c| {
59 if r != sr && c != sc {
60 res[i] = m[(r, c)];
61 i += 1;
62 }
63 }));
64 }
65
66 let s = self.raw_dim();
67 assert!(s[0] == s[1]);
68 let l = s[0];
69 let det = match l {
72 2 | 3 | 4 => {
73 let d = self.det();
74
75 if d.is_zero() {
76 return None;
77 } else {
78 d
79 }
80 },
81 _ => T::zero(),
82 };
83 match l {
84 1 => Some(array![[self[(0, 0)].recip()]]),
85 2 => Some(array![
86 [self[(1, 1)] / det, -self[(0, 1)] / det],
87 [-self[(1, 0)] / det, self[(0, 0)] / det],
88 ]),
89 3 => {
90 let m00 = self[(0, 0)];
91 let m01 = self[(0, 1)];
92 let m02 = self[(0, 2)];
93 let m10 = self[(1, 0)];
94 let m11 = self[(1, 1)];
95 let m12 = self[(1, 2)];
96 let m20 = self[(2, 0)];
97 let m21 = self[(2, 1)];
98 let m22 = self[(2, 2)];
99
100 let x00 = m11 * m22 - m21 * m12;
101 let x01 = m02 * m21 - m01 * m22;
102 let x02 = m01 * m12 - m02 * m11;
103 let x10 = m12 * m20 - m10 * m22;
104 let x11 = m00 * m22 - m02 * m20;
105 let x12 = m10 * m02 - m00 * m12;
106 let x20 = m10 * m21 - m20 * m11;
107 let x21 = m20 * m01 - m00 * m21;
108 let x22 = m00 * m11 - m10 * m01;
109
110 Some(array![
111 [x00 / det, x01 / det, x02 / det],
112 [x10 / det, x11 / det, x12 / det],
113 [x20 / det, x21 / det, x22 / det]
114 ])
115 }
116 4 => {
117 let m00 = self[(0, 0)];
118 let m01 = self[(0, 1)];
119 let m02 = self[(0, 2)];
120 let m03 = self[(0, 3)];
121 let m10 = self[(1, 0)];
122 let m11 = self[(1, 1)];
123 let m12 = self[(1, 2)];
124 let m13 = self[(1, 3)];
125 let m20 = self[(2, 0)];
126 let m21 = self[(2, 1)];
127 let m22 = self[(2, 2)];
128 let m23 = self[(2, 3)];
129 let m30 = self[(3, 0)];
130 let m31 = self[(3, 1)];
131 let m32 = self[(3, 2)];
132 let m33 = self[(3, 3)];
133
134 let x00 = m11 * m22 * m33
135 - m11 * m32 * m23
136 - m12 * m21 * m33
137 + m12 * m31 * m23
138 + m13 * m21 * m32
139 - m13 * m31 * m22;
140
141 let x10 = -m10 * m22 * m33
142 + m10 * m32 * m23
143 + m12 * m20 * m33
144 - m12 * m30 * m23
145 - m13 * m20 * m32
146 + m13 * m30 * m22;
147
148 let x20 = m10 * m21 * m33
149 - m10 * m31 * m23
150 - m11 * m20 * m33
151 + m11 * m30 * m23
152 + m13 * m20 * m31
153 - m13 * m30 * m21;
154
155 let x30 = -m10 * m21 * m32
156 + m10 * m31 * m22
157 + m11 * m20 * m32
158 - m11 * m30 * m22
159 - m12 * m20 * m31
160 + m12 * m30 * m21;
161
162 let x01 = -m01 * m22 * m33
163 + m01 * m32 * m23
164 + m02 * m21 * m33
165 - m02 * m31 * m23
166 - m03 * m21 * m32
167 + m03 * m31 * m22;
168
169 let x11 = m00 * m22 * m33
170 - m00 * m32 * m23
171 - m02 * m20 * m33
172 + m02 * m30 * m23
173 + m03 * m20 * m32
174 - m03 * m30 * m22;
175
176 let x21 = -m00 * m21 * m33
177 + m00 * m31 * m23
178 + m01 * m20 * m33
179 - m01 * m30 * m23
180 - m03 * m20 * m31
181 + m03 * m30 * m21;
182
183 let x31 = m00 * m21 * m32
184 - m00 * m31 * m22
185 - m01 * m20 * m32
186 + m01 * m30 * m22
187 + m02 * m20 * m31
188 - m02 * m30 * m21;
189
190 let x02 = m01 * m12 * m33
191 - m01 * m32 * m13
192 - m02 * m11 * m33
193 + m02 * m31 * m13
194 + m03 * m11 * m32
195 - m03 * m31 * m12;
196
197 let x12 = -m00 * m12 * m33
198 + m00 * m32 * m13
199 + m02 * m10 * m33
200 - m02 * m30 * m13
201 - m03 * m10 * m32
202 + m03 * m30 * m12;
203
204 let x22 = m00 * m11 * m33
205 - m00 * m31 * m13
206 - m01 * m10 * m33
207 + m01 * m30 * m13
208 + m03 * m10 * m31
209 - m03 * m30 * m11;
210
211 let x03 = -m01 * m12 * m23
212 + m01 * m22 * m13
213 + m02 * m11 * m23
214 - m02 * m21 * m13
215 - m03 * m11 * m22
216 + m03 * m21 * m12;
217
218 let x32 = -m00 * m11 * m32
219 + m00 * m31 * m12
220 + m01 * m10 * m32
221 - m01 * m30 * m12
222 - m02 * m10 * m31
223 + m02 * m30 * m11;
224
225 let x13 = m00 * m12 * m23
226 - m00 * m22 * m13
227 - m02 * m10 * m23
228 + m02 * m20 * m13
229 + m03 * m10 * m22
230 - m03 * m20 * m12;
231
232 let x23 = -m00 * m11 * m23
233 + m00 * m21 * m13
234 + m01 * m10 * m23
235 - m01 * m20 * m13
236 - m03 * m10 * m21
237 + m03 * m20 * m11;
238
239 let x33 = m00 * m11 * m22
240 - m00 * m21 * m12
241 - m01 * m10 * m22
242 + m01 * m20 * m12
243 + m02 * m10 * m21
244 - m02 * m20 * m11;
245
246 Some(array![
247 [x00 / det, x01 / det, x02 / det, x03 / det],
248 [x10 / det, x11 / det, x12 / det, x13 / det],
249 [x20 / det, x21 / det, x22 / det, x23 / det],
250 [x30 / det, x31 / det, x32 / det, x33 / det]
251 ])
252 }
253 _ => {
254 if let Some(res) = self.lu_inv() {
255 Some(res)
256 } else {
257 let det = self.det();
258 if !det.is_zero() {
260 let mut cofactors: Array2<T> = Array2::zeros((l, l));
261 let mut res: Vec<T> = vec![T::zero(); (l - 1) * (l - 1)];
262 for r in 0 .. l {
263 for c in 0 .. l {
264 submat(&mut res, self, r, c);
266 let d = determinant(&res, l - 1);
267 cofactors[(c, r)] = if ((r + c) % 2) == 0 { d } else { -d };
268 }
269 }
270 Some(cofactors / det)
271 } else {
272 None
273 }
274 }
275 }
276 }
277 }
278
279 fn inv_diag(&self) -> Option<Self> {
280 let s = self.raw_dim();
281 assert!(s[0] == s[1]);
282
283 let mut res: Array2<T> = Array2::zeros(s);
284
285 for i in 0 .. s[0] {
286 res[(i, i)] = self[(i, i)].recip();
287 }
288
289 Some(res)
290 }
291
292 fn cholesky(&self) -> Self {
293 let s = self.raw_dim();
294 assert!(s[0] == s[1]);
295 let l = s[0];
296 let mut res: Array2<T> = Array2::zeros(s);
297 for i in 0 .. l {
298 for j in 0 ..= i {
299 let mut s = T::zero();
300 for k in 0 .. j {
301 if !res[(i, k)].is_zero() && !res[(j, k)].is_zero() {
302 s += res[(i, k)] * res[(j, k)];
303 }
304 }
305
306 res[(i, j)] = if i == j {
307 (self[(i, i)] - s).sqrt()
308 } else {
309 res[(j, j)].recip() * (self[(i, j)] - s)
310 };
311 }
312 }
313
314 res
315 }
316
317 fn lu(&self) -> Option<(Self, Self, Self)> {
318 fn pivot<T: Float>(A: &Array2<T>) -> Array2<T> {
319 let n = A.raw_dim()[0];
320 let mut P: Array2<T> = Array::eye(n);
321
322 for (idx, col) in A.axis_iter(Axis(1)).enumerate() {
323 let mut mp = idx;
325 for i in idx .. n {
326 if col[mp].abs() < col[i].abs() {
327 mp = i;
328 }
329 }
330 if mp != idx {
332 let (l, r) = P.multi_slice_mut((s![idx, ..], s![mp, ..]));
333
334 Zip::from(l).and(r).for_each(std::mem::swap);
335 }
336 }
337
338 P
339 }
340
341 let d = self.raw_dim();
342 let n = d[0];
343 assert_eq!(n, d[1], "LU decomposition must take a square matrix.");
344
345 let P = pivot(self);
346 let pA = P.dot(self);
347
348 let mut L: Array2<T> = Array::eye(n);
349 let mut U: Array2<T> = Array::zeros((n, n));
350
351 for c in 0 .. n {
352 for r in 0 .. n {
353 let s = U.slice(s![0..r, c]).dot(&L.slice(s![r, 0..r]));
354
355 if s.is_nan() || s.is_infinite() {
356 return None;
357 }
358
359 let pAs = pA[[r, c]] - s;
360
361 if r < c + 1 { U[[r, c]] = pAs;
363 } else { L[[r, c]] = (pAs) / U[[c, c]];
365 }
366 }
367 }
368
369 Some((L, U, P))
370 }
371
372 fn lu_inv(&self) -> Option<Self> {
373 let d = self.raw_dim();
374 let n = d[0];
375
376 assert!(d[0] == d[1]);
377
378 if let Some((l, u, p)) = self.lu() {
379 let lt = linv(&l, n);
380 let ut = uinv(&u, n);
381
382 Some(ut.dot(<).dot(&p))
383 } else {
384 None
385 }
386 }
387
388 fn inv_lt(&self) -> Option<Self> {
389 let n = check_diag(self);
390
391 if n == 0 {
392 None
393 } else {
394 Some(linv(self, n))
395 }
396 }
397
398 fn inv_ut(&self) -> Option<Self> {
399 let n = check_diag(self);
400
401 if n == 0 {
402 None
403 } else {
404 Some(uinv(self, n))
405 }
406 }
407}
408
409fn determinant<T>(vm: &[T], l: usize) -> T
410where
411 T: Copy + Debug + Zero + One + Mul + PartialEq + Sub<Output=T> + Neg<Output=T> + Sum<T>,
412{
413 match l {
417 1 => vm[0],
418 2 => vm[0] * vm[3] - vm[2] * vm[1],
419 3 => {
420 let m00 = vm[0];
421 let m01 = vm[1];
422 let m02 = vm[2];
423 let m10 = vm[3];
424 let m11 = vm[4];
425 let m12 = vm[5];
426 let m20 = vm[6];
427 let m21 = vm[7];
428 let m22 = vm[8];
429
430 m00 * (m11 * m22 - m12 * m21) - m01
431 * (m10 * m22 - m12 * m20) + m02
432 * (m10 * m21 - m11 * m20)
433 }
434 4 => {
435 let m00 = vm[0];
436 let m01 = vm[1];
437 let m02 = vm[2];
438 let m03 = vm[3];
439 let m10 = vm[4];
440 let m11 = vm[5];
441 let m12 = vm[6];
442 let m13 = vm[7];
443 let m20 = vm[8];
444 let m21 = vm[9];
445 let m22 = vm[10];
446 let m23 = vm[11];
447 let m30 = vm[12];
448 let m31 = vm[13];
449 let m32 = vm[14];
450 let m33 = vm[15];
451
452 m00 * m21 * m32 * m13 -
454 m00 * m21 * m33 * m12 -
455 m00 * m22 * m31 * m13 +
456 m00 * m22 * m33 * m11 +
457 m00 * m23 * m31 * m12 -
458 m00 * m23 * m32 * m11 -
459 m21 * m30 * m02 * m13 +
460 m21 * m30 * m03 * m12 -
461 m21 * m32 * m03 * m10 +
462 m21 * m33 * m02 * m10 +
463 m22 * m30 * m01 * m13 -
464 m22 * m30 * m03 * m11 +
465 m22 * m31 * m03 * m10 -
466 m22 * m33 * m01 * m10 -
467 m23 * m30 * m01 * m12 +
468 m23 * m30 * m02 * m11 -
469 m23 * m31 * m02 * m10 +
470 m23 * m32 * m01 * m10 +
471 m31 * m02 * m13 * m20 -
472 m31 * m03 * m12 * m20 -
473 m32 * m01 * m13 * m20 +
474 m32 * m03 * m11 * m20 +
475 m33 * m01 * m12 * m20 -
476 m33 * m02 * m11 * m20
477 }
478 5 => {
479 let m00 = vm[0];
480 let m01 = vm[1];
481 let m02 = vm[2];
482 let m03 = vm[3];
483 let m04 = vm[4];
484 let m10 = vm[5];
485 let m11 = vm[6];
486 let m12 = vm[7];
487 let m13 = vm[8];
488 let m14 = vm[9];
489 let m20 = vm[10];
490 let m21 = vm[11];
491 let m22 = vm[12];
492 let m23 = vm[13];
493 let m24 = vm[14];
494 let m30 = vm[15];
495 let m31 = vm[16];
496 let m32 = vm[17];
497 let m33 = vm[18];
498 let m34 = vm[19];
499 let m40 = vm[20];
500 let m41 = vm[21];
501 let m42 = vm[22];
502 let m43 = vm[23];
503 let m44 = vm[24];
504
505 m00 * m11 * m22 * m33 * m44
507 - m00 * m11 * m22 * m34 * m43
508 - m00 * m11 * m23 * m32 * m44
509 + m00 * m11 * m23 * m34 * m42
510 + m00 * m11 * m24 * m32 * m43
511 - m00 * m11 * m24 * m33 * m42
512 - m00 * m12 * m21 * m33 * m44
513 + m00 * m12 * m21 * m34 * m43
514 + m00 * m12 * m23 * m31 * m44
515 - m00 * m12 * m23 * m34 * m41
516 - m00 * m12 * m24 * m31 * m43
517 + m00 * m12 * m24 * m33 * m41
518 + m00 * m13 * m21 * m32 * m44
519 - m00 * m13 * m21 * m34 * m42
520 - m00 * m13 * m22 * m31 * m44
521 + m00 * m13 * m22 * m34 * m41
522 + m00 * m13 * m24 * m31 * m42
523 - m00 * m13 * m24 * m32 * m41
524 - m00 * m14 * m21 * m32 * m43
525 + m00 * m14 * m21 * m33 * m42
526 + m00 * m14 * m22 * m31 * m43
527 - m00 * m14 * m22 * m33 * m41
528 - m00 * m14 * m23 * m31 * m42
529 + m00 * m14 * m23 * m32 * m41
530 - m01 * m10 * m22 * m33 * m44
531 + m01 * m10 * m22 * m34 * m43
532 + m01 * m10 * m23 * m32 * m44
533 - m01 * m10 * m23 * m34 * m42
534 - m01 * m10 * m24 * m32 * m43
535 + m01 * m10 * m24 * m33 * m42
536 + m01 * m12 * m20 * m33 * m44
537 - m01 * m12 * m20 * m34 * m43
538 - m01 * m12 * m23 * m30 * m44
539 + m01 * m12 * m23 * m34 * m40
540 + m01 * m12 * m24 * m30 * m43
541 - m01 * m12 * m24 * m33 * m40
542 - m01 * m13 * m20 * m32 * m44
543 + m01 * m13 * m20 * m34 * m42
544 + m01 * m13 * m22 * m30 * m44
545 - m01 * m13 * m22 * m34 * m40
546 - m01 * m13 * m24 * m30 * m42
547 + m01 * m13 * m24 * m32 * m40
548 + m01 * m14 * m20 * m32 * m43
549 - m01 * m14 * m20 * m33 * m42
550 - m01 * m14 * m22 * m30 * m43
551 + m01 * m14 * m22 * m33 * m40
552 + m01 * m14 * m23 * m30 * m42
553 - m01 * m14 * m23 * m32 * m40
554 + m02 * m10 * m21 * m33 * m44
555 - m02 * m10 * m21 * m34 * m43
556 - m02 * m10 * m23 * m31 * m44
557 + m02 * m10 * m23 * m34 * m41
558 + m02 * m10 * m24 * m31 * m43
559 - m02 * m10 * m24 * m33 * m41
560 - m02 * m11 * m20 * m33 * m44
561 + m02 * m11 * m20 * m34 * m43
562 + m02 * m11 * m23 * m30 * m44
563 - m02 * m11 * m23 * m34 * m40
564 - m02 * m11 * m24 * m30 * m43
565 + m02 * m11 * m24 * m33 * m40
566 + m02 * m13 * m20 * m31 * m44
567 - m02 * m13 * m20 * m34 * m41
568 - m02 * m13 * m21 * m30 * m44
569 + m02 * m13 * m21 * m34 * m40
570 + m02 * m13 * m24 * m30 * m41
571 - m02 * m13 * m24 * m31 * m40
572 - m02 * m14 * m20 * m31 * m43
573 + m02 * m14 * m20 * m33 * m41
574 + m02 * m14 * m21 * m30 * m43
575 - m02 * m14 * m21 * m33 * m40
576 - m02 * m14 * m23 * m30 * m41
577 + m02 * m14 * m23 * m31 * m40
578 - m03 * m10 * m21 * m32 * m44
579 + m03 * m10 * m21 * m34 * m42
580 + m03 * m10 * m22 * m31 * m44
581 - m03 * m10 * m22 * m34 * m41
582 - m03 * m10 * m24 * m31 * m42
583 + m03 * m10 * m24 * m32 * m41
584 + m03 * m11 * m20 * m32 * m44
585 - m03 * m11 * m20 * m34 * m42
586 - m03 * m11 * m22 * m30 * m44
587 + m03 * m11 * m22 * m34 * m40
588 + m03 * m11 * m24 * m30 * m42
589 - m03 * m11 * m24 * m32 * m40
590 - m03 * m12 * m20 * m31 * m44
591 + m03 * m12 * m20 * m34 * m41
592 + m03 * m12 * m21 * m30 * m44
593 - m03 * m12 * m21 * m34 * m40
594 - m03 * m12 * m24 * m30 * m41
595 + m03 * m12 * m24 * m31 * m40
596 + m03 * m14 * m20 * m31 * m42
597 - m03 * m14 * m20 * m32 * m41
598 - m03 * m14 * m21 * m30 * m42
599 + m03 * m14 * m21 * m32 * m40
600 + m03 * m14 * m22 * m30 * m41
601 - m03 * m14 * m22 * m31 * m40
602 + m04 * m10 * m21 * m32 * m43
603 - m04 * m10 * m21 * m33 * m42
604 - m04 * m10 * m22 * m31 * m43
605 + m04 * m10 * m22 * m33 * m41
606 + m04 * m10 * m23 * m31 * m42
607 - m04 * m10 * m23 * m32 * m41
608 - m04 * m11 * m20 * m32 * m43
609 + m04 * m11 * m20 * m33 * m42
610 + m04 * m11 * m22 * m30 * m43
611 - m04 * m11 * m22 * m33 * m40
612 - m04 * m11 * m23 * m30 * m42
613 + m04 * m11 * m23 * m32 * m40
614 + m04 * m12 * m20 * m31 * m43
615 - m04 * m12 * m20 * m33 * m41
616 - m04 * m12 * m21 * m30 * m43
617 + m04 * m12 * m21 * m33 * m40
618 + m04 * m12 * m23 * m30 * m41
619 - m04 * m12 * m23 * m31 * m40
620 - m04 * m13 * m20 * m31 * m42
621 + m04 * m13 * m20 * m32 * m41
622 + m04 * m13 * m21 * m30 * m42
623 - m04 * m13 * m21 * m32 * m40
624 - m04 * m13 * m22 * m30 * m41
625 + m04 * m13 * m22 * m31 * m40
626 }
627 6 => {
628 let m00 = vm[0];
630 let m01 = vm[1];
631 let m02 = vm[2];
632 let m03 = vm[3];
633 let m04 = vm[4];
634 let m05 = vm[5];
635 let m10 = vm[6];
636 let m11 = vm[7];
637 let m12 = vm[8];
638 let m13 = vm[9];
639 let m14 = vm[10];
640 let m15 = vm[11];
641 let m20 = vm[12];
642 let m21 = vm[13];
643 let m22 = vm[14];
644 let m23 = vm[15];
645 let m24 = vm[16];
646 let m25 = vm[17];
647 let m30 = vm[18];
648 let m31 = vm[19];
649 let m32 = vm[20];
650 let m33 = vm[21];
651 let m34 = vm[22];
652 let m35 = vm[23];
653 let m40 = vm[24];
654 let m41 = vm[25];
655 let m42 = vm[26];
656 let m43 = vm[27];
657 let m44 = vm[28];
658 let m45 = vm[29];
659 let m50 = vm[30];
660 let m51 = vm[31];
661 let m52 = vm[32];
662 let m53 = vm[33];
663 let m54 = vm[34];
664 let m55 = vm[35];
665
666 m00 * m11 * m22 * m33 * m44 * m55
668 - m00 * m11 * m22 * m33 * m45 * m54
669 - m00 * m11 * m22 * m34 * m43 * m55
670 + m00 * m11 * m22 * m34 * m45 * m53
671 + m00 * m11 * m22 * m35 * m43 * m54
672 - m00 * m11 * m22 * m35 * m44 * m53
673 - m00 * m11 * m23 * m32 * m44 * m55
674 + m00 * m11 * m23 * m32 * m45 * m54
675 + m00 * m11 * m23 * m34 * m42 * m55
676 - m00 * m11 * m23 * m34 * m45 * m52
677 - m00 * m11 * m23 * m35 * m42 * m54
678 + m00 * m11 * m23 * m35 * m44 * m52
679 + m00 * m11 * m24 * m32 * m43 * m55
680 - m00 * m11 * m24 * m32 * m45 * m53
681 - m00 * m11 * m24 * m33 * m42 * m55
682 + m00 * m11 * m24 * m33 * m45 * m52
683 + m00 * m11 * m24 * m35 * m42 * m53
684 - m00 * m11 * m24 * m35 * m43 * m52
685 - m00 * m11 * m25 * m32 * m43 * m54
686 + m00 * m11 * m25 * m32 * m44 * m53
687 + m00 * m11 * m25 * m33 * m42 * m54
688 - m00 * m11 * m25 * m33 * m44 * m52
689 - m00 * m11 * m25 * m34 * m42 * m53
690 + m00 * m11 * m25 * m34 * m43 * m52
691 - m00 * m12 * m21 * m33 * m44 * m55
692 + m00 * m12 * m21 * m33 * m45 * m54
693 + m00 * m12 * m21 * m34 * m43 * m55
694 - m00 * m12 * m21 * m34 * m45 * m53
695 - m00 * m12 * m21 * m35 * m43 * m54
696 + m00 * m12 * m21 * m35 * m44 * m53
697 + m00 * m12 * m23 * m31 * m44 * m55
698 - m00 * m12 * m23 * m31 * m45 * m54
699 - m00 * m12 * m23 * m34 * m41 * m55
700 + m00 * m12 * m23 * m34 * m45 * m51
701 + m00 * m12 * m23 * m35 * m41 * m54
702 - m00 * m12 * m23 * m35 * m44 * m51
703 - m00 * m12 * m24 * m31 * m43 * m55
704 + m00 * m12 * m24 * m31 * m45 * m53
705 + m00 * m12 * m24 * m33 * m41 * m55
706 - m00 * m12 * m24 * m33 * m45 * m51
707 - m00 * m12 * m24 * m35 * m41 * m53
708 + m00 * m12 * m24 * m35 * m43 * m51
709 + m00 * m12 * m25 * m31 * m43 * m54
710 - m00 * m12 * m25 * m31 * m44 * m53
711 - m00 * m12 * m25 * m33 * m41 * m54
712 + m00 * m12 * m25 * m33 * m44 * m51
713 + m00 * m12 * m25 * m34 * m41 * m53
714 - m00 * m12 * m25 * m34 * m43 * m51
715 + m00 * m13 * m21 * m32 * m44 * m55
716 - m00 * m13 * m21 * m32 * m45 * m54
717 - m00 * m13 * m21 * m34 * m42 * m55
718 + m00 * m13 * m21 * m34 * m45 * m52
719 + m00 * m13 * m21 * m35 * m42 * m54
720 - m00 * m13 * m21 * m35 * m44 * m52
721 - m00 * m13 * m22 * m31 * m44 * m55
722 + m00 * m13 * m22 * m31 * m45 * m54
723 + m00 * m13 * m22 * m34 * m41 * m55
724 - m00 * m13 * m22 * m34 * m45 * m51
725 - m00 * m13 * m22 * m35 * m41 * m54
726 + m00 * m13 * m22 * m35 * m44 * m51
727 + m00 * m13 * m24 * m31 * m42 * m55
728 - m00 * m13 * m24 * m31 * m45 * m52
729 - m00 * m13 * m24 * m32 * m41 * m55
730 + m00 * m13 * m24 * m32 * m45 * m51
731 + m00 * m13 * m24 * m35 * m41 * m52
732 - m00 * m13 * m24 * m35 * m42 * m51
733 - m00 * m13 * m25 * m31 * m42 * m54
734 + m00 * m13 * m25 * m31 * m44 * m52
735 + m00 * m13 * m25 * m32 * m41 * m54
736 - m00 * m13 * m25 * m32 * m44 * m51
737 - m00 * m13 * m25 * m34 * m41 * m52
738 + m00 * m13 * m25 * m34 * m42 * m51
739 - m00 * m14 * m21 * m32 * m43 * m55
740 + m00 * m14 * m21 * m32 * m45 * m53
741 + m00 * m14 * m21 * m33 * m42 * m55
742 - m00 * m14 * m21 * m33 * m45 * m52
743 - m00 * m14 * m21 * m35 * m42 * m53
744 + m00 * m14 * m21 * m35 * m43 * m52
745 + m00 * m14 * m22 * m31 * m43 * m55
746 - m00 * m14 * m22 * m31 * m45 * m53
747 - m00 * m14 * m22 * m33 * m41 * m55
748 + m00 * m14 * m22 * m33 * m45 * m51
749 + m00 * m14 * m22 * m35 * m41 * m53
750 - m00 * m14 * m22 * m35 * m43 * m51
751 - m00 * m14 * m23 * m31 * m42 * m55
752 + m00 * m14 * m23 * m31 * m45 * m52
753 + m00 * m14 * m23 * m32 * m41 * m55
754 - m00 * m14 * m23 * m32 * m45 * m51
755 - m00 * m14 * m23 * m35 * m41 * m52
756 + m00 * m14 * m23 * m35 * m42 * m51
757 + m00 * m14 * m25 * m31 * m42 * m53
758 - m00 * m14 * m25 * m31 * m43 * m52
759 - m00 * m14 * m25 * m32 * m41 * m53
760 + m00 * m14 * m25 * m32 * m43 * m51
761 + m00 * m14 * m25 * m33 * m41 * m52
762 - m00 * m14 * m25 * m33 * m42 * m51
763 + m00 * m15 * m21 * m32 * m43 * m54
764 - m00 * m15 * m21 * m32 * m44 * m53
765 - m00 * m15 * m21 * m33 * m42 * m54
766 + m00 * m15 * m21 * m33 * m44 * m52
767 + m00 * m15 * m21 * m34 * m42 * m53
768 - m00 * m15 * m21 * m34 * m43 * m52
769 - m00 * m15 * m22 * m31 * m43 * m54
770 + m00 * m15 * m22 * m31 * m44 * m53
771 + m00 * m15 * m22 * m33 * m41 * m54
772 - m00 * m15 * m22 * m33 * m44 * m51
773 - m00 * m15 * m22 * m34 * m41 * m53
774 + m00 * m15 * m22 * m34 * m43 * m51
775 + m00 * m15 * m23 * m31 * m42 * m54
776 - m00 * m15 * m23 * m31 * m44 * m52
777 - m00 * m15 * m23 * m32 * m41 * m54
778 + m00 * m15 * m23 * m32 * m44 * m51
779 + m00 * m15 * m23 * m34 * m41 * m52
780 - m00 * m15 * m23 * m34 * m42 * m51
781 - m00 * m15 * m24 * m31 * m42 * m53
782 + m00 * m15 * m24 * m31 * m43 * m52
783 + m00 * m15 * m24 * m32 * m41 * m53
784 - m00 * m15 * m24 * m32 * m43 * m51
785 - m00 * m15 * m24 * m33 * m41 * m52
786 + m00 * m15 * m24 * m33 * m42 * m51
787 - m01 * m10 * m22 * m33 * m44 * m55
788 + m01 * m10 * m22 * m33 * m45 * m54
789 + m01 * m10 * m22 * m34 * m43 * m55
790 - m01 * m10 * m22 * m34 * m45 * m53
791 - m01 * m10 * m22 * m35 * m43 * m54
792 + m01 * m10 * m22 * m35 * m44 * m53
793 + m01 * m10 * m23 * m32 * m44 * m55
794 - m01 * m10 * m23 * m32 * m45 * m54
795 - m01 * m10 * m23 * m34 * m42 * m55
796 + m01 * m10 * m23 * m34 * m45 * m52
797 + m01 * m10 * m23 * m35 * m42 * m54
798 - m01 * m10 * m23 * m35 * m44 * m52
799 - m01 * m10 * m24 * m32 * m43 * m55
800 + m01 * m10 * m24 * m32 * m45 * m53
801 + m01 * m10 * m24 * m33 * m42 * m55
802 - m01 * m10 * m24 * m33 * m45 * m52
803 - m01 * m10 * m24 * m35 * m42 * m53
804 + m01 * m10 * m24 * m35 * m43 * m52
805 + m01 * m10 * m25 * m32 * m43 * m54
806 - m01 * m10 * m25 * m32 * m44 * m53
807 - m01 * m10 * m25 * m33 * m42 * m54
808 + m01 * m10 * m25 * m33 * m44 * m52
809 + m01 * m10 * m25 * m34 * m42 * m53
810 - m01 * m10 * m25 * m34 * m43 * m52
811 + m01 * m12 * m20 * m33 * m44 * m55
812 - m01 * m12 * m20 * m33 * m45 * m54
813 - m01 * m12 * m20 * m34 * m43 * m55
814 + m01 * m12 * m20 * m34 * m45 * m53
815 + m01 * m12 * m20 * m35 * m43 * m54
816 - m01 * m12 * m20 * m35 * m44 * m53
817 - m01 * m12 * m23 * m30 * m44 * m55
818 + m01 * m12 * m23 * m30 * m45 * m54
819 + m01 * m12 * m23 * m34 * m40 * m55
820 - m01 * m12 * m23 * m34 * m45 * m50
821 - m01 * m12 * m23 * m35 * m40 * m54
822 + m01 * m12 * m23 * m35 * m44 * m50
823 + m01 * m12 * m24 * m30 * m43 * m55
824 - m01 * m12 * m24 * m30 * m45 * m53
825 - m01 * m12 * m24 * m33 * m40 * m55
826 + m01 * m12 * m24 * m33 * m45 * m50
827 + m01 * m12 * m24 * m35 * m40 * m53
828 - m01 * m12 * m24 * m35 * m43 * m50
829 - m01 * m12 * m25 * m30 * m43 * m54
830 + m01 * m12 * m25 * m30 * m44 * m53
831 + m01 * m12 * m25 * m33 * m40 * m54
832 - m01 * m12 * m25 * m33 * m44 * m50
833 - m01 * m12 * m25 * m34 * m40 * m53
834 + m01 * m12 * m25 * m34 * m43 * m50
835 - m01 * m13 * m20 * m32 * m44 * m55
836 + m01 * m13 * m20 * m32 * m45 * m54
837 + m01 * m13 * m20 * m34 * m42 * m55
838 - m01 * m13 * m20 * m34 * m45 * m52
839 - m01 * m13 * m20 * m35 * m42 * m54
840 + m01 * m13 * m20 * m35 * m44 * m52
841 + m01 * m13 * m22 * m30 * m44 * m55
842 - m01 * m13 * m22 * m30 * m45 * m54
843 - m01 * m13 * m22 * m34 * m40 * m55
844 + m01 * m13 * m22 * m34 * m45 * m50
845 + m01 * m13 * m22 * m35 * m40 * m54
846 - m01 * m13 * m22 * m35 * m44 * m50
847 - m01 * m13 * m24 * m30 * m42 * m55
848 + m01 * m13 * m24 * m30 * m45 * m52
849 + m01 * m13 * m24 * m32 * m40 * m55
850 - m01 * m13 * m24 * m32 * m45 * m50
851 - m01 * m13 * m24 * m35 * m40 * m52
852 + m01 * m13 * m24 * m35 * m42 * m50
853 + m01 * m13 * m25 * m30 * m42 * m54
854 - m01 * m13 * m25 * m30 * m44 * m52
855 - m01 * m13 * m25 * m32 * m40 * m54
856 + m01 * m13 * m25 * m32 * m44 * m50
857 + m01 * m13 * m25 * m34 * m40 * m52
858 - m01 * m13 * m25 * m34 * m42 * m50
859 + m01 * m14 * m20 * m32 * m43 * m55
860 - m01 * m14 * m20 * m32 * m45 * m53
861 - m01 * m14 * m20 * m33 * m42 * m55
862 + m01 * m14 * m20 * m33 * m45 * m52
863 + m01 * m14 * m20 * m35 * m42 * m53
864 - m01 * m14 * m20 * m35 * m43 * m52
865 - m01 * m14 * m22 * m30 * m43 * m55
866 + m01 * m14 * m22 * m30 * m45 * m53
867 + m01 * m14 * m22 * m33 * m40 * m55
868 - m01 * m14 * m22 * m33 * m45 * m50
869 - m01 * m14 * m22 * m35 * m40 * m53
870 + m01 * m14 * m22 * m35 * m43 * m50
871 + m01 * m14 * m23 * m30 * m42 * m55
872 - m01 * m14 * m23 * m30 * m45 * m52
873 - m01 * m14 * m23 * m32 * m40 * m55
874 + m01 * m14 * m23 * m32 * m45 * m50
875 + m01 * m14 * m23 * m35 * m40 * m52
876 - m01 * m14 * m23 * m35 * m42 * m50
877 - m01 * m14 * m25 * m30 * m42 * m53
878 + m01 * m14 * m25 * m30 * m43 * m52
879 + m01 * m14 * m25 * m32 * m40 * m53
880 - m01 * m14 * m25 * m32 * m43 * m50
881 - m01 * m14 * m25 * m33 * m40 * m52
882 + m01 * m14 * m25 * m33 * m42 * m50
883 - m01 * m15 * m20 * m32 * m43 * m54
884 + m01 * m15 * m20 * m32 * m44 * m53
885 + m01 * m15 * m20 * m33 * m42 * m54
886 - m01 * m15 * m20 * m33 * m44 * m52
887 - m01 * m15 * m20 * m34 * m42 * m53
888 + m01 * m15 * m20 * m34 * m43 * m52
889 + m01 * m15 * m22 * m30 * m43 * m54
890 - m01 * m15 * m22 * m30 * m44 * m53
891 - m01 * m15 * m22 * m33 * m40 * m54
892 + m01 * m15 * m22 * m33 * m44 * m50
893 + m01 * m15 * m22 * m34 * m40 * m53
894 - m01 * m15 * m22 * m34 * m43 * m50
895 - m01 * m15 * m23 * m30 * m42 * m54
896 + m01 * m15 * m23 * m30 * m44 * m52
897 + m01 * m15 * m23 * m32 * m40 * m54
898 - m01 * m15 * m23 * m32 * m44 * m50
899 - m01 * m15 * m23 * m34 * m40 * m52
900 + m01 * m15 * m23 * m34 * m42 * m50
901 + m01 * m15 * m24 * m30 * m42 * m53
902 - m01 * m15 * m24 * m30 * m43 * m52
903 - m01 * m15 * m24 * m32 * m40 * m53
904 + m01 * m15 * m24 * m32 * m43 * m50
905 + m01 * m15 * m24 * m33 * m40 * m52
906 - m01 * m15 * m24 * m33 * m42 * m50
907 + m02 * m10 * m21 * m33 * m44 * m55
908 - m02 * m10 * m21 * m33 * m45 * m54
909 - m02 * m10 * m21 * m34 * m43 * m55
910 + m02 * m10 * m21 * m34 * m45 * m53
911 + m02 * m10 * m21 * m35 * m43 * m54
912 - m02 * m10 * m21 * m35 * m44 * m53
913 - m02 * m10 * m23 * m31 * m44 * m55
914 + m02 * m10 * m23 * m31 * m45 * m54
915 + m02 * m10 * m23 * m34 * m41 * m55
916 - m02 * m10 * m23 * m34 * m45 * m51
917 - m02 * m10 * m23 * m35 * m41 * m54
918 + m02 * m10 * m23 * m35 * m44 * m51
919 + m02 * m10 * m24 * m31 * m43 * m55
920 - m02 * m10 * m24 * m31 * m45 * m53
921 - m02 * m10 * m24 * m33 * m41 * m55
922 + m02 * m10 * m24 * m33 * m45 * m51
923 + m02 * m10 * m24 * m35 * m41 * m53
924 - m02 * m10 * m24 * m35 * m43 * m51
925 - m02 * m10 * m25 * m31 * m43 * m54
926 + m02 * m10 * m25 * m31 * m44 * m53
927 + m02 * m10 * m25 * m33 * m41 * m54
928 - m02 * m10 * m25 * m33 * m44 * m51
929 - m02 * m10 * m25 * m34 * m41 * m53
930 + m02 * m10 * m25 * m34 * m43 * m51
931 - m02 * m11 * m20 * m33 * m44 * m55
932 + m02 * m11 * m20 * m33 * m45 * m54
933 + m02 * m11 * m20 * m34 * m43 * m55
934 - m02 * m11 * m20 * m34 * m45 * m53
935 - m02 * m11 * m20 * m35 * m43 * m54
936 + m02 * m11 * m20 * m35 * m44 * m53
937 + m02 * m11 * m23 * m30 * m44 * m55
938 - m02 * m11 * m23 * m30 * m45 * m54
939 - m02 * m11 * m23 * m34 * m40 * m55
940 + m02 * m11 * m23 * m34 * m45 * m50
941 + m02 * m11 * m23 * m35 * m40 * m54
942 - m02 * m11 * m23 * m35 * m44 * m50
943 - m02 * m11 * m24 * m30 * m43 * m55
944 + m02 * m11 * m24 * m30 * m45 * m53
945 + m02 * m11 * m24 * m33 * m40 * m55
946 - m02 * m11 * m24 * m33 * m45 * m50
947 - m02 * m11 * m24 * m35 * m40 * m53
948 + m02 * m11 * m24 * m35 * m43 * m50
949 + m02 * m11 * m25 * m30 * m43 * m54
950 - m02 * m11 * m25 * m30 * m44 * m53
951 - m02 * m11 * m25 * m33 * m40 * m54
952 + m02 * m11 * m25 * m33 * m44 * m50
953 + m02 * m11 * m25 * m34 * m40 * m53
954 - m02 * m11 * m25 * m34 * m43 * m50
955 + m02 * m13 * m20 * m31 * m44 * m55
956 - m02 * m13 * m20 * m31 * m45 * m54
957 - m02 * m13 * m20 * m34 * m41 * m55
958 + m02 * m13 * m20 * m34 * m45 * m51
959 + m02 * m13 * m20 * m35 * m41 * m54
960 - m02 * m13 * m20 * m35 * m44 * m51
961 - m02 * m13 * m21 * m30 * m44 * m55
962 + m02 * m13 * m21 * m30 * m45 * m54
963 + m02 * m13 * m21 * m34 * m40 * m55
964 - m02 * m13 * m21 * m34 * m45 * m50
965 - m02 * m13 * m21 * m35 * m40 * m54
966 + m02 * m13 * m21 * m35 * m44 * m50
967 + m02 * m13 * m24 * m30 * m41 * m55
968 - m02 * m13 * m24 * m30 * m45 * m51
969 - m02 * m13 * m24 * m31 * m40 * m55
970 + m02 * m13 * m24 * m31 * m45 * m50
971 + m02 * m13 * m24 * m35 * m40 * m51
972 - m02 * m13 * m24 * m35 * m41 * m50
973 - m02 * m13 * m25 * m30 * m41 * m54
974 + m02 * m13 * m25 * m30 * m44 * m51
975 + m02 * m13 * m25 * m31 * m40 * m54
976 - m02 * m13 * m25 * m31 * m44 * m50
977 - m02 * m13 * m25 * m34 * m40 * m51
978 + m02 * m13 * m25 * m34 * m41 * m50
979 - m02 * m14 * m20 * m31 * m43 * m55
980 + m02 * m14 * m20 * m31 * m45 * m53
981 + m02 * m14 * m20 * m33 * m41 * m55
982 - m02 * m14 * m20 * m33 * m45 * m51
983 - m02 * m14 * m20 * m35 * m41 * m53
984 + m02 * m14 * m20 * m35 * m43 * m51
985 + m02 * m14 * m21 * m30 * m43 * m55
986 - m02 * m14 * m21 * m30 * m45 * m53
987 - m02 * m14 * m21 * m33 * m40 * m55
988 + m02 * m14 * m21 * m33 * m45 * m50
989 + m02 * m14 * m21 * m35 * m40 * m53
990 - m02 * m14 * m21 * m35 * m43 * m50
991 - m02 * m14 * m23 * m30 * m41 * m55
992 + m02 * m14 * m23 * m30 * m45 * m51
993 + m02 * m14 * m23 * m31 * m40 * m55
994 - m02 * m14 * m23 * m31 * m45 * m50
995 - m02 * m14 * m23 * m35 * m40 * m51
996 + m02 * m14 * m23 * m35 * m41 * m50
997 + m02 * m14 * m25 * m30 * m41 * m53
998 - m02 * m14 * m25 * m30 * m43 * m51
999 - m02 * m14 * m25 * m31 * m40 * m53
1000 + m02 * m14 * m25 * m31 * m43 * m50
1001 + m02 * m14 * m25 * m33 * m40 * m51
1002 - m02 * m14 * m25 * m33 * m41 * m50
1003 + m02 * m15 * m20 * m31 * m43 * m54
1004 - m02 * m15 * m20 * m31 * m44 * m53
1005 - m02 * m15 * m20 * m33 * m41 * m54
1006 + m02 * m15 * m20 * m33 * m44 * m51
1007 + m02 * m15 * m20 * m34 * m41 * m53
1008 - m02 * m15 * m20 * m34 * m43 * m51
1009 - m02 * m15 * m21 * m30 * m43 * m54
1010 + m02 * m15 * m21 * m30 * m44 * m53
1011 + m02 * m15 * m21 * m33 * m40 * m54
1012 - m02 * m15 * m21 * m33 * m44 * m50
1013 - m02 * m15 * m21 * m34 * m40 * m53
1014 + m02 * m15 * m21 * m34 * m43 * m50
1015 + m02 * m15 * m23 * m30 * m41 * m54
1016 - m02 * m15 * m23 * m30 * m44 * m51
1017 - m02 * m15 * m23 * m31 * m40 * m54
1018 + m02 * m15 * m23 * m31 * m44 * m50
1019 + m02 * m15 * m23 * m34 * m40 * m51
1020 - m02 * m15 * m23 * m34 * m41 * m50
1021 - m02 * m15 * m24 * m30 * m41 * m53
1022 + m02 * m15 * m24 * m30 * m43 * m51
1023 + m02 * m15 * m24 * m31 * m40 * m53
1024 - m02 * m15 * m24 * m31 * m43 * m50
1025 - m02 * m15 * m24 * m33 * m40 * m51
1026 + m02 * m15 * m24 * m33 * m41 * m50
1027 - m03 * m10 * m21 * m32 * m44 * m55
1028 + m03 * m10 * m21 * m32 * m45 * m54
1029 + m03 * m10 * m21 * m34 * m42 * m55
1030 - m03 * m10 * m21 * m34 * m45 * m52
1031 - m03 * m10 * m21 * m35 * m42 * m54
1032 + m03 * m10 * m21 * m35 * m44 * m52
1033 + m03 * m10 * m22 * m31 * m44 * m55
1034 - m03 * m10 * m22 * m31 * m45 * m54
1035 - m03 * m10 * m22 * m34 * m41 * m55
1036 + m03 * m10 * m22 * m34 * m45 * m51
1037 + m03 * m10 * m22 * m35 * m41 * m54
1038 - m03 * m10 * m22 * m35 * m44 * m51
1039 - m03 * m10 * m24 * m31 * m42 * m55
1040 + m03 * m10 * m24 * m31 * m45 * m52
1041 + m03 * m10 * m24 * m32 * m41 * m55
1042 - m03 * m10 * m24 * m32 * m45 * m51
1043 - m03 * m10 * m24 * m35 * m41 * m52
1044 + m03 * m10 * m24 * m35 * m42 * m51
1045 + m03 * m10 * m25 * m31 * m42 * m54
1046 - m03 * m10 * m25 * m31 * m44 * m52
1047 - m03 * m10 * m25 * m32 * m41 * m54
1048 + m03 * m10 * m25 * m32 * m44 * m51
1049 + m03 * m10 * m25 * m34 * m41 * m52
1050 - m03 * m10 * m25 * m34 * m42 * m51
1051 + m03 * m11 * m20 * m32 * m44 * m55
1052 - m03 * m11 * m20 * m32 * m45 * m54
1053 - m03 * m11 * m20 * m34 * m42 * m55
1054 + m03 * m11 * m20 * m34 * m45 * m52
1055 + m03 * m11 * m20 * m35 * m42 * m54
1056 - m03 * m11 * m20 * m35 * m44 * m52
1057 - m03 * m11 * m22 * m30 * m44 * m55
1058 + m03 * m11 * m22 * m30 * m45 * m54
1059 + m03 * m11 * m22 * m34 * m40 * m55
1060 - m03 * m11 * m22 * m34 * m45 * m50
1061 - m03 * m11 * m22 * m35 * m40 * m54
1062 + m03 * m11 * m22 * m35 * m44 * m50
1063 + m03 * m11 * m24 * m30 * m42 * m55
1064 - m03 * m11 * m24 * m30 * m45 * m52
1065 - m03 * m11 * m24 * m32 * m40 * m55
1066 + m03 * m11 * m24 * m32 * m45 * m50
1067 + m03 * m11 * m24 * m35 * m40 * m52
1068 - m03 * m11 * m24 * m35 * m42 * m50
1069 - m03 * m11 * m25 * m30 * m42 * m54
1070 + m03 * m11 * m25 * m30 * m44 * m52
1071 + m03 * m11 * m25 * m32 * m40 * m54
1072 - m03 * m11 * m25 * m32 * m44 * m50
1073 - m03 * m11 * m25 * m34 * m40 * m52
1074 + m03 * m11 * m25 * m34 * m42 * m50
1075 - m03 * m12 * m20 * m31 * m44 * m55
1076 + m03 * m12 * m20 * m31 * m45 * m54
1077 + m03 * m12 * m20 * m34 * m41 * m55
1078 - m03 * m12 * m20 * m34 * m45 * m51
1079 - m03 * m12 * m20 * m35 * m41 * m54
1080 + m03 * m12 * m20 * m35 * m44 * m51
1081 + m03 * m12 * m21 * m30 * m44 * m55
1082 - m03 * m12 * m21 * m30 * m45 * m54
1083 - m03 * m12 * m21 * m34 * m40 * m55
1084 + m03 * m12 * m21 * m34 * m45 * m50
1085 + m03 * m12 * m21 * m35 * m40 * m54
1086 - m03 * m12 * m21 * m35 * m44 * m50
1087 - m03 * m12 * m24 * m30 * m41 * m55
1088 + m03 * m12 * m24 * m30 * m45 * m51
1089 + m03 * m12 * m24 * m31 * m40 * m55
1090 - m03 * m12 * m24 * m31 * m45 * m50
1091 - m03 * m12 * m24 * m35 * m40 * m51
1092 + m03 * m12 * m24 * m35 * m41 * m50
1093 + m03 * m12 * m25 * m30 * m41 * m54
1094 - m03 * m12 * m25 * m30 * m44 * m51
1095 - m03 * m12 * m25 * m31 * m40 * m54
1096 + m03 * m12 * m25 * m31 * m44 * m50
1097 + m03 * m12 * m25 * m34 * m40 * m51
1098 - m03 * m12 * m25 * m34 * m41 * m50
1099 + m03 * m14 * m20 * m31 * m42 * m55
1100 - m03 * m14 * m20 * m31 * m45 * m52
1101 - m03 * m14 * m20 * m32 * m41 * m55
1102 + m03 * m14 * m20 * m32 * m45 * m51
1103 + m03 * m14 * m20 * m35 * m41 * m52
1104 - m03 * m14 * m20 * m35 * m42 * m51
1105 - m03 * m14 * m21 * m30 * m42 * m55
1106 + m03 * m14 * m21 * m30 * m45 * m52
1107 + m03 * m14 * m21 * m32 * m40 * m55
1108 - m03 * m14 * m21 * m32 * m45 * m50
1109 - m03 * m14 * m21 * m35 * m40 * m52
1110 + m03 * m14 * m21 * m35 * m42 * m50
1111 + m03 * m14 * m22 * m30 * m41 * m55
1112 - m03 * m14 * m22 * m30 * m45 * m51
1113 - m03 * m14 * m22 * m31 * m40 * m55
1114 + m03 * m14 * m22 * m31 * m45 * m50
1115 + m03 * m14 * m22 * m35 * m40 * m51
1116 - m03 * m14 * m22 * m35 * m41 * m50
1117 - m03 * m14 * m25 * m30 * m41 * m52
1118 + m03 * m14 * m25 * m30 * m42 * m51
1119 + m03 * m14 * m25 * m31 * m40 * m52
1120 - m03 * m14 * m25 * m31 * m42 * m50
1121 - m03 * m14 * m25 * m32 * m40 * m51
1122 + m03 * m14 * m25 * m32 * m41 * m50
1123 - m03 * m15 * m20 * m31 * m42 * m54
1124 + m03 * m15 * m20 * m31 * m44 * m52
1125 + m03 * m15 * m20 * m32 * m41 * m54
1126 - m03 * m15 * m20 * m32 * m44 * m51
1127 - m03 * m15 * m20 * m34 * m41 * m52
1128 + m03 * m15 * m20 * m34 * m42 * m51
1129 + m03 * m15 * m21 * m30 * m42 * m54
1130 - m03 * m15 * m21 * m30 * m44 * m52
1131 - m03 * m15 * m21 * m32 * m40 * m54
1132 + m03 * m15 * m21 * m32 * m44 * m50
1133 + m03 * m15 * m21 * m34 * m40 * m52
1134 - m03 * m15 * m21 * m34 * m42 * m50
1135 - m03 * m15 * m22 * m30 * m41 * m54
1136 + m03 * m15 * m22 * m30 * m44 * m51
1137 + m03 * m15 * m22 * m31 * m40 * m54
1138 - m03 * m15 * m22 * m31 * m44 * m50
1139 - m03 * m15 * m22 * m34 * m40 * m51
1140 + m03 * m15 * m22 * m34 * m41 * m50
1141 + m03 * m15 * m24 * m30 * m41 * m52
1142 - m03 * m15 * m24 * m30 * m42 * m51
1143 - m03 * m15 * m24 * m31 * m40 * m52
1144 + m03 * m15 * m24 * m31 * m42 * m50
1145 + m03 * m15 * m24 * m32 * m40 * m51
1146 - m03 * m15 * m24 * m32 * m41 * m50
1147 + m04 * m10 * m21 * m32 * m43 * m55
1148 - m04 * m10 * m21 * m32 * m45 * m53
1149 - m04 * m10 * m21 * m33 * m42 * m55
1150 + m04 * m10 * m21 * m33 * m45 * m52
1151 + m04 * m10 * m21 * m35 * m42 * m53
1152 - m04 * m10 * m21 * m35 * m43 * m52
1153 - m04 * m10 * m22 * m31 * m43 * m55
1154 + m04 * m10 * m22 * m31 * m45 * m53
1155 + m04 * m10 * m22 * m33 * m41 * m55
1156 - m04 * m10 * m22 * m33 * m45 * m51
1157 - m04 * m10 * m22 * m35 * m41 * m53
1158 + m04 * m10 * m22 * m35 * m43 * m51
1159 + m04 * m10 * m23 * m31 * m42 * m55
1160 - m04 * m10 * m23 * m31 * m45 * m52
1161 - m04 * m10 * m23 * m32 * m41 * m55
1162 + m04 * m10 * m23 * m32 * m45 * m51
1163 + m04 * m10 * m23 * m35 * m41 * m52
1164 - m04 * m10 * m23 * m35 * m42 * m51
1165 - m04 * m10 * m25 * m31 * m42 * m53
1166 + m04 * m10 * m25 * m31 * m43 * m52
1167 + m04 * m10 * m25 * m32 * m41 * m53
1168 - m04 * m10 * m25 * m32 * m43 * m51
1169 - m04 * m10 * m25 * m33 * m41 * m52
1170 + m04 * m10 * m25 * m33 * m42 * m51
1171 - m04 * m11 * m20 * m32 * m43 * m55
1172 + m04 * m11 * m20 * m32 * m45 * m53
1173 + m04 * m11 * m20 * m33 * m42 * m55
1174 - m04 * m11 * m20 * m33 * m45 * m52
1175 - m04 * m11 * m20 * m35 * m42 * m53
1176 + m04 * m11 * m20 * m35 * m43 * m52
1177 + m04 * m11 * m22 * m30 * m43 * m55
1178 - m04 * m11 * m22 * m30 * m45 * m53
1179 - m04 * m11 * m22 * m33 * m40 * m55
1180 + m04 * m11 * m22 * m33 * m45 * m50
1181 + m04 * m11 * m22 * m35 * m40 * m53
1182 - m04 * m11 * m22 * m35 * m43 * m50
1183 - m04 * m11 * m23 * m30 * m42 * m55
1184 + m04 * m11 * m23 * m30 * m45 * m52
1185 + m04 * m11 * m23 * m32 * m40 * m55
1186 - m04 * m11 * m23 * m32 * m45 * m50
1187 - m04 * m11 * m23 * m35 * m40 * m52
1188 + m04 * m11 * m23 * m35 * m42 * m50
1189 + m04 * m11 * m25 * m30 * m42 * m53
1190 - m04 * m11 * m25 * m30 * m43 * m52
1191 - m04 * m11 * m25 * m32 * m40 * m53
1192 + m04 * m11 * m25 * m32 * m43 * m50
1193 + m04 * m11 * m25 * m33 * m40 * m52
1194 - m04 * m11 * m25 * m33 * m42 * m50
1195 + m04 * m12 * m20 * m31 * m43 * m55
1196 - m04 * m12 * m20 * m31 * m45 * m53
1197 - m04 * m12 * m20 * m33 * m41 * m55
1198 + m04 * m12 * m20 * m33 * m45 * m51
1199 + m04 * m12 * m20 * m35 * m41 * m53
1200 - m04 * m12 * m20 * m35 * m43 * m51
1201 - m04 * m12 * m21 * m30 * m43 * m55
1202 + m04 * m12 * m21 * m30 * m45 * m53
1203 + m04 * m12 * m21 * m33 * m40 * m55
1204 - m04 * m12 * m21 * m33 * m45 * m50
1205 - m04 * m12 * m21 * m35 * m40 * m53
1206 + m04 * m12 * m21 * m35 * m43 * m50
1207 + m04 * m12 * m23 * m30 * m41 * m55
1208 - m04 * m12 * m23 * m30 * m45 * m51
1209 - m04 * m12 * m23 * m31 * m40 * m55
1210 + m04 * m12 * m23 * m31 * m45 * m50
1211 + m04 * m12 * m23 * m35 * m40 * m51
1212 - m04 * m12 * m23 * m35 * m41 * m50
1213 - m04 * m12 * m25 * m30 * m41 * m53
1214 + m04 * m12 * m25 * m30 * m43 * m51
1215 + m04 * m12 * m25 * m31 * m40 * m53
1216 - m04 * m12 * m25 * m31 * m43 * m50
1217 - m04 * m12 * m25 * m33 * m40 * m51
1218 + m04 * m12 * m25 * m33 * m41 * m50
1219 - m04 * m13 * m20 * m31 * m42 * m55
1220 + m04 * m13 * m20 * m31 * m45 * m52
1221 + m04 * m13 * m20 * m32 * m41 * m55
1222 - m04 * m13 * m20 * m32 * m45 * m51
1223 - m04 * m13 * m20 * m35 * m41 * m52
1224 + m04 * m13 * m20 * m35 * m42 * m51
1225 + m04 * m13 * m21 * m30 * m42 * m55
1226 - m04 * m13 * m21 * m30 * m45 * m52
1227 - m04 * m13 * m21 * m32 * m40 * m55
1228 + m04 * m13 * m21 * m32 * m45 * m50
1229 + m04 * m13 * m21 * m35 * m40 * m52
1230 - m04 * m13 * m21 * m35 * m42 * m50
1231 - m04 * m13 * m22 * m30 * m41 * m55
1232 + m04 * m13 * m22 * m30 * m45 * m51
1233 + m04 * m13 * m22 * m31 * m40 * m55
1234 - m04 * m13 * m22 * m31 * m45 * m50
1235 - m04 * m13 * m22 * m35 * m40 * m51
1236 + m04 * m13 * m22 * m35 * m41 * m50
1237 + m04 * m13 * m25 * m30 * m41 * m52
1238 - m04 * m13 * m25 * m30 * m42 * m51
1239 - m04 * m13 * m25 * m31 * m40 * m52
1240 + m04 * m13 * m25 * m31 * m42 * m50
1241 + m04 * m13 * m25 * m32 * m40 * m51
1242 - m04 * m13 * m25 * m32 * m41 * m50
1243 + m04 * m15 * m20 * m31 * m42 * m53
1244 - m04 * m15 * m20 * m31 * m43 * m52
1245 - m04 * m15 * m20 * m32 * m41 * m53
1246 + m04 * m15 * m20 * m32 * m43 * m51
1247 + m04 * m15 * m20 * m33 * m41 * m52
1248 - m04 * m15 * m20 * m33 * m42 * m51
1249 - m04 * m15 * m21 * m30 * m42 * m53
1250 + m04 * m15 * m21 * m30 * m43 * m52
1251 + m04 * m15 * m21 * m32 * m40 * m53
1252 - m04 * m15 * m21 * m32 * m43 * m50
1253 - m04 * m15 * m21 * m33 * m40 * m52
1254 + m04 * m15 * m21 * m33 * m42 * m50
1255 + m04 * m15 * m22 * m30 * m41 * m53
1256 - m04 * m15 * m22 * m30 * m43 * m51
1257 - m04 * m15 * m22 * m31 * m40 * m53
1258 + m04 * m15 * m22 * m31 * m43 * m50
1259 + m04 * m15 * m22 * m33 * m40 * m51
1260 - m04 * m15 * m22 * m33 * m41 * m50
1261 - m04 * m15 * m23 * m30 * m41 * m52
1262 + m04 * m15 * m23 * m30 * m42 * m51
1263 + m04 * m15 * m23 * m31 * m40 * m52
1264 - m04 * m15 * m23 * m31 * m42 * m50
1265 - m04 * m15 * m23 * m32 * m40 * m51
1266 + m04 * m15 * m23 * m32 * m41 * m50
1267 - m05 * m10 * m21 * m32 * m43 * m54
1268 + m05 * m10 * m21 * m32 * m44 * m53
1269 + m05 * m10 * m21 * m33 * m42 * m54
1270 - m05 * m10 * m21 * m33 * m44 * m52
1271 - m05 * m10 * m21 * m34 * m42 * m53
1272 + m05 * m10 * m21 * m34 * m43 * m52
1273 + m05 * m10 * m22 * m31 * m43 * m54
1274 - m05 * m10 * m22 * m31 * m44 * m53
1275 - m05 * m10 * m22 * m33 * m41 * m54
1276 + m05 * m10 * m22 * m33 * m44 * m51
1277 + m05 * m10 * m22 * m34 * m41 * m53
1278 - m05 * m10 * m22 * m34 * m43 * m51
1279 - m05 * m10 * m23 * m31 * m42 * m54
1280 + m05 * m10 * m23 * m31 * m44 * m52
1281 + m05 * m10 * m23 * m32 * m41 * m54
1282 - m05 * m10 * m23 * m32 * m44 * m51
1283 - m05 * m10 * m23 * m34 * m41 * m52
1284 + m05 * m10 * m23 * m34 * m42 * m51
1285 + m05 * m10 * m24 * m31 * m42 * m53
1286 - m05 * m10 * m24 * m31 * m43 * m52
1287 - m05 * m10 * m24 * m32 * m41 * m53
1288 + m05 * m10 * m24 * m32 * m43 * m51
1289 + m05 * m10 * m24 * m33 * m41 * m52
1290 - m05 * m10 * m24 * m33 * m42 * m51
1291 + m05 * m11 * m20 * m32 * m43 * m54
1292 - m05 * m11 * m20 * m32 * m44 * m53
1293 - m05 * m11 * m20 * m33 * m42 * m54
1294 + m05 * m11 * m20 * m33 * m44 * m52
1295 + m05 * m11 * m20 * m34 * m42 * m53
1296 - m05 * m11 * m20 * m34 * m43 * m52
1297 - m05 * m11 * m22 * m30 * m43 * m54
1298 + m05 * m11 * m22 * m30 * m44 * m53
1299 + m05 * m11 * m22 * m33 * m40 * m54
1300 - m05 * m11 * m22 * m33 * m44 * m50
1301 - m05 * m11 * m22 * m34 * m40 * m53
1302 + m05 * m11 * m22 * m34 * m43 * m50
1303 + m05 * m11 * m23 * m30 * m42 * m54
1304 - m05 * m11 * m23 * m30 * m44 * m52
1305 - m05 * m11 * m23 * m32 * m40 * m54
1306 + m05 * m11 * m23 * m32 * m44 * m50
1307 + m05 * m11 * m23 * m34 * m40 * m52
1308 - m05 * m11 * m23 * m34 * m42 * m50
1309 - m05 * m11 * m24 * m30 * m42 * m53
1310 + m05 * m11 * m24 * m30 * m43 * m52
1311 + m05 * m11 * m24 * m32 * m40 * m53
1312 - m05 * m11 * m24 * m32 * m43 * m50
1313 - m05 * m11 * m24 * m33 * m40 * m52
1314 + m05 * m11 * m24 * m33 * m42 * m50
1315 - m05 * m12 * m20 * m31 * m43 * m54
1316 + m05 * m12 * m20 * m31 * m44 * m53
1317 + m05 * m12 * m20 * m33 * m41 * m54
1318 - m05 * m12 * m20 * m33 * m44 * m51
1319 - m05 * m12 * m20 * m34 * m41 * m53
1320 + m05 * m12 * m20 * m34 * m43 * m51
1321 + m05 * m12 * m21 * m30 * m43 * m54
1322 - m05 * m12 * m21 * m30 * m44 * m53
1323 - m05 * m12 * m21 * m33 * m40 * m54
1324 + m05 * m12 * m21 * m33 * m44 * m50
1325 + m05 * m12 * m21 * m34 * m40 * m53
1326 - m05 * m12 * m21 * m34 * m43 * m50
1327 - m05 * m12 * m23 * m30 * m41 * m54
1328 + m05 * m12 * m23 * m30 * m44 * m51
1329 + m05 * m12 * m23 * m31 * m40 * m54
1330 - m05 * m12 * m23 * m31 * m44 * m50
1331 - m05 * m12 * m23 * m34 * m40 * m51
1332 + m05 * m12 * m23 * m34 * m41 * m50
1333 + m05 * m12 * m24 * m30 * m41 * m53
1334 - m05 * m12 * m24 * m30 * m43 * m51
1335 - m05 * m12 * m24 * m31 * m40 * m53
1336 + m05 * m12 * m24 * m31 * m43 * m50
1337 + m05 * m12 * m24 * m33 * m40 * m51
1338 - m05 * m12 * m24 * m33 * m41 * m50
1339 + m05 * m13 * m20 * m31 * m42 * m54
1340 - m05 * m13 * m20 * m31 * m44 * m52
1341 - m05 * m13 * m20 * m32 * m41 * m54
1342 + m05 * m13 * m20 * m32 * m44 * m51
1343 + m05 * m13 * m20 * m34 * m41 * m52
1344 - m05 * m13 * m20 * m34 * m42 * m51
1345 - m05 * m13 * m21 * m30 * m42 * m54
1346 + m05 * m13 * m21 * m30 * m44 * m52
1347 + m05 * m13 * m21 * m32 * m40 * m54
1348 - m05 * m13 * m21 * m32 * m44 * m50
1349 - m05 * m13 * m21 * m34 * m40 * m52
1350 + m05 * m13 * m21 * m34 * m42 * m50
1351 + m05 * m13 * m22 * m30 * m41 * m54
1352 - m05 * m13 * m22 * m30 * m44 * m51
1353 - m05 * m13 * m22 * m31 * m40 * m54
1354 + m05 * m13 * m22 * m31 * m44 * m50
1355 + m05 * m13 * m22 * m34 * m40 * m51
1356 - m05 * m13 * m22 * m34 * m41 * m50
1357 - m05 * m13 * m24 * m30 * m41 * m52
1358 + m05 * m13 * m24 * m30 * m42 * m51
1359 + m05 * m13 * m24 * m31 * m40 * m52
1360 - m05 * m13 * m24 * m31 * m42 * m50
1361 - m05 * m13 * m24 * m32 * m40 * m51
1362 + m05 * m13 * m24 * m32 * m41 * m50
1363 - m05 * m14 * m20 * m31 * m42 * m53
1364 + m05 * m14 * m20 * m31 * m43 * m52
1365 + m05 * m14 * m20 * m32 * m41 * m53
1366 - m05 * m14 * m20 * m32 * m43 * m51
1367 - m05 * m14 * m20 * m33 * m41 * m52
1368 + m05 * m14 * m20 * m33 * m42 * m51
1369 + m05 * m14 * m21 * m30 * m42 * m53
1370 - m05 * m14 * m21 * m30 * m43 * m52
1371 - m05 * m14 * m21 * m32 * m40 * m53
1372 + m05 * m14 * m21 * m32 * m43 * m50
1373 + m05 * m14 * m21 * m33 * m40 * m52
1374 - m05 * m14 * m21 * m33 * m42 * m50
1375 - m05 * m14 * m22 * m30 * m41 * m53
1376 + m05 * m14 * m22 * m30 * m43 * m51
1377 + m05 * m14 * m22 * m31 * m40 * m53
1378 - m05 * m14 * m22 * m31 * m43 * m50
1379 - m05 * m14 * m22 * m33 * m40 * m51
1380 + m05 * m14 * m22 * m33 * m41 * m50
1381 + m05 * m14 * m23 * m30 * m41 * m52
1382 - m05 * m14 * m23 * m30 * m42 * m51
1383 - m05 * m14 * m23 * m31 * m40 * m52
1384 + m05 * m14 * m23 * m31 * m42 * m50
1385 + m05 * m14 * m23 * m32 * m40 * m51
1386 - m05 * m14 * m23 * m32 * m41 * m50
1387 }
1388 _ =>
1389 {
1391 let l1 = l - 1;
1393 let mut res: Vec<T> = vm.to_vec(); (0 .. l)
1396 .map(|i| {
1397 if !vm[i].is_zero() {
1399 for r in 0 .. l1 {
1400 let d = r * l1;
1401 let s = (r + 1) * l;
1402 for c in 0 .. l1 {
1403 res[d + c] = vm[s + c + usize::from(c >= i)];
1404 }
1405 }
1406 let v = vm[i] * determinant(&res, l1);
1407 if (i % 2) == 0 { v } else { -v }
1408 } else {
1409 T::zero()
1410 }
1411 })
1412 .sum::<T>()
1413 }
1414 }
1415}
1416
1417fn check_diag<T: Float>(m: &Array2<T>) -> usize {
1418 let d = m.raw_dim();
1419 assert!(d[0] == d[1]);
1420
1421 for i in 0 .. d[0] {
1422 if m[(i, i)].is_zero() {
1423 return 0;
1424 }
1425 }
1426
1427 d[0]
1428}
1429
1430fn linv<T: Float>(l: &Array2<T>, n: usize) -> Array2<T> {
1431 let mut m: Array2<T> = Array2::zeros((n, n));
1432
1433 for i in 0 .. n {
1434 m[(i, i)] = l[(i, i)].recip();
1435
1436 for j in 0 .. i {
1437 for k in j .. i {
1438 m[(i, j)] = m[(i, j)] + l[(i, k)] * m[(k, j)];
1439 }
1440
1441 m[(i, j)] = -m[(i, j)] * m[(i, i)];
1442 }
1443 }
1444
1445 m
1446}
1447
1448fn uinv<T: Float>(u: &Array2<T>, n: usize) -> Array2<T> {
1449 linv(&u.t().to_owned(), n).t().to_owned()
1450}
1451
1452#[cfg(test)]
1453mod test_inverse {
1454 use ndarray::prelude::*;
1455 use crate::Inverse;
1456 use num_traits::Float;
1457 use std::fmt::Debug;
1458
1459 const EPS: f64 = 1e-8;
1460
1461 fn to_vec<T: Float>(m: &Array2<T>) -> Vec<T> {
1462 m.indexed_iter().map(|(_, &e)| e).collect()
1463 }
1464
1465 fn compare_vecs<T: Float + Debug>(v1: &Array2<T>, v2: &Array2<T>, epsilon: T) -> bool {
1467
1468 to_vec(v1)
1469 .into_iter()
1470 .zip(to_vec(v2))
1471 .map(|(a, b)| Float::abs(a.abs() - b.abs()) < epsilon)
1472 .all(|b| b)
1473 }
1474
1475 #[test]
1476 fn test_1x1() {
1477 let a: Array2<f64> = array![[2.0]];
1478 let inv = a.inv().expect("Sods Law");
1479 let back = inv.inv();
1480 assert!(compare_vecs(&a, &back.unwrap(), EPS));
1481 let expected = array![[0.5]];
1482 assert!(compare_vecs(&inv, &expected, EPS));
1483 }
1484
1485 #[test]
1486 fn test_2x2() {
1487 let a: Array2<f64> = array![[1.0, 2.0], [3.0, 4.0]];
1488 let inv = a.inv().expect("Sods Law");
1489 let back = inv.inv();
1490 assert!(compare_vecs(&a, &back.unwrap(), EPS));
1491 let expected = array![[-2.0, 1.0], [1.5, -0.5]];
1492 assert!(compare_vecs(&inv, &expected, EPS));
1493 }
1494
1495 #[test]
1496 fn test_3x3() {
1497 let a: Array2<f64> = array![[1.0, 0.0, 3.0], [2.0, 1.0, 6.0], [1.0, 0.0, 9.0]];
1498 let inv = a.inv().expect("Sods Law");
1499 let back = inv.inv();
1500 assert!(compare_vecs(&a, &back.unwrap(), EPS));
1501 let expected = array![[1.5, 0.0, -0.5], [-2.0, 1.0, -0.0], [-0.16666667, 0.0, 0.16666667]];
1502 assert!(compare_vecs(&inv, &expected, EPS));
1503 }
1504
1505 #[test]
1506 fn test_4x4() {
1507 let a: Array2<f64> = array![
1508 [-68.0, 68.0, -16.0, 4.0],
1509 [-36.0, 35.0, -9.0, 3.0],
1510 [48.0, -47.0, 11.0, -3.0],
1511 [64.0, -64.0, 16.0, -4.0]
1512 ];
1513 let inv = a.inv().expect("Sods Law");
1514 let back = inv.inv();
1515 eprintln!("{:?}", to_vec(&inv));
1516 assert!(compare_vecs(&a, &back.unwrap(), EPS));
1518 let expected = array![
1519 [1.25, 0.5, 1.5, 0.5],
1520 [1.5, 0.5, 1.5, 0.75],
1521 [1.5, 0.5, 0.5, 1.5],
1522 [2.0, 2.0, 2.0, 1.75]
1523 ];
1524 assert!(compare_vecs(&inv, &expected, EPS));
1525 }
1526
1527 #[test]
1528 fn test_5x5() {
1529 let a: Array2<f64> = array![
1530 [10.0, 1.0, 7.0, 1.0, 5.0],
1531 [2.0, 4.0, 8.0, 3.0, 2.0],
1532 [5.0, 1.0, 2.0, 9.0, 10.0],
1533 [6.0, 9.0, 9.0, 7.0, 3.0],
1534 [1.0, 8.0, 8.0, 10.0, 5.0],
1535 ];
1536 let inv = a.inv().expect("Sods Law");
1537 let back = inv.inv().unwrap();
1538 assert!(compare_vecs(&a, &back, EPS));
1541 let expected: Array2<f64> = array![
1542 [-11.98, 15.64, 9.32, 10.34, -19.12],
1543 [ 33.62, -44.16, -26.08, -28.46, 53.28],
1544 [ -9.36, 12.48, 7.24, 7.88, -14.84],
1545 [-37.2 , 48.6 , 28.8 , 31.6 , -58.8 ],
1546 [ 37.98, -49.64, -29.32, -32.34, 60.12]
1547 ];
1548 assert!(compare_vecs(&inv, &expected, EPS));
1549 }
1550
1551 #[test]
1552 fn test_5x5b() {
1553 let a: Array2<f64> = array![
1554 [8.12658, 1.02214, 2.67754, 7.00746, 0.16869],
1555 [3.31581, 0.63792, 4.28436, 5.61229, 1.19785],
1556 [8.13028, 6.48437, 6.91838, 2.68860, 6.60971],
1557 [2.85578, 4.40362, 1.13265, 8.32311, 4.41366],
1558 [7.93379, 4.09682, 4.86444, 6.67282, 3.57098]
1559 ];
1560 let inv = a.inv().expect("Sods Law");
1561 let back = inv.inv().unwrap();
1562 assert!(compare_vecs(&a, &back, EPS));
1564 let expected: Array2<f64> = array![
1565 [ 0.50074759, -0.02510542, 0.33877844, 0.09275141, -0.7569348 ],
1566 [-1.58848174, -0.64356647, -1.32943575, -0.57849057, 3.46664012],
1567 [-0.49917931, 0.12378294, -0.30452856, -0.25789124, 0.86447499],
1568 [-0.0513424 , 0.06348965, -0.1164256 , 0.0585164 , 0.12430141],
1569 [ 1.48578933, 0.50685455, 1.40491123, 0.69956389, -3.42524083]
1570 ];
1571 assert!(compare_vecs(&inv, &expected, EPS));
1572 }
1573
1574 #[test]
1575 fn test_6x6() {
1576 let a: Array2<f64> = array![
1577 [1.0, 1.0, 3.0, 4.0, 9.0, 3.0],
1578 [10.0, 10.0, 1.0, 2.0, 2.0, 5.0],
1579 [2.0, 9.0, 6.0, 10.0, 10.0, 9.0],
1580 [10.0, 9.0, 9.0, 7.0, 3.0, 6.0],
1581 [7.0, 6.0, 6.0, 2.0, 9.0, 5.0],
1582 [3.0, 8.0, 1.0, 4.0, 1.0, 5.0]
1583 ];
1584 let inv = a.inv().expect("Sods Law");
1585 eprintln!("{:?}", to_vec(&inv));
1586 let back = inv.inv();
1587 assert!(compare_vecs(&a, &back.unwrap(), EPS));
1589 let expected = array![
1590 [-0.53881418, 0.57793399, 0.6653423 , -0.08374083, -0.16962103,
1591 -1.18215159],
1592 [ 2.16075795, -1.52750611, -2.44070905, 0.44132029, 0.32457213,
1593 3.77017115],
1594 [ 0.21454768, -0.41503667, -0.39425428, 0.14792176, 0.19743276,
1595 0.62102689],
1596 [ 0.70018337, -0.24052567, -0.525978 , 0.20354523, -0.21179707,
1597 0.73471883],
1598 [ 0.85055012, -0.47157702, -0.82793399, 0.1106357 , 0.1146088 ,
1599 1.20415648],
1600 [-3.90709046, 2.46699267, 4.17114914, -0.87041565, -0.31051345,
1601 -6.07579462]
1602 ];
1603 assert!(compare_vecs(&inv, &expected, EPS));
1604 }
1605
1606 #[test]
1607 fn test_7x7() {
1608 let a:Array2<f64> = array![
1609 [4.3552 , 6.25851, 4.12662, 1.93708, 0.21272, 3.25683, 6.53326],
1610 [4.24746, 1.84137, 6.71904, 0.59754, 3.5806 , 3.63597, 5.347 ],
1611 [2.30479, 1.70591, 3.05354, 1.82188, 5.27839, 7.9166 , 2.04607],
1612 [2.40158, 6.38524, 7.90296, 4.69683, 6.63801, 7.32958, 1.45936],
1613 [0.42456, 6.47456, 1.55398, 8.28979, 4.20987, 0.90401, 4.94587],
1614 [5.78903, 1.92032, 6.20261, 5.78543, 1.94331, 8.25178, 7.47273],
1615 [1.44797, 7.41157, 7.69495, 8.90113, 3.05983, 0.41582, 6.42932]];
1616 let inv = a.inv().expect("Sods Law");
1617 let back = inv.inv();
1618 assert!(compare_vecs(&a, &back.unwrap(), EPS));
1619 let expected = array![[-0.0819729 , 0.30947774, -0.90286793, 0.49545542,
1620 0.47316173, 0.30651669, -0.71946285],
1621 [ 0.15235561, -0.07559205, -0.0070691 , 0.06757486,
1622 0.01026404, -0.08019532, -0.01972633],
1623 [-0.03064041, -0.00625996, 0.07810921, -0.01780696,
1624 -0.19716224, -0.03302802, 0.20558502],
1625 [-0.11358415, -0.02543206, -0.24035635, 0.1260263 ,
1626 0.1346972 , 0.16555344, -0.111583 ],
1627 [-0.1016236 , 0.21480081, -0.0698384 , 0.05521421,
1628 0.20158013, -0.05049775, -0.16205801],
1629 [ 0.0760218 , -0.20124886, 0.34308958, -0.12448475,
1630 -0.18988991, -0.02734616, 0.18705112],
1631 [ 0.08020185, -0.02906765, 0.46181332, -0.36087425,
1632 -0.15255756, -0.14045526, 0.31376606]];
1633 assert!(compare_vecs(&inv, &expected, EPS));
1634 }
1635
1636 #[cfg(not(debug_assertions))]
1637 use std::time::{Instant, Duration};
1638
1639 #[test]
1640 #[cfg(not(debug_assertions))]
1641 fn bench_6x6() {
1644 let a: Array2<f64> = array![
1645 [1.0, 1.0, 3.0, 4.0, 9.0, 3.0],
1646 [10.0, 10.0, 1.0, 2.0, 2.0, 5.0],
1647 [2.0, 9.0, 6.0, 10.0, 10.0, 9.0],
1648 [10.0, 9.0, 9.0, 7.0, 3.0, 6.0],
1649 [7.0, 6.0, 6.0, 2.0, 9.0, 5.0],
1650 [3.0, 8.0, 1.0, 4.0, 1.0, 5.0]
1651 ];
1652
1653 let start = Instant::now();
1654 let mut _q = a.inv();
1655 for _ in 0..100000 {
1656 _q = a.inv();
1657 }
1658 assert!(start.elapsed() < Duration::new(2, 0));
1660 }
1661
1662 #[test]
1663 fn test_diag() {
1664 let a: Array2<f64> = array![
1665 [1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
1666 [0.0, 10.0, 0.0, 0.0, 0.0, 0.0],
1667 [0.0, 0.0, 6.0, 0.0, 0.0, 0.0],
1668 [0.0, 0.0, 0.0, 7.0, 0.0, 0.0],
1669 [0.0, 0.0, 0.0, 0.0, 9.0, 0.0],
1670 [0.0, 0.0, 0.0, 0.0, 0.0, 5.0]
1671 ];
1672 let inv = a.inv_diag().expect("Sods Law");
1673 let back = inv.inv_diag();
1674 assert!(compare_vecs(&a, &back.unwrap(), EPS));
1675 let expected = a.inv().expect("Sods Law");
1676 assert!(compare_vecs(&inv, &expected, EPS));
1677 }
1678
1679 #[test]
1680 fn test_cholesky() {
1681 let a = array![[25.0, 15.0, -5.0],
1682 [15.0, 18.0, 0.0],
1683 [-5.0, 0.0, 11.0]];
1684
1685 let res = a.cholesky();
1686
1687 let expected = array![[5.0, 0.0, 0.0],
1688 [3.0, 3.0, 0.0],
1689 [-1.0, 1.0, 3.]];
1690
1691 assert!(compare_vecs(&res, &expected, EPS));
1692
1693 let a = array![[18.0, 22.0, 54.0, 42.0],
1694 [22.0, 70.0, 86.0, 62.0],
1695 [54.0, 86.0, 174.0, 134.0],
1696 [42.0, 62.0, 134.0, 106.0]];
1697
1698 let res = a.cholesky();
1699
1700 let expected = array![[4.242640687119285, 0.0, 0.0, 0.0],
1701 [5.185449728701349, 6.565905201197403, 0.0, 0.0],
1702 [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0.0],
1703 [9.899494936611667, 1.624553864213788, 1.8497110052313648, 1.3926212476456026]];
1704
1705 assert!(compare_vecs(&res, &expected, EPS));
1706
1707 let a = array![[2.0, 0.1], [0.1, 3.0]];
1708
1709 let res = a.cholesky().t().to_owned();
1710
1711 let expected = array![[1.41421356, 0.07071068], [0.0, 1.73060683]];
1712
1713 assert!(compare_vecs(&res, &expected, EPS));
1714
1715 let expected_sq = array![[2.005, 0.12237238], [0.12237238, 2.995]];
1716
1717 assert!(compare_vecs(&res.dot(&res.t()), &expected_sq, EPS));
1718 }
1719
1720 #[test]
1721 fn test_triangular() {
1722 let a = array![[4.242640687119285, 0.0, 0.0, 0.0],
1723 [5.185449728701349, 6.565905201197403, 0.0, 0.0],
1724 [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0.0],
1725 [9.899494936611667, 1.624553864213788, 1.8497110052313648, 1.3926212476456026]];
1726
1727 assert!(compare_vecs(&a.inv().unwrap(), &a.inv_lt().unwrap(), EPS));
1728 assert!(compare_vecs(&a.t().to_owned().inv().unwrap(), &a.t().to_owned().inv_ut().unwrap(), EPS));
1729 }
1730}