hdu6706 huntian oy

hdu
好久没写数论函数题了，上一次写还是在纪中学min25筛的时候了，赶紧来一道补下手感

题面：求

\[\sum_{i=1}^n\sum_{j=1}^igcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1] \]

保证\(n,a,b\leq 10^9,gcd(a,b)=1\)

知道\((i^a-j^a)\)这个玩意因式分解后有\((i-j)\)，不妨大力猜一下\(gcd(i^a-j^a,i^b-j^b)=i-j\)，证明明天再补
那么

\[\begin{aligned} 原式=&\sum_{i=1}^n\sum_{j=1}^i(i-j)[gcd(i,j)=1]\\ =&\sum_{i=1}^n\sum_{j=1}^i(i-j)\sum_{d|gcd(i,j)}\mu(d)\\ =&\sum_{d=1}^n\mu(d)d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^ii-j\\ \end{aligned} \]

后面那个你拿平方和公式算一下应该有\(\sum_{i=1}^n\sum_{j=1}^ii-j=\frac{n^3-n}{6}\),前面那一部分就是\(\mu·id\)，卷上\(id\)之后就是\(\epsilon\),杜教筛直接求

感觉我回忆杜教筛的时间都比前面推导的时间长。。。

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<vector>
#include<math.h>
#include<queue>
#include<set>
#include<map>
#include<unordered_map>
using namespace std;
typedef long long ll;
typedef long double db;
typedef pair<int,int> pii;
const int N=6000000;
const db pi=acos(-1.0);
#define lowbit(x) (x)&(-x)
#define sqr(x) (x)*(x)
#define rep(i,a,b) for (register int i=a;i<=b;i++)
#define per(i,a,b) for (register int i=a;i>=b;i--)
#define fir first
#define sec second
#define mp(a,b) make_pair(a,b)
#define pb(a) push_back(a)
#define maxd 1000000007
#define inv6 166666668
#define eps 1e-8
int pri[N+10],tot=0,mu[N+10];
int sum[N+10];
bool vis[N+10];
unordered_map<int,ll> has;

int read()
{
    int x=0,f=1;char ch=getchar();
    while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
    while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
    return x*f;
}

void sieve()
{
	mu[1]=1;
	rep(i,2,N)
	{
		if (!vis[i]) {pri[++tot]=i;mu[i]=-1;}
		for (int j=1;j<=tot && i*pri[j]<=N;j++)
		{
			vis[i*pri[j]]=1;
			if (i%pri[j]==0) break;
			else mu[i*pri[j]]-=mu[i];
		}
	}
	rep(i,1,N) sum[i]=(sum[i-1]+i*mu[i]+maxd)%maxd;
}

ll calc(int n)
{
	ll ans=1ll*n*n%maxd*n%maxd;
	ans=(ans-n+maxd)%maxd;
	ans=ans*inv6%maxd;
	return ans;
}

ll query(int n)
{
	if (n<=N) return sum[n];
	if (has[n]) return has[n];
	ll ans=1;int l,r;
	for (l=2;l<=n;l=r+1)
	{
		r=n/(n/l);
		ll tmp=(1ll*(r+l)*(r-l+1)/2)%maxd;
		ans=(ans-tmp*query(n/l)%maxd+maxd)%maxd;
	}
	has[n]=ans;
	return ans;
}
		
ll f(int n,int a,int b)
{
	ll ans=0;
	int l,r;
	for (l=1;l<=n;l=r+1)
	{
		r=n/(n/l);
		ll nowsum=(query(r)-query(l-1)+maxd)%maxd;
		ans=(ans+nowsum*calc(n/l)%maxd)%maxd; 
	}
	return ans;
}

int main()
{
	int T=read();
	sieve();
	while (T--)
	{
		int n=read(),a=read(),b=read();
		printf("%lld\n",f(n,a,b));
	}
	return 0;
}