多项式基础操作全家桶

多项式入门教程

基础概念

形如 $F(x)=a_1+a_2x+a_3x^2+...+a_nx^n$ 的东西叫多项式。

然后，这个 $n$ 可以是无穷大的。

其中上面的 $a_i$ 称为多项式的第 $i$ 项系数。

多项式乘法

即各项相乘，假设有两个多项式 $F(x)=\sum_{i=0}^na_ix^i$ 和 $G(x)=\sum_{i=0}^nb_ix^i$

那么：

$$F(x)\times G(x)=\sum_{i=0}^{n+m} (\sum_{x+y=i} a_xb_y)x^i$$

然后我们发现这个式子里有一个部分很有意思，我们专门拿出来：

设两多项式的乘积的第 $i$ 项系数为 $c_i$ 那么有：

$$c_i=\sum_{x+y=i}a_xb_y$$

我们称这个为加法卷积。如果把那个加号换成别的就变成其他的什么卷积了。

那么我们如何计算多项式乘法？你发现可以容易的在 $O(n^2)$ 内计算。但是这是不优的，我们接下来介绍一种 $O(n\log n)$ 内计算多项式乘法的算法，先了解一些前置知识。

复数

什么叫复数？形如 $a+bi$ 的叫复数。其中 $i^2=-1$，$a$ 称为一个复数的实部，$b$ 称为一个复数的虚部。$|a+bi|$ 表示其模长。其共轭复数为 $a-bi$，对于复数 $z$ ，其共轭复数记作 $\overline{z}$。

然后其四则运算为：

$$(a+bi)+(c+di)=(a+c)+(b+d)i\\ (a+bi)\times(c+di)=ac+adi+bci-bd=(ac-bd)+(ad+bc)i$$

然后模长 $|z|=|a+b_i|=\sqrt{z\overline{z} }=\sqrt{a^2+b^2}$

$a+bi\not=0$ 时存在逆元 $\frac{1}{a+bi}=\frac{a-bi}{|a+b_i|}$

复数几何意义

我们把 $x$ 坐标表示实部， $y$ 坐标表示虚部的平面称为复平面。可以发现是数轴的扩展。

可以把复数看做一个在复平面上的向量 $(a,b)$，记其幅角为 $\theta$ ，模长为 $r$ 那么有 $a=r\cos{\theta}$ 且 $b=r\sin{\theta}$，那么其实复数还可表示为 $r(\cos{\theta}+i\sin{\theta})$ 。

然后你会惊奇的发现复数相乘相当于模长相乘，幅角相加。

$$r_1(\cos{\theta_1}+i\sin{\theta_1})\times r_2(\cos{\theta_2}+i\sin{\theta_2})\\ =r_1r_2(\cos{\theta_1}\cos{\theta_2}+i\cos{\theta_1}\sin{\theta_2}+i\sin{\theta_1}\cos{\theta_2}-\sin{\theta_1}\sin{\theta_2})\\ =r_1r_2(\cos{(\theta_1+\theta_2)}+i\sin{(\theta_1+\theta_2)})$$

单位根

我们称满足 $z^n=1$ 的复数 $z$ 为 $n$ 次单位根。

首先由于模长相乘，幅角相加的乘法，我们知道 $n$ 次单位根模长为 $1$，$n\theta=2k\pi$。因此不同的 $n$ 次单位根总共有 $n$ 个，它们将单位圆 $n$ 等分。

一般我们用 $w^i_n$ 表示第 $i+1$ 个 $n$ 次单位根，即 $\theta_n^i=\frac{2i\pi}{n}$ 。

一些性质

性质1

$$w^k_n=w^{k+n}_n$$

性质2

$$w^{ak}_{an}=w^k_n$$

角度所占比相同，模长都为 $1$ ，因此固然相同。

性质3

$$w^i_n\times w^j_n=w^{i+j}_{n}$$

单位根可以看做 $\frac{1}{i}$ 圆周。

性质4

$$w^{k+\frac{n}{2}}_n=-w^k_n$$

旋转了一个半圆，因此相反。

性质5

$$\sum_{i=0}^{n-1}w_n^{ai}=n[a=0]$$

证明：

当 $a=0$ 时，等价于 $n$ 个第 $1$ 个 $n$ 次单位根相加，就是 $n$ 个 $1$ 相加，于是得到 $n$。

当 $a\not=0$，

先了解一下等比数列求和公式

对于 $a,aq,aq^2,...,aq^n$ 求 $S_n=a+aq+aq^2+...+aq^n$

其等价于 $af(q),f(x)=1+x+x^2+...+x^n$

我们可以知道 $f(x)=xf(x)+1-x^{n+1}$ 即 $f(x)=\frac{1-x^{n+1} }{1-x}$

那么 $S_n=af(q)=a\frac{1-q^{n+1} }{1-q}$

好，那么把 $w_n^0$ 代入 $a$ ，$w_n^a$ 代入 $q$ 得到

$$原式=w^0_n\frac{1-(w_n^a)^{n}}{1-w_n^a}=\frac{1-w^0_n}{1-w^a_n}=0$$

Q.E.D.

离散傅里叶变换

简称 DFT (Discrete Fourier Transform)。

其逆变换称为 IDFT 。

所以实际上
我们知道， $n$ 个方程可以确定 $n$ 个未知数。也就是说对于一个 $n$ 次多项式，我们只需要 $n+1$ 个点值就可以求出多项式的系数。点值是什么？将 $x$ 的取值代入多项式得到的值。

然后经过简单验证可以发现，$F,G$ 在 $x$ 处的点值的乘积，等于 $F\times G$ 在 $x$ 处的点值。

因为我们发现多项式在点值上相乘更快，所以我们考虑一种系数表示转点值表示再转回系数表示的做法。

而这个东西的作用和点值表示有一些相同的性质，以及一个更为优秀的性质。因为傅里叶变换的本质是把单位根代入多项式。

回到正题，傅里叶变换。

假设对于一个多项式 $f$ 我们知道其在 $0$ 到 $n-1$ 处的系数。

那么定义其傅里叶变换 $F[f]$ 满足：（两式的单位根的次方的正负号可以交换）

$$F[f](k)=\sum_{i=0}^{n-1}f[i]w_n^{ik}\\ F^{-1}[f](k)=\frac{1}{n}\sum_{i=0}^{n-1}f[i]w_{n}^{-ik}$$

先证明一下这两个互为逆运算，即证明：

$$F^{-1}[F[f] ](k)=f[k]\\ F[F^{-1}[f] ](k)=f[k]$$

证明：

怎么证明？全部展开！

$$F^{-1}[F[f] ](k)\\ =\frac{1}{n}\sum_{i=0}^{n-1}F[f](i)w^{-ik}_n\\ =\frac{1}{n}\sum_{i=0}^{n-1}w^{-ik}_n\sum_{j=0}^{n-1}f[j]w^{ji}_n\\ =\frac{1}{n}\sum_{j=0}^{n-1}f[j]\sum_{i=0}^{n-1}w^{i(j-k)}_n\\ =f[k]$$

