GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 96.2% 128 / 0 / 133
Functions: 100.0% 18 / 0 / 18
Branches: 98.6% 69 / 40 / 110
fft/engines/real.hpp
Line Branch Exec Source
1 #pragma once
2
3 #include <algorithm>
4 #include <cassert>
5 #include <span>
6 #include <utility>
7 #include <vector>
8
9 #include "fft/core.hpp"
10 #include "fft/engine.hpp"
11
12 namespace wala::fft::engines {
13
14 // Convolve real (floating point) values by packing into complex numbers with
15 // a'[t] = a[2t] + i * a[2t+1]
16 // We use conjugate symmetry to untangle/retangle the two.
17 // TODO: Add type bounds?
18 template <typename dbl = double> struct real {
19 using value_type = dbl;
20 static constexpr bool commutative = true;
21 using cnum = cplx<dbl>;
22 using core = fft_core<cnum>;
23
2/2
✓ Branch 67 → 68 taken 2 times.
✓ Branch 84 → 85 taken 6 times.
 768 struct transformed {
24 vector<cnum> v;
25 435 int size() const { return 2 * sz(v); }
26 };
27 using product = transformed;
28 // Precision is caller-managed for this engine (see add), so scale is untracked.
29 static constexpr int unit_scale = 0;
30 template <int A = 0> using transformed_t = transformed;
31 template <int K = 0> using product_t = product;
32
33 428 static int packed_size(int n) { return std::max(n / 2, 1); }
34 4 static void pack(std::span<const dbl> a, std::span<cnum> c) {
35
4/4
✓ Branch 6 → 7 taken 95 times.
✓ Branch 6 → 12 taken 99 times.
✓ Branch 18 → 3 taken 194 times.
✓ Branch 18 → 19 taken 4 times.
 198 for (int i = 0; i < sz(a); i++) (i & 1 ? c[i/2].y : c[i/2].x) = a[i];
36 4 }
37 // Spectrum of the real (odd = false) or imaginary (odd = true) part of the packed
38 // sequence at bitrev entry t, by conjugate symmetry with the entry of w^{-k}.
39 17176 static cnum part(const transformed& f, int t, bool odd) {
40 17176 cnum g = conj(f.v[core::conj_index(t)]);
41
4/4
✓ Branch 12 → 13 taken 8588 times.
✓ Branch 12 → 24 taken 8588 times.
✓ Branch 36 → 37 taken 8588 times.
✓ Branch 36 → 38 taken 8588 times.
 17176 return odd ? (f.v[t] - g) * cnum(0, dbl(-0.5)) : (f.v[t] + g) * cnum(dbl(0.5));
42 }
43 // Given the spectra (s0, s1) of a real sequence x at w_{2mo}^q and w_{2mo}^{q+mo},
44 // the packed-transform entry of x at packed size mo: the even/odd interleaves of x
45 // have spectra (s0 +- s1)/2 (the odd one twisted by w_{2mo}^{-q}).
46 4150 static cnum retangle(cnum s0, cnum s1, int mo, int q) {
47 4150 cnum s = (s0 + s1) * cnum(dbl(0.5));
48 4150 cnum d = (s0 - s1) * cnum(dbl(0.5)) * core::inv_rt[mo + q];
49 4150 return s + cnum(-d.y, d.x);
50 }
51
52 252 static transformed transform(std::span<const dbl> a, int n) {
53 252 assert(sz(a) <= 2 * n);
54 252 transformed r;
55
1/1
✓ Branch 15 → 16 taken 252 times.
 252 r.v.assign(packed_size(n), cnum(0));
56
2/2
✓ Branch 42 → 18 taken 6706 times.
✓ Branch 42 → 43 taken 252 times.
 6958 for (int i = 0; i < sz(a); i++) {
57
2/2
✓ Branch 18 → 19 taken 123 times.
✓ Branch 18 → 20 taken 6583 times.
 6706 int j = i < n ? i : i - n;
58
2/2
✓ Branch 23 → 24 taken 3281 times.
✓ Branch 23 → 33 taken 3425 times.
 6706 ((j & 1) ? r.v[j/2].y : r.v[j/2].x) += a[i];
59 }
60
1/1
✓ Branch 46 → 47 taken 252 times.
 252 core::forward(std::span<cnum>(r.v));
61 252 return r;
62 }
63 51 static void extend_to(transformed& t, int m, std::span<const dbl> coeffs) {
64 51 assert(!(m & (m-1)) && sz(coeffs) <= 2 * m);
65
2/2
✓ Branch 11 → 12 taken 43 times.
✓ Branch 11 → 149 taken 8 times.
 141 if (t.size() >= m) return;
66
3/3
✓ Branch 17 → 18 taken 39 times.
✓ Branch 17 → 46 taken 4 times.
✓ Branch 19 → 20 taken 39 times.
 164 if (t.size() == 0) { t = transform(coeffs, m); return; }
67
1/1
✓ Branch 48 → 49 taken 4 times.
 4 auto buf = buffer_pool<cnum>::get((sz(coeffs) + 1) / 2);
68 16 std::fill(buf.span().begin(), buf.span().end(), cnum(0));
69 4 pack(coeffs, buf.span());
70
2/2
✓ Branch 146 → 98 taken 9 times.
✓ Branch 146 → 147 taken 4 times.
 26 while (t.size() < m) {
71 9 int s = sz(t.v);
72
1/1
✓ Branch 103 → 104 taken 9 times.
 9 t.v.resize(2 * s);
73 // packed coeffs past 2s are zero: they didn't fit in the transform we're a prefix of
74
1/1
✓ Branch 134 → 135 taken 9 times.
 18 core::extend(
75 9 std::span<cnum>(t.v),
76 36 std::span<const cnum>(buf.span()).first(size_t(min(sz(buf.span()), 2 * s)))
77 );
78 }
79 4 }
80 20 static transformed downsample(const transformed& t, int n, bool odd) { return half(t, n, odd); }
81 12 static transformed upsample(const transformed& t, int n, bool odd) { return unhalf(t, n, odd); }
82 // A(-x) negates the odd (imaginary-slot) coefficients, i.e. conjugates the packed
83 // sequence; the transform of a conjugated sequence is the conjugate at w^(-k).
84 4 static transformed negate_arg(const transformed& t, int n) {
85 4 int m = packed_size(n);
86 4 assert(n >= 2 && sz(t.v) >= m);
87
1/1
✓ Branch 16 → 35 taken 4 times.
 4 transformed r; r.v.resize(m);
88
2/2
✓ Branch 35 → 17 taken 43 times.
✓ Branch 35 → 36 taken 4 times.
 47 for (int j = 0; j < m; j++) r.v[j] = conj(t.v[core::conj_index(j)]);
89 4 return r;
90 }
91 // b[k] = a[-k mod n] packs to b'[t] = Re a'[-t] + i Im a'[-t-1], so the packed spectrum
92 // at w = w_m^k is X(w^-1) + i w^-1 Y(w^-1) in terms of the real/imaginary-part spectra.
93 15 static transformed of_reverse(const transformed& t, int n) {
94 15 int m = packed_size(n);
95 15 assert(sz(t.v) >= m);
96
1/1
✓ Branch 15 → 16 taken 15 times.
 15 transformed r; r.v.resize(m);
97
2/2
✓ Branch 16 → 17 taken 6 times.
✓ Branch 16 → 29 taken 9 times.
 15 if (m == 1) { r.v[0] = t.v[0]; return r; }
98
1/1
✓ Branch 29 → 63 taken 9 times.
 9 core::init(m);
99
2/2
✓ Branch 63 → 30 taken 126 times.
✓ Branch 63 → 64 taken 9 times.
 135 for (int j = 0; j < m; j++) {
100 126 int k = core::brev(j, m);
101
2/2
✓ Branch 32 → 33 taken 63 times.
✓ Branch 32 → 38 taken 63 times.
 126 cnum iw = k < m / 2 ? core::inv_rt[m / 2 + k] : -core::inv_rt[k];
102 126 int cj = core::conj_index(j);
103 126 cnum x = part(t, cj, false), y = part(t, cj, true) * iw;
104 126 r.v[j] = x + cnum(-y.y, y.x);
105 }
106 return r;
107 }
108 20 static transformed half(const transformed& f, int n, bool odd) {
109 40 assert(n >= 2 && f.size() >= 2 * n);
110 20 int mo = n / 2;
111 20 core::init(2 * mo);
112
1/1
✓ Branch 17 → 31 taken 20 times.
 20 transformed r; r.v.resize(mo);
113
2/2
✓ Branch 31 → 18 taken 310 times.
✓ Branch 31 → 32 taken 20 times.
 330 for (int u = 0; u < mo; u++) {
114 310 r.v[u] = retangle(part(f, 2*u, odd), part(f, 2*u+1, odd), mo, core::brev(u, mo));
115 }
116 20 return r;
117 }
118 // Inverse of half: the size-n packed transform of the input spread as evens or odds.
119 12 static transformed unhalf(const transformed& f, int n, bool odd) {
120 24 assert(n >= 2 && 2 * f.size() >= n);
121 12 int mo = packed_size(n);
122
1/1
✓ Branch 17 → 18 taken 12 times.
 12 transformed r; r.v.resize(mo);
123
1/2
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 28 taken 12 times.
 12 if (mo == 1) {
124 ✗ r.v[0] = part(f, 0, false);
125 } else {
126 // The spread-as-evens sequence packs to the plain (purely real) input sequence,
127 // so untangle its packed transform into the full spectrum:
128 // A(w) = E(w^2) + w O(w^2) at w = +-w_{2mi}^brev(u, mi).
129 12 int mi = mo / 2;
130
1/1
✓ Branch 28 → 60 taken 12 times.
 12 core::init(2 * mi);
131
2/2
✓ Branch 60 → 29 taken 164 times.
✓ Branch 60 → 61 taken 12 times.
 176 for (int u = 0; u < mi; u++) {
132 164 cnum e = part(f, u, false);
133 164 cnum o = part(f, u, true) * core::rt[mi + core::brev(u, mi)];
134 164 r.v[2*u+0] = e + o;
135 164 r.v[2*u+1] = e - o;
136 }
137 }
138 // spreading as odds instead packs into the imaginary slots, scaling the spectrum by i
139
4/4
✓ Branch 61 → 62 taken 6 times.
✓ Branch 61 → 110 taken 6 times.
✓ Branch 89 → 67 taken 164 times.
✓ Branch 89 → 90 taken 6 times.
 188 if (odd) for (auto& z : r.v) z = cnum(-z.y, z.x);
140 12 return r;
141 }
142 132 static product mul(const transformed& a, const transformed& b, int n) {
143 132 int m = packed_size(n);
144 396 assert(a.size() >= n && b.size() >= n);
145 132 core::init(2 * m);
146
1/1
✓ Branch 23 → 67 taken 132 times.
 132 product p; p.v.resize(m);
147
2/2
✓ Branch 67 → 24 taken 3686 times.
✓ Branch 67 → 68 taken 132 times.
 3818 for (int t = 0; t < m; t++) {
148 3686 int k = core::brev(t, m);
149 3686 cnum w = core::rt[m + k];
150 3686 cnum xa = part(a, t, false), ya = part(a, t, true);
151 3686 cnum xb = part(b, t, false), yb = part(b, t, true);
152 // full spectra at w_{2m}^k and w_{2m}^{k+m} = -w_{2m}^k
153 3686 cnum p0 = (xa + w * ya) * (xb + w * yb);
154 3686 cnum p1 = (xa - w * ya) * (xb - w * yb);
155 3686 p.v[t] = retangle(p0, p1, m, k);
156 }
157 132 return p;
158 }
159 5 static product sq(const transformed& a, int n) { return mul(a, a, n); }
160 8 static product mul2(
161 const transformed& a1, const transformed& b1,
162 const transformed& a2, const transformed& b2,
163 int n
164 ) {
165 8 int m = packed_size(n);
166 40 assert(a1.size() >= n && b1.size() >= n && a2.size() >= n && b2.size() >= n);
167 8 core::init(2 * m);
168
1/1
✓ Branch 35 → 103 taken 8 times.
 8 product p; p.v.resize(m);
169
2/2
✓ Branch 103 → 36 taken 154 times.
✓ Branch 103 → 104 taken 8 times.
 162 for (int t = 0; t < m; t++) {
170 154 int k = core::brev(t, m);
171 154 cnum w = core::rt[m + k];
172 154 cnum xa1 = part(a1, t, false), ya1 = part(a1, t, true);
173 154 cnum xb1 = part(b1, t, false), yb1 = part(b1, t, true);
174 154 cnum xa2 = part(a2, t, false), ya2 = part(a2, t, true);
175 154 cnum xb2 = part(b2, t, false), yb2 = part(b2, t, true);
176 154 cnum p0 = (xa1 + w * ya1) * (xb1 + w * yb1) + (xa2 + w * ya2) * (xb2 + w * yb2);
177 154 cnum p1 = (xa1 - w * ya1) * (xb1 - w * yb1) + (xa2 - w * ya2) * (xb2 - w * yb2);
178 154 p.v[t] = retangle(p0, p1, m, k);
179 }
180 8 return p;
181 }
182 ✗ static product add(product&& a, const product& b) {
183 ✗ assert(a.size() == b.size());
184 ✗ for (int i = 0; i < sz(a.v); i++) a.v[i] = a.v[i] + b.v[i];
185 ✗ return std::move(a);
186 }
187 // sum_k finish(p)[k] * b[k] = Re sum_t p'[t] conj(b'[t]) = Re (1/m) sum_j P[j] conj(B[j])
188 5 static dbl dot(const product& p, const transformed& t, int n) {
189 5 int m = packed_size(n);
190 5 assert(sz(p.v) >= m && sz(t.v) >= m);
191 dbl r = 0;
192
2/2
✓ Branch 41 → 14 taken 44 times.
✓ Branch 41 → 42 taken 5 times.
 49 for (int j = 0; j < m; j++) r += p.v[j].x * t.v[j].x + p.v[j].y * t.v[j].y;
193 5 return r / dbl(m);
194 }
195 152 template <typename Op = assign_op> static void finish(product&& p, std::span<dbl> out, Op op = {}) {
196 152 int m = sz(p.v);
197 152 assert(sz(out) <= 2 * m);
198
2/2
void wala::fft::engines::real<double>::finish<wala::fft::detail::cut_op<double, wala::fft::assign_op> >(wala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, wala::fft::detail::cut_op<double, wala::fft::assign_op>):
✓ Branch 12 → 13 taken 7 times.
void wala::fft::engines::real<double>::finish<wala::fft::assign_op>(wala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, wala::fft::assign_op):
✓ Branch 12 → 13 taken 145 times.
 152 core::inverse(std::span<cnum>(p.v));
199 152 dbl d = dbl(1) / dbl(m);
200
8/8
void wala::fft::engines::real<double>::finish<wala::fft::detail::cut_op<double, wala::fft::assign_op> >(wala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, wala::fft::detail::cut_op<double, wala::fft::assign_op>):
✓ Branch 15 → 16 taken 117 times.
✓ Branch 15 → 23 taken 118 times.
✓ Branch 34 → 15 taken 235 times.
✓ Branch 34 → 35 taken 7 times.
void wala::fft::engines::real<double>::finish<wala::fft::assign_op>(wala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, wala::fft::assign_op):
✓ Branch 15 → 16 taken 3706 times.
✓ Branch 15 → 23 taken 3718 times.
✓ Branch 35 → 15 taken 7424 times.
✓ Branch 35 → 36 taken 145 times.
 7811 for (int i = 0; i < sz(out); i++) op(out[i], (i & 1 ? p.v[i/2].y : p.v[i/2].x) * d);
201 152 }
202 };
203
204 } // namespace wala::fft::engines
205