zoukankan      html  css  js  c++  java
  • 【51NOD 1847】奇怪的数学题(莫比乌斯反演,杜教筛,min_25筛,第二类斯特林数)

    【51NOD 1847】奇怪的数学题(莫比乌斯反演,杜教筛,min_25筛,第二类斯特林数)

    题面

    51NOD

    [sum_{i=1}^nsum_{j=1}^nsgcd(i,j)^k ]

    其中(sgcd)表示次大公约数。

    题解

    明摆着(sgcd)就是在(gcd)的基础上除掉(gcd)的最小因数。
    所以直接枚举(gcd)

    [egin{aligned} ans&=sum_{i=1}^nsum_{j=1}^n sgcd(i,j)^k\ &=sum_{i=1}^nsum_{j=1}^n (frac{gcd(i,j)}{min_p(gcd(i,j))})^k\ &=sum_{d=1}^n (frac{d}{min_p(d)})^ksum_{i=1}^nsum_{j=1}^n [gcd(i,j)=1]\ &=sum_{d=1}^n (frac{d}{min_p(d)})^ksum_{i=1}^{n/d}mu(i)[frac{n}{id}]^2\ &=sum_{T=1}^n [frac{n}{T}]^2sum_{d|T}(frac{d}{min_p(d)})^kmu(frac{T}{d}) end{aligned}]

    好啦好啦,又有一个看起来要求前缀和的狄利克雷卷积啦。
    (displaystyle f(d)=(frac{d}{min_p(d)})^k)
    那么后面那个东西就是((f*mu)(T))
    要求它的前缀和,emmmm,杜教筛。
    看到(mu),emmm,构造(g(x)=1)
    那么令(displaystyle S(n)=sum_{i=1}^n (f*mu)(i))
    那么套杜教筛的式子:

    [g(1)S(n)=sum_{i=1}^n (f*mu*g)(i)-sum_{i=2}^n g(i)S([frac{n}{i}]) ]

    (g(x)=1)带进去就得到了:

    [S(n)=sum_{i=1}^n f(i)-sum_{i=2}^nS([frac{n}{i}]) ]

    行,来筛(f)前缀和,涉及到了最小质因子???(min\_25)筛。
    考虑一下(f)的前缀和是什么东西,首先所有质数的贡献都是(1),这里需要计算质数个数。
    然后就是枚举最小质因子求剩下部分的(k)次幂的和就行了。
    而很有趣的一点就是(f)所要求的所有合数的贡献,恰好就是计算所有质数(k)次方贡献时减去的部分。所以直接加上就好啦。
    然后这里涉及到了自然数幂和的问题。推导戳这里

    [sum_{i=1}^n i^k=sum_{j=0}^kegin{Bmatrix}k\jend{Bmatrix}frac{(n+1)^{underline {j+1}}}{j+1} ]

    这样一来就可以直接算啦。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    #define uint unsigned int
    #define MAX 100000
    uint n,blk;int k;
    bool zs[MAX];
    int pri[MAX],tot;
    uint prik[MAX],sprik[MAX];
    uint w[MAX],g[MAX],h[MAX],f[MAX];
    int id1[MAX],id2[MAX],m;
    uint S[60][60];
    int getid(uint x){return (x<=blk)?id1[x]:id2[n/x];}
    uint fpow(uint a,int b)
    {
    	uint s=1;
    	while(b){if(b&1)s*=a;a*=a;b>>=1;}
    	return s;
    }
    void pre(int n)
    {
    	for(int i=2;i<=n;++i)
    	{
    		if(!zs[i])pri[++tot]=i,prik[tot]=fpow(i,k);
    		for(int j=1;j<=tot&&i*pri[j]<=n;++j)
    		{
    			zs[i*pri[j]]=true;
    			if(i%pri[j]==0)break;
    		}
    	}
    	for(int i=1;i<=tot;++i)sprik[i]=sprik[i-1]+prik[i];
    }
    uint Sumk(int n)
    {
    	uint ret=0;
    	for(int i=1;i<=k;++i)
    	{
    		uint s=1;
    		for(int j=0;j<=i;++j)
    			if((n+1-j)%(i+1))s*=n+1-j;
    			else s*=(n+1-j)/(i+1);
    		ret+=S[k][i]*s;
    	}
    	return ret;
    }
    bool vis[MAX];
    uint M[MAX];
    uint Solve(uint n)
    {
    	if(vis[getid(n)])return M[getid(n)];
    	uint ret=f[getid(n)];
    	for(uint i=2,j;i<=n;i=j+1)
    		j=n/(n/i),ret-=(j-i+1)*Solve(n/i);
    	vis[getid(n)]=true;return M[getid(n)]=ret;
    }
    int main()
    {
    	scanf("%u%d",&n,&k);pre(blk=sqrt(n));S[0][0]=1;
    	for(int i=1;i<=k;++i)
    		for(int j=1;j<=i;++j)
    			S[i][j]=S[i-1][j-1]+j*S[i-1][j];
    	for(uint i=1,j;i<=n;i=j+1)
    	{
    		j=n/(n/i);w[++m]=n/i;
    		g[m]=w[m]-1;h[m]=Sumk(w[m])-1;
    		if(w[m]<=blk)id1[w[m]]=m;
    		else id2[n/w[m]]=m;
    	}
    	for(int j=1;j<=tot;++j)
    		for(int i=1;i<=m&&1u*pri[j]*pri[j]<=w[i];++i)
    		{
    			int k=getid(w[i]/pri[j]);
    			g[i]-=g[k]-(j-1);
    			h[i]-=prik[j]*(h[k]-sprik[j-1]);
    			f[i]+=h[k]-sprik[j-1];
    		}
    	for(int i=1;i<=m;++i)f[i]+=g[i];
    	uint ans=0;
    	for(uint i=1,j,lt=0,nw;i<=n;i=j+1)
    	{
    		j=n/(n/i);uint s=n/i;s*=s;
    		nw=Solve(j);ans+=s*(nw-lt);lt=nw;
    	}
    	printf("%u
    ",ans);
    	return 0;
    }
    
  • 相关阅读:
    网络七层
    微信小程序开发工具 常用快捷键
    BZOJ 1026 windy数 (数位DP)
    BZOJ 1026 windy数 (数位DP)
    CodeForces 55D Beautiful numbers (SPOJ JZPEXT 数位DP)
    CodeForces 55D Beautiful numbers (SPOJ JZPEXT 数位DP)
    HDU 3709 Balanced Number (数位DP)
    HDU 3709 Balanced Number (数位DP)
    UVA 11361 Investigating Div-Sum Property (数位DP)
    UVA 11361 Investigating Div-Sum Property (数位DP)
  • 原文地址:https://www.cnblogs.com/cjyyb/p/10172276.html
Copyright © 2011-2022 走看看