同理另一个得证。

Q.E.D.

然后实际上，你发现：

$$F[f](k)=f(w^{k}_n)\\ F^{-1}[f](k)=\frac{1}{n}f(w^{-k}_n)\\ $$

卷积定理

$$f\times g[x]=\sum_{i+j=x}f[i]g[j] \quad i,j\in[0,n-1]$$

然后有：

$$F[f\times g](x)=(F[f]\times F[g])(x)\\ $$

证明：

$$F[f\times g](x)=\sum_{i=0}^{n-1}(f\times g)[i]w^{ix}\\ =\sum_{i=0}^{n-1}\sum_{a+b=i} f[a]g[b]w^{ix}_n\\ =\sum_{i=0}^{n-1}\sum_{a+b=i} f[a]w^{ax}_ng[b]w^{bx}_n\\ =(\sum_{i=0}^{n-1}f[i]w^{ix}_n)(\sum_{i=0}^{n-1}g[i]w^{ix}_n)\\ =(F[f]\times F[g])(x)$$

也就是傅里叶变换的先后是无关的，所以可以用傅里叶变换的乘积进行逆傅里叶变换得到乘积。

快速傅里叶变换

简称 FFT (Fast Fourier Transform)，用来快速计算傅里叶变换。

为方便处理，我们把 $n$ 取到大于等于自身的最小的 $2$ 的幂次。（不用担心会计算错误，因为不存在的高位系数都是 $0$）。

记 $m=\frac{n}{2}$

设：

$$F_0[f](k)=\sum_{i=0}^{m-1} f[2i]w^{2ik}_n=\sum_{i=0}^{m-1} f[2i]w^{ik}_m\\ F_1[f](k)=\sum_{i=0}^{m-1} f[2i+1]w^{(2i+1)k}_n=w^{k}_n\sum_{i=0}^{m-1} f[2i+1]w^{ik}_m$$

然后有 $F[f](k)=F_0[f](k)+F_1[f](k)$

观察一下可以发现，$F_0[f](k+m)=F_0[f](k)$ 且 $F_1[f](k+m)=w_n^{m}F_1[f](k)$ 且 $w^{k+m}_n=-w^{k}_n$ 。（三个结论暴力拆解易证）

所以 $F_0$ 和 $F_1$ 都只计算前 $m$ 项，这个分治处理，边界是 $m-1=1$，此时两式都是原多项式。

然后复杂度 $T(n)=2T(\frac{n}{2})+O(n)$ 于是总复杂度 $O(n\log n)$

然后再观察一下 $F$ 和 $F^-1$:

$$F[f](k)=\sum_{i=0}^{n-1}f[i]w_n^{-ik}\\ F^{-1}[f](k)=\frac{1}{n}\sum_{i=0}^{n-1}f[i]w_{n}^{ik}$$

如果令 $w'_n=w^{-1}_n$ ：

$$F^{-1}[f](k)=\frac{1}{n}\sum_{i=0}^{n-1}f[i]w_{n}^{ik}$$

发现和 $F$ 长得其实巨像无比，所以我们只需要将 $w'_n$ 作为单位根进行傅里叶变换就可以实现逆傅里叶变换。

代码实现

由于考虑到代码常数的问题，我们就不考虑实现递归版了。因为这个理解和实现起来都简单。这个递归版的实现可以看 这里

为了保证一个能够过题的常数，我们考虑以下的一些优化。

别的没什么好说的，我们重点说一下 FFT 的二进制翻转优化。

先给出结论，递归到最底层时，位置为 $i$ 处理的是在原序列位置为 $rev(i)$ 的地方，其中 $rev(i)$ 表示 $i$ 的二进制翻转。

我们直接把序列放到这个位置便是对的，因为我们只在乎递归边界的值的大小，中间的位置我们不需要在意。

那么为什么递归的最底层满足位置的二进制翻转？

证明：

考虑归纳法证明。

当 $n=1$，有 $rev(0)=0$ 成立。

当 $n=2$，有 $rev(0)=0,rev(1)=1$ 成立。

当 $n=4$， 有 $rev(0)=0,rev(1)=2,rev(2)=1,rev(3)=3$ 成立。

当 $n=8$，有：

0 1 2 3 4 5 6 7 第1层
0 2 4 6|1 3 5 7 第2层
0 4|2 6|1 5|3 7 第3层
0|4|2|6|1|5|3|7 第4层

成立。

假设结论对于 $n=2^k$ 成立，证明对于 $n=2^{k+1}$ 也成立。

有

$$0~1~2~...~2^{k+1}-1\quad 第一层\\ 0~2~...~2^{k+1}-2\mid 1~3~...~2^{k+1}-1\quad 第二层\\ $$

左边的位置全部除去二进制下最低位的 $0$ 得到

$$0~1~...~2^k-1$$

我们知道，对于 $n=2^k$ 结论成立，于是左边重新加回最低位 $0$ ，发现此操作前后二进制翻转不变，因此左半边结论成立。

右半边除去最低位 $1$ 得到相同的式子。

由于对于 $n=2^k$ 成立，右边加回最低位 $1$ ，发现操作前后二进制翻转相差 $2^k$，减去左半边下标 $2^k$ 正好对应右半边的二进制翻转。于是右半边结论成立。

综上结论成立。Q.E.D.

那么我们直接把这位的多项式的值换到二进制翻转的位置计算即可。至于怎么线性求二进制翻转？看 这里

可能你有点疑惑，为什么递归最底层的值不需要操作可以直接用？你代入 $F[f](0)$ 你发现 $w^{0}=1$ ，直接就是原本的多项式。

