fft/poly.hpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #pragma once | ||
| 2 | |||
| 3 | #include <algorithm> | ||
| 4 | #include <cassert> | ||
| 5 | #include <concepts> | ||
| 6 | #include <cstddef> | ||
| 7 | #include <initializer_list> | ||
| 8 | #include <span> | ||
| 9 | #include <type_traits> | ||
| 10 | #include <utility> | ||
| 11 | #include <vector> | ||
| 12 | |||
| 13 | #include "fft/series_core.hpp" | ||
| 14 | |||
| 15 | namespace wala::poly { | ||
| 16 | |||
| 17 | // polynomial class | ||
| 18 | // As above, we represent polynomials by a series::exact containing the coefficients in reverse order. | ||
| 19 | // This representation should be internal-only: | ||
| 20 | // all accesses/constructors use the logical order though: P[k] = [x^k] P. | ||
| 21 | // To use the representation, use rev_series() / from_rev_series() | ||
| 22 |
11/11✓ Branch 9 → 10 taken 38 times.
✓ Branch 11 → 12 taken 38 times.
✓ Branch 213 → 214 taken 6 times.
✓ Branch 308 → 309 taken 1 time.
✓ Branch 318 → 319 taken 6 times.
✓ Branch 334 → 335 taken 6 times.
✓ Branch 351 → 352 taken 1 time.
✓ Branch 455 → 456 taken 6 times.
✓ Branch 509 → 510 taken 1 time.
✓ Branch 626 → 627 taken 1 time.
✓ Branch 777 → 778 taken 1 time.
|
763 | template <fft::engine E> struct vec { |
| 23 | using T = typename E::value_type; | ||
| 24 | using engine_t = E; | ||
| 25 | series::exact<E> c; | ||
| 26 | |||
| 27 | ✗ | vec() = default; | |
| 28 | // zero polynomial with `len` coefficient slots | ||
| 29 | 49 | explicit vec(int len) : c(size_t(len), T{}) {} | |
| 30 | // coefficient (x^0-first) order | ||
| 31 | 1970659 | vec(std::initializer_list<T> coeffs) : c(std::rbegin(coeffs), std::rend(coeffs)) {} | |
| 32 | 40 | explicit vec(std::span<const T> coeffs) : c(coeffs.rbegin(), coeffs.rend()) {} | |
| 33 | |||
| 34 | 104 | const series::exact<E>& rev_series() const { return c; } | |
| 35 | 83 | static vec from_rev_series(series::exact<E> s) { | |
| 36 | 83 | vec r; | |
| 37 | 121 | r.c = std::move(s); | |
| 38 | 36 | return r; | |
| 39 | } | ||
| 40 | |||
| 41 | // This should rarely be used | ||
| 42 | 1 | series::trunc<E> unrev_series(int n) const { | |
| 43 | 3 | series::trunc<E> r(size_t(n), T{}); | |
| 44 |
1/1✓ Branch 38 → 39 taken 1 time.
|
7 | std::copy(begin(), begin() + std::min(n, len()), r.begin()); |
| 45 | 1 | return r; | |
| 46 | } | ||
| 47 | |||
| 48 | // logical (coefficient) order | ||
| 49 | 52 | auto begin() { return c.rbegin(); } | |
| 50 | 52 | auto end() { return c.rend(); } | |
| 51 | 2 | auto begin() const { return c.rbegin(); } | |
| 52 | ✗ | auto end() const { return c.rend(); } | |
| 53 | |||
| 54 |
1/2✗ Branch 8 → 9 not taken.
✓ Branch 8 → 10 taken 11 times.
|
66182077 | int len() const { return c.len(); } |
| 55 | ✗ | int degree() const { return len() - 1; } | |
| 56 |
5/8✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 7210203 times.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 7210241 times.
✓ Branch 16 → 17 taken 7210241 times.
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 644060 times.
✓ Branch 22 → 23 taken 644060 times.
|
22276269 | T& operator[](int i) { return c[len() - 1 - i]; } |
| 57 | 296 | const T& operator[](int i) const { return c[len() - 1 - i]; } | |
| 58 | 1 | T leading() const { return c.front(); } | |
| 59 | // multiply by x^k: appends the new zero constant terms to the storage | ||
| 60 | 1 | void shift(int k = 1) { | |
| 61 |
2/3✓ Branch 12 → 13 taken 1 time.
✗ Branch 12 → 68 not taken.
✓ Branch 36 → 37 taken 1 time.
|
7 | if (len() > 0) c.insert(c.end(), size_t(k), T(0)); |
| 62 | 1 | } | |
| 63 | // grow (zero-filled leading coefficients) or shrink to n coefficients | ||
| 64 | 1 | void resize(int n) { | |
| 65 |
2/3✓ Branch 12 → 13 taken 1 time.
✗ Branch 12 → 78 not taken.
✓ Branch 46 → 47 taken 1 time.
|
8 | if (n >= len()) c.insert(c.begin(), size_t(n - len()), T(0)); |
| 66 | ✗ | else c.erase(c.begin(), c.begin() + (len() - n)); | |
| 67 | 1 | } | |
| 68 | |||
| 69 | 87 | T operator()(const T& x) const { | |
| 70 | 87 | T r{}; | |
| 71 |
2/2✓ Branch 36 → 15 taken 2218 times.
✓ Branch 36 → 37 taken 87 times.
|
2479 | for (const T& v : c) r = r * x + v; |
| 72 | 87 | return r; | |
| 73 | } | ||
| 74 | |||
| 75 | 1 | vec& operator+=(const vec& o) { | |
| 76 |
1/2✗ Branch 22 → 23 not taken.
✓ Branch 22 → 44 taken 1 time.
|
3 | if (o.len() > len()) resize(o.len()); |
| 77 |
2/2✓ Branch 55 → 35 taken 23 times.
✓ Branch 55 → 56 taken 1 time.
|
25 | for (int i = 0; i < o.len(); i++) (*this)[i] += o[i]; |
| 78 | 1 | return *this; | |
| 79 | } | ||
| 80 |
1/1✓ Branch 308 → 309 taken 1 time.
|
2 | friend vec operator+(vec a, const vec& b) { a += b; return a; } |
| 81 | 1 | vec& operator-=(const vec& o) { | |
| 82 |
1/2✓ Branch 22 → 23 taken 1 time.
✗ Branch 22 → 44 not taken.
|
4 | if (o.len() > len()) resize(o.len()); |
| 83 |
2/2✓ Branch 55 → 35 taken 37 times.
✓ Branch 55 → 56 taken 1 time.
|
39 | for (int i = 0; i < o.len(); i++) (*this)[i] -= o[i]; |
| 84 | 1 | return *this; | |
| 85 | } | ||
| 86 |
1/1✓ Branch 351 → 352 taken 1 time.
|
2 | friend vec operator-(vec a, const vec& b) { a -= b; return a; } |
| 87 | 54 | friend bool operator==(const vec& a, const vec& b) { return a.c == b.c; } | |
| 88 | |||
| 89 | ✗ | vec& operator*=(const T& n) { for (T& v : c) v *= n; return *this; } | |
| 90 | ✗ | friend vec operator*(vec a, const T& n) { a *= n; return a; } | |
| 91 | ✗ | friend vec operator*(const T& n, vec a) { a *= n; return a; } | |
| 92 | |||
| 93 | ✗ | vec& operator*=(const vec& o) { return *this = (*this) * o; } | |
| 94 | }; | ||
| 95 | |||
| 96 | // any polynomial representation exposing its reversed coefficient series | ||
| 97 | template <typename P> | ||
| 98 | concept like = requires(const P& p) { | ||
| 99 | typename P::engine_t; | ||
| 100 | { p.len() } -> std::same_as<int>; | ||
| 101 | p.rev_series(); | ||
| 102 | requires series::like<std::remove_cvref_t<decltype(p.rev_series())>>; | ||
| 103 | requires std::remove_cvref_t<decltype(p.rev_series())>::exact_v; | ||
| 104 | }; | ||
| 105 | |||
| 106 | // immutable polynomial carrying the whole-sequence transform of its rev_series | ||
| 107 | template <fft::engine E> | ||
| 108 | 11821616 | struct cached { | |
| 109 | using T = typename E::value_type; | ||
| 110 | using engine_t = E; | ||
| 111 | |||
| 112 | 3940466 | cached() = default; | |
| 113 | // moving coefficients in or out is free: implicit on rvalues, explicit copy otherwise | ||
| 114 | 1970233 | cached(vec<E>&& p) : c(std::move(p.c)) {} | |
| 115 | 3 | explicit cached(const vec<E>& p) : c(p.c) {} | |
| 116 | 36 | operator vec<E>() && { return vec<E>::from_rev_series(std::move(c)); } | |
| 117 | |||
| 118 | 4622697 | const series::cached_exact<E>& rev_series() const { return c; } | |
| 119 | 1970197 | static cached from_rev_series(series::cached_exact<E> s) { | |
| 120 | 1970197 | cached r; | |
| 121 | 1970197 | r.c = std::move(s); | |
| 122 | 206 | return r; | |
| 123 | } | ||
| 124 | |||
| 125 |
1/2✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 2651880 times.
|
2652480 | int len() const { return c.len(); } |
| 126 | ✗ | int degree() const { return len() - 1; } | |
| 127 | 117 | const T& operator[](int i) const { return c[len() - 1 - i]; } | |
| 128 | ✗ | T leading() const { return c[0]; } | |
| 129 | |||
| 130 | 4 | T operator()(const T& x) const { | |
| 131 | 4 | T r{}; | |
| 132 |
2/2✓ Branch 34 → 13 taken 158 times.
✓ Branch 34 → 35 taken 4 times.
|
170 | for (const T& v : c) r = r * x + v; |
| 133 | 4 | return r; | |
| 134 | } | ||
| 135 | |||
| 136 | private: | ||
| 137 | series::cached_exact<E> c; | ||
| 138 | }; | ||
| 139 | |||
| 140 | // rev(a*b) = rev(a)*rev(b); the series product reuses/adopts transforms | ||
| 141 | template <like A, like B> requires fft::same_engine<A, B> | ||
| 142 | 1326079 | cached<typename A::engine_t> operator*(const A& a, const B& b) { | |
| 143 |
3/3wala::poly::cached<wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > >::engine_t> wala::poly::operator*<wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > >, wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&):
✓ Branch 13 → 14 taken 7 times.
wala::poly::cached<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > >::engine_t> wala::poly::operator*<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > >, wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&):
✓ Branch 13 → 14 taken 1 time.
wala::poly::cached<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > >::engine_t> wala::poly::operator*<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > >, wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > const&):
✓ Branch 13 → 14 taken 131 times.
|
2652297 | return cached<typename A::engine_t>::from_rev_series(a.rev_series() * b.rev_series()); |
| 144 | } | ||
| 145 | template <like A> | ||
| 146 | 2 | cached<typename A::engine_t> square(const A& a) { | |
| 147 |
2/2wala::poly::cached<wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > >::engine_t> wala::poly::square<wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&):
✓ Branch 8 → 9 taken 1 time.
wala::poly::cached<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > >::engine_t> wala::poly::square<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > const&):
✓ Branch 8 → 9 taken 1 time.
|
4 | return cached<typename A::engine_t>::from_rev_series(square(a.rev_series())); |
| 148 | } | ||
| 149 | // rev(a*b + c*d) = rev(a)*rev(b) + rev(c)*rev(d) | ||
| 150 | template <like A, like B, like C, like D> | ||
| 151 | requires fft::same_engine<A, B> && fft::same_engine<A, C> && fft::same_engine<A, D> | ||
| 152 | 644116 | cached<typename A::engine_t> multiply_add2( | |
| 153 | const A& a, const B& b, const C& c, const D& d) { | ||
| 154 | return cached<typename A::engine_t>::from_rev_series( | ||
| 155 |
1/1✓ Branch 23 → 24 taken 65 times.
|
1288427 | multiply_add2(a.rev_series(), b.rev_series(), c.rev_series(), d.rev_series())); |
| 156 | } | ||
| 157 | template <like A, like B> requires fft::same_engine<A, B> | ||
| 158 | 5 | bool operator==(const A& a, const B& b) { | |
| 159 | 15 | return a.rev_series() == b.rev_series(); | |
| 160 | } | ||
| 161 | |||
| 162 | // finite-support linear form | ||
| 163 | // These are one side of the pairing <vec P, series::vec S> = [x^0] P(1/x) S(x). | ||
| 164 | // (Strictly speaking, this is actually <>_d where we take polynomials of degree < d.) | ||
| 165 | // The main point of this wrapper is that if we have <*, S> and want <P *, S>, that's a middle product by P. | ||
| 166 | // | ||
| 167 | // TODO: Should we split it apart into <*, S> and <P, *>? | ||
| 168 | // | ||
| 169 | // Some use cases of this pairing: | ||
| 170 | // <P, 1/(1-ax)> = P(a) | ||
| 171 | // if we represent P as a "polynomial" in the differential operator D (x^k = k! D^k): | ||
| 172 | // <P, e^{aD}> = P(a) | ||
| 173 | template <fft::engine E> | ||
| 174 | 7956844 | struct form { | |
| 175 | using T = typename E::value_type; | ||
| 176 | // coeffs of S in <*, S>; always whole-cached: the kernel transform is | ||
| 177 | // what repeated middle products against the same form reuse | ||
| 178 | series::cached_exact<E> c; | ||
| 179 | |||
| 180 | 2652204 | form() = default; | |
| 181 | ✗ | explicit form(int len) : c(series::exact<E>(size_t(len), T{})) {} | |
| 182 | // We don't provide coefficient-list constructors, to avoid ordering confusion. | ||
| 183 | |||
| 184 | 1327387 | const series::cached_exact<E>& rev_series() const { return c; } | |
| 185 | 2652210 | static form from_rev_series(series::cached_exact<E> s) { | |
| 186 | 2652210 | form r; | |
| 187 | 308 | r.c = std::move(s); | |
| 188 | 288 | return r; | |
| 189 | } | ||
| 190 | 29 | static form from_poly(const vec<E>& p) { return from_rev_series(series::cached_exact<E>(p.rev_series())); } | |
| 191 | |||
| 192 |
2/4✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 11 times.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 11 times.
|
2576555 | int len() const { return c.len(); } |
| 193 | |||
| 194 | // Restrict the form's domain: only valid against exact series of length n | ||
| 195 | 17 | form for_length(int n) const { | |
| 196 |
2/3✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 11 times.
✓ Branch 10 → 11 taken 6 times.
|
23 | auto r = series::exact<E>(c); |
| 197 |
5/6✓ Branch 7 → 8 taken 8 times.
✓ Branch 7 → 10 taken 3 times.
✓ Branch 8 → 9 taken 8 times.
✓ Branch 24 → 25 taken 6 times.
✗ Branch 24 → 91 not taken.
✓ Branch 59 → 60 taken 6 times.
|
59 | if (n >= len()) r.insert(r.begin(), size_t(n - len()), T(0)); |
| 198 |
0/1✗ Branch 138 → 139 not taken.
|
3 | else r.erase(r.begin(), r.begin() + (len() - n)); |
| 199 | 40 | return from_rev_series(std::move(r)); | |
| 200 | 17 | } | |
| 201 | |||
| 202 | // the functional p -> p(z) on polynomials of length up to len (weight z^i on [x^i]) | ||
| 203 | 1 | static form polynomial_evaluation(T z, int len) { | |
| 204 | 3 | series::exact<E> k(size_t(len), T{}); | |
| 205 | 1 | T p = T(1); | |
| 206 |
2/2✓ Branch 24 → 19 taken 40 times.
✓ Branch 24 → 25 taken 1 time.
|
41 | for (int i = 0; i < len; i++) { k[i] = p; p *= z; } |
| 207 | 3 | return from_rev_series(std::move(k)); | |
| 208 | 1 | } | |
| 209 | |||
| 210 | template <like P> | ||
| 211 | 4 | T operator()(const P& p) const { | |
| 212 | 12 | assert(p.len() <= len()); | |
| 213 | 4 | T r{}; | |
| 214 |
4/4wala::modnum<998244353> wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::operator()<wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&) const:
✓ Branch 60 → 32 taken 88 times.
✓ Branch 60 → 61 taken 3 times.
wala::modnum<998244353> wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::operator()<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > const&) const:
✓ Branch 66 → 35 taken 39 times.
✓ Branch 66 → 67 taken 1 time.
|
262 | for (int i = 0; i < p.len(); i++) r += c[i] * p[i]; // weights multiply from the left |
| 215 | 4 | return r; | |
| 216 | } | ||
| 217 | |||
| 218 | // <*, S> -> <q x *, S> | ||
| 219 | template <like P> | ||
| 220 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 2651880 times.
|
2652141 | form composed_with(const P& q) const { |
| 221 | 2652924 | assert(q.len() > 0 && q.len() <= len()); | |
| 222 |
2/2wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > > wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::composed_with<wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&) const:
✓ Branch 48 → 49 taken 1 time.
wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > > wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::composed_with<wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > >(wala::poly::cached<wala::fft::engines::ntt<wala::modnum<998244353> > > const&) const:
✓ Branch 54 → 55 taken 260 times.
|
5304804 | return from_rev_series(middle_product(c, q.rev_series())); |
| 223 | } | ||
| 224 | |||
| 225 | // <P, *> -> <P, s x *> | ||
| 226 | template <series::like S> requires std::same_as<typename S::engine_t, E> | ||
| 227 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 11 times.
|
19 | form composed_with(const S& s) const { |
| 228 | 32 | if constexpr (!S::exact_v) assert(s.len() >= len()); | |
| 229 |
2/2wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > > wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::composed_with<wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, false> >(wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, false> const&) const:
✓ Branch 28 → 29 taken 7 times.
wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > > wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::composed_with<wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, true> >(wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, true> const&) const:
✓ Branch 7 → 8 taken 1 time.
|
20 | series::vec<E, S::exact_v> r = c * s; |
| 230 |
3/3wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > > wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::composed_with<wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, false> >(wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, false> const&) const:
✓ Branch 11 → 12 taken 11 times.
✓ Branch 42 → 43 taken 7 times.
wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > > wala::poly::form<wala::fft::engines::ntt<wala::modnum<998244353> > >::composed_with<wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, true> >(wala::series::vec<wala::fft::engines::ntt<wala::modnum<998244353> >, true> const&) const:
✓ Branch 30 → 31 taken 1 time.
|
27 | r.resize(size_t(len())); |
| 231 | 62 | return from_rev_series(series::exact<E>(std::move(r))); | |
| 232 | 19 | } | |
| 233 | }; | ||
| 234 | |||
| 235 | // ==== the D basis ==== | ||
| 236 | |||
| 237 | // to_D_basis writes P in the iterated-integral basis x^k / k! = I^k 1 (I = D^{-1}, integration from 0), scaling [x^k] by k!: | ||
| 238 | // Q = to_D_basis(P) means P = Q(I) 1. | ||
| 239 | // from_D_basis undoes it. | ||
| 240 | // In this basis D is division by the variable (dropping the constant term), | ||
| 241 | // so the shift operator e^{cD} is multiplication by e^{c/x}, i.e. by e^{cx} mod x^len on rev_series. | ||
| 242 | template <fft::engine E> | ||
| 243 | 70 | vec<E> to_D_basis(vec<E> p) { | |
| 244 | using T = typename E::value_type; | ||
| 245 | 70 | T f = 1; | |
| 246 |
5/6✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 7210241 times.
✓ Branch 9 → 3 taken 7210203 times.
✓ Branch 9 → 10 taken 38 times.
✓ Branch 27 → 7 taken 268 times.
✓ Branch 27 → 28 taken 32 times.
|
7210841 | for (int i = 1; i < p.len(); i++) { |
| 247 | 7210739 | f *= i; | |
| 248 | 7210471 | p[i] *= f; | |
| 249 | } | ||
| 250 | 70 | return p; | |
| 251 | } | ||
| 252 | template <fft::engine E> | ||
| 253 | 70 | vec<E> from_D_basis(vec<E> p) { | |
| 254 | using T = typename E::value_type; | ||
| 255 | 70 | T f = 1; | |
| 256 |
5/6✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 7210241 times.
✓ Branch 9 → 3 taken 7210203 times.
✓ Branch 9 → 10 taken 38 times.
✓ Branch 24 → 7 taken 268 times.
✓ Branch 24 → 25 taken 32 times.
|
7210841 | for (int i = 1; i < p.len(); i++) f *= i; |
| 257 | 102 | f = inv(f); | |
| 258 |
4/4✓ Branch 15 → 12 taken 7210203 times.
✓ Branch 15 → 16 taken 38 times.
✓ Branch 54 → 43 taken 268 times.
✓ Branch 54 → 55 taken 32 times.
|
7210541 | for (int i = p.len() - 1; i > 0; i--) { |
| 259 | 7210471 | p[i] *= f; | |
| 260 | 7210739 | f *= i; | |
| 261 | } | ||
| 262 | 70 | return p; | |
| 263 | } | ||
| 264 | |||
| 265 | // P(x + c) = e^{cD} P; length p.len(). | ||
| 266 | // Polynomials only: [x^k] P(x + c) depends on every higher coefficient of P, so there is no trunc version. | ||
| 267 | template <fft::engine E> | ||
| 268 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 38 times.
|
64 | vec<E> taylor_shift(vec<E> p, typename E::value_type c) { |
| 269 | using T = typename E::value_type; | ||
| 270 | 64 | int n = p.len(); | |
| 271 | // e^{cx} mod x^n: [x^i] = c^i / i! | ||
| 272 |
1/2✓ Branch 5 → 6 taken 38 times.
✗ Branch 5 → 8 not taken.
|
116 | series::trunc<E> e(size_t(n), T{}); |
| 273 |
3/4✓ Branch 5 → 6 taken 38 times.
✗ Branch 5 → 8 not taken.
✓ Branch 23 → 24 taken 22 times.
✓ Branch 23 → 45 taken 4 times.
|
86 | if (n > 0) e[0] = T(1); |
| 274 |
4/4✓ Branch 9 → 7 taken 7210203 times.
✓ Branch 9 → 10 taken 38 times.
✓ Branch 46 → 34 taken 222 times.
✓ Branch 46 → 47 taken 26 times.
|
7210489 | for (int i = 1; i < n; i++) e[i] = e[i-1] * c; |
| 275 | 64 | T f = 1; | |
| 276 |
4/4✓ Branch 15 → 11 taken 7210203 times.
✓ Branch 15 → 16 taken 38 times.
✓ Branch 60 → 52 taken 222 times.
✓ Branch 60 → 61 taken 26 times.
|
7210489 | for (int i = 1; i < n; i++) f *= i; |
| 277 | 90 | f = inv(f); | |
| 278 |
4/4✓ Branch 19 → 18 taken 7210203 times.
✓ Branch 19 → 20 taken 38 times.
✓ Branch 81 → 70 taken 222 times.
✓ Branch 81 → 82 taken 26 times.
|
7210489 | for (int i = n - 1; i > 0; i--) { |
| 279 | 7210425 | e[i] *= f; | |
| 280 | 7210647 | f *= i; | |
| 281 | } | ||
| 282 |
2/2✓ Branch 21 → 22 taken 38 times.
✓ Branch 102 → 103 taken 26 times.
|
298 | return from_D_basis(vec<E>::from_rev_series(series::exact<E>(to_D_basis(std::move(p)).rev_series() * e))); |
| 283 | 64 | } | |
| 284 | |||
| 285 | // ==== multipoint evaluation / interpolation ==== | ||
| 286 | |||
| 287 | // Subproduct tree over points a[0:N] | ||
| 288 | // BFS-order tree, each node holds prod (x - a[i]) as a cached vec. | ||
| 289 | template <fft::engine E> | ||
| 290 | 12 | struct subproduct_tree { | |
| 291 | using T = typename E::value_type; | ||
| 292 | int N; | ||
| 293 | std::vector<cached<E>> nodes; | ||
| 294 | |||
| 295 |
1/1✓ Branch 11 → 12 taken 12 times.
|
32 | explicit subproduct_tree(std::span<const T> pts) : N(sz(pts)), nodes(size_t(2) * N) { |
| 296 | 32 | assert(N > 0); | |
| 297 |
4/4✓ Branch 12 → 5 taken 1325960 times.
✓ Branch 12 → 13 taken 20 times.
✓ Branch 90 → 17 taken 142 times.
✓ Branch 90 → 91 taken 12 times.
|
1326134 | for (int i = 0; i < N; i++) { |
| 298 |
2/2✓ Branch 6 → 7 taken 1325960 times.
✓ Branch 29 → 30 taken 142 times.
|
2652914 | nodes[N + i] = vec<E>{-pts[i], T(1)}; |
| 299 | } | ||
| 300 |
4/4✓ Branch 19 → 14 taken 1325940 times.
✓ Branch 19 → 20 taken 20 times.
✓ Branch 123 → 94 taken 130 times.
✓ Branch 123 → 124 taken 12 times.
|
1326102 | for (int i = N - 1; i > 0; i--) { |
| 301 |
2/2✓ Branch 14 → 15 taken 1325940 times.
✓ Branch 105 → 106 taken 130 times.
|
2652270 | nodes[i] = nodes[2*i] * nodes[2*i+1]; |
| 302 | } | ||
| 303 | 32 | } | |
| 304 | |||
| 305 | // number of points under node i | ||
| 306 | ✗ | int size(int i) const { return nodes[i].len() - 1; } | |
| 307 | // prod (x - z_j) over node i's leaves; length size(i) + 1 | ||
| 308 |
2/4✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 11 times.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 9 times.
|
32 | const cached<E>& prod(int i) const { return nodes[i]; } |
| 309 | |||
| 310 | // Computes, for each i, f(product_{j != i} (1 - a[j] x)). Requires f.len() == N. | ||
| 311 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 20 times.
|
32 | std::vector<T> pushdown(form<E> f) const { |
| 312 | 32 | assert(f.len() == N); | |
| 313 |
1/1✓ Branch 22 → 23 taken 12 times.
|
32 | std::vector<form<E>> down(size_t(2) * N); |
| 314 | 32 | down[1] = std::move(f); | |
| 315 |
4/4✓ Branch 21 → 9 taken 1325940 times.
✓ Branch 21 → 22 taken 20 times.
✓ Branch 96 → 33 taken 130 times.
✓ Branch 96 → 97 taken 12 times.
|
1326102 | for (int i = 1; i < N; i++) { |
| 316 | // the form's kernel transform serves both children's middle products | ||
| 317 |
2/2✓ Branch 9 → 10 taken 1325940 times.
✓ Branch 42 → 43 taken 130 times.
|
1326330 | down[2*i+0] = down[i].composed_with(nodes[2*i+1]); |
| 318 |
2/2✓ Branch 13 → 14 taken 1325940 times.
✓ Branch 65 → 66 taken 130 times.
|
1326330 | down[2*i+1] = down[i].composed_with(nodes[2*i+0]); |
| 319 | 2652270 | down[i] = form<E>{}; // done with the parent; free it early | |
| 320 | } | ||
| 321 |
2/2✓ Branch 22 → 23 taken 20 times.
✓ Branch 105 → 106 taken 12 times.
|
44 | std::vector<T> out(size_t(N), T{}); |
| 322 |
4/4✓ Branch 25 → 24 taken 1325960 times.
✓ Branch 25 → 26 taken 20 times.
✓ Branch 131 → 109 taken 142 times.
✓ Branch 131 → 132 taken 12 times.
|
1326276 | for (int i = 0; i < N; i++) out[i] = down[N + i].rev_series()[0]; |
| 323 | 44 | return out; | |
| 324 | 32 | } | |
| 325 | |||
| 326 | // Compute sum_i leaf_vals[i] prod_{j!=i} (x - a[j]) (transpose of pushdown) | ||
| 327 | 15 | cached<E> combine_up(std::span<const T> leaf_vals) const { | |
| 328 | 15 | assert(sz(leaf_vals) == N); | |
| 329 |
1/1✓ Branch 11 → 12 taken 6 times.
|
15 | std::vector<cached<E>> up(size_t(2) * N); |
| 330 |
4/4✓ Branch 12 → 6 taken 644060 times.
✓ Branch 12 → 13 taken 9 times.
✓ Branch 79 → 14 taken 71 times.
✓ Branch 79 → 80 taken 6 times.
|
644146 | for (int i = 0; i < N; i++) { |
| 331 |
2/2✓ Branch 6 → 7 taken 644060 times.
✓ Branch 21 → 22 taken 71 times.
|
1288475 | up[N + i] = vec<E>{leaf_vals[i]}; |
| 332 | } | ||
| 333 |
4/4✓ Branch 25 → 14 taken 644051 times.
✓ Branch 25 → 26 taken 9 times.
✓ Branch 145 → 83 taken 65 times.
✓ Branch 145 → 146 taken 6 times.
|
644131 | for (int i = N - 1; i > 0; i--) { |
| 334 |
2/2✓ Branch 14 → 15 taken 644051 times.
✓ Branch 102 → 103 taken 65 times.
|
644246 | up[i] = multiply_add2(up[2*i+0], nodes[2*i+1], up[2*i+1], nodes[2*i+0]); |
| 335 | 644246 | up[2*i+0] = cached<E>{}; | |
| 336 | 1288297 | up[2*i+1] = cached<E>{}; | |
| 337 | } | ||
| 338 | 27 | return std::move(up[1]); | |
| 339 | 15 | } | |
| 340 | }; | ||
| 341 | |||
| 342 | template <fft::engine E> | ||
| 343 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 11 times.
|
17 | std::vector<typename E::value_type> multipoint( |
| 344 | const vec<E>& p, | ||
| 345 | std::span<const typename E::value_type> pts | ||
| 346 | ) { | ||
| 347 |
2/4std::__debug::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 16 taken 6 times.
std::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 11 times.
|
17 | if (pts.empty()) return {}; |
| 348 | 17 | int N = sz(pts); | |
| 349 |
1/1✓ Branch 18 → 19 taken 6 times.
|
17 | subproduct_tree<E> tree{pts}; |
| 350 |
3/4std::__debug::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✓ Branch 36 → 37 taken 6 times.
std::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 11 times.
✓ Branch 7 → 8 taken 11 times.
|
35 | auto q = series::trunc<E>(tree.prod(1).rev_series()); |
| 351 |
2/2std::__debug::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✓ Branch 47 → 48 taken 6 times.
std::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✓ Branch 10 → 11 taken 11 times.
|
23 | q.resize(p.len()); // inverse precision must cover the form's window |
| 352 |
6/6std::__debug::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✓ Branch 50 → 51 taken 6 times.
✓ Branch 52 → 53 taken 6 times.
✓ Branch 54 → 55 taken 6 times.
std::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✓ Branch 11 → 12 taken 11 times.
✓ Branch 12 → 13 taken 11 times.
✓ Branch 13 → 14 taken 11 times.
|
40 | form<E> f = form<E>::from_poly(p).composed_with(ps_inv(q)); |
| 353 |
4/4std::__debug::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✓ Branch 78 → 79 taken 6 times.
✓ Branch 81 → 82 taken 6 times.
std::vector<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type, std::allocator<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type> > wala::poly::multipoint<wala::fft::engines::ntt<wala::modnum<998244353> > >(wala::poly::vec<wala::fft::engines::ntt<wala::modnum<998244353> > > const&, std::span<wala::fft::engines::ntt<wala::modnum<998244353> >::value_type const, 18446744073709551615ul>):
✓ Branch 17 → 18 taken 11 times.
✓ Branch 18 → 19 taken 11 times.
|
34 | return tree.pushdown(f.for_length(N)); |
| 354 | 29 | } | |
| 355 | |||
| 356 | template <fft::engine E> | ||
| 357 | 15 | vec<E> interpolate( | |
| 358 | std::span<const typename E::value_type> pts, | ||
| 359 | std::span<const typename E::value_type> vals | ||
| 360 | ) { | ||
| 361 | using T = typename E::value_type; | ||
| 362 | 15 | assert(sz(pts) == sz(vals)); | |
| 363 |
2/4✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 9 times.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 20 taken 6 times.
|
15 | if (pts.empty()) return {}; |
| 364 | 15 | int N = sz(pts); | |
| 365 | using ps = series::trunc<E>; | ||
| 366 |
1/1✓ Branch 22 → 23 taken 6 times.
|
15 | subproduct_tree<E> tree{pts}; |
| 367 |
3/4✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 9 times.
✓ Branch 9 → 10 taken 9 times.
✓ Branch 40 → 41 taken 6 times.
|
33 | auto root = ps(tree.prod(1).rev_series()); |
| 368 |
2/2✓ Branch 10 → 11 taken 9 times.
✓ Branch 41 → 42 taken 6 times.
|
15 | root.shrink(N); |
| 369 | |||
| 370 | // We need to evaluate the derivative of the root at each point | ||
| 371 | 644075 | ps deriv_root = root; | |
| 372 |
4/4✓ Branch 16 → 12 taken 644060 times.
✓ Branch 16 → 17 taken 9 times.
✓ Branch 59 → 49 taken 71 times.
✓ Branch 59 → 60 taken 6 times.
|
644146 | for (int i = 0; i < N; i++) { |
| 373 |
1/2✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 644060 times.
|
1288262 | deriv_root[i] *= T(N - i); |
| 374 | } | ||
| 375 |
1/1✓ Branch 91 → 92 taken 6 times.
|
21 | std::vector<T> denoms = tree.pushdown( |
| 376 |
5/5✓ Branch 17 → 18 taken 9 times.
✓ Branch 18 → 19 taken 9 times.
✓ Branch 20 → 21 taken 9 times.
✓ Branch 66 → 67 taken 6 times.
✓ Branch 68 → 69 taken 6 times.
|
69 | form<E>::from_rev_series(series::exact<E>(ps_inv(root) * deriv_root)) |
| 377 | ); | ||
| 378 | |||
| 379 |
2/2✓ Branch 28 → 29 taken 9 times.
✓ Branch 159 → 160 taken 6 times.
|
21 | std::vector<T> leaf_vals(size_t(N), T{}); |
| 380 |
4/4✓ Branch 32 → 30 taken 644060 times.
✓ Branch 32 → 33 taken 9 times.
✓ Branch 176 → 163 taken 71 times.
✓ Branch 176 → 177 taken 6 times.
|
644146 | for (int i = 0; i < N; i++) leaf_vals[i] = vals[i] / denoms[i]; |
| 381 |
2/2✓ Branch 35 → 36 taken 9 times.
✓ Branch 180 → 181 taken 6 times.
|
30 | return tree.combine_up(std::span<const T>(leaf_vals)); |
| 382 | 39 | } | |
| 383 | |||
| 384 | } // namespace wala::poly | ||
| 385 |