4491. Primes in GCD Table

Problem code: PGCD

Johnny has created a table which encodes the results of some operation -- a function of two arguments. But instead of a boring multiplication table of the sort you learn by heart at prep-school, he has created a GCD (greatest common divisor) table! So he now has a table (of height a and width b), indexed from (1,1) to (a,b), and with the value of field (i,j) equal to gcd(i,j). He wants to know how many times he has used prime numbers when writing the table.

Input

First, t ≤ 10, the number of test cases. Each test case consists of two integers, 1 ≤ a,b < 10⁷.

Output

For each test case write one number - the number of prime numbers Johnny wrote in that test case.

Example

Input:

2 10 10 100 100

Output:

30 2791

题解参考：

ans = sigma(p, sigma(d, μ(d) * (n/pd) * (m/pd)))Let s = pd, thenans = sigma(s, sigma(p, μ(s/p) * (n/s) * (m/s)))    = sigma(s, (n/s) * (m/s) * sigma(p, μ(s/p)))Let g(x) = sigma(p, μ(x/p)), thenans = sigma(s, (n/s) * (m/s) * g(s))

如果我们能预处理g(x)的话就能和前面一样分块搞了。这个时候我们多么希望g(x)和μ(x)一样是积性函数。看完题解后，发现有一个不是积性函数，胜似积性函数的性质。由于题解没有给证明，所以就意淫了一个证明。

考虑质数p'，g(p'x) = sigma(p | p'x, μ(p'x/p))。

当x % p' == 0，那么考虑sigma中的变量p的所有取值，它和g(x)中p是相同的。而μ(x)这个函数，如果x有平方因子的话就等于0，因此当p != p'时μ(p'x/p) = 0，当p == p'是，μ(p'x/p) = μ(x)。所以g(p'x) = μ(x)。

当x % p' != 0，同样考虑p，会发现它的取值只比g(x)中的p多出一个p'。同理按照p是否等于p'讨论，可以得到g(p'x) = -g(x) + μ(x)。

因此g(x)可以在线性筛素数的时候算出。剩下的就是前缀和、分块了。

1 /* *********************************************** 2 Author        :kuangbin 3 Created Time  :2013-10-19 22:01:05 4 File Name     :E:\2013ACM\专题学习\数学\莫比乌斯反演\SPOJ_PGCD.cpp 5 ************************************************ */ 6  7 #include 
     
       8 #include 
      
        9 #include 
       
        10 #include 
        
         11 #include 
         
          12 #include 
          
           13 #include 
           
            14 #include 
            15 #include 
             
              16 #include 
              
               17 #include 
               
                18 #include 
                
                 19 using namespace std;20 21 const int MAXN = 10000000;22 bool check[MAXN+10];23 int prime[MAXN+10];24 int mu[MAXN+10];25 int g[MAXN+10];26 int sum[MAXN+10];27 void Moblus()28 {29 memset(check,false,sizeof(check));30 mu[1] = 1;31 int tot = 0;32 for(int i = 2; i <= MAXN; i++)33 {34 if(!check[i])35 {36 prime[tot++] = i;37 mu[i] = -1;38 g[i] = 1;39 }40 for(int j = 0;j < tot;j++)41 {42 if(i * prime[j] > MAXN)break;43 check[i*prime[j]] = true;44 if(i % prime[j] == 0)45 {46 mu[i * prime[j]] = 0;47 g[i * prime[j]] = mu[i];48 break;49 }50 else51 {52 mu[i * prime[j]] = -mu[i];53 g[i * prime[j]] = -g[i] + mu[i];54 }55 }56 }57 sum[0] = 0;58 for(int i = 1;i <= MAXN;i++)59 sum[i] = sum[i-1] + g[i];60 }61 int main()62 {63 //freopen("in.txt","r",stdin);64 //freopen("out.txt","w",stdout);65 Moblus();66 int T;67 int n,m;68 scanf("%d",&T);69 while(T--)70 {71 scanf("%d%d",&n,&m);72 if(n > m)swap(n,m);73 long long ans = 0;74 int last = 0;75 for(int i = 1;i <= n;i = last+1)76 {77 last = min(n/(n/i),m/(m/i));78 ans += (long long)(sum[last] - sum[i-1])*(n/i)*(m/i);79 }80 printf("%lld\n",ans);81 }82 return 0;83 }

转载地址：http://sfdzl.baihongyu.com/

你可能感兴趣的文章