1. 程式人生 > 實用技巧 >【洛谷P4449】於神之怒加強版

【洛谷P4449】於神之怒加強版

題目

題目連結:https://www.luogu.com.cn/problem/P4449
給定 \(n,m,k\),計算

\[\sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k \]

\(10^9+7\) 取模的結果。

思路

菜到真的只會模板題了 /kk。

\[\sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k \]

\[=\sum^{n}_{d=1}(\sum^n_{d|i}\mu(\frac{i}{d})\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i} \rfloor\times d^k) \]

\[=\sum^{n}_{i=1}\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i} \rfloor(\sum_{d|i}\mu(\frac{i}{d})\times d^k) \]

前面整除分塊,後面預處理即可。
時間複雜度 \(O(n\log n+T\sqrt{n})\)

程式碼

#include <bits/stdc++.h>
#define reg register
using namespace std;
typedef long long ll;

const int N=5000010,MOD=1e9+7;
int Q,cnt,n,m,t,mu[N],prm[N],sum[N];
ll ans;
bool v[N];

ll fpow(ll x,ll k)
{
	ll ans=1;
	for (;k;k>>=1,x=x*x%MOD)
		if (k&1) ans=ans*x%MOD;
	return ans;
}

void findprm(int n)
{
	mu[1]=1;
	for (reg int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++cnt]=i,mu[i]=-1;
		for (reg int j=1;j<=cnt;j++)
		{
			if (i>n/prm[j]) break;
			v[i*prm[j]]=1; mu[i*prm[j]]=-mu[i];
			if (i%prm[j]==0)
			{
				mu[i*prm[j]]=0;
				break;
			}
		}
	}
}

int main()
{
	scanf("%d%d",&Q,&t);
	findprm(N-10);
	for (reg int i=1;i<=N-10;i++)
	{
		ll power=fpow(i,t);
		for (reg int j=i;j<=N-10;j+=i)
			sum[j]=(sum[j]+power*mu[j/i])%MOD;
	}
	for (reg int i=1;i<=N-10;i++)
		sum[i]=(sum[i]+sum[i-1])%MOD;
	while (Q--)
	{
		ans=0;
		scanf("%d%d",&n,&m);
		for (reg int l=1,r;l<=min(n,m);l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans=(ans+1LL*(n/l)*(m/l)%MOD*(sum[r]-sum[l-1]))%MOD;
		}
		printf("%lld\n",(ans%MOD+MOD)%MOD);
	}
	return 0;
}