GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 96.0% 96 / 0 / 100
Functions: 100.0% 14 / 0 / 14
Branches: 100.0% 50 / 30 / 80
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 ecnerwala::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 3 times.
 487 struct transformed {
24 vector<cnum> v;
25 305 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 262 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 13104 static cnum part(const transformed& f, int t, bool odd) {
40 13104 cnum g = conj(f.v[core::conj_index(t)]);
41
4/4
✓ Branch 12 → 13 taken 6552 times.
✓ Branch 12 → 24 taken 6552 times.
✓ Branch 36 → 37 taken 6552 times.
✓ Branch 36 → 38 taken 6552 times.
 13104 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 3195 static cnum retangle(cnum s0, cnum s1, int mo, int q) {
47 3195 cnum s = (s0 + s1) * cnum(dbl(0.5));
48 3195 cnum d = (s0 - s1) * cnum(dbl(0.5)) * core::inv_rt[mo + q];
49 3195 return s + cnum(-d.y, d.x);
50 }
51
52 171 static transformed transform(std::span<const dbl> a, int n) {
53 171 assert(sz(a) <= 2 * n);
54 171 transformed r;
55
1/1
✓ Branch 15 → 16 taken 171 times.
 171 r.v.assign(packed_size(n), cnum(0));
56
2/2
✓ Branch 42 → 18 taken 4678 times.
✓ Branch 42 → 43 taken 171 times.
 4849 for (int i = 0; i < sz(a); i++) {
57
2/2
✓ Branch 18 → 19 taken 26 times.
✓ Branch 18 → 20 taken 4652 times.
 4678 int j = i < n ? i : i - n;
58
2/2
✓ Branch 23 → 24 taken 2283 times.
✓ Branch 23 → 33 taken 2395 times.
 4678 ((j & 1) ? r.v[j/2].y : r.v[j/2].x) += a[i];
59 }
60
1/1
✓ Branch 46 → 47 taken 171 times.
 171 core::forward(std::span<cnum>(r.v));
61 171 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 8 static transformed downsample(const transformed& t, int n, bool odd) { return half(t, n, odd); }
81 // A(-x) negates the odd (imaginary-slot) coefficients, i.e. conjugates the packed
82 // sequence; the transform of a conjugated sequence is the conjugate at w^(-k).
83 4 static transformed negate_arg(const transformed& t, int n) {
84 4 int m = packed_size(n);
85 4 assert(n >= 2 && sz(t.v) >= m);
86
1/1
✓ Branch 16 → 35 taken 4 times.
 4 transformed r; r.v.resize(m);
87
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)]);
88 4 return r;
89 }
90 8 static transformed half(const transformed& f, int n, bool odd) {
91 16 assert(n >= 2 && f.size() >= 2 * n);
92 8 int mo = n / 2;
93 8 core::init(2 * mo);
94
1/1
✓ Branch 17 → 31 taken 8 times.
 8 transformed r; r.v.resize(mo);
95
2/2
✓ Branch 31 → 18 taken 146 times.
✓ Branch 31 → 32 taken 8 times.
 154 for (int u = 0; u < mo; u++) {
96 146 r.v[u] = retangle(part(f, 2*u, odd), part(f, 2*u+1, odd), mo, core::brev(u, mo));
97 }
98 8 return r;
99 }
100 79 static product mul(const transformed& a, const transformed& b, int n) {
101 79 int m = packed_size(n);
102 237 assert(a.size() >= n && b.size() >= n);
103 79 core::init(2 * m);
104
1/1
✓ Branch 23 → 67 taken 79 times.
 79 product p; p.v.resize(m);
105
2/2
✓ Branch 67 → 24 taken 2895 times.
✓ Branch 67 → 68 taken 79 times.
 2974 for (int t = 0; t < m; t++) {
106 2895 int k = core::brev(t, m);
107 2895 cnum w = core::rt[m + k];
108 2895 cnum xa = part(a, t, false), ya = part(a, t, true);
109 2895 cnum xb = part(b, t, false), yb = part(b, t, true);
110 // full spectra at w_{2m}^k and w_{2m}^{k+m} = -w_{2m}^k
111 2895 cnum p0 = (xa + w * ya) * (xb + w * yb);
112 2895 cnum p1 = (xa - w * ya) * (xb - w * yb);
113 2895 p.v[t] = retangle(p0, p1, m, k);
114 }
115 79 return p;
116 }
117 5 static product sq(const transformed& a, int n) { return mul(a, a, n); }
118 8 static product mul2(
119 const transformed& a1, const transformed& b1,
120 const transformed& a2, const transformed& b2,
121 int n
122 ) {
123 8 int m = packed_size(n);
124 40 assert(a1.size() >= n && b1.size() >= n && a2.size() >= n && b2.size() >= n);
125 8 core::init(2 * m);
126
1/1
✓ Branch 35 → 103 taken 8 times.
 8 product p; p.v.resize(m);
127
2/2
✓ Branch 103 → 36 taken 154 times.
✓ Branch 103 → 104 taken 8 times.
 162 for (int t = 0; t < m; t++) {
128 154 int k = core::brev(t, m);
129 154 cnum w = core::rt[m + k];
130 154 cnum xa1 = part(a1, t, false), ya1 = part(a1, t, true);
131 154 cnum xb1 = part(b1, t, false), yb1 = part(b1, t, true);
132 154 cnum xa2 = part(a2, t, false), ya2 = part(a2, t, true);
133 154 cnum xb2 = part(b2, t, false), yb2 = part(b2, t, true);
134 154 cnum p0 = (xa1 + w * ya1) * (xb1 + w * yb1) + (xa2 + w * ya2) * (xb2 + w * yb2);
135 154 cnum p1 = (xa1 - w * ya1) * (xb1 - w * yb1) + (xa2 - w * ya2) * (xb2 - w * yb2);
136 154 p.v[t] = retangle(p0, p1, m, k);
137 }
138 8 return p;
139 }
140 ✗ static product add(product&& a, const product& b) {
141 ✗ assert(a.size() == b.size());
142 ✗ for (int i = 0; i < sz(a.v); i++) a.v[i] = a.v[i] + b.v[i];
143 ✗ return std::move(a);
144 }
145 93 template <typename Op = assign_op> static void finish(product&& p, std::span<dbl> out, Op op = {}) {
146 93 int m = sz(p.v);
147 93 assert(sz(out) <= 2 * m);
148
2/2
void ecnerwala::fft::engines::real<double>::finish<ecnerwala::fft::detail::cut_op<double, ecnerwala::fft::assign_op> >(ecnerwala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, ecnerwala::fft::detail::cut_op<double, ecnerwala::fft::assign_op>):
✓ Branch 12 → 13 taken 7 times.
void ecnerwala::fft::engines::real<double>::finish<ecnerwala::fft::assign_op>(ecnerwala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, ecnerwala::fft::assign_op):
✓ Branch 12 → 13 taken 86 times.
 93 core::inverse(std::span<cnum>(p.v));
149 93 dbl d = dbl(1) / dbl(m);
150
8/8
void ecnerwala::fft::engines::real<double>::finish<ecnerwala::fft::detail::cut_op<double, ecnerwala::fft::assign_op> >(ecnerwala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, ecnerwala::fft::detail::cut_op<double, ecnerwala::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 ecnerwala::fft::engines::real<double>::finish<ecnerwala::fft::assign_op>(ecnerwala::fft::engines::real<double>::transformed&&, std::span<double, 18446744073709551615ul>, ecnerwala::fft::assign_op):
✓ Branch 15 → 16 taken 2755 times.
✓ Branch 15 → 23 taken 2763 times.
✓ Branch 35 → 15 taken 5518 times.
✓ Branch 35 → 36 taken 86 times.
 5846 for (int i = 0; i < sz(out); i++) op(out[i], (i & 1 ? p.v[i/2].y : p.v[i/2].x) * d);
151 93 }
152 };
153
154 /* namespace ecnerwala::fft::engines */ }
155