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$。

可是要考慮算重的狀況，手動討論一下就好了。

 1 #include<iostream>
 2 #include<cstring>
 3 #include<cstdio>
 4 #include<cmath>
 5 #define int LL
 6 #define LL long long
 7 using namespace std;
 8 const int mod=1e9+7;
 9 LL n;
10 bool isprime[100010];
11 int prime[100010],cnt,mu[100010];
12 void shai(int n)
13 {    
14     mu[1]=1;
15     for(int i=2;i<=n;i++)isprime[i]=1;
16     for(int i=2;i<=n;i++)
17     {
18         if(isprime[i])mu[i]=-1,prime[++cnt]=i;
19         for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
20         {
21             isprime[prime[j]*i]=0;
22             if(i%prime[j]==0)break;
23             else mu[i*prime[j]]=-mu[i];
24         }
25     }
26 }
27 signed main()
28 {
29     cin>>n;int maxn=sqrt(n);shai(maxn+1);
30     LL ans=0;
31     for(int i=1;i<=maxn;i++)
32     {
33         LL res=0;
34         int m=n/(i*i);
35         for(int a=1;a*a*a<=m;a++)
36         {
37             int maxb=sqrt((1.0*m)/a);
38             for(int b=a;b<=maxb;b++)
39             if(m/(a*b)>=b)
40             {
41                 if(a==b)res=(res+(m/(a*b)-b)*3+1)%mod;
42                 else    res=(res+(m/(a*b)-b)*6+3)%mod;
43             }
44         }
45         ans=(ans+mu[i]*res%mod)%mod;
46     }
47     printf("%lld\n",(n%mod*(n%mod)%mod-ans%mod+mod)%mod);
48 }

View Code