User:DQ9911

多项式全家桶

快速傅立叶变换

定义单位根:${\textstyle w_{n}^{k}=e^{\frac {2k\pi i}{n}}}$ <br>DFT: ${\textstyle y_{k}=\displaystyle \sum _{i=0}^{n-1}a_{i}(w_{n}^{k})^{i}}$ <br>IDFT: ${\textstyle a_{k}=\displaystyle {\frac {1}{n}}\sum _{i=0}^{n-1}y_{i}(w_{n}^{-k})^{i}}$

一个特殊的技巧：

设${\textstyle C(x)=A(x)+iB(x),D(x)=A(x)-iB(x)}$.

则有${\textstyle DFT(C)[i]={\overline {DFT(D)[n-i]}}}$

使用这个技巧可以将3次dft优化到2次dft, 但有精度问题，一定要预处理单位根

typedef double db;
typedef complex<db> cpd;
/* 约提速两倍
 * struct cpd {
 *     db x, y;
 *     cpd(db _x=0,db _y=0){x=_x,y=_y;}
 *     cpd &operator += (const cpd&rhs){return *this={x+rhs.x,y+rhs.y};}
 *     cpd &operator -= (const cpd&rhs){return *this={x-rhs.x,y-rhs.y};}
 *     cpd &operator *= (const cpd&rhs){return *this={x*rhs.x-y*rhs.y,x*rhs.y+y*rhs.x};}
 *     cpd &operator /= (const db&rhs){return *this={x/rhs,y/rhs};}
 *     cpd operator + (const cpd&rhs){cpd ret=*this;return ret+=rhs;}
 *     cpd operator - (const cpd&rhs){cpd ret=*this;return ret-=rhs;}
 *     cpd operator * (const cpd&rhs){cpd ret=*this;return ret*=rhs;}
 *     cpd operator / (const db&rhs){cpd ret=*this;return ret/=rhs;}
 *     friend cpd exp(const cpd&rhs){return cpd(cos(rhs.y),sin(rhs.y))*exp(rhs.x);}
 *     friend cpd conj(const cpd&rhs){return cpd(rhs.x, -rhs.y);}
 *     db real(){return x;}
 *     db imag(){return y;}
 * };
 */
const double pi = acos(-1.0);
const int N = 1 << 18 | 7;
const int L = 19; // lg(n) + 1
vector<cpd> w[2][L];
void fft_init() {
    for(int i = 1; i < L; i++) {
        w[0][i].resize(1 << i), w[1][i].resize(1 << i);
        for(int j = 0; j < 1 << i; j++) {
            w[0][i][j] = exp(cpd(0, j * 2 * pi / (1 << i)));
            w[1][i][j] = conj(w[0][i][j]);
        }
    }
}
void fft(cpd *y, int len, int on) { // on == 0: dfn, on == 1 : idft
    static int r[N], i, j, k, l, nl;
    static cpd u, v;
    l = __builtin_ctz(len);
    if(l != nl) for(i = 0, nl = l; i < len; i++)
        r[i] = (r[i>>1]>>1)|((i&1)<<l-1);
    for(i = 0; i < len; i++) 
        if(i<r[i]) swap(y[i],y[r[i]]);
    for(i = 1, l = 1; i < len; i<<=1, l++)
        for(j = 0; j < len; j += i<<1)
            for(k = j; k < j + i; k++)
                u=y[k], v=y[k+i]*w[on][l][k-j],
                y[k]=u+v, y[k+i]=u-v;
    if(on) for(i=0;i<len;i++) y[i]/=len;
}

任意模数FFT

先拆系数，直接暴力卷积是7次dft
${\textstyle A(x)=\sum _{i=0}^{n}(a_{i}\times 2^{15}+b_{i})x^{i}}$
${\textstyle B(x)=\sum _{i=0}^{n}(c_{i}\times 2^{15}+d_{i})x^{i}}$
只需要求${\textstyle a\times c}$、${\textstyle a\times d}$、${\textstyle b\times c}$、${\textstyle b\times d}$四个卷积
优化方法:
设 ${\textstyle P(x)=\displaystyle \sum _{i=0}^{n}(a_{i}+ib_{i})x^{i}}$、${\textstyle Q(x)=\displaystyle \sum _{i=0}^{n}(c_{i}+id_{i})x^{i}}$
使用上面提到共轭复数的技巧，可将dft次数优化到4次