#include<bits/stdc++.h>
#define ll long long
#define db double 
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define Better_IO true
namespace IO{
    #if Better_IO
        char buf[(1<<20)+3],*p1(buf),*p2(buf);
        const int lim=1<<20;
        inline char gc(){
            if(p1==p2) p2=(p1=buf)+fread(buf,1,lim,stdin);
            return p1==p2?EOF:*p1++;
        }
        #define pc putchar
    #else
        #define gc getchar
        #define pc putchar
    #endif
    inline bool blank(char c){
        return c=='\t' or c=='\n' or c=='\r' or c==' ' or c==EOF;
    }
    inline void gs(char *s){
        char ch=gc();
        while(blank(ch) ) ch=gc();
        while(!blank(ch) ) {*s++=ch;ch=gc();}
        *s=0;
    }
    inline void ps(char *s){
        while(*s!=0) {pc(*s++);}
    }
    template<class T>
    inline void read(T &s){
        s=0;char ch=gc();bool f=0;
        while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
        while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
        if(ch=='.'){
            db p=0.1;ch=gc();
            while('0'<=ch&&ch<='9') {s=s+p*(ch^48);ch=gc();p*=0.1;}
        }
        s=f?-s:s;
    }
    template<class T,class ...A>
    inline void read(T &s,A &...a){
        read(s);read(a...);
    }
};
using IO::gs;
using IO::ps;
using IO::read;
using IO::gc;
const int N=(1<<20)+3;
const db pi=acos(-1);
struct Complex{
    db a,b;
    Complex(){};
    Complex(db _a,db _b){
        a=_a;b=_b;
    }
    inline friend Complex operator + (Complex x,Complex y){
        return {x.a+y.a,x.b+y.b};
    }
    inline friend Complex operator - (Complex x,Complex y){
        return {x.a-y.a,x.b-y.b};
    }
    inline friend Complex operator * (Complex x,Complex y){
        return {x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a};
    }
};
int n,m;
Complex F[N<<1],G[N<<1],w[N];
int rev[N<<1];
inline void getrev(){
    for(int i=0;i<n;++i){
        rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
    }
}
inline void FFT(Complex *f,bool I){
    for(int i=0;i<n;++i){
        if(rev[i]<i)
            std::swap(f[i],f[rev[i] ]);
    }
    for(int p=2;p<=n;p<<=1){
        int len=p>>1;
        w[0]=Complex(1,0);w[1]=Complex(cos(2*pi/p),sin(2*pi/p) );
        if(I) w[1].b*=-1;
        for(int i=2;i<len;++i){
            w[i]=w[i-1]*w[1];
        }
        for(int k=0;k<n;k+=p){
            for(int l=k;l<k+len;++l){
                Complex tmp=f[l+len]*w[l-k];
                f[l+len]=f[l]-tmp;
                f[l]=f[l]+tmp;
            }
        }
    }
}
int main(){
    file(a);
    read(n,m);
    for(int i=0;i<=n;++i){
        read(F[i].a);
    }
    for(int i=0;i<=m;++i){
        read(G[i].a);
    }
    for(m+=n,n=1;n<=m;n<<=1);
    getrev();
    FFT(F,0);FFT(G,0);
    for(int i=0;i<n;++i) F[i]=F[i]*G[i];
    FFT(F,1);
    for(int i=0;i<=m;++i){
        printf("%d ",(int)(F[i].a/n+0.5) );
    }
    return 0;
}

快速数论变换

简称 NTT (Number-theoretic transform)。

很多时候计数问题都是在取模意义下的，而 FFT 使用浮点数，不单单是巨慢无比，更可怕的是基本上用不了了，怎么办？于是考虑选择一个模意义下的玩意来代替单位根。最终大佬们找到了 原根。然后就形成了 NTT。

原根与阶

阶

定义：$a$ 在模 $p$ 意义下的阶为一个最小的整数 $k$ 使得 $a^k\equiv 1\pmod p$。其中 $gcd(a,p)=1$。记作 $\delta_p(a)$。

本质上是幂次的最小循环节。

性质1

根据欧拉定理，我们发现一个数 $a$ 的阶会小于等于 $\varphi(a)$。

性质2

若 $a^n\equiv 1\pmod m$ , 则 $\delta_m(a)|n$

证明：

令 $n=\delta_m(a)q+r\quad (0\le r<\delta_m(a) )$

当 $r>0$ 有 $a^r\equiv a^r(a^{\delta_m(a)})^q\equiv 1\pmod m$

然后就是有一个比 $\delta_m(a)$ 更小的满足阶的性质的东西 $r$ ，于是矛盾。

于是 $r=0$ ，即 $\delta_m(a)|n$ ，得证，Q.E.D.

性质3

对于 $a,b,m$ ，且 $gcd(a,m)=gcd(b,m)=1$ 则：

$\delta_m(ab)=\delta_m(a)\delta_m(b)$ 的充分必要条件是 $gcd(\delta_m(a),\delta_m(b) )=1$

证明：

必要性证明：

由阶的定义可得

$$(ab)^{lcm(\delta_m(a),\delta_m(b) )}\equiv (a^{\delta_m(a)})^{k_1}(b^{\delta_m(b)})^{k_2} \equiv 1 \pmod m$$

由 性质2 可知

$$\delta_m(ab)|lcm(\delta_m(a),\delta_m(b) )$$

由 $\delta_m(ab)=\delta_m(a)\delta_m(b)$ 得：

$$\delta_m(a)\delta_m(b)|lcm(\delta_m(a),\delta_m(b) )$$

由公式 $ab=lcm(a,b)\gcd(a,b)$ 可知：

$gcd(\delta_m(a),\delta_m(b) )=1$

Q.E.D.

充分性证明：

$$1\equiv (ab)^{\delta_m(ab)\delta_m(b)}\equiv a^{\delta_m(ab)} \pmod m$$

可得 $\delta(a)|\delta(ab)\delta_m(b)$ 再结合 $gcd(\delta_m(a),\delta_m(b) )=1$ 得到 $\delta_m(a)|\delta_m(ab)$

同理可得： $\delta_m(b)|\delta_m(ab)$

因此 $\delta_m(a)\delta_m(b)|\delta_m(ab)$

然后类似地可以得到 $\delta_m(ab)|\delta_m(a)\delta_m(b)$

性质4

对于 $\gcd(a,m)=1$，有:

$$\delta_m(a^k)=\frac{\delta_m(a)}{\gcd(\delta(a),k)}$$

证明：

由：

$$(a^k)^{\delta_m(a^k)}\equiv 1 \pmod m$$

可得 $\delta_m(a)|k\delta_m(a^k)$ ，则 $\frac{\delta_m(a) }{\gcd(\delta_m(a),k)}|\delta_m(a^k)$

再由：

$$(a^k)^{\frac{\delta_m(a)}{\gcd(\delta_m(a),k)}}=(a^{\delta_m(a)})^{\frac{k}{\gcd(\delta_m(a),k)}}\equiv 1 \pmod m$$

于是有 $\delta_m(a^k)|\frac{\delta_m(a) }{\gcd(\delta_m(a),k)}$

然后把刚才这个式子与上面得到的式子结合，可得：

$$\delta_m(a^k)=\frac{\delta_m(a) }{\gcd(\delta_m(a),k)}$$

Q.E.D.

原根

定义：一个数 $a$ 在模 $p$ 意义下满足 $\delta_p(a)=\varphi(p)$ ，则称之为 $p$ 的原根。

证明 性质1 前，需要证明两个定理：

拉格朗日定理

对于素数 $p$，对于模 $p$ 意义下的整系数多项式 $f(x)$:

$$f(x)=a_0+a_1x^1+...+a^nx^n(p \not | a_n)$$

其同余方程 $f(x)\equiv 0 \pmod p$ 在模 $p$ 意义下至多有 $n$ 个不同解。

证明：

