【JZOJ4872】太阳神 题解

题目大意

$~~~~~~$求这玩意儿：

$$\sum_{a=1}^n\sum_{b=1}^n~[~lcm(a,b)>n~]$$
$~~~~~~$n<=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<=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);
}