首页 > 其他分享 >FFT 学习笔记

FFT 学习笔记

时间:2024-07-01 21:31:34浏览次数:18  
标签:frac 原根 FFT 笔记 学习 delta alpha 2n equiv

\(\text{FFT}\) 学习笔记

多项式

确定一个多项式,往往只需要知道每一次项前的系数是多少即可。众所周知,一个朴素的多项式往往可以被写成

\[f(x)=\sum_{n\ge 0}a_nx^n \]

的形式,在这种形式下的两个多项式 \(f,g\) 的乘积 \(h\) 往往可以按照

\[h(x)=(f*g)(x)=\sum_{n\ge 0}(\sum_{i=0}^{n}a_ib_{n-i})x^n \]

的方式计算。也就是说,得到的新多项式的每一次项前的系数可以写作

\[c_n=\sum_{i=0}^na_ib_{n-i} \]

这种形式被叫做卷积,其特征是两个相乘的数所对应的下标相加是一个定值。不难发现,如此计算,对于两个次数为 \(n\) 的多项式来说,复杂度是 \(O(n^2)\),然而这非常的慢,这个时候就需要利用 \(\text{FFT}\),即快速傅里叶变换进行求解。

考虑 \(\text{FFT}\) 是如何优化的。对于一个 \(n\) 次多项式来说,只需要给定 \(n+1\) 个点,便可以确定一个唯一的 \(n\) 次多项式。那么我们在两个 \(n\) 次函数上取横坐标相同的 \(n+1\) 个点 \((x_i,y_i)\),那么对于新函数来说,新得到的 \(n+1\) 个点的坐标一定为 \((x_i,{y_1}_i{y_2}_i)\)。不难发现,此时,我们只需要 \(O(n)\) 的时间复杂度即可确定新得到的多项式。那么,现在需要做得就是对于多项式的系数表示和点值表示进行一个快速的转换。

复数

即形如 \(a+bi\) 形式的数,其中 \(a,b\) 为实数,\(i\) 为虚数单位,即 \(\sqrt{-1}\)。

接着引入复平面,也就是以实数为 \(x\) 轴,虚数为 \(y\) 轴建立的平面直角坐标系,一个复数都在复平面上对应一个点。

复数的加减只需要将实数部分和虚数部分相加减即可。对于乘法则有

\[(a+bi)(c+di)=(ac-bd)+(ad+bc)i \]

单位根

对于一个 \(n\) 次方程,一般来说都有 \(n\) 个根,而对于方程

\[x^n=1 \]

的这 \(n\) 个根,对应到复平面上来看,我们以原点为圆心,做一个半径为 \(1\) 的圆,称为单位圆,并将其 \(n\) 等分,与 \(x\) 正半轴相交的点为一个根(即 \(1\))。接着从 \(x\) 轴出发逆时针旋转到第 \(1\) 个点,这个点即为单位根,记作 \(w_n\)。接下来将第 \(i\) 个点记作 \(w_n^i\),不难发现有以下性质:

  1. \(w_{n}^k=w_{n}^{n+k}\)
  2. \(w_n^k=w_{2n}^{2k}\)
  3. \(w_{2n}^{n+k}=-w_{2n}^{k}\)
  4. \((w_n^k)^a=w_n^{ak}\)
  5. \(w_n^k\) 对于 \(k\in[0,n-1]\) 互不相等

对于单位根的计算,可以利用三角函数简单求解,即

\[w_n=\text{cos}(\frac{2\pi}{n})+\text{sin}(\frac{2\pi}{n})i \]

离散傅里叶变换

即 \(\text{DFT}\),是一种将系数表示快速转换成点值表示的变换。

对于一个 \(2n-1\) 次多项式 \(f(x)\):

\[f(x)=a_0+a_1x+a_2x^2+\cdots+a_{2n-1}x^{2n-1} \]

如果我们分别设

\[g(x)=a_0+a_2x+a_4x^2+\cdots\\h(x)=a_1+a_3x+a_5x^2+\cdots \]

那么显然有 \(f(x)=g(x^2)+xh(x^2)\),不妨代入 \(x=w_{2n}^k\),那么有

\[\begin{aligned}f(w_{2n}^k)&=g((w_{2n}^k)^2)+w_{2n}^kh((w_{2n}^k)^2)\\&=g(w_{2n}^{2k})+w_{2n}^{k}h(w_{2n}^{2k})\\&=g(w_n^k)+w_{2n}^{k}h(w_n^k) \end{aligned} \]

同理有

\[\begin{aligned}f(w_{2n}^{n+k})&=g((w_{2n}^{n+k})^2)+w_{2n}^{n+k}h((w_{2n}^{n+k})^2)\\&=g(w_{2n}^{2n+2k})-w_{2n}^{k}h(w_{2n}^{2n+2k})\\&=g(w_{2n}^{2k})-w_{2n}^{k}h(w_{2n}^{2k})\\&=g(w_n^k)-w_{2n}^{k}h(w_n^k) \end{aligned} \]

