常数变易法的证明/P6613 题解
可能更好的阅读体验
题解区的很多大佬都是强解微分方程的，但常数变易法蒟蒻对此有许多不解。
因此蒟蒻在此给出本题常数变易法的证明。
一阶线性微分方程的常数变易法
这里有基础的大佬可以直接跳过。
让我们先考虑 一阶齐次线性微分方程：
显然，我们已经知道了
其中
这个过程被称为分离变量法，相当地简单。
但对于一般的 一阶线性微分方程，多了一个与
为求出该方程的通解，拉格朗日花了整整 11 年的时间（在此佩服大佬的毅力），
最后发现，如果把
只要使
可以看到，如果使
至于
由此得到了 一阶线性微分方程 的通解：
拉格朗日发现，这个通解满足
其实不难理解，
因此之后求解的时候就略去拆开
然后将常数变为函数代入回去，这就是一阶线性微分方程的常数变易法。
本题微分方程的常数变易法
接下来，让我们考虑本题的微分方程，它满足以下形式：
之前的经验告诉我们考虑对应的齐次方程：
其中
推广到非齐次，我们猜想，本题微分方程通解的形式为：
接下来证明为什么这么做是正确的。
同样将
当
于是得
由此得解，时间复杂度
Code
/*
this code is made by warzone
2021.5.22 20:00
*/
#include<stdio.h>
#include<string.h>
typedef unsigned char byte;
typedef unsigned long long ull;
typedef long long ll;
typedef unsigned int word;
struct READ{//快读
char c;
inline READ(){c=getchar();}
template<typename type>
inline READ& operator>>(register type& num){
while('0'>c||c>'9') c=getchar();
for(num=0;'0'<=c&&c<='9';c=getchar())
num=num*10+(c-'0');
return *this;
}
}cin;
class WRITE{//快写
private:
char out[1<<20],*top;
public:
inline WRITE(){top=out;}
inline ~WRITE(){if(top!=out) fwrite(out,1,top-out,stdout);}
inline WRITE& operator <<(char c){
*(top++)=c;
if(top==out+(1<<20))
fwrite(top=out,1,1<<20,stdout);
return *this;
}
inline void outnum(word num){
if(num) outnum(num/10),*this<<(char)(num%10+'0');
}
template<typename type>
inline WRITE& operator <<(type num){
if(num) return outnum(num),*this;
return *this<<'0';
}
}cout;
#define mx 18
#define mx_ 17
word root[1<<mx],inv[1<<mx],realid[1<<mx];
ull size,ni2;
const word mod=998244353;
inline ull pow(register ull a,register word b){
register ull ans=1;
for(;b;b>>=1){
if(b&1) (ans*=a)%=mod;
(a*=a)%=mod;
}
return ans;
}
//预处理单位根，翻转的 id
#define loading() \
register ull num1=pow(3,(mod-1)>>mx); \
register ull num2=pow(num1,mod-2); \
register word head,i,floor; \
root[1<<mx_]=inv[1<<mx_]=1; \
ni2=pow(2,mod-2); \
for(i=1;head=i,i<(1<<mx);++i) \
for(floor=0;floor<mx;++floor,head>>=1) \
realid[i]=realid[i]<<1|(head&1); \
for(i=1;i<(1<<mx_);++i){ \
root[1<<mx_|i]=num1*root[1<<mx_|(i-1)]%mod; \
inv[1<<mx_|i]=num2*inv[1<<mx_|(i-1)]%mod; \
} \
for(i=(1<<mx_)-1;i;--i) \
root[i]=root[i<<1],inv[i]=inv[i<<1];
#define load() \
register ull num1,num2;\
register word head,i,floor;
#define nttfor(size) \
for(floor=1;floor<1u<<(size);floor<<=1) \
for(head=0;head<1u<<(size);head+=floor<<1) \
for(i=0;i<floor;++i)
//floor 变换区间大小
//head 变换区间头指针
//i 变换位置
#define ntt(num,root)( \
num1=(num)[head|i], \
num2=(num)[head|i|floor], \
(num)[head|i]=(num1+((num2*=root[i|floor])%=mod))%mod, \
(num)[head|i|floor]=(num1+mod-num2)%mod)
#define id(size,i) (realid[i]>>(mx-(size)))
#define modx(num,size) memset(num+(1<<(size)),0,4<<(size))
#define set0(num,size) memset(num,0,4<<(size))
#define FOR(size) for(i=0;i<1u<<(size);i++)
#define newton(size) \
ull ni=ni2,size_=1; \
while(size_++,ni=ni*ni2%mod,size_-2<(size))
word in[1<<mx],eax[1<<mx],ebx[1<<mx],ecx[1<<mx],edx[1<<mx];
word out[1<<mx];
inline void _1(word size){//求 eax 的逆元，放入 ebx
ebx[0]=pow(eax[0],mod-2);
load()
newton(size){
FOR(size_-1){
head=id(size_,i);
ecx[head]=eax[i]? mod-eax[i]:0;
edx[head]=ebx[i];
head=id(size_,i+(1<<(size_-1)));
ecx[head]=edx[head]=0;
}
nttfor(size_) ntt(ecx,root),ntt(edx,root);
FOR(size_) ebx[id(size_,i)]=(ull)(ecx[i])*edx[i]%mod;
nttfor(size_) ntt(ebx,inv);
modx(ebx,size_-1);
FOR(size_) ecx[id(size_,i)]=ni*ebx[i]%mod;
ecx[0]=(ecx[0]+2)%mod;
nttfor(size_) ntt(ecx,root);
FOR(size_) ebx[id(size_,i)]=(ull)(ecx[i])*edx[i]%mod;
nttfor(size_) ntt(ebx,inv);
modx(ebx,size_-1);
FOR(size_-1) ebx[i]=ni*ebx[i]%mod;
}
}
inline void ln(word size){//求 eax 的 ln，放入 ebx
_1(size);
load()
FOR(size){
head=id(size+1,i);
ecx[head]=ebx[i];
edx[head]=(ll)(eax[i+1])*(i+1)%mod;
head=id(size+1,i+(1<<size));
ecx[head]=edx[head]=0;
}
nttfor(size+1) ntt(ecx,root),ntt(edx,root);
FOR(size+1) ebx[id(size+1,i)]=(ull)(ecx[i])*edx[i]%mod;
nttfor(size+1) ntt(ebx,inv);
modx(ebx,size);
for(register word i=(1u<<size)-1,ni=pow(1u<<(size+1),mod-2);i>0;i--)
ebx[i]=pow(i,mod-2)*ebx[i-1]%mod*ni%mod;
ebx[0]=0;
}
inline void exp(word size){// 求 in 的逆元，放入 eax
eax[0]=1;
load();
newton(size){
ln(size_-1);
FOR(size_-1){
head=id(size_,i);
ecx[head]=eax[i];
edx[head]=(mod-ebx[i]+in[i])%mod;
head=id(size_,i+(1<<(size_-1)));
ecx[head]=edx[head]=0;
}
edx[0]=(1ull+edx[0])%mod;
nttfor(size_) ntt(ecx,root),ntt(edx,root);
FOR(size_) eax[id(size_,i)]=(ull)(ecx[i])*edx[i]%mod;
nttfor(size_) ntt(eax,inv);
modx(eax,size_-1);
FOR(size_-1) eax[i]=ni*eax[i]%mod;
}
}
int main(){
loading();
word n,size=0;
cin>>n;
for(;1u<<size<=n+1;size++);
for(i=0;i<=n;++i)
cin>>out[id(size+1,i)];
for(i=1;i<=n+1;++i)
cin>>num1,in[i]=pow(i,mod-2)*num1%mod;
exp(size);
FOR(size){
ecx[id(size+1,i)]=eax[i]? mod-eax[i]:0;
ecx[id(size+1,i+(1u<<size))]=0;
}
nttfor(size+1) ntt(out,root),ntt(ecx,root);
FOR(size+1) ebx[id(size+1,i)]=(ull)(out[i])*ecx[i]%mod;
nttfor(size+1) ntt(ebx,inv);
eax[0]=1;
num1=pow(1u<<(size+1),mod-2);
for(i=1;i<1u<<(size);++i)
eax[i]=pow(i,mod-2)*num1%mod*ebx[i-1]%mod;
ln(size);
++in[0];
for(register word i=0;i<=n;i++)
cout<<((ull)(in[i])+mod-ebx[i])%mod<<' ';
return 0;
}