HZOJ 太阳神

所以我刚学反演还没学反演就要做这么一道神仙题……

首先大于n不好求，补集转化。

$ans=n*n-\sum \limits _{i=1}^{n} \sum \limits _{j=1}^{n} \left [ lcm(i,j)\leqslant n\right ] $

那么我们要求：

$\sum \limits _{i=1}^{n} \sum \limits _{j=1}^{n} \left [ \frac{i*j}{gcd(i,j) } \leqslant n \right ]$

枚举d=gcd(i，j),

原式=$\sum \limits _{d=1}^{n} \sum \limits _{i=1}^{n} \sum \limits _{j=1}^{n} \left [ i*j*d\leqslant n ,gcd(i,j)=1 \right ]$

=$\sum \limits _{d=1}^{n} \sum \limits _{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor} \sum \limits _{j=1}^{\left \lfloor \frac{n}{d*i} \right \rfloor} \left [ gcd(i,j)=1 \right ]$

根据莫比乌斯函数的性质：$\sum \limits _{d\mid n}u(d) =\left [ n=1 \right ]$

于是原式=$\sum \limits _{d=1}^{n} \sum \limits _{i=1}^{n} \sum \limits _{j=1}^{n} \sum \limits _{g\mid gcd(i,j)} u(g)$

所以就要反演了？其实就是交换求和的顺序。

个人这步稍难理解（因为我没学过反演），将g提前后相当于求u（g）出现的次数，那么修改g的定义，令${i}'=\frac{i}{g},{j}'=\frac{j}{g}$.

原式=$\sum \limits _{d=1}^{n} \sum \limits _{g=1} u(g) \sum \limits _{{i}'=1}^{\left \lfloor \frac{n}{d*g} \right \rfloor} \sum \limits _{{j}'=1}^{\left \lfloor \frac{n}{d*i*g} \right \rfloor} 1$

将g提前，原式=$\sum \limits _{g=1}^{\sqrt{n}}u(g) \sum \limits _{{i}'=1} \sum \limits _{{j}'=1} \sum \limits _{d=1} \left [ {i}'*{j}'*d\leqslant \frac{n}{g*g} \right ]$

到此式子就推完了，可是看起来还是不是很可做……但是可以发现g是根号n范围内的，u线筛即可，同时枚举g。

不妨设${i}'\leqslant {j}'\leq d$,那么设$m=\frac{n}{g*g}$可以以$O\left ( m^{\frac{1}{3}} \right )$的复杂度枚举i，以$\sqrt{\frac{m}{i}}$的复杂度枚举j，O1算出d的个数，之后乘$A_3^3$。

但是要考虑算重的情况，手动讨论一下就行了。

 #include<iostream>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #define int LL

 #define LL long long

 using namespace std;

 const int mod=1e9+;

 LL n;

 bool isprime[];

 int prime[],cnt,mu[];

 void shai(int n)

 {

     mu[]=;

     for(int i=;i<=n;i++)isprime[i]=;

     for(int i=;i<=n;i++)

     {

         if(isprime[i])mu[i]=-,prime[++cnt]=i;

         for(int j=;j<=cnt&&prime[j]*i<=n;j++)

         {

             isprime[prime[j]*i]=;

             if(i%prime[j]==)break;

             else mu[i*prime[j]]=-mu[i];

         }

     }

 }

 signed main()

 {

     cin>>n;int maxn=sqrt(n);shai(maxn+);

     LL ans=;

     for(int i=;i<=maxn;i++)

     {

         LL res=;

         int m=n/(i*i);

         for(int a=;a*a*a<=m;a++)

         {

             int maxb=sqrt((1.0*m)/a);

             for(int b=a;b<=maxb;b++)

             if(m/(a*b)>=b)

             {

                 if(a==b)res=(res+(m/(a*b)-b)*+)%mod;

                 else    res=(res+(m/(a*b)-b)*+)%mod;

             }

         }

         ans=(ans+mu[i]*res%mod)%mod;

     }

     printf("%lld\n",(n%mod*(n%mod)%mod-ans%mod+mod)%mod);

 }

巴特西

HZOJ 太阳神

最新文章

热门文章