题解 LOJ6485 【LJJ学二项式定理】

题目

由于看到正解的单位根反演过于复杂 (也就是看不懂)

所以自己构造了一个算法,理论上这个算法应该还有成长的空间(可以变得普适性更强)

不知道和单位根反演有没有一样,就发表出来了

反正转载前记得要联系本人,联系方式参考 index


【分析】

题目所求为

\(\displaystyle Ans=[\sum_{k=0}^nC_n^ks^ka_{(k\mod 4)}]\mod 998244353\)

为避免混淆,本文中(除代码)所有 \(i\) 均表明虚数的单位,即 \(i^2=-1\),而代码中的变量名为 omega

考虑二项式定理: \(\displaystyle (a+b)^n=\sum_{k=0}^nC_n^ka^kb^{n-k}\)

我们可以很天真的想象一下,如果那玩意儿不是 \(a_{(k\mod 4)}\) 而是某个数的 \((n-k)\) 次方,那这题就很简单了

于是,前一段时间刚学完 NTT 的我就想到了,可不可以把它构造为这种形式?


对于 \(a_{(k\mod n)}\) 这种周期性的函数,我们可以拆成 \(\omega_n^{0,1,2\dots(n-1)}\)\(n\) 个单位根,关于某个函数 \(f(x)\) 的函数值

其中,\(\omega_n^k\) 均指 FFT 思路中的单位复数根

方便起见,我们这样设:\(\begin{cases}a_0=f(\omega_n^0)\\a_1=f(\omega_n^1)\\\vdots\\a_{n-1}=f(\omega_n^{n-1})\end{cases}\)

等等,这不就像 IDFT 吗?

因此我们直接设 \(f(x)\) 是一个多项式函数好了


当然,这一题我们不需要考虑 IDFT 来解题,毕竟 \(n=4\) ,手推就行了

首先,设 \(f(x)=m_0+m_1x+m_2x+m_3x^3\)

\(\omega_4^0=1,\omega_4^1=i,\omega_4^2=-1,\omega_4^3=-i\)

如果不知道为什么是这样,我在文章末尾补充

代入方程组得\(\begin{cases} a_0=f(1)=m_0+m_1+m_2+m_3 \\\ \\ a_1=f(i)=m_0+m_1i-m_2-m_3i \\\ \\ a_2=f(-1)=m_0-m_1+m_2-m_3 \\\ \\ a_3=f(-i)=m_0-m_1i-m_2+m_3i \end{cases}\)

解得\(\begin{cases} m_0={a_0+a_1+a_2+a_3\over 4} \\\ \\ m_1={(a_0-a_2)+(a_3-a_1)i\over 4} \\\ \\ m_2={a_0-a_1+a_2-a_3\over 4} \\\ \\ m_3={(a_0-a_2)+(a_1-a_3)i\over 4} \end{cases}\)

所以所求等式可以进行变型:

\(\displaystyle\quad Ans\)

\(\displaystyle =\sum_{k=0}^n C_n^ks^ka_{(k\mod 4)}\)

\(\displaystyle =\sum_{k=0}^n C_n^ks^kf(i^k)\)

\(\displaystyle =\sum_{k=0}^n C_n^ks^k[ m_0+m_1i^k+m_2(-1)^k+m_3(-i)^k]\)

展开多项式,每项提取系数 \(m\)

\(\displaystyle =m_0\sum_{k=0}^n C_n^ks^k1^k+m_1\sum_{k=0}^n C_n^ks^ki^k+m_2\sum_{k=0}^n C_n^ks^k(-1)^k+m_3\sum_{k=0}^n C_n^ks^k(-i)^k\)

好像更复杂了

等等,根据二项式定理 \(\displaystyle (a+b)^n=\sum_{k=0}^nC_n^ka^kb^{n-k}\)

应该指数和为 \(n\) ,而不是均为 \(k\) 啊!


别急,我们现在来化:

根据 \(i^{-1}=i^{4-1}=i^3=-i\)

所以有 \(i^k=(-i)^{-k}=(-i)^{-k}\times (-i)^n\times (-i)^{-n}=i^n\times (-i)^{n-k}\)

