【JZOJ4872】太阳神 题解

题目大意

  求这玩意儿：

$$\sum_{a=1}^n\sum_{b=1}^n~[~lcm(a,b)>n~]$$

  $n \le 10^{10}$

【首先】

  显然算出 $lcm<=n$ 的数对的数量，然后用 $n^2$ 去减。以下 $ans$ 表示的都是小于等于的数量。

【20%】n<=2000

  暴力。

【40%】n<=10^7

$$\begin{array}{rcl} ans&=&\sum_{a=1}^n\sum_{b=1}^n~[~\frac{ab}{gcd(a,b)}<=n~]\\ &=&\sum_{d=1}^n\sum_{a'=1}^{\lfloor\frac{n}{d^2}\rfloor}\sum_{b'=1}^{\lfloor\frac{n}{d^2}\rfloor}~[~a'b'd<=n~]\\ &=&\sum_{d=1}^n\sum_{t=1}^{\lfloor\frac{n}{d^2}\rfloor}2^{g(t)}\\ 其中~g(t)~表示~t~的质因子种类数 \end{array}$$

【60%】n<=10^8

\begin{array}{rcl} ans&=&\sum_{d=1}^n \sum_{a=1}^n \sum_{b=1}^n~[~gcd(a,b)==d~,~\frac{ab}{d}<=n~]\\ &=&\sum_{d=1}^n \sum_{a'} \sum_{b'}~[~gcd(a',b')==1~,~a'b'd<=n~]\\ &=&\sum_{D=1}^\sqrt n \mu(D) \sum_d \sum_{a''} \sum_{b''}~[~a''b''d<=\lfloor\frac{n}{D^2}\rfloor~]~\\ &=&\sum_{D=1}^\sqrt n \mu(D) \sum_{d=1}^{\lfloor\frac{n}{D^2}\rfloor} \sum_{a''=1}^{\lfloor\frac{n}{D^2d}\rfloor} \lfloor\frac{n}{D^2da''}\rfloor \end{array}

  根号下分块套分块，时间我不会算，大概是非常不满的三个根号，反正不能过就是了。

【100%】n<=10^10

  设 $m=\lfloor\frac{n}{D^2d}\rfloor$，则上述式子的最后一个sigma实际上是在求 $1$~$m$ 的约数个数和。这个东西在合理范围内是可以预处理的。
  因此对于最后一个sigma，如果 $m \le 10^7$，那么我们直接用预处理的结果，否则再根号算。由于上限下降很快，所以要用根号算的其实很少。

代码

#include<cmath>
#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;

typedef long long LL;

const LL mo=1e9+7, maxn=1e7+5;

LL n;

int p0,p[maxn],miu[maxn],fy[maxn],fy1[maxn];
LL Sfy[maxn];
bool bz[maxn];
void Pre(int n)
{
	miu[1]=fy[1]=1;
	fo(i,2,n)
	{
		if (!bz[i])
		{
			p[++p0]=i;
			miu[i]=-1;
			fy[i]=2, fy1[i]=1;
		}
		fo(j,1,p0)
		{
			if ((LL)i*p[j]>n) break;
			bz[i*p[j]]=1;
			if (i%p[j]==0)
			{
				miu[i*p[j]]=0;
				fy[i*p[j]]=fy[i]/(fy1[i]+1)*(fy1[i]+2);
				fy1[i*p[j]]=fy1[i]+1;
				break;
			} else
			{
				miu[i*p[j]]=-miu[i];
				fy[i*p[j]]=fy[i]*2;
				fy1[i*p[j]]=1;
			}
		}
	}
	fo(i,1,n) Sfy[i]=(Sfy[i-1]+fy[i])%mo;
}

LL cala(LL x)
{
	LL re=0;
	for(LL ai=1, aj; ai<=x; ai=aj+1)
	{
		aj=x/(x/ai);
		re=(re+x/ai*(aj-ai+1)%mo)%mo;
	}
	return re;
}

int main()
{
	Pre(maxn-5);
	
	scanf("%lld",&n);
	LL sqrtn=sqrt(n);
	
	
	LL ans=0;
	for(LL D=1; D<=sqrtn; D++)
	{
		LL maxd=n/D/D, ansD=0;
		for(LL di=1, dj; di<=maxd; di=dj+1)
		{
			dj=maxd/(maxd/di);
			LL maxa=maxd/di;
			
			LL t=(maxa<=maxn-5) ?Sfy[maxa] :cala(maxa) ;
			ansD=(ansD+t*(dj-di+1)%mo)%mo;
		}
		ans=(ans+(ansD*miu[D]%mo+mo)%mo)%mo;
	}
	
	n%=mo;
	printf("%lld\n",(n*n%mo-ans+mo)%mo);
}