SPOJ PGCD 4491. Primes in GCD Table && BZOJ 2820 YY的GCD （莫比乌斯反演）

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

题解参考：

http://quartergeek.com/eight-gcd-problems/

ans = sigma(p, sigma(d, μ(d) * (n/pd) * (m/pd)))

Let s = pd, then

ans = 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)), then

ans = 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 <stdio.h>
 8 #include <string.h>
 9 #include <iostream>
10 #include <algorithm>
11 #include <vector>
12 #include <queue>
13 #include <set>
14 #include <map>
15 #include <string>
16 #include <math.h>
17 #include <stdlib.h>
18 #include <time.h>
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             else
51             {
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 }