多项式运算

现在大概有：

多项式乘法逆
多项式开根
多项式带余除法
多项式牛顿迭代
多项式 \(\exp\)
多项式 \(\ln\)
多项式快速幂

多项式复合函数一类的高端内容后面再补。

很多地方的严谨证明因为水平不足略过，等之后数学水平提高再回来修正。

一些有用的资料：

NTT与多项式全家桶 - command_block 的博客

多项式全家桶 - ckj blog

原根&离散对数相关 - command_block 的博客

你也能懂的微积分 - 知乎

前言

运算定义

这里的多项式运算都可以用加法和乘法（加法卷积）定义。

通常对于无限项的多项式，只考虑多项式的前 \(n\) 项，也就是在 \(\pmod {x^n}\) 意义下运算。

相应地有运算性质：

\(F(x) \bmod x^n + G(x) \bmod x^n \equiv (F(x) + G(x)) \bmod x^n \pmod {x^n}\)
\(F(x) \bmod x^n \cdot G(x) \bmod x^n \equiv (F(x) \cdot G(x)) \bmod x^n \pmod {x^n}\)

所以可以在运算的过程中限定多项式的界，从而避免高次项。

当然这个也可以用来卡常。

循环卷积

普通的卷积是线性卷积，另外还有周期卷积和循环卷积。

循环卷积的定义是：

\((F * G)(x) = \sum\limits_{k = 0}^m x^k \sum\limits_{i + j \bmod m = k} F_i G_j\)

更直观的形式是 \((F * G)(x) = \sum\limits_{i = 0}^{m} \sum\limits_{j = 0}^m F_i G_j x^{(i + j) \bmod m}\)

封装好的 NTT 模板：

void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }

ll qpow(ll base, ll power)
{
    ll res = 1;
    while (power)
    {
        if (power & 1) res = res * base % mod;
        base = base * base % mod;
        power >>= 1;
    }
    return res;
}

void NTT(ll *A, int n)
{
    calc_rev(n);
    for (int i = 1; i < n; i++)
        if (rev[i] > i) swap(A[i], A[rev[i]]);
    for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
    {
        ll wn = qpow(g, (mod - 1) / len);
        wp[0] = 1;
        for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
        for (int l = 0, r = len - 1; r <= n; l += len, r += len)
        {
            int w = 0;
            for (int p = l; p < l + m; p++, w++)
            {
                ll x = A[p], y = wp[w] * A[p + m] % mod;
                A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
            }
        }
    }
}

void INTT(ll *A, int n)
{
    NTT(A, n);
    reverse(A + 1, A + n);
    int inv = qpow(n, mod - 2);
    for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;
}

多项式乘法逆

P4238 【模板】多项式乘法逆

求多项式在模意义下的逆元。

朴素做法

这个虽然有一种复杂度为 \(O(n^2)\) 的递推做法，但是有更优的 \(O(n \log n)\) NTT 倍增做法。

当界为 \(\pmod {x^1}\) （只有常数项）时直接求自然数的逆元就行。

反之，假设已经有 \(R_* (x) \equiv F(x)^{-1} \pmod {x^\frac{n}{2}}\)，考虑倍增求出 \(R(x) \equiv F(x)^{-1} \pmod {x^n}\).

因为 \(R(x) \equiv R_*(x) \pmod {x^{\frac{n}{2}}}\)，所以 \(R(x) - R_*(x) \equiv 0 \pmod {x^{\frac{n}{2}}}\).

考虑对两边平方得到 \((R(x) - R_*(x))^2 \equiv 0 \pmod {x^n}\)

展开得到 \(R(x)^2 - 2 R(x) R_*(x) + R_*(x)^2 \equiv 0 \pmod {x^n}\)

两边同乘 \(F(x)\) 得到 \(R(x) - 2 R_*(x) + R_*(x)^2 F(x) \equiv 0 \pmod {x^n}\)

整理即 \(R(x) \equiv 2 R_*(x) - R_*(x)^2 F(x) \pmod {x^n}\)

时间复杂度为 \(T(n) = T(\frac{n}{2}) + O(n \log n) = O(n \log n)\)，下面的倍增法同理。

优化

NTT 代入的是原根，它具有单位根的性质，所以实际上 NTT 求的是循环卷积。

观察倍增的式子：\(R(x) \equiv 2 R_*(x) - R_*(x)^2 F(x) \pmod {x^n}\).