当 $n=0$ 时，有 $p\not|a_0$ ，因此 $f(x)\equiv 0\pmod p$ 是无解的，定理对于 $n=0$ 成立。

假设定理对于 $\deg f<n$ 成立，考虑证明定理对于 $\deg f=n$ 成立。

反证法，假设存在 $n+1$ 个符合条件的 $x_0,x_1,...,x_n$ 。

现在设 $g(x)$ 满足 $f(x)-f(x_0)=(x-x_0)g(x)$ ，因此可知 $g(x)$ 在模 $p$ 意义下最多是 $n-1$ 次多项式，又因为对于 $i\in[0,n]$ 都满足 $f(x_i)\equiv 0 \pmod p$ ，于是有：

对于 $i\in[1,n]$ 有，

$$(x_i-x_0)g(x_i)=f(x_i)-f(x_0)\equiv 0 \pmod p$$

又因为 $x_i\not \equiv x_0\pmod p$ ，得到 $g(x)\equiv 0 \pmod p$ ，此时 $g(x)\equiv 0\pmod p$ 这个同余方程有 $n$ 个根，与归纳法的假设冲突，矛盾。

因此定理对 $\deg f=n$ 也成立，得证。

Q.E.D.

定理1

对于 $a,b,p$ 满足 $gcd(a,p)=gcd(b,p)=1$ ,存在 $c\in Z$ 使得 $\delta_p(c)=lcm(\delta_p(a),\delta_p(b))$

证明：

对于 $\delta_p(a),\delta_p(b)$ 质因数分解，得到：

$$\delta_p(a)=\prod_{i=0}^sp_i^{n_i} \\ \delta_p(b)=\prod_{i=0}^sp_i^{m_i}$$

令 $k_1$ 为所有 $n_i\le m_i$ 的 $p^{n_i}$ 的乘积，$k_2$ 为所有 $n_i>m_i$ 的 $p_i^{m_i}$ 的乘积。

然后取 $x,y$ 使得 $\delta_p(a)=k_1x,\delta_P(b)=k_2y$。可以发现 $\gcd(x,y)=1$ 且 $lcm(\delta_p(a),\delta_p(b) )=xy$ 。

然后由 阶的性质4 ，$\delta_p(a^{k_1})=\frac{\delta_p(a)}{\gcd(\delta_p(a),k_1)}=x,\delta_p(b^{k_2})=y$ 。

再由 阶的性质3 可得：（条件：$\gcd(\delta_p(a),\delta_p(b))=1$ 且 $\gcd(a,p)=\gcd(b,p)=1$）

$$\delta_p(a^{k_1}b^{b_2})=xy=lcm(\delta_P(a),\delta_p(b) )$$

可得 $c=a^{k_1}b^{k_2}$ 。

Q.E.D.

性质1

对于任意奇素数（非 $2$ 素数）$p$ ，其存在原根。

证明：

由上面的 定理1 ，可以发现存在一个 $g$ 使得 $\delta_p(g)=lcm(\delta_p(1),...,\delta_p(p-1) )$

于是，对于任意 $i\in[1,p-1]$ 有 $\delta_p(i)|\delta_p(g)$ 即 $i^{\delta_p(g)}\equiv 1\pmod p$。也就是说，这些 $i$ 全部都是同余方程 $x^{\delta_p(g)}\equiv 1\pmod p$ 的解。

有 拉格朗日定理 可知，$\delta_p(g)\le p-1$ ，同时由费马小定理，有 $\delta_p(g)\ge p-1$ ，可得 $\delta_p(g)=p-1=\varphi(p)$ 。

所以 $g$ 是模 $p$ 意义下的原根。

Q.E.D.

性质2

对于奇素数 $p$ 和 $a\in N^*$ ，$p^a$ 有原根。

证明：

考虑先证明 $p^2$ 的原根存在。

记 $p$ 的一个原根为 $g$ ，其在模 $p^2$ 的阶为 $d$ ，则有 $d|\varphi(p^2)$ （欧拉定理+阶的性质2），再由 $g^d\equiv 1\pmod p$ 于是有 $g^d\equiv 1\pmod{p^2}$，因此 $\varphi(p)|d$。

由于 $\varphi(p)=(p-1)|d$ 且 $d|p(p-1)=\varphi(p^2)$ ，$p$ 为奇素数，可知 $d=p-1$ 或 $d=p(p-1)$ 。

当 $d=p-1$ ，则 $g$ 为 $p^2$ 的原根。当 $d=p(p-1)$ 时，考虑证明 $(g+p)$ 是 $p^2$ 的原根。

首先 $(g+p)$ 为 $p$ 的原根。同上可知 $g+p$ 在 $p^2$ 的阶为 $(p-1)$ 或 $p(p-1)$ 。由 $g^{p-1}\equiv 1\pmod{p^2}$ 结合二项式定理可得：

$$(g+p)^{p-1}\\ \equiv 1-pg^{p-2}\pmod{p^2}$$

可知 $(g+p)^{p-1}\not \equiv 1\pmod{p^2}$ ,则 $(g+p)^{p(p-1)}\equiv 1\pmod{p^2}$。

然后可以归纳证明：

假设 $g$ 是 $p^a$ 的原根 ，$d$ 是 $g$ 对于模 $p^{a+1}$ 的阶，考虑证明 $g$ 是 $p^{a+1}$ 。于是类似的可以得到 $d=p^{a-1}(p-1)$ 或 $d=p^a(p-1)$。

由 欧拉定理: $g^{p^{a-2}(p-1) }\equiv 1\pmod {p^{a-1}}$ ，即 $g^{p^{a-2}(p-1) }=1+kp^{a-1}$

由于 $g$ 为 $p^a$ 的原根,$d>p^{a-2}(p-1)$，有 $g^{p^{a-2}(p-1)}\not \equiv 1 \pmod{p^a}$

于是又有 $p^a\not |kp^{a-1}$ 即 $\gcd(k,p)=1$

然后二项式定理展开，有：

$$g^{p^{a-1}(p-1)}=(1+kp^{a-1})^p\\ =1+kp^a\not \equiv 1\pmod{p^{a+1} }$$

于是 $d=p^a(p-1)$ , $p^{a+1}$ 存在原根 $g$。

Q.E.D.

性质3

对于奇素数 $p$ 和 $a\in N^*$ ，$2p^a$ 有原根。

证明：

假设对于 $p^a$ 存在原根 $g$。首先我们知道对于 $2p^a$ ，它的原根必然是奇数的。而对于 $p^a$ 的原根 $g$ ，以及 $(g+p^a)$ ，其中至少有一个是奇数，我们记这个奇数为 $b$。

设 $b$ 在模 $2p^a$ 的意义下的阶为 $d$ ，于是有 $d|\varphi(2p^a)=\varphi(p^a)$（阶的性质2）且 $\varphi(p^a)|d$（阶的性质2），即 $d=\varphi(p^a)=\varphi(2p^a)$。

