「余姚中学 2019 联测 Day 2」Under Pressure(多项式类欧几里得,DFT)
程序员文章站
2022-05-22 19:22:42
...
对于第位,我们实际上是要求一个多项式:
我们知道类欧几里得实际上很多时候是需要转化为二维直线下的所有点的点权和。
这个题也可以这样转化:
原式=
前面是逆用了等比数列求和公式。
那么我们可以定义函数
那么
因为这个卷积是长度为的循环卷积,根据卷积定理我们只需要求出次单位根,求出多项式在这个单位根的点值后即可得到多项式。
那么求点值就简单了,上面所有形如,的点值都是可以计算的。(注意等比数列有公比为1这个天坑)
等等,哪来的次单位根?
因为,所以我们随便找一个的质数来当模数,就可以保证精确算出多项式的系数的同时还有次单位根。
这个哪里找呢?
直接从以上循环找即可,因为质数密度是的,所以期望找次即可找到,可以判质数,可以写一个压压惊。
然后再简单找一下原根即可。
时间复杂度
时间复杂度的瓶颈在,请务必循环展开,实测有倍的常数差距。
注:可以使用算法来代替,但是因为这个是模意义下的,所以需要写来避免精度流失(用两个不同模数难以简单切换),个人认为是跑不过循环展开的,欢迎大佬写一发来吊打我。
#include<bits/stdc++.h>
#define maxn 1005
#define LL long long
#define mod 998244353
using namespace std;
int n,K;
int t1 = clock() , ts = 0;
namespace Prime{
int Pow(int b,int k,int p){
int r = 1;
for(;k;k>>=1,b=1ll*b*b%p)
if(k&1)
r=1ll*r*b%p;
return r;
}
bool MR(int x){
static int base[6] = {2 , 3 , 7 , 31 , 61 , 101};
int r=x-1,t=0;
for(;!(r&1);r>>=1,t++);
for(int i=0;i<6;i++){
int nw = Pow(base[i],r,x) , pr = 0;
for(int k=0;k<t;k++){
swap(nw , pr);
if((nw = (1ll * pr * pr % x)) == 1 && pr != 1 && pr != x-1)
return 0;
}
if(nw != 1) return 0;
}
return 1;
}
}
using Prime::MR;
int P,w[maxn]={1},w_1[maxn]={0},r_1[maxn];
int Pow(int b,int k){ int r=1;for(;k;k>>=1,b=1ll*b*b%P) if(k&1) r=1ll*r*b%P; return r; }
void IDFT(int *A,int *B){
static int r[maxn];
for(int i=0;i<K;i++){
int t = i ? w[K-i] : w[0];
static int pw[9]={};
for(int k=pw[0]=1;k<=8;k++) pw[k] = 1ll * pw[k-1] * t % P;
r[i] = 0;
int j = K-1;
for(;j-7>=0;j-=8)
r[i] = (1ll * r[i] * pw[8]
+ 1ll * A[j] * pw[7]
+ 1ll * A[j-1] * pw[6]
+ 1ll * A[j-2] * pw[5]
+ 1ll * A[j-3] * pw[4]
+ 1ll * A[j-4] * pw[3]
+ 1ll * A[j-5] * pw[2]
+ 1ll * A[j-6] * pw[1]
+ 1ll * A[j-7] * pw[0]) % P;
for(;j>=0;j--)
r[i] = (1ll * r[i] * t + A[j]) % P;
}
static int iv = Pow(K , P-2);
for(int i=0;i<K;i++)
B[i] = r[i] * 1ll * iv % P;
}
int tot[5]={};
void F(int n,LL A,LL B,LL C,int p,int q,int *f){
if(n < 0) return;
if(B >= C){
LL b = B / C;
F(n,A,B%C,C,p,q,f);
for(int i=1;i<K;i++){
int iq = 1ll * i * q % K , ip = 1ll * i * p % K , biq = 1ll * b * iq % K;
int x = ip ? r_1[ip] * 1ll * w_1[1ll * ip * (n+1) % K] % P :
n+1,
y = iq % K ? r_1[iq] * 1ll * w_1[biq] % P :
b;
f[i] = (1ll * f[i] * w[biq] + 1ll * x * y) % P;
}
return;
}
if(A >= C){
LL a = A / C;
F(n,A%C,B,C,(p+1ll*a*q)%K,q,f);
for(int i=1;i<K;i++){
int iq = 1ll * i * q % K , ip = 1ll * i * p % K;
if(iq){
int x = (p + 1ll * a * q) % K , xi = 1ll * i * x % K;
x = xi ? w_1[1ll * xi * (n+1ll) % K] * 1ll * r_1[xi] % P : n+1;
int y = ip ? w_1[1ll * ip * (n+1ll) % K] * 1ll * r_1[ip] % P : n + 1;
f[i] = (f[i] + 1ll * (x - y) * r_1[iq]) % P;
}
else{
if(ip == 0)
f[i] = (f[i] + 1ll * a * (n * (n+1ll) / 2 % P)) % P;
else{
f[i] = (f[i] - (1ll * w_1[1ll * ip * n % K] * w[ip] % P * r_1[ip]
- 1ll * n * w[1ll * ip * (n+1ll) % K]) % P * r_1[ip] % P * a
) % P;
}
}
}
return;
}
if(! A){
LL b = B / C + 1;
for(int i=1;i<K;i++){
int iq = 1ll * i * q % K , ip = 1ll * i * p % K;
int x = ip ? r_1[ip] * 1ll * w_1[1ll * ip * (n+1) % K] % P :
n+1,
y = iq ? r_1[iq] * 1ll * w_1[1ll * iq * b % K] % P :
b;
f[i] = 1ll * x * y % P;
}
return;
}
LL T = (1ll*A*n+B)/C , R = (1ll * A * n + B) % C , nT = floor(((n + 1ll) * A - R - 1.0) / C) , np = (K-p)%K , nq = (K-q)%K;
F(nT,C,R,A,nq,np,f);
for(int i=1;i<K;i++){
int nqi = 1ll * nq * i % K , npi = 1ll * np * i % K;
int x = nqi ? w[1ll * (nT+1) * nqi % K] * 1ll * w_1[1ll * (T-nT) * nqi % K] % P * r_1[nqi] % P :
T - nT;
int y = npi % K ? w_1[1ll * (n + 1) * npi % K] * 1ll * r_1[npi] % P : n + 1;
f[i] = 1ll * (f[i] + 1ll * x * y % P) * w[(1ll * p * n + 1ll * q * T) % K * i % K] % P;
}
}
int main(){
scanf("%d%d",&n,&K);
static int fn[maxn],A[maxn],B[maxn],C[maxn];
for(int i=0;i<K;i++) scanf("%d",&fn[i]);
LL mx = 0;
for(int i=1;i<=n;i++) scanf("%d%d%d",&A[i],&B[i],&C[i]),mx=max(mx,A[i]*1ll*B[i]+C[i]);
for(P=ceil(1000000002.0/K)*K+1;!MR(P);P+=K);
vector<int>pr;
int t = P - 1;
for(int i=2;i*i<=t;i++)
if(t % i == 0){
pr.push_back(i);
for(;t % i == 0;t /= i);
}
if(t > 1) pr.push_back(t);
for(int i=2;;i++){
bool flg = 0;
for(int j=0;j<pr.size();j++)
if(Pow(i,(P-1) / pr[j]) == 1){
flg = 1;
break;
}
if(!flg){
w[1] = Pow(i,(P-1) / K);
break;
}
}
for(int i=1;i<K;i++) w[i] = 1ll * w[i-1] * w[1] % P , w_1[i] = w[i] - 1 , r_1[i] = Pow(w_1[i],P-2);
int ans = 0;
for(LL k=1;k<=mx && k>=0;k=k*K){
static int g[maxn]={};
for(int i=0;i<K;i++) g[i] = 1;
for(int i=1;i<=n;i++){
static int f[maxn]={};
memset(f,0,sizeof f);
f[0] = A[i] + 1;
if(A[i] * 1ll * B[i] + C[i] >= k){
F(A[i],B[i],C[i],k,0,1,f);
for(int p=1;p<K;p++)
f[p] = (f[p] * (w[p] - 1ll) % P + A[i] + 1) * w[(K-p) % K] % P;
IDFT(f,f);
}
for(int j=0;j<K;j++){
f[j] = (f[j] + P) % P;
f[j] = (f[j] + (j ? f[j-1] : 0)) % mod;
g[j] = 1ll * g[j] * f[j] % mod;
}
}
for(int i=1;i<K;i++)
ans = (ans + 1ll * (g[i] - g[i-1]) * fn[i]) % mod;
}
printf("%d\n",(ans+mod)%mod);
}
上一篇: OpenCV,计算两幅图像的单应矩阵 完成两个图像的融合
下一篇: smarty模板基础知识