需要用到多项式乘法的项只有 \(R_*(x)^2 F(x)\)，其中 \(R_*(x)\) 的次数是 \(\frac{n}{2}\)，\(F(x)\) 的次数是 \(n\)，并且我们只关心 \(R(x)\) 的后 \(\frac{n}{2}\) 项。

可以考虑先算 \(R_*(x) F(x)\)，此时由 \(R_*(x)\) 的定义得其前 \(\frac{n}{2}\) 项中只有第一项为 \(1\)，其余都为 \(0\). 于是循环卷积只会影响前 \(\frac{n}{2}\) 项。

然后再考虑乘上 \(R_*(x)\)，正好还是影响前 \(\frac{n}{2}\) 项。

void invp(ll *f, ll *r, int n)
{
    int k = 1;
    while (k < n) k <<= 1;
    r[0] = qpow(f[0], mod - 2);
    for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
    {
        for (int i = 0; i < len; i++) Rt[i] = r[i], Ft[i] = f[i];
        NTT(Ft, len), NTT(Rt, len);
        for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
        INTT(Rt, len);
        for (int i = 0; i < m; i++) Rt[i] = 0; Rt[0] = 1;
        for (int i = 0; i < len; i++) Ft[i] = r[i];
        NTT(Ft, len), NTT(Rt, len);
        for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
        INTT(Rt, len);
        for (int i = m; i < len; i++) r[i] = (r[i] * 2ll - Rt[i] + mod) % mod;
    }
    memset(Ft, 0, k * sizeof(ll));
    memset(Rt, 0, k * sizeof(ll));
    for (int i = n; i < k; i++) r[i] = 0;
}

多项式开根

P5205 【模板】多项式开根

还是考虑倍增。

令 \(B(x)^2 \equiv A(x) \pmod {x^n}\)，假设现在已知 \(B_*(x) \equiv A(x) \pmod {x^{\frac{n}{2}}}\)，考虑倍增求出 \(B(x)\).

同求逆有 \(B(x) - B_*(x) \equiv 0 \pmod {x^{\frac{n}{2}}}\).

平方再展开得 \(B(x)^2 - 2 B(x) B_*(x) + B_*(x)^2 \equiv 0 \pmod {x^n}\).

又因为 \(B(x)^2 \equiv A(x) \pmod {x^n}\)，代入得：

\(A(x) - 2 B(x) B_*(x) + B_*(x)^2 \equiv 0 \pmod {x^n}\).

整理得 \(B(x) \equiv \frac{\frac{A(x)}{B_*(x)} + B_*(x)}{2}\)

考虑上全家桶就行。

注意常数项是二次剩余，当 \(A(x)\) 常数项为 \(1\) 时显然 \(B(x)\) 常数项也为 \(1\).

void sqrtp(ll *f, ll *s, int n)
{
    s[0] = mod_sqrt(f[0]);
    if (s[0] == -1) return;
    int k = 1;
    while (k < (n << 1)) k <<= 1;
    ll inv2 = qpow(2, mod - 2, mod);
    for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
    {
        invp(s, rt, m);
        for (int i = 0; i < m; i++) ft[i] = f[i];
        NTT(ft, len), NTT(rt, len);
        for (int i = 0; i < len; i++) rt[i] = rt[i] * ft[i] % mod;
        INTT(rt, len);
        for (int i = 0; i < m; i++) rt[i] = (rt[i] + s[i]) % mod;
        for (int i = 0; i < len; i++) s[i] = rt[i] * inv2 % mod, rt[i] = ft[i] = 0;
    }
    for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
    for (int i = n; i < k; i++) s[i] = 0;
}

多项式带余除法

P4512 【模板】多项式除法

整除的情况比较好做，考虑把余数去掉。

舍弃高次项可以直接用界做，所以考虑把整个多项式反转（对称的项互换系数）。

定义多项式 \(F(x)\) 的反转 \(F_R(x) = x^n F(\frac{1}{x})\).

考虑题目所求：\(F(x) = Q(x) G(x) + R(x)\)，其中 \(F\) 是 \(n\) 次已知多项式，\(G\) 是 \(m\) 次已知多项式，\(Q(x)\) 是 \(n - m\) 次未知多项式，\(Q\) 是 \(m - 1\) 次未知多项式（不够用 \(0\) 补足）。

考虑换元，根据上式有 \(F(\frac{1}{x}) = Q(\frac{1}{x}) G(\frac{1}{x}) + R(\frac{1}{x})\).

两边同乘 \(x^n\) 得 \(x^n F(\frac{1}{x}) = x^n Q(\frac{1}{x}) G(\frac{1}{x}) + x^n R(\frac{1}{x})\)

