GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 68.6% 72 / 0 / 105
Functions: 100.0% 2 / 0 / 2
Branches: 82.7% 115 / 8 / 147
fft/ap_sampled_poly.hpp
Line Branch Exec Source
1 #pragma once
2
3 #include <cassert>
4 #include <span>
5 #include <vector>
6
7 #include "fft/multiply.hpp"
8
9 namespace ecnerwala {
10
11 // A polynomial represented by its values evaluated at an Arithmetic Progression (AP).
12 // TODO: The AP is always assumed to be 0..length-1; store an explicit offset/gap instead?
13 // Maybe not, this is just more convenient.
14 template <fft::engine E>
15 431 struct ap_sampled_poly : public std::vector<typename E::value_type> {
16 using T = typename E::value_type;
17
5/5
✓ Branch 5 → 6 taken 64 times.
✓ Branch 6 → 7 taken 8 times.
✓ Branch 24 → 25 taken 7 times.
✓ Branch 71 → 72 taken 8 times.
✓ Branch 479 → 480 taken 49 times.
 168 using std::vector<T>::vector;
18
19 48298492 int len() const {
20
6/12
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 44 times.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 15 times.
✗ Branch 52 → 53 not taken.
✓ Branch 52 → 54 taken 6587187 times.
✗ Branch 56 → 57 not taken.
✓ Branch 56 → 58 taken 29 times.
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 65 taken 11703295 times.
✗ Branch 85 → 86 not taken.
✓ Branch 85 → 87 taken 11703295 times.
 30007923 return int(this->size());
21 }
22 ✗ int degree() const {
23 ✗ return len() - 1;
24 }
25
26 ✗ ap_sampled_poly& operator += (const ap_sampled_poly& o) {
27 ✗ assert(len() == o.len());
28 ✗ for (int i = 0; i < int(o.size()); i++) {
29 ✗ (*this)[i] += o[i];
30 }
31 ✗ return *this;
32 }
33 ✗ friend ap_sampled_poly operator + (const ap_sampled_poly& a, const ap_sampled_poly& b) {
34 ✗ assert(a.size() == b.size());
35 ✗ ap_sampled_poly r(a.size());
36 ✗ for (int i = 0; i < r.len(); i++) {
37 ✗ r[i] = a[i] + b[i];
38 }
39 ✗ return r;
40 }
41 ✗ ap_sampled_poly& operator -= (const ap_sampled_poly& o) {
42 ✗ assert(len() == o.len());
43 ✗ for (int i = 0; i < int(o.size()); i++) {
44 ✗ (*this)[i] -= o[i];
45 }
46 ✗ return *this;
47 }
48 ✗ friend ap_sampled_poly operator - (const ap_sampled_poly& a, const ap_sampled_poly& b) {
49 ✗ assert(a.size() == b.size());
50 ✗ ap_sampled_poly r(a.size());
51 ✗ for (int i = 0; i < r.len(); i++) {
52 ✗ r[i] = a[i] - b[i];
53 }
54 ✗ return r;
55 }
56
57 220 T eval_at(T k) {
58
3/4
✓ Branch 3 → 4 taken 220 times.
✗ Branch 3 → 21 not taken.
✓ Branch 12 → 13 taken 53 times.
✓ Branch 12 → 21 taken 167 times.
 660 if (0 <= int(k) && int(k) < len()) {
59 53 return (*this)[int(k)];
60 } else {
61 // Just do the lagrange interpolation
62
1/1
✓ Branch 24 → 25 taken 167 times.
 167 std::vector<T> terms(*this);
63 {
64 // Inverse factorial terms
65 167 T v = 1;
66
2/2
✓ Branch 45 → 30 taken 2192 times.
✓ Branch 45 → 46 taken 167 times.
 4718 for (int i = 1; i <= len(); i++) v *= T(i);
67 334 v = inv(v);
68
2/2
✓ Branch 93 → 62 taken 2192 times.
✓ Branch 93 → 94 taken 167 times.
 2526 for (int i = len()-1; i >= 0; i--) {
69 4384 v *= T(i+1);
70 2192 terms[i] *= v;
71
2/2
✓ Branch 73 → 74 taken 1029 times.
✓ Branch 73 → 79 taken 1163 times.
 5413 terms[len()-1-i] *= (i & 1) ? -v : v;
72 }
73 }
74 {
75 // Prefix terms
76 167 T v = 1;
77
2/2
✓ Branch 122 → 100 taken 2192 times.
✓ Branch 122 → 123 taken 167 times.
 4718 for (int i = 0; i < len(); i++) {
78 2192 terms[i] *= v;
79 4384 v *= T(k - i);
80 }
81 }
82 {
83 // Suffix terms
84 167 T v = 1;
85
2/2
✓ Branch 151 → 136 taken 2192 times.
✓ Branch 151 → 152 taken 167 times.
 2526 for (int i = len() - 1; i >= 0; i--) {
86 2192 terms[i] *= v;
87 4384 v *= T(k - i);
88 }
89 }
90 167 T res = 0;
91
2/2
✓ Branch 169 → 158 taken 2192 times.
✓ Branch 169 → 170 taken 167 times.
 2526 for (int i = 0; i < len(); i++) res += terms[i];
92 167 return res;
93 167 }
94 }
95
96 126 ap_sampled_poly eval_range(T k, int osz) {
97
3/3
✓ Branch 2 → 3 taken 26 times.
✓ Branch 2 → 5 taken 44 times.
✓ Branch 2 → 8 taken 56 times.
 126 if (osz == 0) {
98 34 return ap_sampled_poly(osz);
99 }
100
3/4
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 10 taken 44 times.
✓ Branch 15 → 16 taken 7 times.
✓ Branch 15 → 27 taken 49 times.
 156 if (len() == 0) {
101 21 return ap_sampled_poly(osz, T(0));
102 }
103
104 // Check for overlaps. We're checking in linear time to avoid unpacking T, but it should be plenty fast.
105 // If the field is very very small and we wrap around several times, our runtime can be bad...
106 // but then something has already gone wrong, why are you evaluating so many points???
107
4/4
✓ Branch 37 → 11 taken 15377168 times.
✓ Branch 37 → 38 taken 29 times.
✓ Branch 388 → 36 taken 693 times.
✓ Branch 388 → 389 taken 49 times.
 15377988 for (int i = -(len() - 1); i <= osz - 1; i++) {
108
7/8
✓ Branch 11 → 12 taken 9541032 times.
✓ Branch 11 → 13 taken 5836136 times.
✓ Branch 14 → 15 taken 7454549 times.
✓ Branch 14 → 16 taken 7922619 times.
✓ Branch 17 → 18 taken 15 times.
✓ Branch 17 → 36 taken 15377153 times.
✗ Branch 53 → 54 not taken.
✓ Branch 53 → 387 taken 693 times.
 46133583 if (k+i == T(0)) {
109 // everything from [i,i+len()-1) of the output is a match
110
1/2
✓ Branch 18 → 19 taken 15 times.
✗ Branch 55 → 56 not taken.
 15 ap_sampled_poly res; res.reserve(osz);
111
1/2
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 15 times.
 15 int lo = std::max(0, i);
112
1/1
✓ Branch 21 → 22 taken 15 times.
 15 int hi = std::min(i+len(), osz);
113 {
114
1/2
✓ Branch 21 → 22 taken 15 times.
✗ Branch 76 → 77 not taken.
 15 auto pref = eval_range(k, lo);
115
1/2
✓ Branch 22 → 23 taken 15 times.
✗ Branch 97 → 98 not taken.
 15 res.insert(res.end(), pref.begin(), pref.end());
116 }
117
1/2
✓ Branch 24 → 25 taken 15 times.
✗ Branch 194 → 195 not taken.
 15 res.insert(res.end(), this->begin() + (lo - i), this->begin() + (hi - i));
118 {
119
5/7
✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 15 times.
✓ Branch 28 → 29 taken 10 times.
✓ Branch 28 → 30 taken 5 times.
✓ Branch 31 → 32 taken 15 times.
✓ Branch 32 → 33 taken 15 times.
✗ Branch 272 → 273 not taken.
 45 auto suff = eval_range(k + hi, osz - hi);
120
1/2
✓ Branch 32 → 33 taken 15 times.
✗ Branch 294 → 295 not taken.
 15 res.insert(res.end(), suff.begin(), suff.end());
121 15 }
122 15 return res;
123 15 }
124 }
125
126
1/1
✓ Branch 393 → 394 taken 49 times.
 78 std::vector<T> inps(*this);
127 {
128 // Inverse factorial terms
129 78 T v = 1;
130
5/6
✗ Branch 44 → 45 not taken.
✓ Branch 44 → 46 taken 6587216 times.
✓ Branch 46 → 40 taken 6587187 times.
✓ Branch 46 → 47 taken 29 times.
✓ Branch 414 → 399 taken 371 times.
✓ Branch 414 → 415 taken 49 times.
 6588056 for (int i = 1; i <= len(); i++) v *= T(i);
131 127 v = inv(v);
132
4/4
✓ Branch 55 → 49 taken 6587187 times.
✓ Branch 55 → 56 taken 29 times.
✓ Branch 462 → 431 taken 371 times.
✓ Branch 462 → 463 taken 49 times.
 6587685 for (int i = len()-1; i >= 0; i--) {
133
2/2
✓ Branch 49 → 50 taken 3293587 times.
✓ Branch 49 → 52 taken 3293600 times.
 6587929 v *= T(i+1);
134
2/2
✓ Branch 49 → 50 taken 3293587 times.
✓ Branch 49 → 52 taken 3293600 times.
 6587558 inps[i] *= v;
135
4/4
✓ Branch 49 → 50 taken 3293587 times.
✓ Branch 49 → 52 taken 3293600 times.
✓ Branch 442 → 443 taken 168 times.
✓ Branch 442 → 448 taken 203 times.
 16468871 inps[len()-1-i] *= (i & 1) ? -v : v;
136 }
137 }
138
2/2
✓ Branch 58 → 59 taken 29 times.
✓ Branch 473 → 474 taken 49 times.
 127 std::vector<T> inv_offsets(len() + osz - 1);
139
2/3
✓ Branch 59 → 60 taken 29 times.
✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 29 times.
 127 ap_sampled_poly results(osz);
140 {
141
1/2
✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 29 times.
 78 T v = 1;
142
4/4
✓ Branch 74 → 63 taken 11703295 times.
✓ Branch 74 → 75 taken 29 times.
✓ Branch 520 → 493 taken 693 times.
✓ Branch 520 → 521 taken 49 times.
 11704115 for (int i = - (len() - 1); i <= osz - 1; i++) {
143
2/2
✓ Branch 65 → 66 taken 6587158 times.
✓ Branch 65 → 67 taken 5116137 times.
 11704681 inv_offsets[i + (len() - 1)] = v;
144
4/4
✓ Branch 68 → 69 taken 5116137 times.
✓ Branch 68 → 70 taken 6587158 times.
✓ Branch 71 → 72 taken 5116137 times.
✓ Branch 71 → 73 taken 6587158 times.
 23407976 v *= k + i;
145
4/4
✓ Branch 71 → 72 taken 5116137 times.
✓ Branch 71 → 73 taken 6587158 times.
✓ Branch 514 → 515 taken 371 times.
✓ Branch 514 → 519 taken 322 times.
 11703988 if (i >= 0) results[i] = v;
146 }
147 // Assert there's no overlap
148 127 assert(v != T(0));
149 127 v = inv(v);
150
5/6
✗ Branch 90 → 91 not taken.
✓ Branch 90 → 92 taken 11703324 times.
✓ Branch 92 → 79 taken 11703295 times.
✓ Branch 92 → 93 taken 29 times.
✓ Branch 586 → 539 taken 693 times.
✓ Branch 586 → 587 taken 49 times.
 11704808 for (int i = osz - 1; i >= -(len() - 1); i--) {
151
2/2
✓ Branch 79 → 80 taken 6587158 times.
✓ Branch 79 → 81 taken 5116137 times.
 11704681 inv_offsets[i + (len() - 1)] *= v;
152
3/4
✓ Branch 82 → 83 taken 5116137 times.
✓ Branch 82 → 84 taken 6587158 times.
✗ Branch 85 → 86 not taken.
✓ Branch 85 → 87 taken 11703295 times.
 23407976 v *= k + i;
153
4/4
✓ Branch 87 → 88 taken 5116137 times.
✓ Branch 87 → 89 taken 6587158 times.
✓ Branch 567 → 568 taken 371 times.
✓ Branch 567 → 578 taken 322 times.
 11704681 if (i + (len() - 1) <= osz - 1) {
154 5116879 results[i + (len() - 1)] *= v;
155 }
156 }
157 }
158
4/6
✗ Branch 95 → 96 not taken.
✓ Branch 95 → 97 taken 29 times.
✓ Branch 97 → 98 taken 29 times.
✗ Branch 98 → 99 not taken.
✓ Branch 98 → 100 taken 29 times.
✓ Branch 593 → 594 taken 49 times.
 78 std::vector<T> prod = fft::middle_product<E>(inv_offsets, inps);
159 78 assert(int(prod.size()) == osz);
160
4/4
✓ Branch 103 → 102 taken 5116137 times.
✓ Branch 103 → 104 taken 29 times.
✓ Branch 605 → 599 taken 371 times.
✓ Branch 605 → 606 taken 49 times.
 5116586 for (int i = 0; i < osz; i++) results[i] *= prod[i];
161 78 return results;
162 225 }
163
164 ✗ void extend_right() {
165 ✗ this->push_back(eval_at(T(len())));
166 }
167 ✗ void extend_left() {
168 ✗ this->insert(this->begin(), eval_at(T(-1)));
169 }
170
171 ✗ [[nodiscard]] ap_sampled_poly prefix_sum_inclusive() const {
172 ✗ ap_sampled_poly r = *this;
173 ✗ r.extend_right();
174 ✗ for (int i = 1; i < r.len(); i++) r[i] += r[i-1];
175 ✗ return r;
176 }
177 };
178
179 /* namespace ecnerwala */ }
180