linalg/char_poly.hpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #pragma once | ||
| 2 | |||
| 3 | #include <vector> | ||
| 4 | #include <bitset> | ||
| 5 | #include <cassert> | ||
| 6 | |||
| 7 | namespace wala { | ||
| 8 | |||
| 9 | // Compute the characteristic polynomial of a square matrix A over some field. | ||
| 10 | // Not numerically stable at all. | ||
| 11 | // Takes argument by value, use std::move if you can. | ||
| 12 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 28 times.
|
28 | template <typename num> std::vector<num> charPoly(std::vector<std::vector<num>> A) { |
| 13 | 28 | int N = int(A.size()); | |
| 14 |
1/1✓ Branch 4 → 5 taken 28 times.
|
28 | std::vector<num> res; res.reserve(N+1); |
| 15 |
1/1✓ Branch 5 → 6 taken 28 times.
|
28 | res.push_back(num(1)); |
| 16 |
2/2✓ Branch 44 → 7 taken 9590 times.
✓ Branch 44 → 45 taken 28 times.
|
9618 | for (int i = 0, deg = 0; i < N; i++) { |
| 17 | 9590 | auto& Ai = A[i]; | |
| 18 | |||
| 19 | 9590 | int c = i+1; | |
| 20 |
4/4✓ Branch 8 → 9 taken 384437 times.
✓ Branch 8 → 11 taken 1696 times.
✓ Branch 9 → 10 taken 376543 times.
✓ Branch 9 → 11 taken 7894 times.
|
386133 | while (c < N && Ai[c] == num(0)) c++; |
| 21 |
2/2✓ Branch 11 → 12 taken 1696 times.
✓ Branch 11 → 23 taken 7894 times.
|
9590 | if (c == N) { |
| 22 |
1/1✓ Branch 12 → 13 taken 1696 times.
|
1696 | res.resize(i+2, num(0)); |
| 23 |
2/2✓ Branch 21 → 14 taken 458172 times.
✓ Branch 21 → 22 taken 1696 times.
|
459868 | for (int x = deg; x >= 0; x--) { |
| 24 | 458172 | num v = res[x]; | |
| 25 |
2/2✓ Branch 19 → 15 taken 1104371 times.
✓ Branch 19 → 20 taken 458172 times.
|
1562543 | for (int y = x+1, z = i; z >= deg; z--, y++) { |
| 26 |
2/2✓ Branch 15 → 16 taken 494484 times.
✓ Branch 15 → 17 taken 609887 times.
|
2208742 | res[y] -= v * Ai[z]; |
| 27 | } | ||
| 28 | } | ||
| 29 | 1696 | deg = i+1; | |
| 30 | 1696 | continue; | |
| 31 | 1696 | } | |
| 32 | |||
| 33 | 7894 | num vc = Ai[c]; | |
| 34 | 7894 | num ivc = inv(vc); | |
| 35 | |||
| 36 | 7894 | Ai[c] = Ai[i+1]; | |
| 37 | 7894 | Ai[i+1] = 0; | |
| 38 | |||
| 39 | 7894 | std::swap(A[i+1], A[c]); | |
| 40 | 7894 | auto& Ai1 = A[i+1]; | |
| 41 |
2/2✓ Branch 26 → 25 taken 3271757 times.
✓ Branch 26 → 39 taken 7894 times.
|
3279651 | for (int k = deg; k < N; k++) { |
| 42 | 3271757 | Ai1[k] *= vc; | |
| 43 | } | ||
| 44 | |||
| 45 |
2/2✓ Branch 39 → 27 taken 1995001 times.
✓ Branch 39 → 41 taken 7894 times.
|
2002895 | for (int k = i+1; k < N; k++) { |
| 46 | 1995001 | auto& Ak = A[k]; | |
| 47 | { | ||
| 48 | 1995001 | auto& x = Ak[i+1]; | |
| 49 | 1995001 | auto& y = Ak[c]; | |
| 50 | 1995001 | num tmp = y; | |
| 51 | 1995001 | y = x; | |
| 52 | 1995001 | x = tmp * ivc; | |
| 53 | } | ||
| 54 | { | ||
| 55 | 1995001 | num v = Ak[i+1]; | |
| 56 |
2/2✓ Branch 32 → 28 taken 894607571 times.
✓ Branch 32 → 33 taken 1995001 times.
|
896602572 | for (int j = deg; j < N; j++) { |
| 57 |
2/2✓ Branch 28 → 29 taken 446296308 times.
✓ Branch 28 → 30 taken 448311263 times.
|
1789215142 | Ak[j] -= v * Ai[j]; |
| 58 | } | ||
| 59 | } | ||
| 60 |
2/2✓ Branch 33 → 34 taken 1987107 times.
✓ Branch 33 → 38 taken 7894 times.
|
1995001 | if (k > i+1) { |
| 61 | 1987107 | num v = Ai[k]; | |
| 62 |
2/2✓ Branch 37 → 35 taken 891335814 times.
✓ Branch 37 → 38 taken 1987107 times.
|
893322921 | for (int j = deg; j < N; j++) { |
| 63 | 891335814 | Ai1[j] += v * Ak[j]; | |
| 64 | } | ||
| 65 | } | ||
| 66 | } | ||
| 67 | |||
| 68 |
2/2✓ Branch 41 → 40 taken 1276756 times.
✓ Branch 41 → 42 taken 7894 times.
|
1284650 | for (int k = deg; k <= i; k++) { |
| 69 | 1276756 | Ai1[k+1] += Ai[k]; | |
| 70 | } | ||
| 71 | } | ||
| 72 | 28 | reverse(res.begin(), res.end()); | |
| 73 | 28 | return res; | |
| 74 | } | ||
| 75 | |||
| 76 | // Compute the characteristic polynomial of a square matrix A over F2. | ||
| 77 | // Takes argument by value, use std::move if you can. | ||
| 78 | // Note that MAXS must be at least N+1 | ||
| 79 | ✗ | template <std::size_t MAXS> std::bitset<MAXS> charPoly(std::vector<std::bitset<MAXS>> A) { | |
| 80 | using bs = std::bitset<MAXS>; | ||
| 81 | ✗ | int N = int(A.size()); | |
| 82 | ✗ | assert(MAXS >= N+1); | |
| 83 | ✗ | bs ans; ans[0] = 1; | |
| 84 | ✗ | int deg = 0; | |
| 85 | ✗ | for (int i = 0; i < N; i++) { | |
| 86 | { | ||
| 87 | ✗ | int j = int(A[i]._Find_next(i)); | |
| 88 | ✗ | if (j >= N) { | |
| 89 | ✗ | bs nans; | |
| 90 | ✗ | for (; deg <= i; ans <<= 1, deg++) { | |
| 91 | ✗ | if (A[i][deg]) nans ^= ans; | |
| 92 | } | ||
| 93 | ✗ | ans ^= nans; | |
| 94 | ✗ | continue; | |
| 95 | } | ||
| 96 | ✗ | if (j != i+1) { | |
| 97 | ✗ | swap(A[j], A[i+1]); | |
| 98 | ✗ | for (auto& a : A) { | |
| 99 | ✗ | bool tmp = a[j]; | |
| 100 | ✗ | a[j] = a[i+1]; | |
| 101 | ✗ | a[i+1] = tmp; | |
| 102 | } | ||
| 103 | } | ||
| 104 | } | ||
| 105 | ✗ | assert(A[i][i+1]); | |
| 106 | ✗ | bs msk = A[i]; msk.flip(i+1); | |
| 107 | ✗ | for (int k = 0; k < N; k++) { | |
| 108 | ✗ | if (msk[k]) A[i+1] ^= A[k]; | |
| 109 | } | ||
| 110 | ✗ | for (auto& a : A) { | |
| 111 | ✗ | if (a[i+1]) a ^= msk; | |
| 112 | } | ||
| 113 | } | ||
| 114 | ✗ | return ans; | |
| 115 | } | ||
| 116 | |||
| 117 | } // namespace wala | ||
| 118 |