又因为 \(x^n F(\frac{1}{x}) = F_R(x), x^n Q(\frac{1}{x}) G(\frac{1}{x}) = Q_R(x) G_R(x), x^n R(\frac{1}{x}) = x^{n - m + 1} R_R(x)\)，代入原式得：

\(F_R(x) = Q_R(x) G_R(x) + x^{n - m + 1} R_R(x)\)

注意到 \(R(x)\) 是 \(n - m\) 次多项式，不妨给上式加上界：

\(F_R(x) = Q_R(x) G_R(x) + x^{n - m + 1} R_R(x) \pmod {x^{n - m + 1}}\)

于是可以求出 \(Q_R(x)\)，从而求出 \(Q(x)\) 和 \(R(x)\).

void divp(ll *f, ll *g, ll *quo, ll *rem, int n, int m)
{
    int len = n - m + 1;
    for (int i = 0; i < min(len, m); i++) rt[i] = g[m - i - 1];
    invp(rt, ft, len);
    for (int i = 0; i < len; i++) rt[i] = f[n - i - 1];
    int k = 1;
    while (k < (len << 1)) k <<= 1;
    NTT(rt, k), NTT(ft, k);
    for (int i = 0; i < k; i++) rt[i] = rt[i] * ft[i] % mod;
    INTT(rt, k);
    for (int i = 0; i < len; i++) quo[i] = rt[len - i - 1];
    for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
    for (int i = 0; i < m; i++) rt[i] = g[i];
    for (int i = 0; i < len; i++) ft[i] = quo[i];
    k = 1;
    while (k < (m + len)) k <<= 1;
    NTT(rt, k), NTT(ft, k);
    for (int i = 0; i < k; i++) rem[i] = rt[i] * ft[i] % mod;
    INTT(rem, k);
    for (int i = 0; i < m - 1; i++) rem[i] = (f[i] - rem[i] + mod) % mod;
    for (int i = m - 1; i < k; i++) rem[i] = 0;
    for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
}

多项式牛顿迭代

需要简单的微积分知识，只需要知道求导和积分的算法即可。

首先定义多项式的求导：\(F^{\prime}(x) = \sum\limits_{i = 1}^n a_{i - 1} i x^i\)

void diffp(ll *f, ll *der, int n)
{
    for (int i = 1; i < n; i++) der[i - 1] = f[i] * i % mod;
    der[n - 1] = 0;
}

以及多项式的积分：\(\int F(x) dx = \sum\limits_{i = 0}^{n - 1} \frac{a_i x^{i + 1}}{i}\)

void intep(ll *f, ll *inte, int n)
{
    for (int i = 1; i < n; i++)
        if (!inv[i]) inv[i] = (i == 1 ? 1 : inv[mod % i] * (mod - mod / i) % mod);
    for (int i = 1; i < n; i++) inte[i] = f[i - 1] * inv[i] % mod;
    inte[0] = 0;
}

也就是对每一项分别求导再求和，符合求导和积分的经典形式。

同时多项式的求导和积分也是互为逆运算。

牛顿迭代：已知函数 \(G(F(x)) = 0\)，求 \(F(x) \bmod x^n\)，这里 \(G(x)\) 通常为根据需要构造的函数。

结论：令 \(F_*(x) \equiv F(x) \pmod {x^{\frac{n}{2}}}\)，\(F(x) \equiv F_*(x) - \frac{G(F_*(x))}{G^{\prime}(F_*(x))} \pmod {x^n}\).

具体证明见多项式计数杂谈 - command_block 的博客

多项式对数函数

P4725 【模板】多项式对数函数（多项式 ln）

求：\(\ln(A(x)) \equiv B(x) \pmod {x^n}\).

两边求导得 \(ln^{\prime}(A(x)) A^{\prime}(x) \equiv B^{\prime}(x) \pmod {x^n}\)

又因为 \(ln^{\prime}(A(x)) = \frac{1}{A(x)}\)，所以得：

\(\frac{A^{\prime}(x)}{A(x)} \equiv B^{\prime}(x) \pmod {x^n}\)

重新积分得 \(B(x) \equiv \frac{A^{\prime}(x)}{A(x)} dx\)

void lnp(ll *f, ll *ln, int n)
{
    diffp(f, ft, n), invp(f, rt, n);
    int k = 1;
    while (k < (n << 1)) k <<= 1;
    NTT(ft, k), NTT(rt, k);
    for (int i = 0; i < k; i++) ft[i] = ft[i] * rt[i] % mod;
    INTT(ft, k);
    intep(ft, ln, n);
    for (int i = 0; i < k; i++) ft[i] = rt[i] = 0;
}