const int M = (1<<15)-1;
const int mod = 1e6 + 3;
void conv(int *a, int *b, int *c, int l1, int l2) {
    static cpd dft[4], x[4][N*2];
    static int i, j, d[4];
    int len = l1 + l2 - 1;
    while(len & (len - 1)) len += len & -len;
    for(i = 0; i < len; i++) {
        x[0][i] = i<l1? cpd(a[i] & M, a[i] >> 15):0;
        x[1][i] = i<l2? cpd(b[i] & M, b[i] >> 15):0;
    }
    fft(x[0], len, 0), fft(x[1], len, 0);
    const cpd half = cpd(0.5, 0);
    for(i = j = 0; i < len; i++, j = len - i) {
        dft[0] = x[0][i] * x[1][i] * half;
        dft[1] = conj(x[0][j]) * x[1][i] * half;
        dft[2] = x[0][i] * conj(x[1][j]) * half;
        dft[3] = conj(x[0][j]) * conj(x[1][j]) * half;
        x[2][i] = dft[0] + dft[1];
        x[3][i] = dft[2] - dft[3];
    }
    fft(x[2], len, 1), fft(x[3], len, 1);
    for(i = 0; i < len; i++) {
        d[0] = ll(x[2][i].x + 0.5) % mod;
        d[1] = ll(x[2][i].y + 0.5) % mod;
        d[2] = ll(x[3][i].y + 0.5) % mod;
        d[3] = ll(x[3][i].x + 0.5) % mod;
        c[i] = ((1ll<<30) * d[3] % mod + d[0]) % mod;
        c[i] = (c[i] + (1ll<<15) * (d[1] + d[2])) % mod;
    }
}

快速数论变换

G为模P意义下的一个原根，此时用${\textstyle G^{\frac {k(P-1)}{n}}}$来代替${\textstyle w_{n}^{k}}$

常用NTT模数及其原根：${\textstyle P=r\cdot 2^{k}+1}$,${\textstyle N}$不能超过${\textstyle 2^{k}}$

NTT的卷积结果为循环卷积, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "http://localhost:6011/en.wikipedia.org/v1/":): {\textstyle c[i] = \displaystyle \sum_{(j+k)\equiv (mod \space n)} a[j] \times b[k]}

P           r   k   G
104857601   25  22  3
167772161   5   25  3
469762049   7   26  3
998244353   119 23  3
1004535809  479 21  3
2013265921  15  27  31
2281701377  17  27  3

typedef long long ll;
const int N = 1 << 18 | 7; 
const int P = 479 << 21 | 1;
inline int sub(int x, int y) {
    return (x-=y)<0?x+P:x;
}
inline int add(int x, int y) {
    return (x+=y)>=P?x-P:x;
}
template<int mod>
int power_mod(int a, int b) {
    int ret = 1;
    for(a %= mod; b; b>>=1,a=(ll)a*a%mod)
        if(b&1) ret = (ll)ret*a%mod;
    return ret;
}
template<int G, int P>
struct NTT {
    NTT() { ntt_init();}
    int w[19][2];
    void ntt_init() {
        for(int i = 1; i < 19; i++) {
            w[i][0] = power_mod<P>(G, P-1>>i);
            w[i][1] = power_mod<P>(w[i][0], P-2);
        }
    }
    void ntt(int *y, int len, int on) {
        static int r[N], nl, ww, wn, u, v;
        int i, j, k, l = __builtin_ctz(len) - 1;
        if(nl != len) {
            for(i = 0, nl = len; i < len; i++)
                r[i] = (r[i>>1]>>1)|(i&1)<<l;
        }
        for(i = 0; i < len; i++)
            if(i < r[i]) swap(y[i], y[r[i]]);
        for(i = 1, l = 1; i < len; i <<= 1, l++)
            for(j = 0, wn = w[l][on]; j < len; j+=i<<1)
                for(k = j, ww = 1; k < j + i; k++, ww = (ll)ww*wn%P)
                    u = y[k], v = (ll)y[k+i]*ww%P,
                    y[k] = add(u,v), y[k+i] = sub(u,v);
        if(on) {
            int invl = power_mod<P>(len, P - 2);
            for(i = 0; i < len; i++) y[i] = (ll)y[i]*invl%P;
        }
    }
};
NTT<3, P> ntt;

