Added a FFTPolynomial file #87
Conversation
|
TBH I don't have a strong opinion on this. It would be nice to join this together with the existing Polynomial.h; there's probably some shared code, and it would make the code foot print of the existing code less, which is nice since it's basically never used. |
|
I was more thinking that after I finish writing this we get rid of the old Polynomial class entirely. I don't think there's much point for it. |
|
There's not much in the old Polynomial class, it's not documented, etc. I made a new file simply so things don't break while I write the new version and that I could submit some smaller PR's instead of a huge 100+line file with like...15 functions. |
|
That's fair. Not that I would mind a large PR; it gives a nice impression of how we imagine things looking in the end. Also, note that PolyRoots.h depends on Polynomial.h and will need some small changes if it's removed. |
|
Ok then, if you insist :^) |
|
It seems to make no impact on performance (perhaps even a negative impact). |
|
I've reorganized the functions into 5 files (with their dependencies by the side (including transitive dependencies). I've tried to balance minimizing number of extraneous dependencies, giving enough space for documentation, and grouping together things by theme.
|
content/numerical/PolynomialBase.h
Outdated
| #include "../number-theory/ModularArithmetic.h" | ||
| #include "FastFourierTransform.h" | ||
| #include "FastFourierTransformMod.h" | ||
| // #include "NumberTheoreticTransform.h" |
There was a problem hiding this comment.
Nit: NTT will be the common case, if that matters to you?
There was a problem hiding this comment.
This is currently commented out for a stupid reason - NumberTheoreticTransform.h doesn't compile properly on its own (due to redefinitions of const mod). It is good to know that NTT is the common case though.
content/numerical/PolynomialPow.h
Outdated
| if (a.empty()) return {0}; | ||
| poly b(sz(a) + 1); | ||
| b[1] = num(1); | ||
| rep(i, 2, sz(b)) b[i] = b[mod%i]*Mod(-mod/i+mod); |
There was a problem hiding this comment.
Doesn't work for complex, probably worth noting.
There was a problem hiding this comment.
Is it even worth noting stuff that doesn't work with complex numbers? Or is it reasonable to assume that every problem involving polynomials/series will be done on a finite field?
There was a problem hiding this comment.
Idk, the field of reals seems like it could be useful sometimes.
| Cd z[] = {1, polar(1.0, M_PI / k)}; | ||
| rep(i, k, 2 * k) rt[i] = Cd(rt[i / 2]) * z[i & 1]; | ||
| } | ||
| } | ||
| rep(i,0,n) if (i < rev[i]) swap(a[i], a[rev[i]]); |
There was a problem hiding this comment.
You can merge this loop with the generating loop if you want, would keep things more compact and (marginally) improve cache locality.
| * Date: 2017-05-10 | ||
| * License: CC0 | ||
| * Source: Wikipedia | ||
| * Author: chilli, Andrew He, Adamant |
There was a problem hiding this comment.
Much of my work comes from THU, probably should give them credit too.
|
I removed the
|
|
If you want, here's my modinv: |
|
Also, empirically, modinv is somewhere between 1x and 2x faster than modular exponentiation. |
Chillee
changed the title
Chillee
changed the title
| ll v; | ||
| Mod() : v(0) {} | ||
| Mod(ll vv) : v(vv % mod) {} | ||
| Mod operator+(Mod b) { return Mod((v + b.v) % mod); } |
| Mod operator*(Mod b) { return Mod((x * b.x) % mod); } | ||
| ll v; | ||
| Mod() : v(0) {} | ||
| Mod(ll vv) : v(vv % mod) {} |
There was a problem hiding this comment.
ll vv = 0 to get rid of the default constructor
| ll x, y, g = euclid(a.x, mod, x, y); | ||
| assert(g == 1); return Mod((x + mod) % mod); | ||
| } | ||
| Mod invert(Mod a) { return a^(mod-2); } |
There was a problem hiding this comment.
assumes the modulo is prime, which is probably worth a comment. Also, inline this into operator/. Or maybe keep this method with a comment about it being only for composite moduli.
| @@ -0,0 +1,36 @@ | |||
| /** | |||
There was a problem hiding this comment.
(need to merge with master to get rid of these fft changes from the diff)
| typedef vector<num> poly; | ||
| poly &operator+=(poly &a, const poly &b) { | ||
| a.resize(max(sz(a), sz(b))); | ||
| rep(i, 0, sz(b)) a[i] = a[i] + b[i]; |
There was a problem hiding this comment.
| rep(i, 0, sz(b)) a[i] = a[i] + b[i]; | |
| rep(i, 0, sz(b)) a[i] = a[i] + b[i]; | |
| rep(i,0,sz(b)) a[i] = a[i] + b[i]; |
(same in other places)
| mine::poly a; | ||
| MIT::poly am; | ||
| for (int i = 0; i < sz; i++) { | ||
| int val = rand(); |
There was a problem hiding this comment.
is this mod a really large number? if so, it would make sense to test zeroes as well, since they are special. It would also make sense to lower the field size to make zeroes happen more naturally. (and also maybe using some fancier finite field to make sure we don't make too many assumptions about it)
|
|
||
| } | ||
|
|
||
| const int NUMITERS=10; |
There was a problem hiding this comment.
way too small. maybe make this non-const and set it in main
| template <class A, class B> void testPow(string name, A f1, B f2, int mxSz = 5, int mxPref=5) { | ||
| for (int it = 0; it < NUMITERS; it++) { | ||
| auto a = genVec((rand() % mxSz) + 1); | ||
| int pref = rand()%mxSz; |
| } | ||
| cout<<endl; | ||
| } | ||
| signed main() { |
There was a problem hiding this comment.
put this in a special benchmarking main, and have a main that just does testing
| auto duration = chrono::duration_cast<chrono::milliseconds>(end - begin).count(); | ||
| cerr << duration << "ms elapsed [" << label << "]" << endl; | ||
| } | ||
| }; |
There was a problem hiding this comment.
move this to utilities/bench.h and replace the nonsense I put in there
This PR is mainly to get a conversation started around how we should be representing our polynomials. Looking around at MIT's implementation and Adamant's implementation, there seem to be 2 primary design choices.
First, how do we represent our polynomials. MIT does it with a
typedef vector<num> poly. Adamant does it with a struct. I'm leaning towards doing it with avector<num>, simply due to having slightly less boilerplate. The primary advantage of doing it in a struct is that we can then define a functions as methods of the struct instead of global properties. For example, I think doingpoly.inv().resize(K)is more natural thanresize(inv(poly), K).Second, how do we make our polynomials work across the 3 "fields" we support. Namely, doubles, ints with NTT, and ints with FFTMOD. MIT does this by wrapping int modulo operations in a struct. I think we can do this too (probably with
ModularArithmetic.hand whatever API changes are necessary). I think with proper enough encapsulation, we should only need to change out one function, namely the*=operator.I've started implementing the other operations, but I wanted to submit a PR with the basics so we could settle on a structure.
Edit: closes #60.