Easy Math Problem[hdu6607]

题目链接

链接

题意

求 $\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)^k\text{lcm}(i,j)[\gcd(i,j)\in \mathbb{P}]\pmod{10^9+7}$.

$n\leq 10^{10}$

题解

化简

化简一波：

$$\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)^k\text{lcm}(i,j)[\gcd(i,j)\in \mathbb{P}]\=\sum_{d=1}^nd^{k-1}[d\in \mathbb{P}]\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)==d]\=\sum_{d=1}^nd^{k+1}[d\in \mathbb{P}]\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij[\gcd(i,j)==1]$$

设 $f(n)=\sum_{i=1}^n\sum_{j=1}^{n}ij[\gcd(i,j)==1]$。

则会有：$f(n)=2\sum_{i=1}^ni\sum_{j=1}^{i-1}j[\gcd(i,j)==1]+\sum_{i=1}^ni^2[\gcd(i,i)==1]$

由经典恒等式：$\sum_{j=1}^{n-1}j[\gcd(n,j)==1]=\frac{n\varphi(n)-[n==1]}{2}$ 可得：

$$f(n)=2(\sum_{i=1}^ni\frac{i\varphi(i)-[i==1]}{2})-1=\sum_{i=1}^{n}i^2\varphi(i)$$

故原式可化简为 ：$\sum_{d=1}^nd^{k+1}[d\in \mathbb{P}]f(\left \lfloor \frac{n}{d} \right \rfloor)$

显然整除分块。现在变成求两部分。

Part 1

对于所有的 $\left \lfloor \frac{n}{i} \right \rfloor$，求 $f(\left \lfloor \frac{n}{i} \right \rfloor)$。

显然杜教筛。

设函数 $g(i)=i^2$，则有：

$$(f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})=\sum_{d|n}d^2\varphi(i)\frac{n^2}{d^2}=n^2\sum_{d|n}\varphi(d)=n^3$$

设 $S(i)=\sum_{i=1}^nf(i)$，那么就有：$S(n)g(1)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S(\left \lfloor \frac{n}{i} \right \rfloor)$

即：$S(n)=\sum_{i=1}^ni^3-\sum_{i=1}^ni^2S(\left \lfloor \frac{n}{i} \right \rfloor)$。

这样就可以求出来了。

时间复杂度 $O(n^{\frac{2}{3}})$。

Part2

对于所有的 $\left \lfloor \frac{n}{i} \right \rfloor$，求 $\sum_{i=1}^ni^{k+1}[i\in \mathbb{P}]$。

Min25筛模板题。

求 $g(n,0)$ 的时候需要求自然数的 $k+1$ 次幂和，拉格朗日插值即可。

时间复杂度 $O(k\sqrt{n}+\frac{n^{\frac{3}{4}}}{\log n})$

程序

