0%

「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

    小调整

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
#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;
}