当然,同理也有 \(1^k=1^{n-k},(-1)^k=(-1)^n\times (-1)^{n-k},(-i)^k=(-i)^n\times i^{n-k}\)

代入就有

\(\displaystyle =m_0\sum_{k=0}^n C_n^ks^k1^{n-k}+m_1\sum_{k=0}^n C_n^ks^k(-i)^{n-k}i^n+m_2\sum_{k=0}^n C_n^ks^k(-1)^{n-k}(-1)^n+m_3\sum_{k=0}^n C_n^ks^ki^{n-k}(-i)^n\)

每一项都提出个常数项:

\(\displaystyle =m_0\sum_{k=0}^n C_n^ks^k1^{n-k}+i^nm_1\sum_{k=0}^n C_n^ks^k(-i)^{n-k}+(-1)^nm_2\sum_{k=0}^n C_n^ks^k(-1)^{n-k}+(-i)^nm_3\sum_{k=0}^n C_n^ks^ki^{n-k}\)

好的,现在就是二项式定理的形式了,直接化简:

\(\displaystyle =m_0(s+1)^n+i^nm_1(s-i)^n+(-1)^nm_2(s-1)^n+(-i)^nm_3(s+i)^n\)

代入上面的 \(m_{0,1,2,3}\) 定义式,这题应该就出来了,复杂度 \(O(T\log n)\)


当然,这里不需要手写一个复数类来实现,那样太复杂了。我们只需要考虑 \(i\) 在模 \(998244353\) 意义下的正整数就行了

由定义得到

\(i^2\equiv(-1)\equiv998244352(\mod 998244353)\)

我们可以算出来,\(i\equiv 911660635(\mod 998244353)\) (计算方法我放在后面)

也就是说,上面的所有 \(i\) 都用 \(911660635\) 代替掉就行了


【代码】

那本蒟蒻就放我 码风极丑的 代码了

#include<cstdio>
#include<algorithm>
using namespace std;
#define f(a,b,c,d) for(register int a=b,c=d;a<=c;a++)
#define g(a,b,c,d) for(register int a=b,c=d;a>=c;a--)
#define LOCAL
typedef int i32;
typedef unsigned int u32;
typedef long long int i64;
typedef unsigned long long int u64;
const i64 Mod=998244353;
const i64 inv4=748683265;//4的逆元
const i64 omega=911660635;

namespace HABIT{//读入输出优化
    #ifdef LOCAL
        inline char gc() { return getchar(); }
    #else
        inline char gc() {
            static char s[1<<20|1]={0},*p1=s,*p2=s;
            return (p1==p2)&&(p2=(p1=s)+fread(s,1,1<<20,stdin),p1==p2)?EOF:*(p1++);
        }
    #endif
    inline i64 read(){
        register i64 ans=0;register char c=gc();register bool neg=0;
        while(c<48||c>57) neg^=!(c^'-'),c=gc();
        while(c>=48&&c<=57) ans=(ans<<3)+(ans<<1)+(c^48),c=gc();
        return neg?-ans:ans;
    }//才不会告诉你们这里我忘了开long long,爆了三次

    char Output_Ans[1<<20|1],*Output_Cur=Output_Ans;
    inline void output() { Output_Cur-=fwrite(Output_Ans,1,Output_Cur-Output_Ans,stdout); }
    inline void print(char c){
        if(Output_Cur-Output_Ans+1>>20) output();
        *(Output_Cur++)=c;
    }
    inline void print(char *s) { while(*s) print( *(s++) ); }
    inline void print(u64 x){
        if(!x) { print('0'); return ; }
        char buf[30]={0},*p=buf+28;
        while(x) *(p--)=x%10+48,x/=10;
        print(p+1);
    }
    inline void print(u32 x) { print( (u64)x ); }
    inline void print(i64 x){
        if(x<0) print('-'),x=-x;
        print( (u64)x );
    }
    inline void print(i32 x){
        if(x<0) print('-'),x=-x;
        print( (u64)x );
    }
}
using namespace HABIT;