#include <set>
#include <map>
#include <queue>
#include <stack>
#include <deque>
#include <cmath>
#include <bitset>
#include <cstdio>
#include <vector>
#include <string>
#include <complex>
#include <utility>
#include <cstring>
#include <iostream>
#include <assert.h>
#include <algorithm>
//#include <unordered_set>
//#include <unordered_map>
using namespace std;
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define ll long long
#define ull unsigned ll
#define db double
#define lb long db
#define pb push_back
#define com complex<db>
#define mp(x,y) make_pair((x),(y))
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define cle(x) memset(x,0,sizeof(x))
#define lowbit(x) ((x)&(-(x)))
#define bit(x,i) (((x)>>(i))&1)
#define fo(i,j,k) for(int i=(j),end_i=(k);i<=end_i;i++)
#define fd(i,j,k) for(int i=(j),end_i=(k);i>=end_i;i--)
inline ll read()
{
	ll x=0; char ch=getchar(); bool f=0;
	for(;ch<'0'||ch>'9';ch=getchar()) if(ch=='-') f=1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+(ch^48);
	return f?-x:x;
}
#define CASET fo(___,1,read())
const ll mod=1e9+7;
inline ll Add(ll x,ll y){x+=y; return (x<mod)?x:x-mod;}
inline ll Dec(ll x,ll y){x-=y%mod; return (x<0)?x+mod:x;}
inline ll Mul(ll x,ll y){return x*y%mod;}
inline ll Pow(ll x,ll y)
{
	y%=(mod-1);
	ll ans=1;
	for(;y;y>>=1,x=x*x%mod) if(y&1) ans=ans*x%mod;
	return ans;
}
const ll inv4=Pow(4,mod-2);
const ll inv6=Pow(6,mod-2);
const int N=4100000;
const int M=500010;
const int K=111;
ll pri[N>>3],phi[N],sg[N];
int tot;
bool vis[N];
void init_prime(int n)
{
	phi[1]=1;
	ll t;
	fo(i,2,n)
	{
		if(!vis[i]) {pri[++tot]=i; phi[i]=i-1;}
		for(int j=1;j<=tot&&(t=(ll)pri[j]*i)<=n;j++)
		{
			vis[t]=1;
			if(i%pri[j]==0) {phi[t]=phi[i]*pri[j]; break;}
			phi[t]=phi[i]*(pri[j]-1);
		}
	}
	fo(i,1,n) sg[i]=Add(sg[i-1],1ll*i*i%mod*phi[i]%mod);
}
ll fac[K],inv[K];
void init_fac(int n)
{
	fac[0]=1;
	fo(i,1,n) fac[i]=fac[i-1]*i%mod;
	inv[n]=Pow(fac[n],mod-2);
	fd(i,n,1) inv[i-1]=inv[i]*i%mod;
}
ll n,ans,Sqr; int k;
map<ll,ll> ma;
inline ll S3(ll n)
{
	n%=mod;
	n=n*(n+1)/2%mod;
	return n*n%mod;
}
inline ll S2(ll n)
{
	n%=mod;
	return n*(n+1)%mod*(n*2+1)%mod*inv6%mod;
}
inline ll S(ll n)
{
	if(n<N) return sg[n];
	if(ma[n]!=0) return ma[n];
	ll ans=S3(n);
	for(ll i=2,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ans=Dec(ans,S(n/i)%mod*(S2(j)-S2(i-1)+mod)%mod);
	}
	return ma[n]=ans;
}
ll pre[K],suf[K],y[K];
inline ll Lagrange(ll n,int k,ll sum=0)
{
	n%=mod; pre[0]=1; suf[k+1]=1;
	fo(i,1,k) pre[i]=Mul(pre[i-1],(n-i+mod)%mod);
	fd(i,k,1) suf[i]=Mul(suf[i+1],(n-i+mod)%mod);
	fo(i,1,k) sum=Add(sum,Mul(pre[i-1],suf[i+1])*Mul(Mul(y[i],inv[i-1]),Mul(inv[k-i],((k-i)&1)?mod-1:1))%mod);
	return sum;
}
int m;
int id1[M],id2[M];
ll g[M],w[M];
void Min25()
{
	m=0;
	fo(i,1,k+2) y[i]=Add(y[i-1],Pow(i,k));
	for(ll i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i); w[++m]=n/i;
		if(w[m]<=Sqr) id1[w[m]]=m; else id2[n/w[m]]=m;
		g[m]=Lagrange(w[m],k+2)-1;
	}
	ll qm=0;
	g[m+1]=0;
	fo(i,1,tot)
	{
		ll tmp=pri[i]*pri[i],p=Pow(pri[i],k),t;
		for(int j=1;j<=m&&tmp<=w[j];j++)
		{
			t=w[j]/pri[i];
			t=(t<=Sqr)?id1[t]:id2[n/t];
			g[j]=Dec(g[j],p*(g[t]-qm+mod)%mod);
		}
		qm=Add(qm,p);
	}
}
int main()
{
	init_prime(N-1);
	init_fac(K-5);
	CASET
	{
		n=read(); k=read()+1;
		Sqr=sqrt(n); ans=0; ma.clear();
		Min25();
		for(ll i=1,j,l=m,h;i<=n;i=j+1,l--)
		{
			j=n/(n/i);
			h=(g[l]-g[l+1]+mod)%mod;
			ans=Add(ans,Mul(S(n/i),h));
		}
		printf("%lld\n",ans%mod);
	}
	return 0;
}