题解 LOJ-6485 【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\) 的值
