teamnote history merge

2025fall/source/Math/Polynomial.cpp

@@ -0,0 +1,149 @@
+namespace GMS {
+template<ll mod, ll w>
+struct Qring : public vl {
+    using poly = Qring<mod, w>;
+    Qring()                 : vl(1, 0)      {}
+    Qring(ll c)             : vl(1, c%mod)  {}
+    Qring(ll c, int n)      : vl(n, c%mod)  {}
+    Qring(const vl& cp)     : vl(cp)        {for(auto &i:*this) i%=mod;}
+    
+    ll& operator[](ll idx) {
+        if((unsigned)idx < size()) return vl::operator[](idx);
+        this->resize(idx+1); return vl::operator[](idx);
+    }
+    ll operator[](ll idx) const {
+        if((unsigned)idx < size()) return vl::operator[](idx);
+        return 0LL;
+    }
+    void adjust() { while(size() > 1 and back() == 0) pop_back(); }
+    void adjust(int n){resize(n, 0);}
+    ll operator()(ll x) {
+        x %= mod; ll ret = 0;
+        for(auto it=rbegin(); it!=rend(); it++)
+            ret = (ret*x+*it)%mod;
+        return ret;
+    }
+    friend poly operator%(const poly& A, int B){ // remainder by x^B
+        poly ret(A); ret.resize(B, 0);
+        return ret;
+    }
+    friend poly operator+(const poly& A, const poly& B) {
+        int n = max(A.size(), B.size()); poly ret(0, n);
+        fors(i, 0, n-1) ret[i] = A[i]+B[i];
+        for(auto&i:ret) if(i >= mod) i -= mod;
+        return ret.adjust(), ret;
+    }
+    friend poly operator-(const poly& A) {
+        int n = A.size(); poly ret(0, n);
+        fors(i, 0, n-1) ret[i] = A[i]?mod-A[i]:0;
+        return ret;
+    }
+    friend poly operator-(const poly& A, const poly& B) {
+        int n = max(A.size(), B.size());
+        poly ret(0, n);
+        fors(i, 0, n-1) ret[i] = (mod+A[i]-B[i])%mod;
+        return ret.adjust(), ret;
+    }
+    friend poly operator*(ll x, const poly& B) {
+        poly ret(B); x %= mod;
+        for(auto &i : ret) i = (i*x)%mod;
+        return ret.adjust(), ret;
+    }
+    friend poly operator*(const poly& A, const poly& B) {
+        poly ret = conv<mod, w>(A, B);
+        // ACL : poly ret = atcoder::convolution<mod>(A, B);
+        return ret.adjust(), ret;
+    }
+    friend poly inv(const poly& A, int t) { assert(A[0] != 0);
+        poly g = pow<mod>(A[0], mod-2); int st=1;
+        while(st < t){st*=2; g = (-A%st*g%st+2)*g%st;}
+        return g.adjust(t), g;
+    }
+    friend poly diff(const poly& A) {
+        int n = A.size(); poly ret(0, n-1);
+        fors(i, 0, n-2) ret[i] = (i+1)*A[i+1]%mod;
+        return ret;
+    }
+    friend poly inte(const poly& A) {
+        static vector<ll> inv(1, 1);
+        int n = A.size(); poly ret(0, n+1);
+        inv.resize(max((int)inv.size(), n+1));
+        forr(i, n) if(inv[i] == 0) inv[i] = pow<mod>(i, mod-2);
+        forr(i, n) ret[i] = inv[i]*A[i-1]%mod;
+        return ret;
+    }
+    friend poly log(const poly& A, int t) { assert(A[0] == 1);
+        return inte(diff(A) * inv(A, t)%t)%t;
+    }    
+    friend poly exp(const poly& A, int t) { assert(A[0] == 0);
+        poly g = 1; int st = 1;
+        while(st < t) {st*=2; g = (A%st-log(g, st)+1)*g%st;}
+        return g.adjust(), g;
+    }
+    friend poly pow(const poly& A, ll b, ll t) {
+        poly ret(A); ret.adjust();
+        if(ret.size() == 1) {
+            ret[0] = pow<mod>(ret[0], b);
+            return ret.adjust(t), ret;
+        }
+        ll idx = 0; while(ret[idx] == 0) idx++;
+        if((__int128_t) idx * b >= t) return poly(0, t);
+        
+        ll c = ret[idx]; ll ic = pow<mod>(ret[idx], mod-2); poly g;
+        int n = ret.size();
+        fors(i, idx, n-1) g[i-idx] = ret[i]*ic%mod;
+        g.resize(t-idx*b);
+        
+        g = exp(b * log(g, t-idx*b), t-idx*b);
+        c = pow<mod>(c, b);
+        
+        ret = poly(0, t); fors(i, idx*b, t-1) ret[i] = g[i-idx*b] * c % mod;
+        
+        return ret;
+    }
+    
+    //Only just Polynomial, not Qring
+    void rev() { reverse(begin(), end()); }
+    friend poly div_quot(poly F, poly G) {
+        F.adjust(); G.adjust();
+        ll df = F.size(), dg = G.size();
+        if(df < dg) return poly(0);
+        F.rev(); G.rev();
+        F = F%(df-dg+1)*inv(G, df-dg+1)%(df-dg+1); F.rev();
+        return F;
+    }
+    friend poly div_rem(poly F, poly G) {return F-G*div_quot(F, G);}
+    friend poly shift(const poly& F, ll c) {
+        ll n = F.size(); c %= mod;
+        poly A(0, n); ll fac = 1;
+        fors(i, 0, n-1) A[i] = F[i]*fac%mod, fac = fac*(i+1)%mod;
+        A.rev();
+        
+        poly C(1, n); fors(i, 1, n-1) C[i] = C[i-1]*c%mod;
+        ll facc = fac = pow<mod>(fac, mod-2)*n%mod;
+        fore(i, n-1, 0) C[i] = C[i]*fac%mod, fac = fac*i%mod;
+        
+        poly B = C*A; B.resize(n); B.rev();
+        fore(i, n-1, 0) B[i] = B[i]*facc%mod, facc = facc*i%mod;
+        return B;
+    }
+    friend void calcG(vector<poly>& G, int i, int l, int r, const vl& p) {
+        if(l == r){ll g = p[l]?mod-p[l]:0; G[i] = vl({g, 1}); return;}
+        int mid = (l+r)/2;
+        calcG(G, i*2, l, mid, p); calcG(G, i*2+1, mid+1, r, p);
+        G[i] = G[i*2]*G[i*2+1];
+    }
+    friend void eval(const vector<poly>& G, int i, int l, int r, poly&& F, vl& ret) {
+        if(l == r){ret[l] = F[0];return;}
+        int mid=(l+r)/2;
+        eval(G, i*2, l, mid, div_rem(F, G[i*2]), ret);
+        eval(G, i*2+1, mid+1, r, div_rem(F, G[i*2+1]), ret);
+    }
+    friend vl multipoint_eval(const poly& A, const vl& B) {
+        int m = B.size();
+        vector<poly> G(4*m); calcG(G, 1, 0, m-1, B);
+        vl ret(m, 0); eval(G, 1, 0, m-1, div_rem(A, G[1]), ret);
+        return ret;
+    }
+};
+}  // namespace GMS