因此当我们求出 \(g(w_n^k)\) 和 \(h(w_n^k)\) 后便可以得到 \(f(w_{2n}^k)\) 和 \(f(w_{2n}^{n+k})\)。

由此我们考虑对于多项式 \(f\) 的 \(2n\) 个点分别代入 \(x=w_{2n}^i\),那么每一个值均可以通过交换系数后求解 \(g,h\),接着通过求解出的 \(g,h\) 更新所求的 \(f\)。为了让合并容易进行,不妨将次数补齐到 \(2^k-1\) 次。复杂度是 \(O(n\log n)\) 的。

但是这样写常数很大,我们考虑模拟交换系数的过程:

\[0,1,2,3,4,5,6,7\\ 0,2,4,6,1,3,5,7\\ 0,4,2,6,1,5,3,7 \]

我们考虑这些数的二进制表示:

\[000,001,010,011,100,101,110,111\\ 000,010,100,110,001,011,101,111\\ 000,100,010,110,001,101,011,111 \]

不难发现最初系数的位置和最终系数的位置满足下标进行了一次二进制下的翻转。

于是我们不难联想到与处理出位置是如何变化的,并存储在 \(r\) 数组中,接着考虑合并时本质上是将 \(2n\) 个位置中的 \(k,n+k\) 合并至 \(k,n+k\),因此可以直接进行覆盖。这样我们就完成了常数上的优化。

代码

void DFT(int n){
    Set(n-1);
    for(int i=0;i<n;i++){
        if(i<r[i])swap(a[i],a[r[i]]);//按照预处理交换系数
    }
    for(int len=1;len<n;len<<=1){
        Complex step(cos(Pi/len),sin(Pi/len));//单位根
        for(int l=0;l<n;l+=(len<<1)){
            Complex w(1,0);
            for(int k=0;k<len;k++,w=w*step){
                Complex A=a[l+k],B=w*a[l+len+k];
                a[l+k]=A+B;
                a[l+len+k]=A-B;
            }
        }
    }
    return ;
}

逆离散傅里叶变换

即 \(\text{IDFT}\),也就是将函数从点值表示转换成系数表示的变换。

考虑 \(\text{DFT}\) 在线性代数上的表示,即为