多项式指数函数

P4726 【模板】多项式指数函数（多项式 exp）

求 \(B(x) \equiv \exp(A(x)) \pmod {x^n}\).

考虑牛顿迭代。设 \(G(F(x)) = \ln(F(x)) - A(x)\)，那么 \(G(B(x)) \equiv 0 \pmod {x^n}\).

根据牛顿迭代结论有 \(B(x) \equiv B_*(x) - \frac{G(B_*(x))}{\frac{1}{B_*(x)}} = B_*(x) - G(B_*(x)) B_*(x) \pmod {x^n}\).

整理得 \(B(x) \equiv (1 - \ln(B_*(x) + A(x)) B_*(x) \pmod {x^n}\).

继续考虑倍增求。

边界为 \(B_0 = 1\)（\(\exp(0) = 1\)）

void expp(ll *f, ll *exp, int n)
{
    int k = 1;
    while (k < n) k <<= 1;
    exp[0] = 1;
    for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
    {
        for (int i = 0; i < len; i++) Fn[i] = exp[i];
        lnp(Fn, Rn, len);
        for (int i = 0; i < len; i++) Rn[i] = (f[i] - Rn[i] + mod) % mod;
        Rn[0] = (Rn[0] + 1) % mod;
        NTT(Fn, len << 1), NTT(Rn, len << 1);
        for (int i = 0; i < (len << 1); i++) Fn[i] = Fn[i] * Rn[i] % mod;
        INTT(Fn, len << 1);
        for (int i = 0; i < len; i++) exp[i] = Fn[i];
    }
    for (int i = 0; i < (k << 1); i++) Fn[i] = Rn[i] = 0;
    for (int i = n; i < k; i++) exp[i] = 0;
}

多项式快速幂

P5245 【模板】多项式快速幂

当 \(A_0 = 1\) 时，根据 \((A(x))^k = \exp(k \ln(A(x)))\) 直接做。

void powp(ll *f, ll *pw, int n, int k)
{
    lnp(f, fn, n);
    for (int i = 0; i < n; i++) fn[i] = fn[i] * k % mod;
    expp(fn, pw, n);
    for (int i = 0; i < n; i++) fn[i] = 0;
}

P5273 【模板】多项式幂函数（加强版）

\(A_0 \neq 1\) 时求不到 \(\exp\)，考虑把 \(A_0\) 化成 \(1\).

直接给整个多项式乘上 \(\frac{1}{A_0}\)，最后再乘回 \(A_0^k\).

如果遇到系数为 \(0\) 的项也没关系，直接把整个多项式移一下 /xk

不过因为是多项式乘法，最后 \(0\) 的个数会变成原来的 \(k\) 倍，可以直接特判。

void powp2(ll *f, ll *pw, int n, int k, int k2)
{
    int pos = 0;
    while ((pos < n) && (f[pos] == 0)) pos++;
    if (1ll * pos * k >= n) return;
    int len0 = pos * k, len = n - len0;
    ll inv = qpow(f[pos], mod - 2, mod), val = qpow(f[pos], k2, mod);
    for (int i = 0; i < len; i++) fn[i] = f[pos + i] * inv % mod;
    lnp(fn, rn, len);
    for (int i = 0; i < len; i++) rn[i] = rn[i] * k % mod, fn[i] = 0;
    expp(rn, fn, len);
    for (int i = 0; i < len0; i++) pw[i] = 0;
    for (int i = len0; i < n; i++) pw[i] = fn[i - len0] * val % mod;
    for (int i = 0; i < len; i++) fn[i] = rn[i] = 0;
}

全家桶模板：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <map>
#include <algorithm>
using namespace std;

typedef long long ll;

const int sz = 3e6 + 5;
const int mod = 998244353;
const int g = 3;

int n;
int rev[sz];
ll F[sz], G[sz], inv[sz], wp[sz];
ll Ft[sz], Rt[sz], ft[sz], rt[sz], Fn[sz], Rn[sz], fn[sz], rn[sz];
ll quo[sz], rem[sz];
map<ll, ll> mp;

void readk(int &k, int &k2)
{
    k = k2 = 0;
    char ch = getchar();
    while ((ch < '0') || (ch > '9')) ch = getchar();
    while ((ch >= '0') && (ch <= '9'))
    {
        ll tmp = 10ll * k + ch - '0', tmp2 = 10ll * k2 + ch - '0';
        k = tmp % mod, k2 = tmp2 % (mod - 1);
        if (tmp >= mod) k += mod;
        if (tmp2 >= mod - 1) k2 += mod - 1;
        ch = getchar();
    }
}

void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }

