您的位置: 首页  >  IT文章  >  [原创题] 雾

[原创题] 雾

分类: IT文章 • 2023-11-22 20:26:55

题目：http://162.105.80.126/contest/%E3%80%90%E5%BC%B1%E7%9C%81%E8%83%A1%E7%AD%96%E3%80%91Round%20%233/%E9%9B%BE

感谢绍兴一中三位神犇前来捧(nue)场，感谢山东神犇。。。。。。

让我们来看一下题目，

$$Ans = sum_{i=1}^{n}sum_{j=1}^{i-1}prod_{k=1}^{min(i-j,j)}[sum_{l=1}^k a(l)]^{varphi(j-k)+varphi(i-j-k)}prod_{i=min(i-j,j)+1}^{max(i-j,j)}f(i,j,k)$$
若$i-j>j$ , $f(i,j,k)=[sum_{l=1}^ka(l)]^{varphi(i-j-k)}$

否则 $f(i,j,k)=[sum_{l=1}^ka(l)]^{varphi(j-k)}$。

注意到$sum_{l=1}^k a(l)$和其他的式子并没有关系，所以可以预先前缀和，下面我们所说的$a(i)$都是前缀和后的$a(i)$。
$$Ans = sum_{i=1}^{n}sum_{j=1}^{i-1}prod_{k=1}^{min(i-j,j)}a(k)^{varphi(j-k)+varphi(i-j-k)}prod_{i=min(i-j,j)+1}^{max(i-j,j)}f(i,j,k)$$

$10 points$： 按照定义$O(n^3)$

显然这个式子是难以优化的，我们考虑加以变形，注意到$f(i,j,k)$函数的定义恰好和前面的连乘的内容相似，我们可以将$f(i,j,k)$代入化简，然后将左侧的连乘拆开为
$$prod_{k=1}^{min(i-j,j)}a(k)^{varphi(j-k)}cdot a(k)^{varphi(i-j-k)}$$
然后做变形
$$prod_{k=1}^{min(i-j,j)}a(k)^{varphi(j-k)} prod_{k=1}^{min(i-j,j)}a(k)^{varphi(i-j-k)}$$
然后发现 $f(i,j,k)$ 其实可以接到左侧式子或者右侧式子（取决于$i-j$和$j$的大小关系）,整理得到：
$$Ans = sum_{i=1}^{n}sum_{j=1}^{i-1}prod_{k=1}^{j}a(k)^{varphi(j-k)}prod_{i=1}^{i-j}a(k)^{varphi(i-j-k)}$$

然后，仔细观察后面那一大坨连乘和前面的求和无关。
设$C(n) = prod_{i=1}^{n}a(i)^{varphi(n-i)}$
然后发现原式变成了：
$$Ans = sum_{i=1}^{n}sum_{j=1}^{i-1}{C(i)cdot C(i-j)}$$
(注意这里的卷积形式，要在首项补$0$,保证$j$是从$i-1$结束)

$30 points$:
假如你不会$FFT$也不会$NTT$那么你可以$O(n^2)$预处理，计算时直接用。
总复杂度：$O(n^2)$

标准的卷积形式，可以用$NTT$ $O(nlogn)$ 解决。
然后回来看$C(n)$，左右用原根代换：
$$ gn^{c'(n)} = prod_{i=1}^{n}gn^{a'(i) cdot {varphi(n-i)}}$$
求出 $C'(n)$ 之后我们就可以 $O(nlogn)$ 快速幂求出 $C(n)$ 了。
而注意到其指数恰为卷积的形式，但是根据费马小定理，应该取$mod-1$的余，模数不太对，不能用 $NTT%$，而 $FFT$ 又炸精度。

考虑别的方法，由于数据中的指数是很小的，最大的情况下数字可以是$4 cdot 115000000000$，取一个超级大的质数，我是$1899956092796929$，原根为$7$，然后就可以了吗，朴素乘法会溢出，快速乘或者用$long double$的大数相乘取余函数。（然而快速乘多$logn$；$long double$乘有概率挂，所以不用是上策）

此外，确定一个数字是原根的多少次方时可以用$BSGS$ $O(nsqrt{n})$注意这里的n小于等于$5cdot 10^4$.

总效率$O(nlog^2n)$还是不能通过所有的数据呀，预计$40$~$50$分。

$50points$:
最后，我们发现我们可以选两个$10^9$级的满足条件的质数，最后$CRT$合并。
然后就可以$50$了（此处应$orz$ $yjp$）

$70points$:
继续改进我们的算法，可以发现更本不用$BSGS$，完全可以预处理出所有的值，用一个$map$或$hash$存下来，取时$O(logn)$查询，总体复杂度：$O(nlogn)$。

$100points$:
可以发现，之所以$70$分，并不是因为复杂度，而是$NTT$太慢了，注意到指数很小，所以可以$FFT$，或者卡卡常数即可得到全分。

注意精度。

#include <cstdio>
#include <cmath>
#include <ctime>
#include <algorithm>
#include <map>

#define N 400010
#define LL unsigned long long
#define mod 1004535809LL
#define gn 3
#define mod2 998244353LL

using namespace std;

const double pi=acos(-1.0);

