zoukankan      html  css  js  c++  java
  • BZOJ 4407 于神之怒加强版

    题目链接:于神之怒加强版

      这个式子还是很妙的,只是我已经思维僵化了

    egin{aligned}
     &sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k \
    =&sum_{g=1}^ng^ksum_{i=1}^{lfloor frac{n}{g} floor}sum_{j=1}^{lfloor frac{m}{g} floor}sum_{d|i,d|j}mu(d) \
    =&sum_{g=1}^ng^ksum_{d=1}^{lfloor frac{n}{g} floor}mu(d)lfloor frac{n}{dg} floorlfloor frac{m}{dg} floor \
    =&sum_{g=1}^ng^ksum_{g|x}^{n}mu(frac{x}{g}) lfloor frac{n}{x} floorlfloor frac{m}{x} floor \
    =&sum_{x=1}^n lfloor frac{n}{x} floorlfloor frac{m}{x} floor  sum_{g|x} g^k mu(frac{x}{g})
    end{aligned}

      然后我们可以令(f(x)=sum_{g|x} g^k mu(frac{x}{g})),显然(f(x))是(id^k)与(mu)的狄利克雷卷积,所以(f(x))也是一个积性函数

      由于(n)的范围有(5 imes 10^6),所以我们来考虑一下如何筛这个(f)函数,也就是考虑在(i)是(p)的倍数的情况下计算(f(ip))((p)为质数)

      令(i=p^kq(qperp p)),然后分情况讨论一下:

      当(q e 1)时,由于(i)的最小质因子为(p),那么(p^{k+1}<i),所以(f(frac{i}{p^k}) imes f(p^{k+1}))就是所求;

      当(q=1)时,由于(f(p^x)=p^{kx}-p^{kx-k}),所以(f(i) imes p^k)即为所求,(p^k)则可以通过预处理所有质数的(k)次方解决。

      最后再分块计算就好了。

      下面贴代码:

    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cmath>
    #define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
    #define maxn 5000010
    #define mod 1000000007
    
    using namespace std;
    typedef long long llg;
    
    int T,k,n,m,pr[maxn],lp;
    int g[maxn],c[maxn];
    bool vis[maxn];
    llg f[maxn];
    
    void gi(llg &x){if(x>=mod) x%=mod;}
    llg mi(llg a,int b){
    	llg s=1;
    	while(b){
    		if(b&1) s=s*a,gi(s);
    		a=a*a,gi(a); b>>=1;
    	}
    	return s;
    }
    
    int main(){
    	File("a");
    	scanf("%d %d",&T,&k); f[1]=c[1]=1;
    	for(int i=2;i<maxn;i++){
    		if(!vis[i]) pr[++lp]=i,g[i]=i,c[i]=mi(i,k),f[i]=c[i]-1;
    		for(int j=1;i*pr[j]<maxn;j++){
    			if(i*pr[j]==31823){
    				int aa;
    				aa++;
    			}
    			vis[i*pr[j]]=1;
    			if(i%pr[j]) g[i*pr[j]]=pr[j],f[i*pr[j]]=f[i]*f[pr[j]]%mod;
    			else{
    				f[i*pr[j]]=g[i]!=i?f[i/g[i]]*f[g[i]*pr[j]]:f[i]*c[pr[j]];
    				gi(f[i*pr[j]]); g[i*pr[j]]=g[i]*pr[j]; break;
    			}
    		}
    	}
    	for(int i=2;i<maxn;i++) f[i]+=f[i-1],gi(f[i]);
    	while(T--){
    		scanf("%d %d",&n,&m);
    		if(n>m) n^=m^=n^=m;
    		llg ans=0;
    		for(int i=1,nt;i<=n;i=nt+1){
    			nt=min(n/(n/i),m/(m/i));
    			ans+=(f[nt]-f[i-1]+mod)*(n/i)%mod*(m/i)%mod;
    			if(ans>=mod) ans%=mod;
    		}
    		printf("%lld
    ",ans);
    	}
    	return 0;
    }
    
  • 相关阅读:
    topcoder srm 681 div1
    topcoder srm 683 div1
    topcoder srm 684 div1
    topcoder srm 715 div1
    topcoder srm 685 div1
    topcoder srm 687 div1
    topcoder srm 688 div1
    topcoder srm 689 div1
    topcoder srm 686 div1
    topcoder srm 690 div1 -3
  • 原文地址:https://www.cnblogs.com/lcf-2000/p/6921584.html
Copyright © 2011-2022 走看看