ll qpow(ll base, ll power, ll mod)
{
    ll res = 1;
    while (power)
    {
        if (power & 1) res = res * base % mod;
        base = base * base % mod;
        power >>= 1;
    }
    return res;
}

ll gcd(ll a, ll b) { return (!b ? a : gcd(b, a % b)); }

ll bsgs(ll a, ll n, ll p, ll ad)
{
    ll m = ceil(sqrt(p));
    for (ll i = 0, t = n; i <= m; i++, t = t * a % p) mp[t] = i;
    for (ll i = 0, pw = qpow(a, m, p), t = ad; i <= m; i++, t = t * pw % p)
        if (mp.count(t) && (i * m - mp[t] >= 0)) return i * m - mp[t];
    return -1;
}

ll exbsgs(ll a, ll n, ll p)
{
    a %= p, n %= p;
    if ((n == 1) || (p == 1)) return 0;
    if (a == 0) return (n == 0 ? 1 : -1);
    ll cnt = 0, ad = 1, d;
    while ((d = gcd(a, p)) != 1)
    {
        if (n % d != 0) return -1;
        cnt++, n /= d, p /= d, ad = ad * (a / d) % p;
        if (ad == n) return cnt;
    }
    ll ans = bsgs(a, n, p, ad);
    if (ans == -1) return -1;
    return ans + cnt;
}

ll mod_sqrt(ll x)
{
    ll base = exbsgs(g, x, mod);
    if ((base == -1) || (base & 1)) return -1;
    return qpow(g, base / 2, mod);
}

void NTT(ll *A, int n)
{
    calc_rev(n);
    for (int i = 1; i < n; i++)
        if (rev[i] > i) swap(A[i], A[rev[i]]);
    for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
    {
        ll wn = qpow(g, (mod - 1) / len, mod);
        wp[0] = 1;
        for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
        for (int l = 0, r = len - 1; r <= n; l += len, r += len)
        {
            int w = 0;
            for (int p = l; p < l + m; p++, w++)
            {
                ll x = A[p], y = wp[w] * A[p + m] % mod;
                A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
            }
        }
    }
}

void INTT(ll *A, int n)
{
    NTT(A, n);
    reverse(A + 1, A + n);
    int inv = qpow(n, mod - 2, mod);
    for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;
}

void invp(ll *f, ll *r, int n)
{
    int k = 1;
    while (k < n) k <<= 1;
    r[0] = qpow(f[0], mod - 2, mod);
    for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
    {
        for (int i = 0; i < len; i++) Rt[i] = r[i], Ft[i] = f[i];
        NTT(Ft, len), NTT(Rt, len);
        for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
        INTT(Rt, len);
        for (int i = 0; i < m; i++) Rt[i] = 0; Rt[0] = 1;
        for (int i = 0; i < len; i++) Ft[i] = r[i];
        NTT(Ft, len), NTT(Rt, len);
        for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
        INTT(Rt, len);
        for (int i = m; i < len; i++) r[i] = (r[i] * 2ll - Rt[i] + mod) % mod;
    }
    memset(Ft, 0, k * sizeof(ll));
    memset(Rt, 0, k * sizeof(ll));
    for (int i = n; i < k; i++) r[i] = 0;
}

void divp(ll *f, ll *g, ll *quo, ll *rem, int n, int m)
{
    int len = n - m + 1;
    for (int i = 0; i < min(len, m); i++) rt[i] = g[m - i - 1];
    invp(rt, ft, len);
    for (int i = 0; i < len; i++) rt[i] = f[n - i - 1];
    int k = 1;
    while (k < (len << 1)) k <<= 1;
    NTT(rt, k), NTT(ft, k);
    for (int i = 0; i < k; i++) rt[i] = rt[i] * ft[i] % mod;
    INTT(rt, k);
    for (int i = 0; i < len; i++) quo[i] = rt[len - i - 1];
    for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
    for (int i = 0; i < m; i++) rt[i] = g[i];
    for (int i = 0; i < len; i++) ft[i] = quo[i];
    k = 1;
    while (k < (m + len)) k <<= 1;
    NTT(rt, k), NTT(ft, k);
    for (int i = 0; i < k; i++) rem[i] = rt[i] * ft[i] % mod;
    INTT(rem, k);
    for (int i = 0; i < m - 1; i++) rem[i] = (f[i] - rem[i] + mod) % mod;
    for (int i = m - 1; i < k; i++) rem[i] = 0;
    for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
}

