Added divide, inverse, mod, derive, and integr

Author: Chillee

content/number-theory/ModularArithmetic.h

@@ -13,6 +13,7 @@
 const ll mod = 17; // change to something else
 struct Mod {
 	ll x;
+	Mod():x(0) {}
 	Mod(ll xx) : x(xx) {}
 	Mod operator+(Mod b) { return Mod((x + b.x) % mod); }
 	Mod operator-(Mod b) { return Mod((x - b.x + mod) % mod); }
@@ -27,4 +28,5 @@ struct Mod {
 		Mod r = *this ^ (e / 2); r = r * r;
 		return e&1 ? *this * r : r;
 	}
+	operator ll(){ return x; }
 };

content/numerical/FFTPolynomial.h

@@ -5,18 +5,24 @@
  */
 #pragma once
 
+#include "../number-theory/ModularArithmetic.h"
 #include "FastFourierTransform.h"
+#include "FastFourierTransformMod.h"
 
-typedef double num;
+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)));
+    return vector<Mod>(all(res));
+}
 poly &operator+=(poly &a, poly &b) {
     a.resize(max(sz(a), sz(b)));
-    rep(i,0,sz(b)) a[i] += b[i];
+    rep(i,0,sz(b)) a[i] = a[i] + b[i];
     return a;
 }
 poly &operator-=(poly &a, poly b) {
     a.resize(max(sz(a), sz(b)));
-    rep(i,0,sz(b)) a[i] -= b[i];
+    rep(i,0,sz(b)) a[i] = a[i] - b[i];
     return a;
 }
 poly &operator*=(poly &a, poly &b) { return a = conv(a, b); }
@@ -26,3 +32,44 @@ poly &operator*=(poly &a, poly &b) { return a = conv(a, b); }
         return c oe b; \
     }
 OP(*, *=) OP(+, +=) OP(-, -=);
+poly modK(poly a, int k) { return {a.begin(), a.begin() + min(k, sz(a))}; }
+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));
+    return modK(B, sz(A));
+}
+poly &operator/=(poly &a, poly &b) {
+    if (sz(a) < sz(b))
+        return a = {};
+    int s = sz(a) - sz(b) + 1;
+    reverse(all(a)), reverse(all(b));
+    a.resize(s), b.resize(s);
+    a = a * inverse(b);
+    a.resize(s), reverse(all(a));
+    return a;
+}
+OP(/, /=)
+poly &operator%=(poly &a, poly &b) {
+    if (sz(a) < sz(b))
+        return a;
+    poly c = (a / b) * b;
+    a.resize(sz(b) - 1);
+    rep(i, 0, sz(a)) a[i] = a[i] - c[i];
+    return a;
+}
+OP(%, %=)
+poly deriv(poly a) {
+    if (a.empty())
+        return {};
+    poly b(sz(a) - 1);
+    rep(i, 1, sz(a)) b[i - 1] = a[i] * num(i);
+    return b;
+}
+poly integr(poly a) {
+    if (a.empty())
+        return {};
+    poly b(sz(a) + 1);
+    rep(i, 0, sz(a) + 1) b[i + 1] = a[i] / num(i + 1);
+    return b;
+}

content/numerical/FastFourierTransformMod.h

@@ -18,8 +18,8 @@ template <int M> vl convMod(const vl &a, const vl &b) {
 	vl res(sz(a) + sz(b) - 1);
 	int B=32-__builtin_clz(sz(res)), n = 1<<B, cut=int(sqrt(M));
 	vector<C> L(n), R(n), outs(n), outl(n), rt;
-	rep(i,0,sz(a)) L[i] = C(a[i] / cut, a[i] % cut);
-	rep(i,0,sz(b)) R[i] = C(b[i] / cut, b[i] % cut);
+	rep(i,0,sz(a)) L[i] = Cd(a[i] / cut, a[i] % cut);
+	rep(i,0,sz(b)) R[i] = Cd(b[i] / cut, b[i] % cut);
 	fft(L, n, B, rt), fft(R, n, B, rt);
 	rep(i,0,n) {
 		int j = -i & (n - 1);