多项式求逆

Failed to parse (unknown function "\space"): {\textstyle A(x)B(x) = 1 (mod \space x^n)}
设Failed to parse (unknown function "\space"): {\textstyle A(x)B_0(x) = 1 (mod \space x^n), A(x)B(x) = 1 (mod \space x^{2n})}

则有Failed to parse (unknown function "\space"): {\textstyle A(x)(B(x) - B_0(x)) = 0 (mod \space x^n)} ，平方得到Failed to parse (unknown function "\space"): {\textstyle B^2(x) - 2B_0(x)B(x) + B_0^2(x) = 0 (mod \space x^{2n})}

两边同时乘${\textstyle A(x)}$, Failed to parse (unknown function "\space"): {\textstyle B(x) = B_0(x)(2 - A(x)B_0(x)) (mod \space x^{2n})} , 使用牛顿迭代法复杂度为${\textstyle O(nlogn)}$

void invpoly(int *A, int *B, int len) {
    static int x[N], i;
    if(len == 1) return B[0] = 1, void(B[1] = 0);
    invpoly(A, B, len>>1);
    for(i = 0; i < len; i++) x[i] = A[i], x[i+len]=B[i+len]=0;
    ntt.ntt(x, len<<1, 0), ntt.ntt(B, len<<1, 0);
    for(i = 0; i < 2 * len; i++)
        B[i] = B[i] * (2 + (ll)(P-B[i])*x[i]%P)%P;
    ntt.ntt(B, len<<1, 1);
    for(i = len; i < 2 * len; i++) B[i] = 0;
}

多项式取ln

设${\textstyle B(x)=ln(A(x))}$,两边求导${\textstyle B'(x)={\frac {A'(x)}{A(x)}}}$,求${\textstyle A(x)}$的逆元即可, 复杂度${\textstyle O(nlogn)}$

void logpoly(int *A, int *B, int len) {
    static int x[N], i;
    invpoly(A, x, len);
    for(i = 0; i < len; i++) B[i] = A[i+1]*(i+1ll)%P;
    for(i = len-1; i < 2 * len; i++) B[i] = 0;
    ntt.ntt(x, len<<1, 0), ntt.ntt(B, len<<1, 0);
    for(i = 0; i < 2 * len; i++)
        B[i] = (ll)x[i]*B[i]%P;
    ntt.ntt(B, len<<1, 1);
    for(i = len; i < 2 * len; i++) B[i] = 0;
    for(i = len-1; i; i--) B[i] = (ll)B[i-1]*inv[i]%P;
    B[0] = 0;
}

多项式exp

cdq分治求法:

设${\textstyle B(x)=e^{A(x)}}$, 则${\textstyle B'(x)=A'(x)B(x)}$

设${\textstyle B(x)=\displaystyle \sum _{i=0}^{n}b_{i}x^{i}}$, ${\textstyle A(x)=\displaystyle \sum _{i=0}^{n}a_{i}x^{i}}$.

则有${\textstyle b_{k}=\displaystyle {\frac {1}{k}}\sum _{i=1}^{k}a_{i}b_{k-i}}$

下面做法: 设${\textstyle g(f(x))=lnB(x)-A(x)=0}$,

利用在${\textstyle f_{0}(x)}$出的泰勒展开式得出模${\textstyle x^{n}}$与模${\textstyle x^{2n}}$时${\textstyle B(x)}$的递推式，再牛顿迭代.

复杂度: ${\textstyle O(nlogn)}$（常数堪比${\textstyle nlog^{2}n}$慎用）

void exppoly(int *A, int *B, int len) {
    static int x[N], i;
    if(len == 1) return B[0] = 1, void(B[1] = 0);
    exppoly(A, B, len>>1);
    logpoly(B, x, len); x[0] = sub(x[0], 1);
    for(i = 0; i < len; i++) x[i] = sub(A[i], x[i]);
    for(i = len; i < 2 * len; i++) x[i] = 0;
    ntt.ntt(B, len<<1, 0), ntt.ntt(x, len<<1, 0);
    for(i = 0; i < 2 * len; i++)
        B[i] = (ll)B[i] * x[i] % P;
    ntt.ntt(B, len<<1, 1);
    for(i = len; i < 2 * len; i++) B[i] = 0;
}

