bzoj 3512: DZY Loves Math IV

题目描述

今有 $n,m\in {\mathbb{N}}^{\ast }$，求$$\sum _{i=1}^n\sum _{j=1}^m\phi (ij)\mathrm{m}\mathrm{o}\mathrm{d}(10^9+7)$$
的值，$n\le 10^5,m\le 10^9$。

Solution

设$$\begin{gathered} S(n,m)=\sum _{i=1}^m\phi (ni) \\ n=\prod p_i^{a_i} \\ p=\prod p_i^{a_i-1} \\ q=\prod p_i \end{gathered}$$则$$\begin{array}{rl}n& =p\times q,\\ S(n,m)& =\sum _{i=1}^m\phi (ni)\\ & =p\sum _{i=1}^m\phi (qi)\\ & =p\sum _{i=1}^m\phi (i)\phi (\frac{q}{\gcd (i,q)})\gcd (i,q)\end{array}$$
因为 $\phi \ast I=id$，所以$$\begin{array}{rl}S(n,m)& =p\sum _{i=1}^m\phi (i)\phi (\frac{q}{\gcd (i,q)})\sum _{d|\gcd (i,q)}\phi (d)\\ & =p\sum _{i=1}^m\phi (i)\sum _{d|i,d|q}\phi (\frac{q}{d})\\ & =p\sum _{d|q}\phi (\frac{q}{d})\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\phi (di)\\ & =p\sum _{d|q}\phi (\frac{q}{d})S(d,\lfloor \frac{m}{d}\rfloor )\end{array}$$
贴上我的常数大代码。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<map>
#include<ctime>

#define reg register
const int MAXN=5000010;
const int MOD=1e9+7;
//const int mod=1e6+7;
typedef long long ll;

bool vis[MAXN+10];
int p[MAXN+10];
int phi[MAXN+10];
ll S_phi[MAXN+10];
int len=0;
ll n,m;
int mn[MAXN+10];
//int ck=0,cnt=0;
//std::map<ll,std::map<ll,ll> >mp;
//int hash[mod+10];
std::map<ll,ll>mp_phi;

void init(){
	memset(vis,1,sizeof(vis));vis[1]=0;
	phi[1]=1;
	for(reg int i=2;i<=MAXN;++i){
		if(vis[i]){
			p[++len]=i;
			phi[i]=i-1;mn[i]=i;
		}
		for(reg int j=1;(j<=len)&&(i*p[j]<=MAXN);++j){
			vis[i*p[j]]=0;mn[i*p[j]]=p[j];
			if(i%p[j])
				phi[i*p[j]]=phi[i]*(p[j]-1);
			else{
				phi[i*p[j]]=phi[i]*p[j];
				break;
			}
		}
	}
	S_phi[0]=0;
//	for(reg int i=1;i<=15;++i) printf("%d ",mn[i]);
	for(reg int i=1;i<=MAXN;++i){
		phi[i]%=MOD;
		S_phi[i]=(S_phi[i-1]+phi[i])%MOD;
	}
}
ll Sphi(ll x){
	if(x<MAXN) return S_phi[x];
//	if(hash[x%mod]) return hash[x%mod];
	if(mp_phi[x]) return mp_phi[x];
	x%=MOD;
	int s=x*(x+1)/2%MOD,sum=0;
	for(reg int l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		sum=(sum+(r-l+1)*Sphi(x/l)%MOD)%MOD;
	}
//	return hash[x%mod]=((s-sum+MOD)%MOD);
	return mp_phi[x]=((s-sum+MOD)%MOD);
}
ll S(ll x,ll y){
	if(!x||!y) return 0;
	if(x==1) return Sphi(y);
	if(y==1) return (Sphi(x)-Sphi(x-1)+MOD)%MOD;
//	if(mp[x][y]) return mp[x][y];
//	++cnt;
//	printf("%lld %lld\n",x,y);
	int xx=x,P,Q=1,l=0,pr[10];
//	int st=clock();
	while(x>1){
		int mx=mn[x];
		Q*=mx;
		while(x%mx==0) x/=mx;
		pr[++l]=mx;
	}
	P=(int)(xx/Q);
	ll sum=0,d;
	for(reg int i=0;i<(1<<l);++i){
		d=1;
		for(reg int j=1;j<=l;++j)
			if(i&(1<<(j-1))) d*=pr[j];
//		if(Q%d==0){
			sum=(sum+phi[Q/d]*S(d,y/d)%MOD)%MOD;
//			sum=(sum+phi[d]*S(Q/d,y/(Q/d))%MOD)%MOD;
//		}
	}
	return P*sum%MOD;
}
int main(){
	init();
	scanf("%lld%lld",&n,&m);
	ll ans=0;
	for(reg int i=1;i<=n;++i)
		ans=(ans+S(i,m))%MOD;
	printf("%lld",ans);
//    printf("\n%d %d",ck,cnt);
}

AC

#include<iostream>
#include<cstdio>
#include<vector>
#include<map>
#include<ctime>
using namespace std;
#define MOD 1000000007
#define MAX 2000000
int n,m,ans;
int pri[MAX],tot,phi[MAX],mn[MAX];
bool zs[MAX];
int ck=0,cnt=0;
void pre(int n)
{
    phi[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!zs[i])pri[++tot]=i,phi[i]=i-1,mn[i]=i;
        for(int j=1;i*pri[j]<=n&&j<=tot;++j)
        {
            zs[i*pri[j]]=true;mn[i*pri[j]]=pri[j];
            if(i%pri[j])phi[i*pri[j]]=phi[i]*phi[pri[j]];
            else{phi[i*pri[j]]=phi[i]*pri[j];break;}
        }
    }
    for(int i=1;i<=15;++i) printf("%d ",mn[i]);
    for(int i=1;i<=n;++i)phi[i]=(phi[i-1]+phi[i])%MOD;
}
map<int,int> M,S[100010];
int Phi(int n)
{
    if(n<MAX)return phi[n];
    if(M[n])return M[n];
    int ret=1ll*n*(n+1)/2%MOD;
    for(int i=2,j;i<=n;i=j+1)
    {
        j=n/(n/i);
        ret=(ret+MOD-1ll*(j-i+1)*Phi(n/i)%MOD)%MOD;
    }
    return M[n]=ret;
}
int Solve(int n,int m)
{
    if(!m)return 0;
    if(n==1)return Phi(m);
    if(m==1)return (Phi(n)-Phi(n-1)+MOD)%MOD;
    if(S[n][m])return S[n][m];
	++cnt;
    int ret=0;
    vector<int> fac;int p=1,q=1,N=n;
    int st=clock();
    while(N>1)
    {
        int x=mn[N];q*=x;N/=x;fac.push_back(x);
        while(N%x==0)p*=x,N/=x;
    }
    ck+=clock()-st;
    for(int i=0,l=fac.size();i<1<<l;++i)
    {
        int d=1;
        for(int j=0;j<l;++j)if(i&(1<<j))d*=fac[j];
        ret=(ret+1ll*(Phi(q/d)-Phi(q/d-1)+MOD)*Solve(d,m/d))%MOD;
    }
    return S[n][m]=1ll*ret*p%MOD;
}
int main()
{
    scanf("%d%d",&n,&m);pre(MAX);
    for(int i=1;i<=n;++i)ans=(ans+Solve(i,m))%MOD;
    printf("%d\n",ans);
    printf("\n%d %d",ck,cnt);
    return 0;
}