所以 $b$ 为模 $2p^a$ 的原根。

Q.E.D.

性质4

对于 $p\not = 2,4,p^c,2p^c$ ，其不存在原根。

证明：

分三种情况来讨论：

1.对于 $p=2^a,a\ge 3$

2.对于 $p$ 为 $2^a$ 与一个奇素数的乘积且 $a\ge 2$

3.对于 $p$ 为 $2^a$ 与多个奇素数的乘积且 $a\ge 1$

我们先处理第 $2,3$ 种情况。这两种情况其实可以概括成：

对于 $p=rs$ ，有 $r,s>2$ 且 $\gcd(r,s)=1$：

由 $\gcd(r,s)=1$ 可知 $\varphi(p)=\varphi(r)\varphi(s)$ ，然后由 $r,s>2$ 得到 $4|\varphi(p)$ （因为此时两者的欧拉函数都为偶数）。

记 $d=\frac{1}{2}\varphi(p)$ ，那么 $\varphi(r),\varphi(s)$ 都是 $d$ 的因子。于是对于任意 $\gcd(a,p)=1$，有 $a^d\equiv 1\pmod r$ 和 $a^d\equiv 1\pmod s$ ，则有 $a^d\equiv 1\pmod p$ 。发现 $d<\varphi(p)$，于是没有原根，

而第 $1$ 种情况：

由于互质的限制，我们只需考虑奇数。

因此对于任意奇数 $b=2k+1$ ，有：

$$b^{2^{a-2} }=(2k+1)^{2^{a-1} }\\ \equiv 1+\binom{2^{a-2} }{1}2k+\binom{2^{a-2} }{2}(2k)^2\\ \equiv 1+2^{a-1}(k+(2^{a-2}-1)k^2)\\ \equiv 1\pmod{2^a}$$

解释一下：

其中第二步，由于从第 $4$ 项开始，都具有 $2^{a-1}(2^{a-2}-2)$ 的组合，具有 $2^a$ 的因子，于是可以删去。

其中第三步 $k+(2^{a-2}-1)k^2$ 是偶数。（因为两项奇偶性相同）

于是 $a^{2^{a-1}}\not \equiv 1\pmod{2^a}$，没有原根。

Q.E.D.

性质5:原根存在定理

形如 $1,2,4,p^c,2p^c$ （$p$ 为奇质数）的数才存在原根。

证明：

对于 $p=1,2,4$ 的情况，分别有原根 $1,1,3$。

对于 $p^c$ 我们在 性质2 中已经证明，对于 $2p^c$ 我们在 性质3 中得到证明。

对于 $p\not=2,4,p^c,2p^c$ 的，我们在 性质4 中证明。

综上，该定理成立。

Q.E.D.

性质6:原根的数量级

可以证明，如果 $m$ 有原根，那么其最小原根 $g\le m^{0.25}$ ，证明略（我不道啊）。

所以找最小原根是可以暴力找的。怎么找呢？

性质7:原根判定定理

对于 $m\ge 3,\gcd(a,m)=1$ ，则 $a$ 是模 $m$ 意义下的原根的充要条件是：对于 $\varphi(m)$ 的每个质因数 $p$，都有 $a^{\frac{\varphi(m)}{p}}\not \equiv 1\pmod m$。

证明：

考虑反证，假设存在一个非模 $m$ 的原根 $a$ 使得存在一个质因数 $p$ 满足 $a^{\frac{\varphi(m)}{p} }\not\equiv 1\pmod m$ 。

于是存在 $t<\varphi(m)$ 使得 $a^t\equiv 1\pmod m$。

由 裴蜀定理 可知，存在 $x,y$ 使得 $ty=\varphi(m)x+\gcd(\varphi(m),t)$，于是有 $a^{ty}\equiv 1\equiv a^{\varphi(m)x+\gcd(\varphi(m),t)} \equiv a^{\gcd(t,\varphi(m))} \pmod m$ ，中间的转换根据欧拉定理。

则 $\gcd(t,\varphi(m))|\varphi(m)$ ，且 $\gcd(t,\varphi(m)) \le t<\varphi(m)$。

然后可知存在 $\varphi(m)$ 的质因数 $p$ 使得 $\gcd(t,\varphi(m) )|\frac{\varphi(m)}{p}$ 于是有 $a^{\frac{\varphi(m)}{p}}\equiv 1\pmod m$ ，矛盾。

Q.E.D.

求出原根

有了上面的性质的支撑，我们就可以计算原根了。我们直接暴力枚举然后用原根判定定理检验，其他原根可由这个原根推出。具体来说，$g^i,i\in[0,\varphi(p)]$ 在模 $p$ 意义下互不相同，我们只需要检验其阶是否为 $\varphi(p)$ 即可（凭借 阶的性质4 可以转化为 $i$ 与 $\varphi(p)$ 是否互质，其揭露了原根循环群的本质）。

复杂度 $O(p^{0.25}\varphi(p) \log \varphi(p)+\sqrt p)$ 。

给出参考代码：