\[\begin{bmatrix}y_0\\y_1\\y_2\\\vdots\\y_{2n-1}\end{bmatrix}=\begin{bmatrix}w_{2n}^0&w_{2n}^0&w_{2n}^0&\cdots&w_{2n}^{0}\\w_{2n}^0&w_{2n}^1&w_{2n}^2&\cdots&w_{2n}^{n-1}\\w_{2n}^0&w_{2n}^2&w_{2n}^{4}&\cdots&w_{2n}^{2(n-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\w_{2n}^0&w_{2n}^{n-1}&w_{2n}^{2(n-1)}&\cdots&w_{2n}^{(n-1)^2}\end{bmatrix}\begin{bmatrix}a_0\\a_1\\a_2\\\vdots\\a_{2n-1}\end{bmatrix} \]

记点值矩阵为 \(E\),中间的矩阵为 \(D\),系数矩阵为 \(V\),则 \(E=DV\),那么有 \(V=ED^{-1}\),也就是我们只需要求出 \(D\) 矩阵的逆矩阵即可。

考虑构造 \(A_{ij}=D_{ij}^{-1}\),计算 \(AD\):

\[\begin{aligned}(AD)_{ij}&=\sum_{k=0}^{2n-1}A_{ik}D_{kj}\\&=\sum_{k=0}^{2n-1}w_{2n}^{-ik}w_{2n}^{kj}\\&=\sum_{k=0}^{2n-1}w_{2n}^{k(j-i)}\\&= \left\{\begin{array}{l}n&(i=j)\\\sum_{k=0}^{2n-1}(w_{2n}^{j-i})^k=\frac{(w_{2n}^{j-i})^{2n}-1}{w_{2n}^{j-i}-1}=\frac{(w_{2n}^{2n})^{j-i}-1}{w_{2n}^{j-i}-1}=0&(i\ne j)\end{array}\right.\end{aligned} \]

不难发现 \(\frac{AD}{n}=I\),所以有 \(D^{-1}=\frac{A}{n}\)。更具体的,因为 \(A_{ij}=D_{ij}^{-1}\),所以每个位置代入的数就变成了 \(w_{2n}^{-i}\),与 \(w_{2n}^{i}\) 关于 \(x\) 轴对称,也就是虚数部分相反。只需要对于 \(\text{DFT}\) 函数传于一个参数 \(\text{type}\) 表示如何转换即可。

代码

void DFT(int n,int type){
    Set(n-1);
    for(int i=0;i<n;i++){
        if(i<r[i])swap(a[i],a[r[i]]);
    }
    for(int len=1;len<n;len<<=1){
        Complex step(cos(Pi/len),type*sin(Pi/len));//是否需要取反
        for(int l=0;l<n;l+=(len<<1)){
            Complex w(1,0);
            for(int k=0;k<len;k++,w=w*step){
                Complex A=a[l+k],B=w*a[l+len+k];
                a[l+k]=A+B;
                a[l+len+k]=A-B;
            }
        }
    }
    if(type==-1){
        for(int i=0;i<n;i++)a[i].a/=n;//除以 n 才是正确结果
    }
    return ;
}

下面可以给出一个带封装的总代码:

const int N=2097152;//注意需要取到最大次数向上的最大二次幂数
const double Pi=acos(-1.0);//处理 2pi

int n,m;
int r[N];//预处理数组
struct Complex{
    double a,b;
    Complex(double a=0,double b=0):a(a),b(b){}
    Complex operator +(Complex x){return {a+x.a,b+x.b};}
    Complex operator -(Complex x){return {a-x.a,b-x.b};}
    Complex operator *(Complex x){return {a*x.a-b*x.b,a*x.b+b*x.a};}
};//手写复数
struct Polynomial{
    int len;
    vector<Complex>a;
    Complex& operator [](int x){return a[x];}
    void Set(int len){this->len=len,a.resize(len+1);}//节省空间的做法
    void DFT(int n,int type){
        Set(n-1);
        for(int i=0;i<n;i++){
            if(i<r[i])swap(a[i],a[r[i]]);
        }
        for(int len=1;len<n;len<<=1){
            Complex step(cos(Pi/len),type*sin(Pi/len));
            for(int l=0;l<n;l+=(len<<1)){
                Complex w(1,0);
                for(int k=0;k<len;k++,w=w*step){
                    Complex A=a[l+k],B=w*a[l+len+k];
                    a[l+k]=A+B;
                    a[l+len+k]=A-B;
                }
            }
        }
        if(type==-1){
            for(int i=0;i<n;i++)a[i].a/=n;
        }
        return ;
    }
    Polynomial operator *(Polynomial x){
        Polynomial y=*this,z;
        int len=x.len+y.len;
        int n=1;
        while(n<=len)n<<=1;//向上取到二次幂
        z.Set(n-1);
        for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
        x.DFT(n,1);
        y.DFT(n,1);
        for(int i=0;i<n;i++)z[i]=x[i]*y[i];
        z.DFT(n,-1);
        z.Set(len);
        return z;
    }
};
Polynomial f,g;

int main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cin>>n>>m;
    f.Set(n);
    g.Set(m);
    for(int i=0;i<=n;i++)cin>>f[i].a;
    for(int i=0;i<=m;i++)cin>>g[i].a;
    f=f*g;
    for(int i=0;i<=f.len;i++)cout<<(int)(f[i].a+0.5)<<" ";//精度问题
    return 0;
}

原根

这是一个又臭又长的内容。

拉格朗日定理:设 \(p\) 为素数,对于模 \(p\) 意义下的整系数多项式

\[f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_0(p\nmid a_n) \]

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

证明:使用归纳法,对于 \(n=0\) 时,由于 \(p\nmid a_0\),所以 \(f(x)\equiv 0\pmod p\) 无解,定理对 \(n=0\) 的多项式 \(f(0)\) 都成立。

若命题对于 \(n< x\) 均成立,采用反证法,假设当 \(n=x\) 时有至少 \(x+1\) 个解 \(x_0,x_1,\cdots,x_n\)。不妨设 \(f(x)-f(x_0)=(x-x_0)g(x)\),则 \(g(x)\) 在模 \(p\) 意义下是一个至多 \(n-1\) 次的多项式。因为当 \(x=x_1,x_2,\cdots,x_n\) 时,\(f(x)\equiv 0\pmod p\),所以

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

从而 \(g(x)\equiv 0\pmod p\) 至少有 \(n\) 个根,与归纳假设矛盾。

所以定理对 \(n\) 次多项式也成立,得证。

:由欧拉定理可知,对 \(a\in \mathbb{Z},m\in\mathbb{N}^*\),若 \(\gcd(a,m)=1\),则

\[a^{\varphi(m)}\equiv 1\pmod m \]

因此满足同余式 \(a^n\equiv 1\pmod m\) 的最小正整数 \(n\) 存在,这个 \(n\) 称作 \(a\) 模 \(m\) 的阶,记作 \(\delta_m(a)\)。

原根:设 \(m\in\mathbb{N}^*,a\in\mathbb{Z}\)。若 \(\gcd(a,m)=1\),且 \(\delta_m(a)=\varphi(m)\),则称 \(a\) 为模 \(m\) 的原根。

关于阶有一个较为显然的性质:

若 \(a^n\equiv 1\pmod m\),则 \(\delta_m(a)\mid n\)。

证明:对 \(n\) 除以 \(\delta_m(a)\) 作带余除法,设

\[n=\delta_m(a)q+r,0\le r<\delta_m(a) \]

若 \(r>0\),则

\(a^r\equiv a^r(a^{\delta_m(a)})^q\equiv a^n\equiv 1\pmod m\)

这与 \(\delta_m(a)\) 的最小性矛盾。故 \(r=0\),即 \(\delta_m(a)\mid n\)。

关于阶还有两个重要的性质:

性质 \(1\):设 \(m\in\mathbb{N}^*,a,b\in\mathbb{Z},\gcd(a,m)=\gcd(b,m)=1\),则

\[\delta_m(ab)=\delta_m(a)\delta_m(b) \]

的充要条件是 \(\gcd(\delta_m(a),\delta_m(b))=1\)。

证明:必要性:有 \(a^{\delta_m(a)}\equiv 1\pmod m\) 及 \(b^{\delta_m(b)}\equiv 1\pmod m\),可知

\[(ab)^{\text{lcm}(\delta_m(a),\delta_m(b))}\equiv 1\pmod m \]

由前面所述的性质,有

\[\delta_m(ab)\mid \text{lcm}(\delta_m(a),\delta_m(b)) \]

有由于 \(\delta_m(ab)=\delta_m(a)\delta_m(b)\),所以

\[\delta_m(a)\delta_m(b)\mid \text{lcm}(\delta_m(a),\delta_m(b)) \]

所以 \(\gcd(\delta_m(a),\delta_m(b))=1\)。

充分性:由 \((ab)^{\delta_m(ab)}\equiv 1\pmod m\) 可知

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

故 \(\delta_m(a)\mid\delta_m(ab)\delta_m(b)\)。结合 \(\gcd(\delta_m(a),\delta_m(b))=1\) 即得

\[\delta_m(a)\mid\delta_m(ab) \]

对称地,可得

\[\delta_m(b)\mid\delta_m(ab) \]

所以

\[\delta_m(a)\delta_m(b)\mid\delta_m(ab) \]

另一方面,有

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

\[\delta_m(ab)\mid\delta_m(a)\delta_m(b) \]

综合即可得到 \(\delta_m(ab)=\delta_m(a)\delta_m(b)\)。

性质 \(2\):设 \(k\in\mathbb{N},m\in\mathbb{N}^*,a\in\mathbb{Z},\gcd(a,m)=1\),则

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

证明:注意到

\[a^{k\delta_m(a^k)}=(a^k)^\delta_m(a^k)\equiv 1\pmod m\\ \Rightarrow \delta_m(a)\mid k\delta_m(a^k)\\ \Rightarrow \frac{\delta_m(a)}{\gcd(\delta_m(a),k)}\mid\delta_m(a^k) \]

另一方面,由 \(a^{\delta_m(a)}\equiv 1\pmod m\) 可知

\[(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)\mid\frac{\delta_m(a)}{\gcd(\delta_m(a),k)} \]

综合可得 \(\delta_m(a)=\frac{\delta_m(a)}{\gcd(\delta_m(a),k)}\)。

接下来讨论什么样的 \(m\) 存在原根。

首先,对于 \(m=2,4\),显然存在原根。

定理 \(1\):对于奇素数 \(p\),\(p\) 存在原根。

证明:先证明一个引理:

引理:设 \(a\) 与 \(b\) 是与 \(p\) 互素的两个整数,则存在 \(c\in\mathbb{Z}\) 使得 \(\delta_p(c)=\text{lcm}(\delta_p(a),\delta_p(b))\)。

证明:记 \(r=\delta_m(a),t=\delta_m(b)\),设它们的标准分解为(只要求 \(\max(\alpha_i,\beta_i)>0\))

\[r=\prod_{i=1}^{s}p_i^{\alpha_i},t=\prod_{i=1}^{s}p_i^{\beta_i} \]

令 \(l\) 为所有 \(\alpha_i\le\beta_i\) 的 \(p_i^{\alpha_i}\) 的乘积,\(m\) 为所有 \(\beta_i<\alpha_i\) 的 \(p_i^{\beta_i}\) 的乘积。记\(r=lx,t=my\),则 \(\gcd(x,y)=1,\text{lcm}(r,t)=xy\)。

由性质 \(2\),\(\delta_p(a^l)=x,\delta_p(b^m)=y\),则由性质 \(1\),\(\delta_p(a^lb^m)=xy=\text{lcm}(\delta_p(a),\delta_p(b))\),取 \(c=a^lb^m\) 即可。

对 \(1\) ~ \(p-1\) 依次两两使用引理,可知存在 \(g\in\mathbb{Z}\) 使得

\[\delta_p(g)=\text{lcm}(\delta_p(1),\delta_p(2),\cdots,\delta_p(p-1)) \]

这表明 \(\delta_p(j)\mid \delta_p(g)\),所以 \(j=1,2,\cdots,p-1\) 都是同余方程

\[x^{\delta_p(g)}\equiv 1\pmod p \]

的根。有拉格朗日定理可知,\(\delta_p(g)\ge p-1\)。又由费马小定理,易知 \(\delta_p(g)\le p-1\),所以 \(\delta_p(g)=p-1=\varphi(p)\)。所以 \(g\) 为 \(p\) 的原根。

定理 \(2\):对于奇素数 \(p\),\(\alpha\in\mathbb{N}^*\),\(p^{\alpha}\) 有原根。

证明:先证明一个引理。

引理:存在模 \(p\) 的原根 \(g\),使得 \(g^{p-1}\not\equiv 1\pmod {p^2}\)。

证明:事实上,仍取模 \(p\) 的原根 \(g\),若 \(g\) 不满足条件,我们认定 \(g+p\) 满足条件。易知 \(g+p\) 也是模 \(p\) 的原根。

于是有

\[\begin{aligned}(g+p)^{p-1}&\equiv g^{p-1}+(p-1)pg^{p-2}\\&\equiv 1-pg^{p-2}\\&\not\equiv1\pmod {p^2}\end{aligned} \]

接着,我们证明若 \(g\) 是一个满足引理条件的原根,则对于任意 \(\alpha\in\mathbb{N}^*\),\(g\) 是模 \(p^{\alpha}\) 的原根。

首先证明下面的结论:对于任意 \(\beta\in\mathbb{N}^*\),都可设

\[g^{\varphi(p^\beta)}=1+p^{\beta}k_\beta \]

不难发现当 \(\beta=1\) 时结论成立,现设上式对 \(\beta\) 成立,则

\[\begin{aligned}g^{\varphi(p^{\beta+1})}&=(g^{\varphi(p^\beta)})^p\\&=(1+p^\beta k_\beta)^p\\&\equiv 1+p^{\beta+1}k_\beta\end{aligned} \]

结合 \(p\nmid k_\beta\) 可知结论成立。

所以命题对 \(\beta\in\mathbb{N}^*\) 都成立。

其次,记 \(\delta=\delta_{p^\alpha}(g)\),则有欧拉定理,可知 \(\delta\mid p^{\alpha-1}(p-1)\)。

而由 \(g\) 为模 \(p\) 的原根,及 \(g^\delta\equiv1\pmod {p^\alpha}\) 可知 \((p-1)\mid \delta\)。

所以可设 \(\delta=p^{\beta-1}(p-1)(1\le\beta\le\alpha)\)。利用之前的结论可知

\[g^{\varphi(p^\beta)}\not\equiv1\pmod {p^{\beta+1}}\Rightarrow g^\delta\not\equiv1\pmod {p^{\beta+1}} \]

结合 \(g^\delta\equiv 1\pmod {p^\alpha}\) 可知 \(\beta\ge\alpha\)。所以 \(\beta=\alpha\),即

\[\delta_{p^\alpha}(g)=p^{\alpha-1}(p-1)=\varphi(p^\alpha) \]

从而 \(g\) 是 \(p^\alpha\) 的原根。

定理 \(3\):对于奇素数 \(p\),\(\alpha\in\mathbb{N}^*\),\(2p^\alpha\) 有原根。

证明:设 \(g\) 是模 \(p^\alpha\) 的原根,则 \(g+p^\alpha\) 也是模 \(p^\alpha\) 的原根。

在 \(g\) 和 \(g+p^\alpha\) 中有一个是奇数,设其为 \(G\),则 \(\gcd(G,2p^\alpha)=1\)。

由欧拉定理,\(\delta_{2p^\alpha}(G)\mid\varphi(2p^\alpha)\)。

而 \(G^{\delta_{2p^\alpha}(G)}\equiv 1\pmod{2p^\alpha}\),故 \(G^{\delta_{2p^\alpha}(G)}\equiv 1\pmod{p^\alpha}\)。

利用 \(G\) 为模 \(p^\alpha\) 的原根可知 \(\varphi(p^\alpha)\mid\delta_{2p^\alpha}(G)\)。

结合 \(\varphi(p^\alpha)=\varphi(2p^\alpha)\) 可知 \(G\) 为模 \(2p^\alpha\) 的原根。

定理 \(4\):对于 \(m\not\in\{2,4,p^\alpha,2p^\alpha\}\),则对于任意 \(a\in\mathbb{Z}\),都有 \(\delta_m{a}<\varphi(m)\),即模 \(m\) 的原根不存在。

证明:对于 \(m=2^\alpha,\alpha\in\mathbb{N}^{*},\alpha\ge 3\),则对于任意奇数 \(a=2k+1\) 有

\[\begin{aligned}a^{2^{\alpha-2}}&=(2k+1)^{2^{\alpha-2}}\\&\equiv 1+2^{\alpha-2}(2k)+2^{\alpha-1}(2^{\alpha-2}-1)k^2\\&\equiv 1+2^{\alpha-1}(k+(2^{\alpha-2}-1)k)\\&\equiv1\pmod{2^\alpha} \end{aligned} \]

若 \(m\) 不是 \(2\) 的次幂,则设 \(m=rt\) 满足 \(2<r<t,\gcd(r,t)=1\)。则由欧拉定理可知

\[a^{\varphi(r)}\equiv1\pmod{r},a^{\varphi(t)}\equiv 1\pmod{t} \]

当 \(n>2\) 时,\(\varphi(n)\) 为偶数,所以

\[a^{\frac{\varphi(r)\varphi(t)}{2}}\equiv1\pmod{rt} \]

所以 \(\delta_m(a)\le\frac{\varphi(r)\varphi(t)}{2}=\frac{\varphi(rt)}{2}=\frac{\varphi(m)}{2}<\varphi(m)\)。


上述定理阐述了原根的充要条件,即除了 \(2,4,p^\alpha,2p^\alpha\) 次方以外,其他的数都没有原根。于是我们可以与处理素数即其幂次,\(O(1)\) 判断是否有原根。

如果一个数有原根,可以先求出其最小原根。事实上,\(m\) 的最小原根是不多于 \(m^{\frac{1}{4}}\) 级别的。根据这个结论,我们可以直接从小到大枚举。根据原根的定义,若 \(g\) 为模 \(m\) 的原根,则对于 \(\varphi(m)\) 的任意素因数 \(p\),必有

\[g^{\frac{\varphi(m)}{p}}\not\equiv 1\pmod m \]

同时,只要满足如上条件的 \(g\) 一定是原根。于是我们可以与处理出 \(\varphi(m)\) 的所有质因数,通过快速幂来判断。

假如找到了一个原根 \(g\),不难证明对于所有 \(\gcd(x,\varphi(m))=1\) 的 \(x\),\(g^x\) 均为原根,所以模 \(m\) 的原根有 \(\varphi(\varphi(m))\) 个,所以我们可以在找到 \(g\) 后再枚举找出所有 \(m\) 的原根。

快速数论变换

即 \(\text{NTT}\),考虑到 \(\text{FFT}\) 存在精度问题,对于模意义下的数就需要用到 \(\text{NTT}\) 了。

考虑到此时单位根不再使用,我们考虑模意义下能够代替单位根的东西,也就是上面讲到的原根。考虑原根如何能够满足代替单位根进行运算,考虑一些性质。

  1. \(w_n^k\) 对于 \(k\in[0,n-1]\) 互不相等。

    若模数 \(p\) 是一个质数,那么对于原根 \(g\),\(g^k\) 对于 \(k\in[0,p-1]\) 在模 \(p\) 意义下互不相等,那么 \((g^{\frac{p-1}{n}})^k\) 对于 \(k\in[0,n-1]\) 在模 \(p\) 意义下互不相等。

  2. \(w_n^k=w_n^{n+k}\)。

    \((g^{\frac{p-1}{n}})^{n+k}=g^{p-1}(g^{\frac{p-1}{n}})^k\equiv(g^{\frac{p-1}{n}})^k\pmod p\)

  3. \(w_n^k=w_{2n}^{2k}\)。

    \((g^{\frac{p-1}{n}})^k=g^{\frac{k(p-1)}{n}}=g^{\frac{2k(p-1)}{2n}}=(g^{\frac{p-1}{2n}})^{2k}\pmod p\)

  4. \(w_{2n}^{n+k}=-w_{2n}^{k}\)。

    \(\frac{(g^{\frac{p-1}{2n}})^{n+k}}{(g^{\frac{p-1}{2n}})^k}=(g^{\frac{p-1}{2n}})^n=g^{\frac{p-1}{2}}\equiv -1\pmod p\),所以 \((g^{\frac{p-1}{2n}})^{n+k}\equiv-(g^{\frac{p-1}{2n}})^k\)。

  5. \((w_n^k)^a=w_n^{ak}\)。

    \(((g^{\frac{p-1}{n}})^k)^a=(g^{\frac{p-1}{n}})^{ak}\)。

因此,不难看出,单位根即为 \(g^{\frac{p-1}{n}}\),同时,这也对 \(p\) 作出了相应的限制,假设运算过程中多项式的次数最多到达 \(2^t\),那么 \(p=a2^\alpha+1,\alpha\ge t\)。下面是一些常见的 \(\text{NTT}\) 模数。

\[998244353=119\times 2^{23}+1,g=3\\ 469762049=7\times 2^{26}+1,g=3\\ 1004535809=479\times 2^{21}+1,g=3 \]

那么当题目给出的模数不是符合要求的 \(\text{NTT}\) 模数时应该怎么办呢——大骂出题人。

代码

#define ll long long
const int N=1e7+5,M=1e5+5,R=262144,P=1004535809,G=3,GI=334845270;

int n,m,t,s;
int w[M],r[R];
ll ans;
ll f[N],g[N];
namespace Math{
    ll QuickPow(ll a,int b,const int p){
        ll res=1;
        while(b>0){
            if((b&1)==1)res=res*a%p;
            a=a*a%p;
            b>>=1;
        }
        return res;
    }
    ll Inv(int x,int p){
        return QuickPow(x,p-2,p);
    }
};
struct Polynomial{
    int len;
    vector<ll>a;
    ll& operator [](int x){return a[x];}
    void Set(int len){this->len=len,a.resize(len+1);}
    void NTT(int n,int type){
        Set(n-1);
        for(int i=0;i<n;i++){
            if(i<r[i])swap(a[i],a[r[i]]);
        }
        for(int len=1;len<n;len<<=1){
            ll step=Math::QuickPow(type==1?G:GI,(P-1)/(len<<1),P);
            for(int l=0;l<n;l+=(len<<1)){
                ll w=1;
                for(int k=0;k<len;k++,w=w*step%P){
                    ll A=a[l+k],B=w*a[l+len+k]%P;
                    a[l+k]=(A+B)%P;
                    a[l+len+k]=(A-B)%P;
                }
            }
        }
        if(type==-1){
            ll inv=Math::Inv(n,P);
            for(int i=0;i<n;i++)a[i]=a[i]*inv%P;
        }
        return ;
    }
    Polynomial operator *(Polynomial x){
        Polynomial y=*this,z;
        int len=x.len+y.len;
        int n=1;
        while(n<=len)n<<=1;
        z.Set(n-1);
        for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
        x.NTT(n,1);
        y.NTT(n,1);
        for(int i=0;i<n;i++)z[i]=x[i]*y[i]%P;
        z.NTT(n,-1);
        z.Set(len);
        return z;
    }
}A,B;

int main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cin>>n>>m;
    A.Set(n);
    B.Set(m);
    for(int i=0;i<=n;i++)cin>>A[i];
    for(int i=0;i<=m;i++)cin>>B[i];
    A=A*B;
    for(int i=0;i<=A.len;i++)cout<<A[i]<<" ";
    return 0;
}

常见模型

事实上,所有的模型都离不开以下两种情况:

\[c_k=\sum_{i=0}^{k}a_ib_{k-i}\\ c_k=\sum_{i=0}^{k}a_ib_{k+i} \]

上面的就是朴素的卷积,而下面的称作差卷积,即两个数下标的差为定值。对于这种情况,我们通常选择讲一个序列倒转,即构造 \(b'_i=b_{n-i}\),那么有

\[c_k=\sum_{i=0}^{k}a_ib'_{n-(k+i)}=\sum_{i=0}^{k}a_ib'_{n-k-i} \]

于是发现和始终为 \(n-k\),可以求解。

字符串匹配

虽然很不想承认,但是 \(\text{FFT}\) 确实适用于字符串匹配的情景,考虑如何维护两个串是否相同,相同则意味着每一位都相同,因此我们构造

\[c_k=\sum_{i=0}^{m-1}(a_{k+i}-b_{i})^2 \]

不难发现,当 \(c_k=0\) 时,意味着从模式串的第 \(k\) 位开始能够成功匹配,那么进行一些化简:

\[c_k=\sum_{i=0}^{m-1}a_{k+i}^2-2\sum_{i=0}^{m-1}a_{k+i}b_{i}+\sum_{i=0}^{m-1}b_i^2 \]

很显然,除了中间那一项,其他的项均可以通过预处理 \(O(1)\) 得出,而中间显然是一个差卷积形式,因此构造 \(b'_{i}=b_{m-1-i}\),所以有

\[c_k=\sum_{i=0}^{m-1}a_{k+i}^2+b_i^2-2\sum_{i=0}^{m-1}a_{k+i}b'_{m-1-i} \]

\(\text{FFT}\) 求解即可,复杂度是 \(O(n\log n)\) 的。而对于存在通配符的匹配,我们认为通配符对应值为 \(0\),这意味着只要有通配符那么这一项就认为相等,因此构造

\[c_k=\sum_{i=0}^{m-1}a_{k+i}b_i(a_{k+i}-b{i})^2 \]

化简可得\

\[\begin{aligned}c_k&=\sum_{i=0}^{m-1}a_{k+i}^3b_i+ \sum_{i=0}^{m-1}a_{k+i}b_i^3-2\sum_{i=0}^{m-1}a_{k+i}^2b_i^2\\&=\sum_{i=0}^{m-1}a_{k+i}^3b_{m-1-i}+\sum_{i=0}^{m-1}a_{k+i}b_{m-1-i}^3-2\sum_{i=0}^{m-1}a_{k+i}^2b_{m-1-i}^2 \end{aligned} \]

不难看出只需要进行 \(3\) 次 \(\text{FFT}\) 即可。

背包问题

类似于求解给出 \(n\) 个数,选出 \(k\) 个数(可重复),求解和在某一范围内的方案数。不难想出构造 \(a_i\) 的值为有多少个值为 \(i\) 的数,所求即为询问范围内 \((a^k)_i\) 的和。需要注意,为了保证时间和空间的正确性,每次处理完后,对于值超过范围的数应当忽视,因为他们一定不符合要求。但是复杂度为 \(O(kn\log n)\) 依旧不优,考虑二进制优化,那么复杂度优化为 \(O(n\log n\log k)\)。

标签:frac,原根,FFT,笔记,学习,delta,alpha,2n,equiv
From: https://www.cnblogs.com/DycBlog/p/18278891

相关文章

  • 【吴恩达机器学习-week2】可选实验:特征工程和多项式回归【Feature Engineering and Po
    支持我的工作......
  • educoder 机器学习 --- kNN算法
    第一关:#encoding=utf8importnumpyasnpfromcollectionsimportCounterclasskNNClassifier(object):def__init__(self,k):'''初始化函数:paramk:kNN算法中的k'''self.k=k#用来......
  • webdav协议及我的笔记方案(私有部署)
    背景用markdown用于文章写作,有几年时间了,不是很喜欢折腾,主要就是在电脑上写,用的笔记软件就是typora。由于里面有很多工作相关的,以及个人资料相关的(包含了各种账号、密码啥的),所以不敢往各种云服务上放,还是想着数据由自己来管着。自己管数据的话,就是数据存储到哪里的问题,有很多朋......
  • 机器学习(四)——Lasso线性回归预测构建分类模型(matlab)
    Lasso线性回归(LeastAbsoluteShrinkageandSelectionOperator)是一种能够进行特征选择和正则化的线性回归方法。其重要的思想是L1正则化:其基本原理为在损失函数中加上模型权重系数的绝对值,要想让模型的拟合效果比较好,就要使损失函数尽可能的小,因此这样会使很多权重变为0或者权重......
  • 自动驾驶新篇章:基于大模型的协作驾驶与终身学习框架
    自动驾驶技术受到了学术界和工业界的广泛关注,但当前的自动驾驶系统大多基于数据驱动的方法,存在可解释性、泛化能力和持续学习能力方面的显著不足。而且单车自动驾驶系统缺乏与其他车辆协作和协商的能力,这对于提高驾驶安全性和效率至关重要。交通路口的场景:其中车1(veh1)和......
  • ESP32-点亮TFT2.4电阻触摸屏 学习笔记
    1、下载好arduinoIDE开发软件IDE(IntegratedDevelopmentEnvironment),译为集成开发环境,相当于编辑器编译器加连接器+其他。ArduinoIDE就是Arduino团队提供的一款专门为Arduino设计的编程软件,使用它,我们便能将程序从代码上传至Arduino主板。去官网下载:Software|Arduino,也......
  • 深度学习之激活函数
    激活函数的公式根据不同的函数类型而有所不同。以下是一些常见的激活函数及其数学公式:Sigmoid函数:公式:f(x)=特性:输出范围在0到1之间,常用于二分类问题,将输出转换为概率值。但存在梯度消失问题,尤其在输入值较大或较小时。Tanh函数(双曲正切函数):公式:f(x)=特性:输出范围在-1......
  • 从项目中学习Bus-Off的快慢恢复
    0前言        说到Bus-Off,大家应该都不陌生,使用VH6501干扰仪进行测试的文章在网上数不胜数,但是一般大家都是教怎么去干扰,但是说如何去看快慢恢复以及对快慢恢复做出解释比较少,因此本文以实践的视角来讲解Bus-Off的快慢恢复。1VH6501实现Bus-Off   首先需要......
  • Python学习笔记(二)
    目录while循环语句while循环的嵌套应用补充知识for循环函数猜数字游戏作业while循环语句练习:while循环的嵌套应用注:结束时i=i-1补充知识九九乘法表实例for循环练习:例如解决方案:再外部定义一个i=0for循环的嵌套break和continue的应用练......
  • Python学习笔记(一)
    目录 什么是变量​编辑数据类型转换语句标识符​编辑运算符​编辑字符串格式化数据输入python判断语句if语句​编辑if  elif  else常用的值类型:凡是被双引号包裹起来的都是字符串print输出,可以同时输出多个内容,用逗号隔开单行注释#  多行注释"""......