多项式运算 (求逆/ln/exp等)
多项式运算 (求逆/ln/exp等)
(latest updated on 2021.02.23)
前置知识NTT
所有操作均在对\(P=\text{998244353}\)取模下进行
代码在最下面，由于板子实在有一点长，所以。。。
下文中\(\pmod {x^n}\)表示求出了多项式的前\(n\)项
\([x^i]F(x)\)表示\(F(x)\)第\(i\)项的系数
每个小问题的模板题都可以在洛谷上找到
1.多项式求乘法逆
(为什么叫做乘法逆?因为还有求复合逆和模逆元的)
求 \(G(x)\equiv \frac{1}{F(x)} \pmod {x^n}\)
形象化的理解就是\(F(x)\cdot G(x) \pmod {x^n}\)只有第一项是\(1\)，其他项都是\(0\)
这个由于是第一个操作，很多人还并不是很能理解多项式操作到底是什么东西，所以讲多一点
Part1 O(\(n^2\))
为了便于理解这个问题，先考虑一个最简单的模拟
\([x^i]F\cdot G(x)=\sum [x^j]F(x)[x^{i-j}]G(x)\)
第一项\([x^0]G(x)=\frac{1}{[x^0]F(0)} \pmod P\)，因此求逆的前提条件是\([x^0]F(x)\ne 0\)
考虑从\(1\)到\(n-1\)依次求出每一项，先从前面的项中得到所有\(j>0\)的和\(Sum\)，然后带入\(j=0\)时知道
Part2 O(\(n\log^2n\))
上面这个式子是一个类似\(dp\)转移的东西，可以直接分治NTT优化掉
Part3 \(O(n\log n)\)
考虑倍增求解，设已经求出了
其中递归边界是\(n=1\)时,\([x^0]G(x)=\frac{1}{[x^0]F(0)} \pmod P\)，因此求逆的前提条件是\([x^0]F(x)\ne 0\)
则
我们对于\(H(x)-G(x)\)平方，结果的前\(n\)项不可能由两个\(\ge \frac{n}{2}\)的项相乘得到，而前\(\frac{n}{2}\)项都是\(0\)，所以
所以通过平方可以扩大模数，这很常用
展开平方的式子
两边乘上\(F(x)\)
带入这个式子倍增求解即可
分析复杂度，每次有一个\(H(x)^2F(x)\)，可以通过\(NTT\)求出，倍增过程中访问的长度是\(O(n+\frac{n}{2}+\frac{n}{4}...)=O(n)\)
所以总复杂度就是\(O(n\log n)\)
2.多项式开根号
求\(G(x)^2\equiv F(x) \pmod {x^n}\)
同样的，递归求解，边界是\([x^0]=\sqrt{[x^0]F(x)} \pmod P\)
可以发现我们需要求二次剩余。。。但是一般题目保证了\([x^0]F(x)\in\{0,1\}\)
设已经求出\(H(x)^2\equiv F(x) \pmod{ x^{\lceil \frac{n}{2} \rceil}}\)
带入这个式子倍增求解即可
复杂度为\(O(n\log n)\)，由于需要求逆，实际可能会比较难写
3.多项式求\(\ln\)
对$\begin{aligned} G(x)\equiv \ln F(x) \pmod {x^n} \end{aligned} $ 两边求导，注意这里是复合函数求导!!!
\(\begin{aligned} G'(x)\equiv F'(x)\frac{1}{F(x)} \pmod {x^n}\end{aligned}\)
求出\(G'(x)\)，然后求原函数即可
通常保证\([x^0]F(x)=1\)，否则不好求\(\ln 1\)，所以求出原函数后首项为0
复杂度为\(O(n\log n)\)
4.多项式求exp
多项式求\(\text{exp}\)即求\(G(x)=e^{F(x)} \mod x^n\)
多项式求\(\text{exp}\)常见的解法有两种
CDQ分治+\(\text{NTT}\)
要求\(G(x)=e^{F(x)}\)
式子两边求导(右边要复合函数求导)，\(G'(x)=F'(x) e^{F(x)}\)
也就是说，\(G'(x)=F'(x)G(x)\)
两边同时积分得到\(\begin{aligned} G(x)=\int{F'(x)G(x)}\end{aligned}\)
我们知道，$ [x^i] \begin{aligned}\int H(x) =\frac{ [x^{i-1}]H(x)}{i}\end{aligned} $
带入上面的式子得到\(\displaystyle i\cdot [x^i]G(x)= \sum_{j=0}^{i-1}[x^j]F'(x)\cdot [x^{i-1-j}]G(x)\)
那么对于这个式子，直接使用分治NTT求解，其复杂为\(O(n\log n)\)
牛顿迭代
这是一种渐进意义上更优的做法，但实际在\(10^6\)以下几乎不可能更快，而且代码难写
但是不管平时用不用，牛顿迭代的知识学习一下肯定是最好的
把题目转化为，对于函数\(f(G)=\ln G-F\)
求出在\(\mod x^n\)意义下的零点
其中\(f(x)=\ln x-c\)
考虑迭代求解，设已经求出\(H(x)=e^{F(x)}\pmod {x^{\frac{n}{2}}}\)
边界条件是\([x^0]H(x)=e^{[x^0]F(x)}\)(由于没有办法求\(e^x\)在模意义下的值，所以通常必须要满足\([x^0]F(x)=0\))
带入牛顿迭代的结果
每次求\(\ln\) 复杂度和\(\text{NTT}\)相同，所以总复杂度为\(O(n\log n)\)
事实上这个还有优化的余地，就是在求\(\ln\)的时候，多项式逆的部分可以同步倍增求出，不需要每次都倍增一下(但是好像效果并不是特别明显)
5.多项式\(k\)次幂
\(G(x)\equiv F(x)^k\pmod {x^n}\)
\(\ln G(x)=k \ln F(x) \pmod {x^n}\)
求出\(\ln G(x)\)之后，\(\exp\)回来即可
由于要求\(\ln\)，所以这样求的条件是\([x^0]F(x)=1\) (可以通过平移和系数变换来调整为\([x^0]F(x)=1\))
很显然这个方法对于开根号也是适用的
复杂度\(O(n\log n)\)
6.多项式带余除法
问题：给定\(F(x),G(x)\)，其次数为\(n,m,n>m\)
求\(F(x)=G(x)P(x)+R(x)\)，其中\(P(x)\)次数为\(n-m\)，\(R(x)\)次数为\(m-1\)
考虑先求解\(P(x)\)，下面引入一种翻转运算
\(F^R(x)=x^nF(\frac{1}{x})\)，即将\(F(x)\)的系数翻转排列
用\(\frac{1}{x}\)带入问题的式子，得到
\(\displaystyle F(\frac{1}{x})=G(\frac{1}{x})P(\frac{1}{x})+R(x)\)
\(\displaystyle x^nF(\frac{1}{x})=x^mG(\frac{1}{x})\cdot x^{n-m}P(\frac{1}{x})+x^nR(x)\)
\(\displaystyle F^R(x)=G^R(x)\cdot P^R(x)+x^{n-m+1}R^R(x)\)
要求的\(P(x)\)是\(n-m\)次的，所以\(R^R(x)\cdot x^{n-m+1}\)并没有贡献
此时可以认为\(\displaystyle P^R(x)=\frac{F^R(x)}{G^R(x)}\)，求逆即可得到
得到\(P(x)\)之后，带入\(R(x)=F(x)-G(x)P(x)\)即可
以上是基本运算，如果不想继续吸多项式请直接跳到最下面的代码
所有的操作均用$\text{vector} $来实现，主要是为了理清思路，并且清零问题上会比较容易解决，同时对于每次计算完多项式的长度的要求会显得更加严格
实际在UOJ/Luogu上会非常慢，在LOJ上不错
稍微整理了一下，没怎么卡过常，所以应该还是比较可读的
代码总览(请使用C++11,O2编译运行)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
typedef unsigned long long ull;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }
char IO;
template <class T=int> T rd(){
T s=0; int f=0;
while(!isdigit(IO=getchar())) if(IO=='-') f=1;
do s=(s<<1)+(s<<3)+(IO^'0');
while(isdigit(IO=getchar()));
return f?-s:s;
}
const int N=1<<17,P=998244353;
int n,k;
ll qpow(ll x,ll k=P-2) {
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
/*
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int l=0;l<n;l+=i*2) {
int e=1;
for(int j=l;j<l+i;++j,e=1ll*e*w%P) {
int t=1ll*a[j+i]*e%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
int e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=1;j<i;++j) e[j]=1ll*e[j-1]*w%P;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
int e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
int w[N];
void Init(int N){
w[N>>1]=1;
int t=qpow(3,(P-1)/N);
rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
drep(i,(N>>1)-1,1) w[i]=w[i<<1];
}
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int *e=w+i;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
*/
/*
namespace MTT{
const double PI=acos((double)-1);
int rev[N];
struct Cp{
double x,y;
Cp(){ ; }
Cp(double _x,double _y): x(_x),y(_y){ }
inline Cp operator + (const Cp &t) const { return (Cp){x+t.x,y+t.y}; }
inline Cp operator - (const Cp &t) const { return (Cp){x-t.x,y-t.y}; }
inline Cp operator * (const Cp &t) const { return (Cp){x*t.x-y*t.y,x*t.y+y*t.x}; }
}A[N],B[N],C[N],w[N/2];
#define E(x) ll(x+0.5)%P
void FFT(int n,Cp *a,int f){
rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
w[0]=Cp(1,0);
for(reg int i=1;i<n;i<<=1) {
Cp t=Cp(cos(PI/i),f*sin(PI/i));
for(reg int j=i-2;j>=0;j-=2) w[j+1]=t*(w[j]=w[j>>1]);
// 上面提到的最优板子
for(reg int l=0;l<n;l+=2*i) {
for(reg int j=l;j<l+i;j++) {
Cp t=a[j+i]*w[j-l];
a[j+i]=a[j]-t;
a[j]=a[j]+t;
}
}
}
if(f==-1) rep(i,0,n-1) a[i].x/=n,a[i].y/=n;
}
void Multiply(int n,int m,int *a,int *b,int *res,int P){
// [0,n-1]*[0,m-1]->[0,n+m-2]
int S=(1<<15)-1;
int R=1,cc=-1;
while(R<=n+m-1) R<<=1,cc++;
rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
rep(i,0,n-1) A[i]=Cp((a[i]&S),(a[i]>>15));
rep(i,0,m-1) B[i]=Cp((b[i]&S),(b[i]>>15));
rep(i,n,R-1) A[i]=Cp(0,0);
rep(i,m,R-1) B[i]=Cp(0,0);
FFT(R,A,1),FFT(R,B,1);
rep(i,0,R-1) {
int j=(R-i)%R;
C[i]=Cp((A[i].x+A[j].x)/2,(A[i].y-A[j].y)/2)*B[i];
B[i]=Cp((A[i].y+A[j].y)/2,(A[j].x-A[i].x)/2)*B[i];
}
FFT(R,C,-1),FFT(R,B,-1);
rep(i,0,n+m-2) {
ll a=E(C[i].x),b=E(C[i].y),c=E(B[i].x),d=E(B[i].y);
res[i]=(a+((b+c)<<15)+(d<<30))%P;
}
}
#undef E
}
*/
namespace Polynomial{
typedef vector <int> Poly;
void Show(Poly a,int k=0){
if(!k){ for(int i:a) printf("%d ",i); puts(""); }
else for(int i:a) printf("%d\n",i);
}
int rev[N],w[N];
int Inv[N+1],Fac[N+1],FInv[N+1];
void Init_w() {
int t=qpow(3,(P-1)/N);
w[N>>1]=1;
rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
drep(i,(N>>1)-1,1) w[i]=w[i<<1];
Inv[0]=Inv[1]=Fac[0]=Fac[1]=FInv[0]=FInv[1]=1;
rep(i,2,N) {
Inv[i]=1ll*(P-P/i)*Inv[P%i]%P;
FInv[i]=1ll*FInv[i-1]*Inv[i]%P;
Fac[i]=1ll*Fac[i-1]*i%P;
}
}
int Init(int n){
int R=1,c=-1;
while(R<n) R<<=1,c++;
rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
return R;
}
#define NTTVersion1
#ifdef NTTVersion1
void NTT(int n,Poly &a,int f){
rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int *e=w+i;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,Mod2(a[j+i]);
a[j]+=t,Mod1(a[j]);
}
}
}
if(f==-1) {
reverse(a.begin()+1,a.begin()+n);
ll base=Inv[n];
rep(i,0,n-1) a[i]=a[i]*base%P;
}
}
void NTT(int n,int *a,int f){
rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int *e=w+i;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,Mod2(a[j+i]);
a[j]+=t,Mod1(a[j]);
}
}
}
if(f==-1) {
reverse(a+1,a+n);
ll base=Inv[n];
rep(i,0,n-1) a[i]=a[i]*base%P;
}
}
#else
void NTT(int n,Poly &a,int f){
static int e[N>>1];
rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int t=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*t*(e[j]=e[j>>1])%P;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,Mod2(a[j+i]);
a[j]+=t,Mod1(a[j]);
}
}
}
if(f==-1) {
ll base=Inv[n];
rep(i,0,n-1) a[i]=a[i]*base%P;
}
}
void NTT(int n,int *a,int f){
static int e[N>>1];
rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int t=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*t*(e[j]=e[j>>1])%P;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,Mod2(a[j+i]);
a[j]+=t,Mod1(a[j]);
}
}
}
if(f==-1) {
ll base=Inv[n];
rep(i,0,n-1) a[i]=a[i]*base%P;
}
}
#endif
Poly operator * (Poly a,Poly b){
int n=a.size()+b.size()-1;
int R=Init(n);
a.resize(R),b.resize(R);
NTT(R,a,1),NTT(R,b,1);
rep(i,0,R-1) a[i]=1ll*a[i]*b[i]%P;
NTT(R,a,-1);
a.resize(n);
return a;
}
Poly operator + (Poly a,Poly b) {
int n=max(a.size(),b.size());
a.resize(n),b.resize(n);
rep(i,0,n-1) a[i]+=b[i],Mod1(a[i]);
return a;
}
Poly operator - (Poly a,Poly b) {
int n=max(a.size(),b.size());
a.resize(n),b.resize(n);
rep(i,0,n-1) a[i]-=b[i],Mod2(a[i]);
return a;
}
Poly Poly_Inv(Poly a) { // 多项式乘法逆，注意这里求出的是前a.size()项
int n=a.size();
if(n==1) return Poly{(int)qpow(a[0],P-2)};
Poly b=a; b.resize((n+1)/2); b=Poly_Inv(b);
int R=Init(n<<1);
a.resize(R),b.resize(R);
NTT(R,a,1),NTT(R,b,1);
rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
NTT(R,a,-1);
a.resize(n);
return a;
}
Poly operator / (Poly a,Poly b){ // 多项式带余除法
reverse(a.begin(),a.end()),reverse(b.begin(),b.end());
int n=a.size(),m=b.size();
a.resize(n-m+1),b.resize(n-m+1),b=Poly_Inv(b);
a=a*b,a.resize(n-m+1);
reverse(a.begin(),a.end());
return a;
}
Poly operator % (Poly a,Poly b) { // 多项式取模
int n=b.size()-1;
if((int)a.size()<=n) return a;
Poly t=a/b;
if((int)t.size()>n) t.resize(n);
t=t*b; t.resize(n); a.resize(n);
return a-t;
}
int Quad(int a,int k=0) { // 二次剩余(不是原根法),用于求Sqrt
if(a<=1) return a;
ll x;
while(1) {
x=1ll*rand()*rand()%P;
if(qpow((x*x-a+P)%P,(P-1)/2)!=1) break;
}
ll w=(x*x-a+P)%P;
Pii res=mp(1,0),t=mp(x,1);
auto Mul=[&](Pii a,Pii b){
int x=(1ll*a.first*b.first+1ll*a.second*b.second%P*w)%P,y=(1ll*a.first*b.second+1ll*a.second*b.first)%P;
return mp(x,y);
};
int d=(P+1)/2;
while(d) {
if(d&1) res=Mul(res,t);
t=Mul(t,t);
d>>=1;
}
ll r=(res.first%P+P)%P;
if(k) r=min(r,(P-r)%P);
return r;
}
Poly Sqrt(Poly a){ // 多项式开根号
int n=a.size();
if(n==1) return Poly{Quad(a[0],1)};
Poly b=a; b.resize((n+1)/2),b=Sqrt(b),b.resize(n);
Poly c=Poly_Inv(b);
int R=Init(n*2);
a.resize(R),c.resize(R);
NTT(R,a,1),NTT(R,c,1);
rep(i,0,R-1) a[i]=1ll*a[i]*c[i]%P;
NTT(R,a,-1);
a.resize(n);
rep(i,0,n-1) a[i]=1ll*(P+1)/2*(a[i]+b[i])%P;
return a;
}
Poly Deri(Poly a){ //求导
rep(i,1,a.size()-1) a[i-1]=1ll*i*a[i]%P;
a.pop_back();
return a;
}
Poly IDeri(Poly a) { //原函数
a.pb(0);
drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*Inv[i]%P;
a[0]=0;
return a;
}
Poly Ln(Poly a){ // 多项式求Ln
int n=a.size();
a=Poly_Inv(a)*Deri(a),a.resize(n-1);
return IDeri(a);
}
Poly Exp(Poly a){ // 多项式Exp
int n=a.size();
if(n==1) return Poly{1};
Poly b=a; b.resize((n+1)/2),b=Exp(b); b.resize(n);
Poly c=Ln(b);
rep(i,0,n-1) c[i]=a[i]-c[i],Mod2(c[i]);
c[0]++,b=b*c;
b.resize(n);
return b;
}
void Exp_Solve(Poly &A,Poly &B,int l,int r){
static int X[N],Y[N];
if(l==r) {
B[l]=1ll*B[l]*Inv[l]%P;
return;
}
int mid=(l+r)>>1;
Exp_Solve(A,B,l,mid);
int R=Init(r-l+2);
rep(i,0,R) X[i]=Y[i]=0;
rep(i,l,mid) X[i-l]=B[i];
rep(i,0,r-l-1) Y[i+1]=A[i];
NTT(R,X,1),NTT(R,Y,1);
rep(i,0,R-1) X[i]=1ll*X[i]*Y[i]%P;
NTT(R,X,-1);
rep(i,mid+1,r) B[i]+=X[i-l],Mod1(B[i]);
Exp_Solve(A,B,mid+1,r);
}
Poly CDQ_Exp(Poly F){
int n=F.size(); F=Deri(F);
Poly A(n);
A[0]=1;
Exp_Solve(F,A,0,n-1);
return A;
}
Poly Pow(Poly x,int k) { // 多项式k次幂
x=Ln(x);
rep(i,0,x.size()-1) x[i]=1ll*x[i]*k%P;
return Exp(x);
}
Poly EvaluateTemp[N<<1];
void EvaluateSolve1(Poly &a,int l,int r,int p=1){
if(l==r) { EvaluateTemp[p]=Poly{P-a[l],1}; return; }
int mid=(l+r)>>1;
EvaluateSolve1(a,l,mid,p<<1),EvaluateSolve1(a,mid+1,r,p<<1|1);
EvaluateTemp[p]=EvaluateTemp[p<<1]*EvaluateTemp[p<<1|1];
}
void EvaluateSolve2(Poly &res,Poly F,int l,int r,int p=1){
if(l==r){ res[l]=F[0]; return; }
int mid=(l+r)>>1;
EvaluateSolve2(res,F%EvaluateTemp[p<<1],l,mid,p<<1);
EvaluateSolve2(res,F%EvaluateTemp[p<<1|1],mid+1,r,p<<1|1);
}
Poly Evaluate(Poly a,Poly b,int flag=1){ // 多项式多点求值
Poly res(b.size());
if(flag) EvaluateSolve1(b,0,b.size()-1);
EvaluateSolve2(res,a,0,b.size()-1);
return res;
}
Poly InterpolationSolve(Poly &T,int l,int r,int p=1){
if(l==r) return Poly{T[l]};
int mid=(l+r)>>1;
return InterpolationSolve(T,l,mid,p<<1)*EvaluateTemp[p<<1|1]+InterpolationSolve(T,mid+1,r,p<<1|1)*EvaluateTemp[p<<1];
}
Poly Interpolation(Poly X,Poly Y){ // 多项式快速插值
int n=X.size();
EvaluateSolve1(X,0,n-1);
Poly T=Evaluate(Deri(EvaluateTemp[1]),X,0);
rep(i,0,n-1) T[i]=Y[i]*qpow(T[i])%P;
return InterpolationSolve(T,0,n-1);
}
void FFPTrans(Poly &a,int f){ // FFP<->EGF
int n=a.size();
Poly b(n);
if(f==1) rep(i,0,n-1) b[i]=FInv[i];
else rep(i,0,n-1) b[i]=(i&1)?P-FInv[i]:FInv[i];
a=a*b; a.resize(n);
}
Poly FFPMul(Poly a,Poly b){ // FFP卷积
int n=a.size()+b.size()-1;
a.resize(n),b.resize(n);
FFPTrans(a,1),FFPTrans(b,1);
rep(i,0,n-1) a[i]=1ll*a[i]*b[i]%P*Fac[i]%P;
FFPTrans(a,-1);
return a;
}
Poly PolyToFFP(Poly F){ // 多项式转FFP
int n=F.size();
Poly G(n);
rep(i,0,n-1) G[i]=i;
G=Evaluate(F,G);
rep(i,0,n-1) F[i]=1ll*G[i]*FInv[i]%P;
FFPTrans(F,-1);
return F;
}
Poly FFPToPoly(Poly F){ // FFP转多项式
FFPTrans(F,1);
int n=F.size(); Poly X(n);
rep(i,0,n-1) X[i]=i,F[i]=1ll*F[i]*Fac[i]%P;
EvaluateSolve1(X,0,n-1);
rep(i,0,n-1) {
F[i]=1ll*F[i]*FInv[i]%P*FInv[n-i-1]%P;
if((n-i-1)&1) F[i]=(P-F[i])%P;
}
return InterpolationSolve(F,0,n-1);
}
}
using namespace Polynomial;
Poly Lag(int n,Poly X,Poly Y){
Poly T(n+1),R(n+1),A(n+1);
T[0]=1;
rep(i,0,n) drep(j,i+1,0) T[j]=(1ll*T[j]*(P-X[i])+(j?T[j-1]:0))%P;
rep(i,0,n) {
ll t=1;
rep(j,0,n) if(i!=j) t=t*(X[i]-X[j]+P)%P;
t=qpow(t)*Y[i]%P,R[n+1]=T[n+1];
drep(j,n,0) A[j]=(A[j]+t*R[j+1])%P,R[j]=(T[j]+1ll*R[j+1]*X[i]%P+P)%P;
}
return A;
}
int main(){
int n=rd();
Init_w();
Poly F(n);
rep(i,0,n-1) F[i]=rd();
Show(CDQ_Exp(F));
}