#include<bits/stdc++.h>
#define ll long long
#define db double
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define sky fflush(stdout)
#define Better_IO true
namespace IO{
	#if Better_IO==true
		char buf[(1<<20)+3],*p1(buf),*p2(buf);
		const int lim=1<<20;
		inline char gc(){
			if(p1==p2) p2=(p1=buf)+fread(buf,1,lim,stdin);
			return p1==p2?EOF:*p1++;
		}
		#define pc putchar
	#else
		#define gc getchar
		#define pc putchar
	#endif
	inline bool blank(const char &c){
		return c==' ' or c=='\n' or c=='\t' or c=='\r' or c==EOF;
	}
	inline void gs(char *s){
		char ch=gc();
		while(blank(ch) ) {ch=gc();}
		while(!blank(ch) ) {*s++=ch;ch=gc();}
		*s=0;
	}
	inline void gs(std::string &s){
		s="";char ch=gc();s+='#';
		while(blank(ch) ) {ch=gc();}
		while(!blank(ch) ) {s+=ch;ch=gc();}
	}
	inline void ps(char *s){
		while(*s!=0) pc(*s++);
	}
	inline void ps(const std::string &s){
		for(auto it:s) 
			if(it!='#') pc(it);
	}
	template<class T>
	inline void read(T &s){
		s=0;char ch=gc();bool f=0;
		while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
		while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
		if(ch=='.'){
			db p=0.1;ch=gc();
			while('0'<=ch&&ch<='9') {s=s+p*(ch^48);p*=0.1;ch=gc();}
		}
		s=f?-s:s;
	}
	template<class T,class ...A>
	inline void read(T &s,A &...a){
		read(s);read(a...);
	}
};
using IO::read;
using IO::gs;
using IO::ps;
const int N=2e6+3;
int n,d;
inline int inc(int x,int y){
	return (x+=y)>=n?x-n:x;
}
inline int dec(int x,int y){
	return (x-=y)<0?x+n:x;
}
inline int mul(int x,int y){
	return 1ll*x*y%n;
}
inline int qpow(int a,int b){
	int ans=1;
	while(b){
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}
std::vector<int>pri;
db phi;
inline void get_phi(int x){
	phi=x;
	for(auto y:pri){
		phi*=(y-1);phi/=y;
	}
}
inline bool check(int x){
	for(auto y:pri){
		if(qpow(x,phi/y)==1)
			return 0;
	}
	return 1;
}
int gcd(int x,int y){
	if(y==0) return x;
	return gcd(y,x%y);
}
inline void answer(int x){
	static int ans[N],top;
	top=0;
	for(int i=0;i<phi;++i){
		if(gcd(i,phi)==1){
			ans[++top]=qpow(x,i);
		}
	}
	std::sort(ans+1,ans+1+top);
	printf("%d\n",top);
	for(int i=d;i<=top;i+=d){
		printf("%d ",ans[i]);
	}
}
int main(){
	file(a);
	int T;read(T);
	while(T--){
		read(n,d);
		if(n==2){
			printf("1\n");
			if(d==1) printf("1\n");
			else printf("\n");
			continue;
		}
		int x=n;
		pri.clear();
		for(int i=2;i*i<=x;++i){
			if(x%i==0){
				while(x%i==0) x/=i;
				pri.push_back(i);
			}
		}
		if(x!=1){
			pri.push_back(x);
		}
		get_phi(n);
		x=phi;
		pri.clear();
		for(int i=2;i*i<=x;++i){
			if(x%i==0){
				while(x%i==0) x/=i;
				pri.push_back(i);
			}
		}
		if(x!=1){
			pri.push_back(x);
		}
		bool f=0;
		for(int i=1;i<=n;++i){
			if(check(i) and gcd(i,n)==1){
				answer(i);f=1;
				break;
			}
		}
		if(!f) printf("0\n");
		pc('\n');
	}
	return 0;
}

再论快速数论变换

回到 NTT 本身，考虑使用原根来代替单位根。

首先，我们看看在 FFT 中使用了那些单位根的性质：

1.$\omega_n^k=(\omega_n^1)^k$

因此只需求出一个单位根 $\omega_n^1$。

2.$\omega_n^i,i\in[0,n)$ 互不相同

保证可以得到 $n$ 个不同点值。

3.$\omega_n^k=\omega_n^{k\bmod n}$

有了这条就然后需要推出：

$$\omega_n^{k+\frac{n}{2} }=-\omega_n^{k}\\ \sum_{i=0}^{n-1}w_n^{ai}=n[a=0]$$

展现出一个循环，或者说圆。

4.$\omega_{na}^{ka}=\omega_n^k$

等分性质。

我们看看怎么使得原根满足这几个性质。

其实是考虑一个模意义下阶为 $n$ 的一个元素 $g$。

发现对于模素数 $p$，其原根 $g$ 的阶为 $p-1$。

于是有 $\delta_p(g^k)=\frac{\delta_m(g)}{\gcd(k,\delta_p(g) )}=\frac{p-1}{\gcd(k,p-1)}$

我们发现，当 $n|(p-1)$ 时， $\delta_p(g^{\frac{p-1}{n} })=n$ ，此时 $g^{\frac{p-1}{n}}$ 满足条件 1,2,3。

至于为什么满足 3 的推理：

由 $(g^{\frac{n}{2} })^2=g^n\equiv 1\pmod p$

由 $g^i,i\in[1,n]$ 互不相同，得到 $g^{\frac{n}{2}}\equiv -1\pmod p$

得 $-g^{k+\frac{n}{2} }\equiv g^k\pmod p$ 。、

Q.E.D.

至于另一个推理则同理单位根的证明了。

4 也是易证的: $(g^{\frac{p-1}{kn}})^k=g^{\frac{p-1}{n}}$ 。

于是 $n|(p-1)$ 时，$g^{\frac{n}{2} }$ 是可以代替单位根存在的。

我们发现 $p=998244353$ 时， $p-1=2^{23}\times 7\times 17$。而我们的 $n$ 是 $2$ 的幂次，于是总是整除的。我们可以把这种模数称为 NTT模数。

接下来直接替换即可。顺带一提, $998244353$ 的一个原根为 $3$。

至于 IDFT 是由其逆元来弄的，把原根换成其逆元即可。

代码实现

#include<bits/stdc++.h>
#define ll long long
#define db double
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define sky fflush(stdout)
#define Better_IO true
namespace IO{
	#if Better_IO==true
		char buf[(1<<20)+3],*p1(buf),*p2(buf);
		const int lim=1<<20;
		inline char gc(){
			if(p1==p2) p2=(p1=buf)+fread(buf,1,lim,stdin);
			return p1==p2?EOF:*p1++;
		}
		#define pc putchar
	#else
		#define gc getchar
		#define pc putchar
	#endif
	inline bool blank(const char &c){
		return c==' ' or c=='\n' or c=='\t' or c=='\r' or c==EOF;
	}
	inline void gs(char *s){
		char ch=gc();
		while(blank(ch) ) {ch=gc();}
		while(!blank(ch) ) {*s++=ch;ch=gc();}
		*s=0;
	}
	inline void gs(std::string &s){
		s="";char ch=gc();s+='#';
		while(blank(ch) ) {ch=gc();}
		while(!blank(ch) ) {s+=ch;ch=gc();}
	}
	inline void ps(char *s){
		while(*s!=0) pc(*s++);
	}
	inline void ps(const std::string &s){
		for(auto it:s) 
			if(it!='#') pc(it);
	}
	template<class T>
	inline void read(T &s){
		s=0;char ch=gc();bool f=0;
		while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
		while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
		if(ch=='.'){
			db p=0.1;ch=gc();
			while('0'<=ch&&ch<='9') {s=s+p*(ch^48);p*=0.1;ch=gc();}
		}
		s=f?-s:s;
	}
	template<class T,class ...A>
	inline void read(T &s,A &...a){
		read(s);read(a...);
	}
};
using IO::read;
using IO::gs;
using IO::ps;
const int N=(1<<20)+3;
const int mods=998244353;
inline int inc(int x,int y){
	return (x+=y)>=mods?x-mods:x;
}
inline int dec(int x,int y){
	return (x-=y)<0?x+mods:x;
}
inline int mul(int x,int y){
	return 1ll*x*y%mods;
}
inline int qpow(int a,int b){
	int res=1;
	while(b){
		if(b&1) res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
int n,m;
int g=3,invg=qpow(g,mods-2);
int F[N<<1],G[N<<1],w[N];
int rev[N<<1];
inline void get_rev(){
	for(int i=1;i<n;++i){
		rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
	}
}
inline void NTT(int *f,bool I){
	for(int i=0;i<n;++i){
		if(i<rev[i]) std::swap(f[i],f[rev[i] ]);
	}
	for(int p=2;p<=n;p<<=1){
		int len=p>>1;
		w[0]=1;w[1]=qpow(I?invg:g,(mods-1)/p);
		for(int i=2;i<len;++i){
			w[i]=mul(w[i-1],w[1]);
		}
		for(int k=0;k<n;k+=p){
			for(int l=k;l<k+len;++l){
				int tmp=mul(f[l+len],w[l-k]);
				f[l+len]=dec(f[l],tmp);
				f[l]=inc(f[l],tmp);
			}
		}
	}
}
int main(){
	file(a);
	read(n,m);
	for(int i=0;i<=n;++i){
		read(F[i]);
	}
	for(int i=0;i<=m;++i){
		read(G[i]);
	}
	for(m+=n,n=1;n<=m;n<<=1);
	get_rev();
	NTT(F,0);NTT(G,0);
	for(int i=0;i<n;++i) F[i]=mul(F[i],G[i]);
	NTT(F,1);
	int invn=qpow(n,mods-2);
	for(int i=0;i<=m;++i){
		printf("%d ",mul(F[i],invn) );
	}
	return 0;
}

多项式插值

拉格朗日插值

对于 $n$ 次多项式 $f(x)$ ，我们知道其 $n+1$ 个点值 $(x_i,y_i)$，于是有：

$$f(x)=\sum_{i=1}^{n+1}y_i\prod_{j\not=i} \frac{x-x_j}{x_i-x_j}$$

由性质 $f(x)\equiv f(a) \pmod{(x-a)}$ 结合中国剩余定理易得。

计算的复杂度为 $O(n^2)$。

如果点值连续，可通过小转换得到：

$$f(x)=\sum_{i=1}^{n+1}y_i\frac{\frac{\prod_{j=1}^{n+1}(x-j)}{(x-i)} }{(-1)^{n+1-i}(i-1)!(n+1-i)!}$$

通过预处理可以 $O(n)$ 计算。

多项式快速插值

快速沃尔什变换

快速沃尔什变换, 即 FWT(Fast Walsh Transform).

其意义是位运算卷积，也就是将我们熟知的加法卷积的加法换为位运算。即：

$$c_k=\sum_{i\oplus j = k}a_ib_j$$

其中我们用 $\oplus$ 来代指某种位运算。

我们期望的做法同样是将多项式转为点值进行运算后在转回来。

考虑构造一组 $FWT[f(x)]$ 的变化，形如：

$$FWT[f](i)=\sum_{j=0}^{n-1} c(i,j)f(j)$$

其中 $c(i,j)$ 表示式子的系数。

我们知道， $A\times B=C$ ，则 $FWT[A]\times FWT[B]=FWT[C]$ 。

则：

$$FWT[A](i)*FWT[B](i)=FWT[C](i)\\ \sum_{j=0}^{n-1} c(i,j)A(j)\sum_{j=0}^{n-1} c(i,j)B(j)=\sum_{j=0}^{n-1} c(i,j)C(j)\\ \sum_{j=0}^{n-1}\sum_{k=0}^{n-1}c(i,j)c(i,k)A[j]B[k]=\sum_{j=0}^{n-1} c(i,j)C(j)\\ \sum_{j=0}^{n-1}\sum_{k=0}^{n-1}c(i,j)c(i,k)A[j]B[k]=\sum_{j=0}^{n-1} c(i,j)\sum_{a\oplus b=j} A(a)B(b)\\ \sum_{j=0}^{n-1}\sum_{k=0}^{n-1}c(i,j)c(i,k)A[j]B[k]=\sum_{a=0}^{n-1}\sum_{b=0}^{n-1}A(a)B(b)c(i,a\oplus b)\\\\ $$

得到 $c(i,j)c(i,k)=c(i,j\oplus k)$ 。然后不难发现 $c(i,x)$ 可以把 $x$ 按照当前二进制操作拆分然后相乘。

我们考虑如下的做法：

$$FWT[A](i)=\sum_{j=0}^{n-1}c(i,j)A(j)\\ =\sum_{j=0}^{\frac{n}{2}-1}c(i,j)A(j)+ \sum_{j=\frac{n}{2}}^{n-1}c(i,j)A(j)$$

发现左边没有 $n$ 的最高位 $0$ ，而右边有。

记 $x'$ 表示 $x$ 在二进制上去除最高位，$x_i$ 表示 $x$ 二进制上从高到底的第 $i$ 位。

得到：

$$c(i_0,0) \sum_{j=0}^{\frac{n}{2}-1}c(i',j')A(j)+c(i_0,1) \sum_{j=\frac{n}{2}}^{n-1}c(i',j')A(j)$$

发现一直这么拆是 $O(\log n)$ ，合并两式复杂度为 $O(n)$ ，总复杂度 $O(n\log n)$。

发现 $c$ 是个矩阵，我们对其求逆然后再 FWT ，就可以做到 IFWT 。

然后我们要用的 $c$ 矩阵只有 $2\times 2$。

接下来为 $3$ 种位运算构造矩阵。构造方法是把式子全列出来，然后手动枚举找可行情况，过于繁琐，于是直接给出构造和细节。

Or

$$\left [ \begin{matrix} 1\quad 1\\ 1 \quad 0 \end{matrix} \right ]$$

或

$$\left [ \begin{matrix} 1\quad 0\\ 1 \quad 1 \end{matrix} \right ]$$

一般采用第二种，因为其同时具有子集卷积的性质。也是记忆技巧。

$$\left [ \begin{matrix} 1\quad 0\\ 1 \quad 1 \end{matrix} \right ] \Rightarrow c(i,j)=[i\And j=j]$$

其逆矩阵为：

$$\left [ \begin{matrix} 1\quad 0\\ -1 \quad 1 \end{matrix} \right ]$$

And

$$\left [ \begin{matrix} 0\quad 1\\ 1 \quad 1 \end{matrix} \right ]$$

或

$$\left [ \begin{matrix} 1\quad 1\\ 0 \quad 1 \end{matrix} \right ]$$

采用第二种，其逆矩阵为：

$$\left [ \begin{matrix} 1\quad -1\\ 0 \quad 1 \end{matrix} \right ]$$

记忆技巧，为 $Or$ 矩阵的转置。

具有超集卷积的性质。

Xor

$$\left [ \begin{matrix} 1\quad 1\\ -1 \quad 1 \end{matrix} \right ]$$

或

$$\left [ \begin{matrix} 1\quad 1\\ 1 \quad -1 \end{matrix} \right ]$$

仍然采用第二种。

记忆技巧，$c(i,j)=(-1)^{i\&j}$。

其逆矩阵为：

$$\left [ \begin{matrix} \frac{1}{2}\quad \frac{1}{2}\\ \frac{1}{2} \quad -\frac{1}{2} \end{matrix} \right ]$$

【模板】快速沃尔什变换 (FWT)

#include<bits/stdc++.h>
#define ll long long
#define db double
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define sky fflush(stdout)
#define gc getchar
#define pc putchar
namespace IO{
	template<class T>
	inline void read(T &s){
		s=0;char ch=gc();bool f=0;
		while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
		while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
		if(ch=='.'){
			T p=0.1;ch=gc();
			while('0'<=ch&&ch<='9') {s=s+p*(ch^48);p/=10;ch=gc();}
		}
		s=f?-s:s;
	}
	template<class T,class ...A>
	inline void read(T &s,A &...a){
		read(s);read(a...);
	}
};
using IO::read;
const int mod=998244353;
inline int inc(int x,int y){
	return (x+=y)>=mod?x-mod:x;
}
inline int dec(int x,int y){
	return (x-=y)<0?x+mod:x;
}
inline int mul(int x,int y){
	return 1ll*x*y%mod;
}
inline int qpow(int a,int b){
	int res=1;
	while(b){
		if(b&1) res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
const int N=(1<<17)+3;
int _n,n;
int a[N],b[N],ans1[N],ans2[N],ans3[N];
int oa[N],ob[N];
const int inv2=qpow(2,mod-2);
const int 
	Cor[2][2]={{1,0},{1,1} },ICor[2][2]={{1,0},{mod-1,1} },
	Cand[2][2]={{1,1},{0,1} },ICand[2][2]={{1,mod-1},{0,1} },
	Cxor[2][2]={{1,1},{1,mod-1} },ICxor[2][2]={{inv2,inv2},{inv2,dec(mod,inv2)} };
inline void FWT(int *f,const int c[2][2]){
	for(int p=2;p<=n;p<<=1){
		int len=(p>>1);
		for(int k=0;k<n;k+=p){
			for(int l=k;l<k+len;++l){
				int tl=f[l],tr=f[l+len];
				f[l]=inc(mul(c[0][0],tl),mul(c[0][1],tr) );
				f[l+len]=inc(mul(c[1][0],tl),mul(c[1][1],tr) );
			}
		}
	}
}
int main(){
	file(a);
	read(_n);n=(1<<_n);
	for(int i=0;i<n;++i){
		read(a[i]);
		oa[i]=a[i];
	}
	for(int i=0;i<n;++i){
		read(b[i]);
		ob[i]=b[i];
	}
	FWT(a,Cor);FWT(b,Cor);
	for(int i=0;i<n;++i) ans1[i]=mul(a[i],b[i]);
	FWT(ans1,ICor);
	for(int i=0;i<n;++i){
		a[i]=oa[i];
		b[i]=ob[i];
	}
	FWT(a,Cand);FWT(b,Cand);
	for(int i=0;i<n;++i) ans2[i]=mul(a[i],b[i]);
	FWT(ans2,ICand);
	for(int i=0;i<n;++i){
		a[i]=oa[i];
		b[i]=ob[i];
	}
	FWT(a,Cxor);FWT(b,Cxor);
	for(int i=0;i<n;++i) ans3[i]=mul(a[i],b[i]);
	FWT(ans3,ICxor);
	for(int i=0;i<n;++i){
		printf("%d ",ans1[i]);
	}pc('\n');
	for(int i=0;i<n;++i){
		printf("%d ",ans2[i]);
	}pc('\n');
	for(int i=0;i<n;++i){
		printf("%d ",ans3[i]);
	}pc('\n');
	return 0;
}

FWT本质的探索

把数看成很多 $01$ 维，那么 $Or$ 运算等价于每一维度取 $max$ ， $And$ 运算等价于每一维取 $min$，而 $Xor$ 就是每维加起来 $\bmod 2$。

然后考虑扩展到 $k$ 进制数，那么 $Or$ 运算等价于每一维度取 $max$ ， $And$ 运算等价于每一维取 $min$，而 $Xor$ 就是每维加起来 $\bmod k$。

Or

考虑 $Or$ 运算，发现矩阵里面只会有 $0,1$ 。然后考虑同一行，对于 $y_1<y_2$ 有 $c(x,y_1)c(x,y_2)=c(x,max\{y_1,y_2\})=c(x,y_2)$ 。也就是说，同一行内的构成必定是先 $1$ 后 $0$ 。然后为了保证矩阵可逆，其需要满秩，那么每行互不相同即可。我们采用如下矩阵：（假设 $k=4$）

$$\left [ \begin{matrix} 1\quad 0\quad 0\quad 0 \\ 1\quad 1\quad 0\quad 0 \\ 1\quad 1\quad 1\quad 0 \\ 1\quad 1\quad 1\quad 1 \\ \end{matrix} \right ]$$

及其逆：

$$\left [ \begin{matrix} 1\quad 0\quad 0\quad 0 \\ -1\quad 1\quad 0\quad 0 \\ 0\quad -1\quad 1\quad 0 \\ 0\quad 0\quad -1\quad 1 \\ \end{matrix} \right ]$$

就是每个 $1$ 下面带个 $-1$ 。

两个矩阵可以看做先（对各维子集）做前缀和，再做差分。

套用 FWT ，复杂度 $O(nk\log_k n)$。

直接计算高维前缀和复杂度 $O(n\log_k n)$（每次 $k$ 进制下把一位挖掉），下同。

And

和 $Or$ 正好相反，矩阵也恰好是为 $Or$ 矩阵的转置。

$$\left [ \begin{matrix} 1\quad 1\quad 1\quad 1 \\ 0\quad 1\quad 1\quad 1 \\ 0\quad 0\quad 1\quad 1 \\ 0\quad 0\quad 0\quad 1 \\ \end{matrix} \right ]$$

及其逆：

$$\left [ \begin{matrix} 1\quad -1\quad 0\quad 0 \\ 0\quad 1\quad -1\quad 0 \\ 0\quad 0\quad 1\quad -1 \\ 0\quad 0\quad 0\quad 1 \\ \end{matrix} \right ]$$

两个矩阵可以看做先（对各维超集）做后缀和，再做差分。

Xor

多项式求逆

考虑问题 $\frac{f(x)}{g(x)} = a(x)$，求解 $a(x)$。

朴素递推求法

发现，举例假设 $g(x) = 1 - x ^ 2$ :

$$\frac{f(x)}{g(x)} = a(x)$$

$$f(x) = a(x)(1 - x ^ 2)$$

$$f[i] = a[i] - a[i - 2]$$

$$a[i] = f[i] + a[i - 2]$$

由此，上项由下项递推而来，可以递推求解。

设 $g(x)$ 项数为 $d$, $f(x)$ 项数为 $n$，则复杂度 $O(nd)$。

参考文献

原根存在定理证明

NTT与多项式全家桶