Added authors, a fuzz test, and interp

Author: Chillee

content/numerical/FFTPolynomial.h

@@ -1,5 +1,5 @@
 /**
- * Author: chilli
+ * Author: chilli, Andrew He, Adamant
  * Date: 2019-04-27
  * Description: A FFT based Polynomial class.
  */
@@ -8,12 +8,13 @@
 #include "../number-theory/ModularArithmetic.h"
 #include "FastFourierTransform.h"
 #include "FastFourierTransformMod.h"
+// #include "NumberTheoreticTransform.h"
 
 typedef Mod num;
 typedef vector<num> poly;
 vector<Mod> conv(vector<Mod> a, vector<Mod> b) {
-    // auto res = convMod<mod>(vl(all(a)), vl(all(b)));
-    auto res = conv(vl(all(a)), vl(all(b)));
+    auto res = convMod<mod>(vl(all(a)), vl(all(b)));
+    // auto res = conv(vl(all(a)), vl(all(b)));
     return vector<Mod>(all(res));
 }
 poly &operator+=(poly &a, const poly &b) {
@@ -39,12 +40,25 @@ poly operator*(poly a, const num b) {
     }
 OP(*, *=) OP(+, +=) OP(-, -=);
 poly modK(poly a, int k) { return {a.begin(), a.begin() + min(k, sz(a))}; }
+// Currently there's two of them - the second is the original one (simply following the formula), the first one is a version Adamant says is faster.
+// I haven't been able to replicate the difference in performance, however.
 poly inverse(poly A) {
     poly B = poly({num(1) / A[0]});
-    while (sz(B) < sz(A))
-        B = modK(B * (poly({num(2)}) - A * B), 2 * sz(B));
+    while (sz(B) < sz(A)){
+        poly C = B*modK(A, 2*sz(B));
+        C = poly(C.begin()+sz(B), C.end());
+        C = modK(B*C, sz(B));
+        C.insert(C.begin(), sz(B), 0);
+        B -= C;
+    }
     return modK(B, sz(A));
 }
+// poly inverse(poly A) {
+//     poly B = poly({num(1) / A[0]});
+//     while (sz(B) < sz(A))
+//         B = modK(B * (poly({num(2)}) - modK(A, 2*sz(B)) * B), 2 * sz(B));
+//     return modK(B, sz(A));
+// }
 poly &operator/=(poly &a, poly b) {
     if (sz(a) < sz(b))
         return a = {};
@@ -73,10 +87,11 @@ poly deriv(poly a) {
     return b;
 }
 poly integr(poly a) {
-    if (a.empty())
-        return {0};
+    if (a.empty()) return {0};
     poly b(sz(a) + 1);
-    rep(i, 1, sz(b)) b[i] = a[i - 1] / num(i);
+    b[1] = num(1);
+    rep(i, 2, sz(b)) b[i] = b[mod%i]*Mod(-mod/i+mod);
+    rep(i, 1 ,sz(b)) b[i] = a[i-1] * b[i];
     return b;
 }
 poly log(poly a) { return modK(integr(deriv(a) * inverse(a)), sz(a)); }
@@ -102,4 +117,33 @@ poly pow(poly a, ll m) {
     auto res =  exp(log(a) * num(m)) * (j ^ m);
     res.insert(res.begin(), p*m, 0);
     return modK(res, n);
+}
+
+vector<num> eval(const poly &a, const vector<num> &x) {
+    int n = sz(x);
+    if (!n) return {};
+    vector<poly> up(2 * n);
+    rep(i, 0, n) up[i + n] = poly({num(0) - x[i], 1});
+    for (int i = n - 1; i > 0; i--)
+        up[i] = up[2 * i] * up[2 * i + 1];
+    vector<poly> down(2 * n);
+    down[1] = a % up[1];
+    rep(i, 2, 2 * n) down[i] = down[i / 2] % up[i];
+    vector<num> y(n);
+    rep(i, 0, n) y[i] = down[i + n][0];
+    return y;
+}
+
+poly interp(vector<num> x, vector<num> y) {
+    int n = sz(x);
+    vector<poly> up(n * 2);
+    rep(i, 0, n) up[i + n] = poly({num(0) - x[i], num(1)});
+    for (int i = n - 1; i > 0; i--)
+        up[i] = up[2 * i] * up[2 * i + 1];
+    vector<num> a = eval(deriv(up[1]), x);
+    vector<poly> down(2 * n);
+    rep(i, 0, n) down[i + n] = poly({y[i] * (num(1) / a[i])});
+    for (int i = n - 1; i > 0; i--)
+        down[i] = down[i * 2] * up[i * 2 + 1] + down[i * 2 + 1] * up[i * 2];
+    return down[1];
 }