莫比乌斯反演前缀和除法分块优化（hdu4746）

C - Mophues

 HDU - 4746 

As we know, any positive integer C ( C >= 2 ) can be written as the multiply of some prime numbers: 
    C = p1×p2× p3× ... × pk 
which p1, p2 ... pk are all prime numbers.For example, if C = 24, then: 
    24 = 2 × 2 × 2 × 3 
    here, p1 = p2 = p3 = 2, p4 = 3, k = 4 

Given two integers P and C. if k<=P( k is the number of C's prime factors), we call C a lucky number of P. 

Now, XXX needs to count the number of pairs (a, b), which 1<=a<=n , 1<=b<=m, and gcd(a,b) is a lucky number of a given P ( "gcd" means "greatest common divisor"). 

Please note that we define 1 as lucky number of any non-negative integers because 1 has no prime factor.

Input

The first line of input is an integer Q meaning that there are Q test cases. 
Then Q lines follow, each line is a test case and each test case contains three non-negative numbers: n, m and P (n, m, P <= 5×10 5. Q <=5000).

Output

For each test case, print the number of pairs (a, b), which 1<=a<=n , 1<=b<=m, and gcd(a,b) is a lucky number of P.

Sample Input

2
10 10 0
10 10 1

Sample Output

63
93

第一发t：

题目的大体意思就是让你求1<=i<=n和1<=j<=m范围内gcd（i,j）的素因子个数小于p的（i,j）个数对。我们可以通过莫比乌斯反演在O(n)的复杂度计算i和j内gcd为某个值的具体个数。再求莫比乌斯函数的时候我们也可以在数据范围内所有的i值对应的因子个数预处理出来。

我猜测我可以枚举gcd（i，j），判断gcd（i，j）的素因子个数。然后把素因子个数小于p的累加。我还算了一下时间复杂度。是个调和级数求和乘以min（n,m）。大致为n*log（n）的时间复杂度。考虑t后如果数据特别强就会t，果不其然t了。

        for(int i=1;i<=minn;i++)
         {
             int nn=n,mm=m;
             nn/=i;
             mm/=i;
             ll ans=0;
             int minnn=min(nn,mm);
         for(int j=1;j<=minnn;j++)
            ans+=(ll)mu[j]*(nn/j)*(mm/j);
            if(a[i]<=p)
            res+=ans;
        }

然后我又换了一种思路，我可以去计算每个值对应的范围里gcd的因子个数的莫比乌斯值的对应贡献值。最后对这些值求一个前缀和。代表在在这个值范围内，小于给定素因子个数值的全部莫比乌斯贡献值。

    设A(d)：gcd(a, b)=d的有多少种

    设B(j): gcd(a, b)是j的倍数的有多少种，易知B(j) = (n/j)*(m/j)

   A(1) = μ(1)*B(1) + μ(2)*B(2) + μ(3)*B(3) + ... + μ(p1*p2...)*B(p1*p2...)

    A(2) = μ(1)*B(1*2) + μ(2)*B(2*2) + μ(3)*B(3*2) + ... + μ(p1*p2..)*B(p1*p2..*2)

     ...

     A(d) = μ(1)*B(1*d) + μ(2)*B(2*d) + μ(3)*B(3*d) + ... + μ(p1*p2..)*B(p1*p2..*d)

 

     ans = A(1)+A(2)+...+A(d) = F(1)*B(1) + F(2)*B(2) + ... + F(p1*p2..)*B(p1*p2..)

     于是可以枚举公约数i{表示A(i)}，利用筛法找出i的倍数j，i对B(j)的贡献系数为：F(j)+=μ(j/i)

     总之，求出B(j)的总贡献系数F(j)即可得答案：F(1)*B(1)+F(2)*B(2)+...+F(n)*B(n)

T=5*10^3,n最大为5*10^5.还有预处理的过程还是有t的可能。

为了优化我们对最后一步进行除法分块。

求B[x]的过程就是    (n/x)*(m/x);

B[x+1]  为  (n/(x+1))*(m/(x+1));

我们可以发现有很多  n/x  ==n/(x+1)的情况，这时我们.顺理当然想到除法分块。

代码如下：

#include<bits/stdc++.h>
typedef long long ll;
const int maxn = 5*1e5+5;
using namespace std;
int isprime[maxn],mu[maxn],prime[maxn],n,m,p,cnt,a[maxn],f[maxn][20];
void mobius(){
  cnt=0;
  mu[1]=1;
  a[1]=0;
  for(int i=2;i<maxn;i++)
  {
      if(!isprime[i]){
        prime[cnt++]=i;
        a[i]=1;
        mu[i]=-1;
      }
      for(int j=0;j<cnt&&i*prime[j]<maxn;j++)
      {
          isprime[i*prime[j]]=1;
          a[i*prime[j]]=a[i]+1;
          if(i%prime[j])
             mu[i*prime[j]]=-mu[i];
  
          else
           {
            mu[i*prime[j]]=0;
            break;
           }

      }
  }
     for(int i = 1; i < maxn; i++){
        for(int j = i; j < maxn; j += i){
            f[j][a[i]] += mu[j / i];
        }
    }
    for(int i = 1; i < maxn; i++){
        for(int j = 1; j < 20; j++){
            f[i][j] +=  f[i][j - 1] ;
        }
    }
    for(int i = 1; i < maxn; i++){
        for(int j = 0; j < 20; j++){
            f[i][j] += f[i - 1][j];
        }
    }  
}
int main()
{
    ll res,j;
    int t;
    mobius();
     scanf("%d",&t);
      while(t--)
      {
        scanf("%d %d %d",&n,&m,&p);
        res=0;
        int minn = min(n,m);
        if(p>19)
        p=19;
        for(int i=1;i<=minn;i=j+1)
        {
           j = min(n/(n/i),m/(m/i));
              res+=(n/j)*1ll*(m/j)*(f[j][p]-f[i-1][p]);
        } 
         printf("%lld\n",res);
     }
    return 0;
}