Skip to content

Commit 18b711b

Browse files
committed
Stage laurent.hpp
1 parent a2c3281 commit 18b711b

File tree

1 file changed

+178
-32
lines changed

1 file changed

+178
-32
lines changed

cp-algo/math/laurent.hpp

Lines changed: 178 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -6,31 +6,62 @@
66
#include <vector>
77
#include <cassert>
88
#include <algorithm>
9+
#include <iterator>
910
#include <type_traits>
1011
#include <bit>
12+
#include <ranges>
1113

1214
#include "cvector.hpp"
1315
#include "convolution.hpp"
1416

1517
namespace cp_algo::math {
18+
// Evaluation configuration using three boolean flags
19+
struct eval_options {
20+
// If true: evaluate only k and store it; if false: cache all up to k
21+
bool direct = false;
22+
// If true: lazy mode (can only use dependees up to current index); if false: eager
23+
bool lazy = false;
24+
// If true: ask dependees to precache up to needed indices before use
25+
bool precache_arguments = false;
26+
};
1627
// Base provider interface for lazy coefficient evaluation
1728
template<typename T>
1829
struct provider {
1930
mutable big_vector<T> cache;
31+
mutable int cache_valid = 0; // Number of valid cached coefficients from cache_offset
2032
mutable int cache_offset = 0; // Index of first cached coefficient
2133
mutable bool initialized = false;
2234
mutable bool all_cached = false; // True if all non-zero coeffs are cached
35+
// Ensure provider has initialized offset and state
36+
void initialize_if_needed() const {
37+
if(!initialized) {
38+
cache_offset = offset();
39+
initialized = true;
40+
}
41+
}
2342

2443
virtual ~provider() = default;
2544
virtual int offset() const { return 0; }
2645

2746
// Returns true if this provider requires lazy evaluation (coefficients must be
2847
// computed in order). False means dependencies can be bulk-cached for FFT.
29-
// Examples: multiply needs lazy eval, add/subtract/negate/scale don't.
3048
virtual bool needs_lazy_eval() const { return false; }
3149

3250
// Compute k-th coefficient lazily without caching
3351
virtual T coeff_lazy(int k) const = 0;
52+
// Compute k-th coefficient directly interface (no cache mutation)
53+
// Returns cached value if already available, otherwise delegates to impl.
54+
T coeff_direct(int k) const {
55+
initialize_if_needed();
56+
int idx = k - cache_offset;
57+
if(idx < 0) return T(0);
58+
if(idx < cache_valid) return cache[idx];
59+
return coeff_direct_impl(k);
60+
}
61+
62+
protected:
63+
// Implementation hook for direct computation without caching
64+
virtual T coeff_direct_impl(int k) const = 0;
3465

3566
// Double the number of known coefficients (or cache all for finite series)
3667
// Default: use coeff_lazy, but can be overridden for efficiency (e.g., FFT)
@@ -40,16 +71,42 @@ namespace cp_algo::math {
4071

4172
cache.resize(new_size);
4273
for(int i = old_size; i < new_size; i++) {
43-
cache[i] = coeff_lazy(cache_offset + i);
74+
cache[i] = coeff(cache_offset + i);
75+
}
76+
cache_valid = new_size;
77+
}
78+
79+
// Ensure the cache contains coefficients up to index n (inclusive)
80+
// Uses provider's preferred strategy: sequential for lazy providers, doubling otherwise.
81+
void upcache(int n) const {
82+
if(all_cached) return;
83+
initialize_if_needed();
84+
int idx = n - cache_offset;
85+
if(idx < 0) return;
86+
87+
if(needs_lazy_eval()) {
88+
while(idx >= cache_valid && !all_cached) {
89+
int next_k = cache_offset + cache_valid;
90+
if((int)cache.size() <= cache_valid) cache.resize(cache_valid + 1);
91+
cache[cache_valid] = coeff_lazy(next_k);
92+
cache_valid++;
93+
}
94+
} else {
95+
while(idx >= cache_valid && !all_cached) {
96+
double_up();
97+
}
4498
}
4599
}
46100

47101
// Get coefficient with caching and doubling (default implementation)
48102
virtual T coeff(int k) const {
49-
if(!initialized) {
50-
cache_offset = offset();
51-
initialized = true;
52-
}
103+
static constexpr eval_options defaults{};
104+
return coeff(k, defaults);
105+
}
106+
107+
// Get coefficient with explicit evaluation options
108+
virtual T coeff(int k, eval_options const& opt) const {
109+
initialize_if_needed();
53110

54111
int idx = k - cache_offset;
55112
if(idx < 0) {
@@ -60,31 +117,26 @@ namespace cp_algo::math {
60117
if(all_cached && idx >= (int)cache.size()) {
61118
return T(0);
62119
}
63-
64-
if(needs_lazy_eval()) {
65-
// Sequentially extend cache to the requested index
66-
while(idx >= (int)cache.size() && !all_cached) {
67-
int next_k = cache_offset + (int)cache.size();
68-
cache.push_back(coeff_lazy(next_k));
69-
}
120+
121+
if(opt.direct) {
122+
// Direct evaluation: compute only k and do not mutate cache
123+
if(idx < cache_valid) return cache[idx];
124+
return coeff_direct(k);
70125
} else {
71-
// Extend cache by doubling until we have enough
72-
while(idx >= (int)cache.size() && !all_cached) {
73-
double_up();
126+
// Expand cache honoring evaluation mode and provider's lazy requirement
127+
if(opt.lazy || needs_lazy_eval()) {
128+
upcache(k);
129+
} else {
130+
upcache(k);
74131
}
75132
}
76133

77-
if(idx < (int)cache.size()) {
134+
if(idx < cache_valid) {
78135
return cache[idx];
79136
}
80137

81138
return T(0);
82139
}
83-
84-
// Alias for backwards compatibility
85-
T get(int k) const {
86-
return coeff(k);
87-
}
88140
};
89141

90142
// Constant provider - returns a single coefficient at position offset
@@ -102,6 +154,9 @@ namespace cp_algo::math {
102154
T coeff_lazy(int k) const override {
103155
return k == offset ? value : T(0);
104156
}
157+
T coeff_direct_impl(int k) const override {
158+
return k == offset ? value : T(0);
159+
}
105160

106161
T coeff(int k) const override {
107162
return coeff_lazy(k);
@@ -111,27 +166,24 @@ namespace cp_algo::math {
111166
// Polynomial provider - wraps a vector of coefficients
112167
template<typename T>
113168
struct polynomial_provider : provider<T> {
114-
polynomial_provider(big_vector<T> coeffs, int offset = 0) {
169+
polynomial_provider(auto &&coeffs, int offset = 0) {
115170
// Find first and last non-zero coefficients
116171
auto non_zero = [](const T& x) { return x != T(0); };
117172
auto first = std::ranges::find_if(coeffs, non_zero);
118-
auto last = std::ranges::find_if(coeffs | std::views::reverse, non_zero);
173+
auto last = std::ranges::find_last_if(coeffs, non_zero);
119174

120175
// Extract non-zero range
121176
if(first != coeffs.end()) {
122177
int leading_zeros = first - coeffs.begin();
123-
int trailing_zeros = last - coeffs.rbegin();
124-
coeffs = big_vector<T>(first, coeffs.end() - trailing_zeros);
178+
int trailing_zeros = static_cast<int>(std::ranges::ssize(last)) - 1;
179+
this->cache = big_vector<T>(first, coeffs.end() - trailing_zeros);
125180
offset += leading_zeros;
126-
} else {
127-
// All zeros
128-
coeffs.clear();
129181
}
130182

131183
// Initialize cache directly with the coefficients
132-
this->cache = std::move(coeffs);
133184
this->cache_offset = offset;
134185
this->initialized = true;
186+
this->cache_valid = (int)this->cache.size();
135187
this->all_cached = true;
136188
}
137189

@@ -141,15 +193,22 @@ namespace cp_algo::math {
141193

142194
T coeff_lazy(int k) const override {
143195
int idx = k - this->cache_offset;
144-
if(idx < 0 || idx >= (int)this->cache.size()) {
196+
if(idx < 0 || idx >= this->cache_valid) {
145197
return T(0);
146198
}
147199
return this->cache[idx];
148200
}
201+
T coeff_direct_impl(int k) const override {
202+
return coeff_lazy(k);
203+
}
149204

150205
T coeff(int k) const override {
151206
return coeff_lazy(k);
152207
}
208+
209+
T coeff(int k, eval_options const& /*opt*/) const override {
210+
return coeff_lazy(k);
211+
}
153212
};
154213

155214
// Base class for unary operations
@@ -160,6 +219,10 @@ namespace cp_algo::math {
160219
unary_provider(std::shared_ptr<provider<T>> operand)
161220
: operand(std::move(operand)) {}
162221

222+
bool needs_lazy_eval() const override {
223+
return operand->needs_lazy_eval();
224+
}
225+
163226
virtual T transform(T const& a) const = 0;
164227

165228
int offset() const override {
@@ -169,10 +232,17 @@ namespace cp_algo::math {
169232
T coeff_lazy(int k) const override {
170233
return transform(operand->coeff_lazy(k));
171234
}
235+
T coeff_direct_impl(int k) const override {
236+
return transform(operand->coeff_direct(k));
237+
}
172238

173239
T coeff(int k) const {
174240
return transform(operand->coeff(k));
175241
}
242+
243+
T coeff(int k, eval_options const& opt) const override {
244+
return transform(operand->coeff(k, opt));
245+
}
176246
};
177247

178248
// Base class for binary operations
@@ -183,6 +253,10 @@ namespace cp_algo::math {
183253
binary_provider(std::shared_ptr<provider<T>> lhs, std::shared_ptr<provider<T>> rhs)
184254
: lhs(std::move(lhs)), rhs(std::move(rhs)) {}
185255

256+
bool needs_lazy_eval() const override {
257+
return lhs->needs_lazy_eval() || rhs->needs_lazy_eval();
258+
}
259+
186260
virtual T combine(T const& a, T const& b) const = 0;
187261

188262
int offset() const override {
@@ -192,10 +266,17 @@ namespace cp_algo::math {
192266
T coeff_lazy(int k) const override {
193267
return combine(lhs->coeff_lazy(k), rhs->coeff_lazy(k));
194268
}
269+
T coeff_direct_impl(int k) const override {
270+
return combine(lhs->coeff_direct(k), rhs->coeff_direct(k));
271+
}
195272

196273
T coeff(int k) const {
197274
return combine(lhs->coeff(k), rhs->coeff(k));
198275
}
276+
277+
T coeff(int k, eval_options const& opt) const override {
278+
return combine(lhs->coeff(k, opt), rhs->coeff(k, opt));
279+
}
199280
};
200281

201282
// Addition provider
@@ -272,6 +353,63 @@ namespace cp_algo::math {
272353
}
273354
return result;
274355
}
356+
T coeff_direct_impl(int k) const override {
357+
int n = k - offset();
358+
if(n < 0) return T(0);
359+
T result = T(0);
360+
for(int j = 0; j <= n; j++) {
361+
int i_l = lhs->offset() + j;
362+
int i_r = rhs->offset() + (n - j);
363+
auto a = lhs->coeff(i_l);
364+
auto b = rhs->coeff(i_r);
365+
result += a * b;
366+
}
367+
return result;
368+
}
369+
370+
T coeff(int k, eval_options const& opt) const override {
371+
int n = k - offset();
372+
if(n < 0) return T(0);
373+
374+
if(opt.direct) {
375+
T result = T(0);
376+
for(int j = 0; j <= n; j++) {
377+
int i_l = lhs->offset() + j;
378+
int i_r = rhs->offset() + (n - j);
379+
if(opt.precache_arguments && !opt.lazy) {
380+
lhs->coeff(i_l);
381+
rhs->coeff(i_r);
382+
}
383+
auto a = (opt.lazy || lhs->needs_lazy_eval())
384+
? lhs->coeff_lazy(i_l) : lhs->coeff(i_l);
385+
auto b = (opt.lazy || rhs->needs_lazy_eval())
386+
? rhs->coeff_lazy(i_r) : rhs->coeff(i_r);
387+
result += a * b;
388+
}
389+
return result;
390+
}
391+
392+
int idx = k - this->cache_offset;
393+
if(idx < 0) return T(0);
394+
while(idx >= (int)this->cache.size() && !this->all_cached) {
395+
if(opt.lazy || needs_lazy_eval()) {
396+
int next_k = this->cache_offset + (int)this->cache.size();
397+
if(opt.precache_arguments && !opt.lazy) {
398+
int n2 = next_k - offset();
399+
int lhs_need = lhs->offset() + n2;
400+
int rhs_need = rhs->offset() + n2;
401+
lhs->upcache(lhs_need);
402+
rhs->upcache(rhs_need);
403+
}
404+
if((int)this->cache.size() <= this->cache_valid) this->cache.resize(this->cache_valid + 1);
405+
this->cache[this->cache_valid] = coeff_lazy(next_k);
406+
this->cache_valid++;
407+
} else {
408+
double_up();
409+
}
410+
}
411+
return idx < this->cache_valid ? this->cache[idx] : T(0);
412+
}
275413

276414
void double_up() const override {
277415
int old_size = this->cache.size();
@@ -280,7 +418,9 @@ namespace cp_algo::math {
280418
// Lazy path: compute the next coefficient sequentially
281419
if(needs_lazy_eval()) {
282420
int k = this->cache_offset + old_size;
283-
this->cache.push_back(coeff_lazy(k));
421+
if((int)this->cache.size() <= this->cache_valid) this->cache.resize(this->cache_valid + 1);
422+
this->cache[this->cache_valid] = coeff_lazy(k);
423+
this->cache_valid++;
284424
return;
285425
}
286426

@@ -302,6 +442,7 @@ namespace cp_algo::math {
302442
for(int i = old_size; i < new_size && i < (int)la.size(); i++) {
303443
this->cache[i] = la[i];
304444
}
445+
this->cache_valid = new_size;
305446

306447
// If both operands are fully cached and we reached their total length, mark as done
307448
if(lhs->all_cached && rhs->all_cached) {
@@ -371,6 +512,11 @@ namespace cp_algo::math {
371512
laurent& operator*=(T const& scalar) {
372513
return *this = *this * scalar;
373514
}
515+
516+
// Coefficient with explicit evaluation options
517+
T coeff(int k, eval_options const& opt) const {
518+
return impl->coeff(k, opt);
519+
}
374520
};
375521

376522
template<typename T>

0 commit comments

Comments
 (0)