void sqrtp(ll *f, ll *s, int n)
{
    s[0] = mod_sqrt(f[0]);
    if (s[0] == -1) return;
    int k = 1;
    while (k < (n << 1)) k <<= 1;
    ll inv2 = qpow(2, mod - 2, mod);
    for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
    {
        invp(s, rt, m);
        for (int i = 0; i < m; i++) ft[i] = f[i];
        NTT(ft, len), NTT(rt, len);
        for (int i = 0; i < len; i++) rt[i] = rt[i] * ft[i] % mod;
        INTT(rt, len);
        for (int i = 0; i < m; i++) rt[i] = (rt[i] + s[i]) % mod;
        for (int i = 0; i < len; i++) s[i] = rt[i] * inv2 % mod, rt[i] = ft[i] = 0;
    }
    for (int i = 0; i < k; i++) rt[i] = ft[i] = 0;
    for (int i = n; i < k; i++) s[i] = 0;
}

void diffp(ll *f, ll *der, int n)
{
    for (int i = 1; i < n; i++) der[i - 1] = f[i] * i % mod;
    der[n - 1] = 0;
}

void intep(ll *f, ll *inte, int n)
{
    for (int i = 1; i < n; i++)
        if (!inv[i]) inv[i] = (i == 1 ? 1 : inv[mod % i] * (mod - mod / i) % mod);
    for (int i = 1; i < n; i++) inte[i] = f[i - 1] * inv[i] % mod;
    inte[0] = 0;
}

void lnp(ll *f, ll *ln, int n)
{
    diffp(f, ft, n), invp(f, rt, n);
    int k = 1;
    while (k < (n << 1)) k <<= 1;
    NTT(ft, k), NTT(rt, k);
    for (int i = 0; i < k; i++) ft[i] = ft[i] * rt[i] % mod;
    INTT(ft, k);
    intep(ft, ln, n);
    for (int i = 0; i < k; i++) ft[i] = rt[i] = 0;
}

void expp(ll *f, ll *exp, int n)
{
    int k = 1;
    while (k < n) k <<= 1;
    exp[0] = 1;
    for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
    {
        for (int i = 0; i < len; i++) Fn[i] = exp[i];
        lnp(Fn, Rn, len);
        for (int i = 0; i < len; i++) Rn[i] = (f[i] - Rn[i] + mod) % mod;
        Rn[0] = (Rn[0] + 1) % mod;
        NTT(Fn, len << 1), NTT(Rn, len << 1);
        for (int i = 0; i < (len << 1); i++) Fn[i] = Fn[i] * Rn[i] % mod;
        INTT(Fn, len << 1);
        for (int i = 0; i < len; i++) exp[i] = Fn[i];
    }
    for (int i = 0; i < (k << 1); i++) Fn[i] = Rn[i] = 0;
    for (int i = n; i < k; i++) exp[i] = 0;
}

void powp(ll *f, ll *pw, int n, int k)
{
    lnp(f, fn, n);
    for (int i = 0; i < n; i++) fn[i] = fn[i] * k % mod;
    expp(fn, pw, n);
    for (int i = 0; i < n; i++) fn[i] = 0;
}

void powp2(ll *f, ll *pw, int n, int k, int k2)
{
    int pos = 0;
    while ((pos < n) && (f[pos] == 0)) pos++;
    if (1ll * pos * k >= n) return;
    int len0 = pos * k, len = n - len0;
    ll inv = qpow(f[pos], mod - 2, mod), val = qpow(f[pos], k2, mod);
    for (int i = 0; i < len; i++) fn[i] = f[pos + i] * inv % mod;
    lnp(fn, rn, len);
    for (int i = 0; i < len; i++) rn[i] = rn[i] * k % mod, fn[i] = 0;
    expp(rn, fn, len);
    for (int i = 0; i < len0; i++) pw[i] = 0;
    for (int i = len0; i < n; i++) pw[i] = fn[i - len0] * val % mod;
    for (int i = 0; i < len; i++) fn[i] = rn[i] = 0;
}

int main()
{
    return 0;
}

posted @ 2023-01-24 20:20 kymru 阅读(84) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

Loading

黄粱酒尚温，孰醒孰梦深