• [原创题] 雾


    题目: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
  • 相关阅读:
    《Python》进程收尾线程初识
    《Python》进程之间的通信(IPC)、进程之间的数据共享、进程池
    L02-RHEL6.5环境中安装JDK1.8
    L01-RHEL6.5中部署NTP(ntp server + client)
    P01-Python中列表的复制问题
    数据库模式(三级模式+两级映射)
    事务的四大性质:ACID
    JAVA_接口_默认方法&静态方法
    2018年最新Java面试题及答案整理
    Socket通信原理
  • 原文地址:https://www.cnblogs.com/lawyer/p/4564231.html
Copyright © 2020-2023  润新知