【LOJ#150】挑战多项式
常用多项式模板整合
题目大意：
给定 $n$ 次多项式 $F(x)$ 以及整数 $k$，求多项式 $G(x)$ 满足：
在模 $998244353$ 意义下进行。
保证 $F(0)$ 在模意义下有二次剩余。
$\sqrt{F(x)}$ 取常数项较小的解。
题解：
这题整合了一些常用的多项式模板。
多项式求导、积分
直接对每项分别求导、积分即可。
多项式求逆
给定 $A(x)$，求 $B(x)$ 满足 $A(x)B(x)\equiv 1\pmod{x^n}$。
假设我们已经知道了 $A(x)B_0(x)\equiv 1\pmod{x^n}$，要求 $B(x)$ 满足 $A(x)B(x)\equiv 1\pmod{x^{2n}}$。
当 $n=1$ 时直接求逆计算，然后倍增即可。
时间复杂度 $O(n\log n)$。
多项式对数函数
给定 $A(x)$，求 $B(x)$ 满足 $B(x)\equiv\ln A(x)\pmod{x^n}$。保证 $A(0)=1$。
以下等号都指模 $x^n$ 意义下相等。
考虑求导，则 $B’(x)=(\ln A(x))’$。
根据链式法则，$B’(x)=\ln’(A(x))\cdot A’(x)=\frac{A’(x)}{A(x)}$。
多项式求逆即可。
时间复杂度 $O(n\log n)$。
多项式指数函数
给定 $A(x)$，求 $B(x)$ 满足 $B(x)\equiv e^{A(x)}\pmod{x^n}$。保证 $A(0)=0$。
以下等号都指模 $x^n$ 意义下相等。
考虑求导，则 $B’(x)=(e^{A(x)})’$。
根据链式法则，$B’(x)=e^{A(x)}A’(x)=A’(x)B(x)$。
使用 CDQ 分治 FFT 进行计算即可。
时间复杂度 $O(n\log^2 n)$，正常写法（指不刻意卡常）下比传统的 $O(n\log n)$ 牛顿迭代做法要快。
多项式幂函数
给定 $A(x),k$，求 $B(x)$ 满足 $B(x)\equiv A(x)^k\pmod{x^n}$。保证 $A(0)\neq 0$。
由于 $\exp(k\cdot \ln A(x))=A(x)^k$，使用直接使用多项式对指函数即可。
但是 $A(x)$ 的常数项可能不为 $1$，我们可以先除掉 $A(0)$，然后最后乘一个 $A(0)^k$ 即可。
时间复杂度同指数函数。
多项式开根
给定 $A(x)$，求 $B(x)$ 满足 $B(x)^2\equiv A(x)\pmod{x^n}$。
假设我们已经知道了 $B_0(x)^2\equiv A(x)\pmod{x^n}$，要求 $B(x)$ 满足 $B(x)^2\equiv A(x)\pmod{x^{2n}}$。
当 $n=1$ 时直接求二次剩余，然后倍增即可。中间嵌套多项式求逆。
时间复杂度 $O(n\log n)$。
二次剩余
给定 $n$，求 $b=\sqrt n$ 在模 $p$ 意义下的值。保证 $p$ 是质数。
以下等号都指模 $p$ 意义相等。
根据费马小定理，$a^{p-1}=1$。
所以 $b^{p-1}=1$，即 $n^{\frac{p-1}2}=1$。
所以如果 $n^{\frac{p-1}2}=-1$，则 $n$ 在模 $p$ 意义下没有二次剩余。
考虑找到一个 $a$ 满足 $a^2-n$ 在模 $p$ 意义下没有二次剩余（即 $(a^2-n)^{\frac{p-1}2}=-1$），令 $w^2=a^2-n$。
可以得到以下几个定理：
- $w^{p-1}=-1$，$w^p=-w$。
- $(a+w)^p=a^p+w^p=a-w$，前面一步可以用二项式展开来证。
那么 $(a+w)^{p+1}=(a+w)(a-w)=a^2-w^2=n$，所以 $b=(a+w)^{\frac{p+1}2}$。
具体实现时可以使用维护复数的方法，将数表示为 $a+bw$ 的形式，然后进行快速幂即可。
对于这个 $a$，暴力枚举即可，期望 $2$ 次能找到符合条件的数。
Code：
//Challenge polynomial
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
const int N=262144,md=998244353,_2=(md+1)/2;
namespace poly{
int W[N],rev[N],lim,_iv,fac[N],iv[N],Inv[N];
inline int pow(int a,int b){
int res=1;
for(;b;b>>=1,a=(LL)a*a%md)if(b&1)res=(LL)res*a%md;
return res;
}
inline void start(){
for(int i=1;i<N;i<<=1){
const int g=pow(3,(md-1)/(i<<1));
W[i]=1;
for(int j=1;j<i;++j)W[i+j]=(LL)W[i+j-1]*g%md;
}
for(int i=*fac=Inv[1]=1;i<N;++i)fac[i]=(LL)fac[i-1]*i%md;
iv[N-1]=pow(fac[N-1],md-2);
for(int i=N-2;~i;--i)iv[i]=(i+1LL)*iv[i+1]%md;
for(int i=2;i<N;++i)Inv[i]=(LL)(md-md/i)*Inv[md%i]%md;
}
inline void init(int n){
int l=-1;
for(lim=1;lim<n;lim<<=1)++l;
for(int i=1;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
_iv=pow(lim,md-2);
}
void DFT(int*b){
static unsigned long long a[N];
for(int i=0;i<lim;++i)a[rev[i]]=b[i];
for(int i=1;i<lim;i<<=1)
for(int j=0;j<lim;j+=i<<1)
for(int k=0;k<i;++k){
const unsigned long long x=a[j+k],y=a[j+k+i]*W[i+k]%md;
a[j+k]+=y,a[j+k+i]=x+md-y;
}
for(int i=0;i<lim;++i)b[i]=a[i]%md;
}
void IDFT(int*a){
DFT(a),reverse(a+1,a+lim);
for(int i=0;i<lim;++i)a[i]=(LL)a[i]*_iv%md;
}
void INV(int*a,int*b,int n){
if(n==1)*b=pow(*a,md-2);else{
const int len=(n+1)/2;
INV(a,b,len);
init(n+len);
static int A[N],B[N];
for(int i=0;i<n;++i)A[i]=a[i];
for(int i=n;i<lim;++i)A[i]=0;
for(int i=0;i<len;++i)B[i]=b[i];
for(int i=len;i<lim;++i)B[i]=0;
DFT(A),DFT(B);
for(int i=0;i<lim;++i)B[i]=B[i]*(2-(LL)A[i]*B[i]%md+md)%md;
IDFT(B);
for(int i=len;i<n;++i)b[i]=B[i];
}
}
void DER(int*a,int*b,int n){
for(int i=1;i<n;++i)
b[i-1]=(LL)a[i]*i%md;
b[n-1]=0;
}
void INT(int*a,int*b,int n){
for(int i=n-2;~i;--i)
b[i+1]=a[i]*(LL)Inv[i+1]%md;
b[0]=0;
}
void LN(int*a,int*b,int n){
static int c[N],d[N];
DER(a,c,n),INV(a,d,n);
init(n<<1);
for(int i=n;i<lim;++i)c[i]=d[i]=0;
DFT(c),DFT(d);
for(int i=0;i<lim;++i)c[i]=(LL)c[i]*d[i]%md;
IDFT(c);
INT(c,b,n);
for(int i=0;i<lim;++i)c[i]=d[i]=0;
}
void solve(int l,int r,int*a,int*b){
if(l==r)return;
const int mid=(l+r)>>1;
solve(l,mid,a,b);
static int A[N],B[N];
init(r-l+2);
for(int i=0;i<lim;++i)A[i]=B[i]=0;
for(int i=l;i<=mid;++i)A[i-l]=b[i];
for(int i=0;i<=r-l+1;++i)B[i]=a[i];
DFT(A),DFT(B);
for(int i=0;i<lim;++i)A[i]=(LL)A[i]*B[i]%md;
IDFT(A);
for(int i=mid+1;i<=r;++i)
b[i]=(b[i]+(LL)A[i-l-1]*Inv[i])%md;
solve(mid+1,r,a,b);
}
void EXP(int*a,int*b,int n){
static int c[N];
DER(a,c,n);
for(int i=0;i<n;++i)b[i]=0;
b[0]=1;
solve(0,n-1,c,b);
}
void POW(int*a,int*b,int k,int n){
static int d[N],c[N];
int s=a[0],t=pow(s,md-2);
for(int i=0;i<n;++i)d[i]=(LL)a[i]*t%md;
LN(d,c,n);
for(int i=0;i<n;++i)c[i]=(LL)c[i]*k%md;
EXP(c,b,n);
s=pow(s,k);
for(int i=0;i<n;++i)b[i]=(LL)b[i]*s%md;
for(int i=0;i<lim;++i)c[i]=d[i]=0;
}
namespace Cipolla{
int w2;
struct cpx{
int a,b;
inline cpx operator*(const cpx&rhs)const{return(cpx){((LL)a*rhs.a+(LL)b*rhs.b%md*w2)%md,((LL)a*rhs.b+(LL)b*rhs.a)%md};}
};
inline cpx pow(cpx a,int b){
cpx res=(cpx){1,0};
for(;b;b>>=1,a=a*a)if(b&1)res=res*a;
return res;
}
int _sqrt(int n){
for(int a=1;;++a)if(poly::pow((int)(((LL)a*a-n+md)%md),(md-1)/2)!=1){
w2=((LL)a*a-n+md)%md;
int res=pow({a,1},_2).a;
return res>md-res?md-res:res;
}
}
}
void SQRT(int*a,int*b,int n){
if(n==1)*b=Cipolla::_sqrt(*a);else{
SQRT(a,b,(n+1)/2);
static int A[N],B[N],C[N];
INV(b,C,n);
init(n<<1);
for(int i=0;i<n;++i)B[i]=b[i],A[i]=a[i];
for(int i=n;i<lim;++i)B[i]=C[i]=A[i]=0;
DFT(A),DFT(B),DFT(C);
for(int i=0;i<lim;++i)A[i]=(LL)(A[i]+(LL)B[i]*B[i])%md*_2%md*C[i]%md;
IDFT(A);
for(int i=(n+1)/2;i<n;++i)b[i]=A[i];
}
}
}
int F[N],n,k,A[N],B[N];
int main(){
poly::start();
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
cin>>n>>k;++n;
for(int i=0;i<n;++i)cin>>F[i];
poly::SQRT(F,A,n);
poly::INV(A,B,n);
poly::INT(B,A,n);
poly::EXP(A,B,n);
for(int i=1;i<n;++i)
A[i]=(F[i]-B[i]+md)%md;
A[0]=1;
poly::LN(A,B,n);
++B[0];B[0]%=md;
poly::POW(B,A,k,n);
poly::DER(A,B,n);
for(int i=0;i<n-1;++i)
cout<<B[i]<<' ';
return 0;
}