「LOJ 150」挑战多项式

LOJ #150

题意

给定 \(n\) 次多项式 \(F(x)\)，求 \(G(x)\) 满足

\[ G(x)\equiv \left(\left(1+\ln\left(2+F(x)-F(0)-\exp\left(\int\frac{1}{\sqrt{F(x)}}\right)\right)\right)^k\right)' \pmod{x^n} \]

保证常数项是模 \(998244353\) 的二次剩余。

注意 \(\pm\sqrt{F(x)}\) 均为合法解，你只需要输出 \(\sqrt{F(x)}\)，舍去 \(-\sqrt{F(x)}\)，我们认为两个解中常数项较小的解为 \(\sqrt{F(x)}\)。

所有运算在模 \(998244353\) 意义下进行。

---

做法

参考多项式一些基础的操作和二次剩余Cipolla's algorithm

---

代码

存板子为主，里面还有多点求值和快速插值

定期不更新

常数应该还看得过去.

\(vector\)很方便啊

• update 2018-12-28 10:40:33

  用了定义的模意义的整数，改了 NTT 写法，优化了常数

• update 2019-1-15 12:15:02

  关于这里面 NTT 的写法，unsigned long long 大约是 \(18\times 10^{18}\) 的范围，由于模意义下数字几乎是随机的，可以认为每层会增加 \(\frac{P}{2}\)，其中 \(P\) 是 \(10^9\) 左右的模数。

  而一般要做的至多 \(21\) 层，或者常数大的操作本身只有 \(2^{18}\)，卡不爆 unsigned long long

  所以可以认为是对的

  如果有能卡掉的请务必告知

• update 2019-1-15 12:22:27

  关于这个乘法小范围的特判，没搞清楚 vector::size() 是什么类型，不过肯定是 unsigned 的，怕爆 int 还是转 unsigned long long 的靠谱一点吧。

• update 2019-1-15 12:27:21

  之前的代码里的插值和求值都是错的.

  由于只需要用到余数我调用了 PolyDiv(*, *, t, t)，之前单独改除法的时候忘了这个问题就把特判改了改，于是全错了.

  现在修复了

• update 2019-1-17 11:06:52

  小调整

#include<cstdio>
#include<algorithm>
#include<ctype.h>
#include<string.h>
#include<math.h>
#include<vector>

using namespace std;
#define ull unsigned long long

inline char read() {
	static const int IN_LEN = 1000000;
	static char buf[IN_LEN], *s, *t;
	return (s == t ? t = (s = buf) + fread(buf, 1, IN_LEN, stdin), (s == t ? -1 : *s++) : *s++);
}
template<class T>
inline void read(T &x) {
	static bool iosig;
	static char c;
	for (iosig = false, c = read(); !isdigit(c); c = read()) {
		if (c == '-') iosig = true;
		if (c == -1) return;
	}
	for (x = 0; isdigit(c); c = read()) x = ((x + (x << 2)) << 1) + (c ^ '0');
	if (iosig) x = -x;
}
const int OUT_LEN = 10000000;
char obuf[OUT_LEN], *ooh = obuf;
inline void print(char c) {
	if (ooh == obuf + OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh = obuf;
	*ooh++ = c;
}
template<class T>
inline void print(T x) {
	static int buf[30], cnt;
	if (x == 0) print('0');
	else {
		if (x < 0) print('-'), x = -x;
		for (cnt = 0; x; x /= 10) buf[++cnt] = x % 10 + 48;
		while (cnt) print((char)buf[cnt--]);
	}
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }

const unsigned N = 1<<18, P = 998244353;
struct Z{
	unsigned x;
	Z(const unsigned _x=0):x(_x){}
	inline Z operator +(const Z &rhs)const{ return x+rhs.x<P?x+rhs.x:x+rhs.x-P;}
	inline Z operator -(const Z &rhs)const{ return x<rhs.x?x-rhs.x+P:x-rhs.x;}
	inline Z operator -()const{ return x?P-x:0;}
	inline Z operator *(const Z &rhs)const{ return static_cast<ull>(x)*rhs.x%P;}
	inline Z operator +=(const Z &rhs){ return x=x+rhs.x<P?x+rhs.x:x+rhs.x-P, *this;}
	inline Z operator -=(const Z &rhs){ return x=x<rhs.x?x-rhs.x+P:x-rhs.x, *this;}
	inline Z operator *=(const Z &rhs){ return x=static_cast<ull>(x)*rhs.x%P, *this;}
};
int n, k;
vector<Z> f;

namespace Poly{

	Z w[N];// for DFT
	Z Inv[N];// for Integral
	vector<Z> ans;// for Evaluate()
	vector<vector<Z>> p;// for Evaluate() & Interpolate()

	inline Z Pow(Z x, int y=P-2){
		Z ans=1;
		for(; y; y>>=1, x=x*x) if(y&1) ans=ans*x;
		return ans;
	}
	inline void Init(){
		for(unsigned i=1; i<N; i<<=1){
			w[i]=1;
			Z t=Pow(3, (P-1)/i/2);
			for(unsigned j=1; j<i; ++j) w[i+j]=w[i+j-1]*t;
		}
		Inv[1]=1;
		for(unsigned i=2; i<N; ++i) Inv[i]=Inv[P%i]*(P-P/i);
	}
	inline pair<Z,Z> Mul(pair<Z,Z> x, pair<Z,Z> y, Z f){
		return make_pair((x.first*y.first+x.second*y.second*f),(x.second*y.first+x.first*y.second));
	}
	inline Z Quadratic_residue(Z a){
		if(a.x<=1) return a.x;
		if(Pow(a, (P-1)/2).x!=1) return -1;
		Z x, f;
		do x=(((ull)rand()<<15)^rand())%(a.x-1)+1; while(Pow(f=x*x-a, (P-1)/2).x==1);
		pair<Z,Z> ans=make_pair(1, 0), t=make_pair(x, 1);
		for(int i=(P+1)/2; i; i>>=1, t=Mul(t, t, f)) if(i&1) ans=Mul(ans, t, f);
		return min(ans.first.x, P-ans.first.x);
	}
	inline int Get(int x){ int n=1; while(n<=x) n<<=1; return n;}
	inline void DFT(vector<Z> &f, int n){
		static ull F[N];
		if((int)f.size()!=n) f.resize(n);
		for(int i=0, j=0; i<n; ++i){
			F[i]=f[j].x;
			for(int k=n>>1; (j^=k)<k; k>>=1);
		}
		for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){
			Z *W=w+i;
			ull *F0=F+j, *F1=F+j+i;
			for(int k=j; k<j+i; ++k, ++W, ++F0, ++F1){
				ull t=(*F1)*(W->x)%P;
				(*F1)=*F0+P-t, (*F0)+=t;
			}
		}
		for(int i=0; i<n; ++i) f[i]=F[i]%P;
	}
	inline void IDFT(vector<Z> &f, int n){
		f.resize(n), reverse(f.begin()+1, f.end());
		DFT(f, n);
		Z I=Pow(n);
		for(int i=0; i<n; ++i) f[i]=f[i]*I;
	}
	inline vector<Z> operator +(const vector<Z> &f, const vector<Z> &g){
		vector<Z> ans=f;
		for(unsigned i=0; i<f.size(); ++i) ans[i]+=g[i];
		return ans;
	}
	inline vector<Z> operator *(const vector<Z> &f, const vector<Z> &g){
		if((ull)f.size()*g.size()<=1000){
			vector<Z> ans;
			ans.resize(f.size()+g.size()-1);
			for(unsigned i=0; i<f.size(); ++i) for(unsigned j=0; j<g.size(); ++j)
				ans[i+j]+=f[i]*g[j];
			return ans;
		}
		static vector<Z> F, G;
		F=f, G=g;
		int p=Get(f.size()+g.size()-2);
		DFT(F, p), DFT(G, p);
		for(int i=0; i<p; ++i) F[i]*=G[i];
		IDFT(F, p);
		return F.resize(f.size()+g.size()-1), F;
	}
	vector<Z> &PolyInv(const vector<Z> &f, int n=-1){
		if(n==-1) n=f.size();
		if(n==1){
			static vector<Z> ans;
			return ans.clear(), ans.push_back(Pow(f[0])), ans;
		}
		vector<Z> &ans=PolyInv(f, (n+1)/2);
		vector<Z> tmp(&f[0], &f[0]+n);
		int p=Get(n*2-2);
		DFT(tmp, p), DFT(ans, p);
		for(int i=0; i<p; ++i) ans[i]=((Z)2-ans[i]*tmp[i])*ans[i];
		IDFT(ans, p);
		return ans.resize(n), ans;
	}
	// a=d*b+r
	inline void PolyDiv(const vector<Z> &a, const vector<Z> &b, vector<Z> &d, vector<Z> &r){
		if(b.size()>a.size()) return d.clear(), (void)(r=a);

		vector<Z> A=a, B=b, iB;
		int n=a.size(), m=b.size();
		reverse(A.begin(), A.end()), reverse(B.begin(), B.end());
		B.resize(n-m+1), iB=PolyInv(B, n-m+1);
		d=A*iB;
		d.resize(n-m+1), reverse(d.begin(), d.end());

		r=b*d, r.resize(m-1);
		for(int i=0; i<m-1; ++i) r[i]=a[i]-r[i];
	}
	inline vector<Z> Derivative(const vector<Z> &a){
		vector<Z> ans(a.size()-1);
		for(unsigned i=1; i<a.size(); ++i) ans[i-1]=a[i]*i;
		return ans;
	}
	inline vector<Z> Integral(const vector<Z> &a){
		vector<Z> ans(a.size()+1);
		for(unsigned i=0; i<a.size(); ++i) ans[i+1]=a[i]*Inv[i+1];
		return ans;
	}
	inline vector<Z> PolyLn(const vector<Z> &f){
		vector<Z> ans=Derivative(f)*PolyInv(f);
		ans.resize(f.size()-1);
		return Integral(ans);
	}
	vector<Z> PolyExp(const vector<Z> &f, int n=-1){
		if(n==-1) n=f.size();
		if(n==1) return {1};
		vector<Z> ans=PolyExp(f, (n+1)/2), tmp;
		ans.resize(n), tmp=PolyLn(ans);
		for(Z &i:tmp) i=-i;
		++tmp[0].x;
		ans=ans*(tmp+f);
		return ans.resize(n), ans;
	}
	vector<Z> PolySqrt(const vector<Z> &f, int n=-1){
		if(n==-1) n=f.size();
		if(n==1) return {Quadratic_residue(f[0])};// !
		vector<Z> ans=PolySqrt(f, (n+1)/2), tmp(&f[0], &f[0]+n);
		ans.resize(n), ans=ans+tmp*PolyInv(ans);
		for(Z &i:ans) i.x=(i.x&1?(i.x+P)/2:i.x/2);
		return ans;
	}
	inline vector<Z> PolyPower(const vector<Z> &f, int k){
		vector<Z> ans=PolyLn(f);
		for(Z &i:ans) i=i*k;
		return PolyExp(ans);
	}
	void Evaluate_Interpolate_Init(int l, int r, int t, const vector<Z> &a){
		if(l==r) return p[t].clear(), p[t].push_back(-a[l]), p[t].push_back(1);
		int mid=(l+r)/2, k=t<<1;
		Evaluate_Interpolate_Init(l, mid, k, a), Evaluate_Interpolate_Init(mid+1, r, k|1, a);
		p[t]=p[k]*p[k|1];
	}
	void Evaluate(int l, int r, int t, const vector<Z> &f, const vector<Z> &a){
		if(r-l+1<=512){
			for(int i=l; i<=r; ++i){
				Z x=0, a1=a[i], a2=a[i]*a[i], a3=a[i]*a2, a4=a[i]*a3, a5=a[i]*a4, a6=a[i]*a5, a7=a[i]*a6, a8=a[i]*a7;
				int j=f.size();
				while(j>=8) x=x*a8+f[j-1]*a7+f[j-2]*a6+f[j-3]*a5+f[j-4]*a4+f[j-5]*a3+f[j-6]*a2+f[j-7]*a1+f[j-8], j-=8;
				while(j--) x=x*a[i]+f[j];
				ans.push_back(x);
			}
			return;
		}
		vector<Z> tmp;
		PolyDiv(f, p[t], tmp, tmp);
		Evaluate(l, (l+r)/2, t<<1, tmp, a), Evaluate((l+r)/2+1, r, t<<1|1, tmp, a);
	}
	inline vector<Z> Evaluate(const vector<Z> &f, const vector<Z> &a, int flag=-1){
		if(flag==-1) p.resize(a.size()<<2), Evaluate_Interpolate_Init(0, a.size()-1, 1, a);
		return ans.clear(), Evaluate(0, a.size()-1, 1, f, a), ans;
	}
	vector<Z> Interpolate(int l, int r, int t, const vector<Z> &x, const vector<Z> &f){
		if(l==r) return {f[l]};
		int mid=(l+r)/2, k=t<<1;
		return Interpolate(l, mid, k, x, f)*p[k|1]+Interpolate(mid+1, r, k|1, x, f)*p[k];
	}
	inline vector<Z> Interpolate(const vector<Z> &x, const vector<Z> &y){
		int n=x.size();
		p.resize(n<<2), Evaluate_Interpolate_Init(0, n-1, 1, x);
		vector<Z> f=Evaluate(Derivative(p[1]), x, 0);
		for(int i=0; i<n; ++i) f[i]=y[i]*Pow(f[i]);
		return Interpolate(0, n-1, 1, x, f);
	}
	inline vector<Z> Solve(const vector<Z> &f, int k){
		vector<Z> ans=PolyExp(Integral(PolyInv(PolySqrt(f))));
		for(unsigned i=1; i<f.size(); ++i) ans[i]=f[i]-ans[i];
		ans=PolyLn(ans);
		++ans[0].x;
		return Derivative(PolyPower(ans, k));
	}
}

int main() {
	Poly::Init();
	read(n), read(k), f.resize(n+1);
	for(int i=0; i<=n; ++i) read(f[i].x);
	f=Poly::Solve(f, k), f.resize(n);
	while(!f[f.size()-1].x) f.erase(f.end()-1);
	for(Z i:f) print(i.x), print(' ');
	return flush(), 0;
}