多项式求幂

${\textstyle f^{k}(x)=e^{klnf(x)}}$

先求${\textstyle lnf(x)}$, 再求${\textstyle e^{kf(x)}}$即可. 细节略.

多项式除法

${\textstyle A(x)=B(x)C(x)+D(x)}$,设A, B, C, D分别为Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "http://localhost:6011/en.wikipedia.org/v1/":): {\textstyle n、m, n - m, m - 1} 阶多项式

${\textstyle x^{n}A({\frac {1}{x}})=x^{m}B({\frac {1}{x}})x^{n-m}C({\frac {1}{x}})+x^{m-1}D({\frac {1}{x}})x^{n-m+1}}$

即${\textstyle rev(A(x))=rev(B(x))rev(C(x))+x^{n-m+1}rev(D(x))}$

在模${\textstyle x^{n-m+1}}$意义下计算${\textstyle B(x)}$逆元即可

可持久化字典树

用于求解异或最大值问题

通过记录每个节点插入次数实现可持久化

struct Trie {
    int cnt[N * 31], c[N * 31][2], tot;
    void init() { 
        tot = 0; 
        c[0][0]=c[0][1]=0;
    }
    int newnode() {
        tot++;
        c[tot][0] = c[tot][1] = 0;
        return tot;
    }
    void insert(int last, int &now, int x) { 
    /* last表示上个版本的根节点，now表示修改后版本的根节点 */
        int tmp = now = newnode();
        for(int i = 29; ~i; i--) {
            int k = x>>i&1;
            c[tmp][k] = newnode(), c[tmp][k^1] = c[last][k^1];
            tmp = c[tmp][k], last = c[last][k];
            cnt[tmp] = cnt[last] + 1;
        }
    }
    int query(int l, int r, int x) { //查询区间异或最大值
        for(int i = 29; ~i; i--) {
            int k = ~x>>i&1;
            if(cnt[c[r][k]] > cnt[c[l][k]])
                x ^= k << i, r = c[r][k], l = c[l][k];
            else x ^= (k^1)<<i, r = c[r][k^1], l = c[l][k^1];
        }
        return x;
    }
} per_trie;

拉格朗日插值

给一个${\textstyle n}$阶多项式的${\textstyle n+1}$个点${\textstyle (x_{0},f(x_{0})),(x_{1},f(x_{1})),(x_{2},f(x_{2})),...,(x_{n},f(x_{n}))}$，且${\textstyle x_{i}\neq x_{j}}$

现给一个${\textstyle x}$的值，要求计算${\textstyle f(x)}$的值。

构造多项式如下

${\textstyle f(x)=\displaystyle \sum _{i=0}^{n}f(x_{i})\prod _{j\neq i}{\frac {x-x_{j}}{x_{i}-x_{j}}}}$

只需${\textstyle O(n^{2})}$预处理, ${\textstyle O(n)}$求值即可.

void Lagrange(int *x, int *y, int *a, int n) {
    for(int i = 0; i <= n; i++) {
        a[i] = 1;
        for(int j = 0; j <= n; j++)
            if(j != i) a[i]=(ll)a[i]*(x[i]-x[j]+mod)%mod;
        a[i] = (ll)power_mod(a[i], mod-2)*y[i]%mod;
    }
}
int calc(int *a, int *x, int *y, int n, int v) {
    static int t1[N], invi[N]; // invi[i]:inv(v-x[i])
    for(int i=0; i<=n; i++)
        if(v==x[i]) return y[i];
/* 这段可以不写下面直接求逆元，复杂度多个lg，实际影响不大
    t1[0] = (v-x[0]+mod)%mod;
    for(int i = 1; i <= n; i++)
        t1[i] = (ll)t1[i-1]*(v-x[i]+mod)%mod;
    invi[n] = power_mod(t1[n], mod-2);
    for(int i = n-1; ~i; i--)
        invi[i] = (ll)invi[i+1]*(v-x[i+1]+mod)%mod;
    for(int i = 1; i <= n; i++)
        invi[i] = (ll)invi[i]*t1[i-1]%mod;
*/
    int ret = 0;
    for(int i = 0; i <= n; i++)
        ret = (ret + (ll)t1[n]*invi[i]%mod*a[i])%mod;
    return ret;
}