inline i64 fpow(i64 a,i64 x){//快速幂
    i64 ans=1; if(a<0) a+=Mod;
    for(;x;x>>=1,a=a*a%Mod) if(x&1) ans=ans*a%Mod;
    return ans;
}

inline i64 ans(){
    i64 d_N=read(),d_S=read()%Mod;
    i64 d_A0=read(),d_A1=read(),d_A2=read(),d_A3=read();

    i64 d_M0=(d_A0+d_A1+d_A2+d_A3)%Mod;
    i64 d_M1=( (d_A3-d_A1)*omega%Mod+d_A0-d_A2+Mod+Mod)%Mod;
    i64 d_M2=(d_A0-d_A1+d_A2-d_A3+Mod+Mod)%Mod;
    i64 d_M3=( (d_A1-d_A3)*omega%Mod+d_A0-d_A2+Mod+Mod)%Mod;
    //1/4最后再来算

    if(d_N&1){
        d_M2=Mod-d_M2;
        d_M1=d_M1*omega%Mod;
        d_M3=d_M3*omega%Mod;
        if( (d_N&3)>>1 ) d_M1=Mod-d_M1;
        else d_M3=Mod-d_M3;
    }
    else if( (d_N&3)>>1 ){
        d_M3=Mod-d_M3;
        d_M1=Mod-d_M1;
    }

    d_N%=(Mod-1);
    d_M0=d_M0*fpow(d_S+1,d_N)%Mod;
    d_M1=d_M1*fpow(d_S-omega,d_N)%Mod;
    d_M2=d_M2*fpow(d_S-1,d_N)%Mod;
    d_M3=d_M3*fpow(d_S+omega,d_N)%Mod;
    return (d_M0+d_M1+d_M2+d_M3)*inv4%Mod;
}

int main(){
    f(i,1,I,read()) print( ans() ),print('\n');
    output();
    return 0;
}

跑得飞快,目前是总榜第3

最后安利一下 本蒟蒻的博客


【补充说明-单位根】

考虑到复数的运算性质,如果两个复数相乘,模长相乘,幅角相加

我们要保证这个函数的周期性,只能设其模长为 \(1\)

而根据欧拉恒等式 \(e^{i\theta}=\cos\theta+i\sin\theta\) 模长即为 \(1\)

而周期为 \(n\) 所以令 \(\omega_n^k=e^{{2\pi\over n}k\dot i}\)

所以 \(\omega_n^k=\cos{2k\pi\over n}+i\sin{2k\pi\over n}\)

所以代入可以得到 \(\omega_4^0=1,\omega_4^1=i,\omega_4^2=-1,\omega_4^3=-i\)


【补充说明-同余意义下的开方】

第一种方法,这一题题目中的模数是确定的,直接暴力跑结果就行了,理论上 \(100s\) 内能出来

第二种方法,用原根

根据欧拉定理, \(gcd(a,m)=1\Rightarrow a^{\varphi(m)}\equiv 1(\mod m)\)

根据 \(998244353\) 是质数,且 \(998244353=2^{23}\times 7\times 17+1\)

避免眼花,后面这个数字写为 \(p\)

所以,首先有 \(\forall a<p\) 的正整数 \(a\) ,都有 \(a^{p-1}\equiv 1(\mod p)\)

当然,这并不说明只有指数那么大的时候,才会是 \(1\) 。但肯定说明,令 \(a^t\equiv 1(\mod p)\) 成立的最小 \(t\) ,一定是 \((p-1)\) 的因数

所以我们遍历过去,找到最小的 \(r\) 使得 \(r^{p-1\over 2}\equiv1,r^{p-1\over 7}\equiv 1,r^{p-1\over 17}\equiv 1(\mod p)\) 均不成立,那就肯定使得上述 \(t=p-1\) 本身了

对于满足上一段所述条件的,我们称之为原根 \(r\)。我算了一下,最小的 \(r\)\(3\)

所以显然, \(i^4\equiv 1\equiv r^{p-1\over 4}(\mod p)\)

快速幂即得 \(i\) 的值

原文地址:https://www.cnblogs.com/JustinRochester/p/12234247.html