struct Ex{
    double r, i;
    Ex(){}
    Ex(double r, double i):r(r), i(i){}
    friend Ex operator+(const Ex &a,const Ex &b){return Ex(a.r+b.r,a.i+b.i);}
    friend Ex operator-(const Ex &a,const Ex &b){return Ex(a.r-b.r,a.i-b.i);}
    friend Ex operator*(const Ex &a,const Ex &b){return Ex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
}a0[N],b0[N],c0[N];

map<LL,int> MT;

inline LL qmult(LL a,LL b,LL P){
    LL ans=0;
    for(;b;b>>=1,a=(a<<1)%P)
        if(b&1) ans=(ans+a)%P;
    return ans;
}

inline void swap(LL &a,LL &b){
    if(a^b) a^=b,b^=a,a^=b;
}

inline LL qpow(LL x,LL n,LL P){
    LL ans=1;
    for(;n;n>>=1,x=x*x%P)
        if(n&1) ans=ans*x%P;
    return ans;
}

int R[N];

inline LL Mod(LL x,LL P){
    if(x>=P) x-=P;
    if(x>=P) x-=P;
    if(x>=P) x%=P;
    return x;
}

LL wt[N];

inline void fnt(LL *x,int n,int t,LL P){
    for(int i=0;i<n;i++) if(i<R[i]) swap(x[i],x[R[i]]);
    for(int m=1;m<n;m<<=1){
        LL wn=qpow(gn,(P-1)/(m<<1),P);
        if(t==-1) wn=qpow(wn,P-2,P);
        wt[0]=1;
        for(int i=1;i<m;i++) wt[i]=Mod(wt[i-1]*wn,P);
        for(int k=0;k<n;k+=(m<<1))
            for(int i=0;i<m;i++){
                LL &A=x[i+m+k];
                LL &B=x[i+k],t=Mod(A*wt[i],P);
                A=Mod(B-t+P,P); B=Mod(B+t,P);
            }
    }
    if(t==-1){
        LL tmp=qpow(n,P-2,P);
        for(int i=0;i<n;i++)
            x[i]=x[i]*tmp%P;
    }
}

inline void fft(Ex x[],int n,int t){
    for(int i=0;i<n;i++) if(i<R[i]) swap(x[i],x[R[i]]);
    for(int m=1;m<n;m<<=1){
        Ex wn(cos(pi/m*t),sin(pi/m*t));
        for(int k=0;k<n;k+=(m<<1)){
            Ex wt(1,0);
            for(int i=0;i<m;i++,wt=wt*wn){
                Ex &A=x[k+m+i],&B=x[k+i],t=wt*A;
                A=B-t; B=B+t;
            }
        }
    }
    if(t==-1) for(int i=0;i<=n;i++) x[i].r/=n;
}

LL a1[N],b1[N],c1[N];

inline void multpoly2(LL *A,LL *B,LL *C,int n){
    int nt,L=0;
    for(nt=1;nt<=n+n;nt<<=1) L++;
    for(int i=0;i<nt;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    for(int i=0;i<nt;i++) a0[i].r=b0[i].r=c0[i].r=0;
    for(int i=0;i<nt;i++) a0[i].i=b0[i].i=c0[i].i=0;
    for(int i=0;i<n;i++) a0[i].r=A[i],b0[i].r=B[i];
    fft(a0,nt,1);
    fft(b0,nt,1);
    for(int i=0;i<nt;i++) c0[i]=a0[i]*b0[i];
    fft(c0,nt,-1);
    for(int i=0;i<n;i++) C[i]=(c0[i].r+0.5);
}

inline void multpoly1(LL *A,LL *B,LL *C,int n,LL P){
    int nt,L=0;
    for(nt=1;nt<=n+n;nt<<=1) L++;
    for(int i=0;i<nt;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    for(int i=0;i<nt;i++) a1[i]=b1[i]=c1[i]=0;
    for(int i=0;i<n;i++) a1[i]=A[i],b1[i]=B[i];
    fnt(a1,nt,1,P);
    fnt(b1,nt,1,P);
    for(int i=0;i<nt;i++) c1[i]=a1[i]*b1[i]%P;
    fnt(c1,nt,-1,P);
    for(int i=0;i<n;i++) C[i]=c1[i];
}

inline LL CRT(LL a1,LL a2){
    LL ans=0,M=mod*(LL)mod2;
    ans+=qmult(qmult(a1,mod2,M),qpow(mod2,mod-2,mod),M);
    ans+=qmult(qmult(a2,mod,M),qpow(mod,mod2-2,mod2),M);
    return (ans%M+M)%M;
}

int n,phi[N];
LL ans=0,a[N],b[N];
LL c[N],ap[N];
double ti;

int main(){
    scanf("%d",&n);
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(phi[i]) continue;
        for(int j=i;j<=n;j+=i){
            if(!phi[j]) phi[j]=j;
            phi[j]=phi[j]/i*(i-1);
        }
    }
    for(int i=0;i<=100000;i++)
        MT[qpow(gn,i,mod)]=i;
    for(int i=0;i<n;i++) scanf("%lld",&a[i]);
    for(int i=1;i<n;i++) a[i]=(a[i]+a[i-1])%mod;
    for(int i=0;i<n;i++) ap[i]=MT[a[i]];
    for(int i=0;i<n;i++) b[i]=phi[i]%23;
    multpoly2(ap,b,c,n);
    for(int i=0;i<n;i++) c[i]=c[i]%(mod-1);
    for(int i=0;i<n;i++) c[i]=qpow(gn,c[i],mod);
    for(int i=n-1;~i;i--) c[i+1]=c[i],c[i]=0;
    multpoly1(c,c,c,n+1,mod);
    for(int i=0;i<=n;i++) ans=(ans+c[i])%mod;
    printf("%lld
",ans);
    return 0;
}
View Code

相关推荐