HDU4746 Mophues (Mobius Inversion)

HDU4746 Mophues (Mobius Inversion)
MophuesTime Limit: 2000/1000 MS (Java/Others) Memory Limit: 327670/327670 K (Java/Others)
Total Submission (s): 647 Accepted Submission (s): 263


Problem DescriptionAs 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.
InputThe 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 × 105. Q <= 5000 ).
OutputFor 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

210 10 010 10 1

Sample Output

6393

Example: 5000 groups. Ask you how many logarithm GCD prime factor numbers in [1, n] and [1, m] are less than p.
Idea: first consider a relatively simple version: [1, A] and [1, B]. The number of pairs satisfies GCD <= d. First, define two functions: a (, b, d) indicates the logarithm of GCD (a, B) = d, B (a, B, d) indicates GCD (a, B) is the logarithm of a multiple of d. B (a, B, d) = (A/d) * (B/d) according to the principle of rejection: a (a, B, d) = (1) * B (a, B, d) + (-1) * B (a, B, 2 * d) + (-1) * B (a, d, b, 3 * d) + (0) * B (a, B, 4 * d) + (-1) * B (a, B, 5 * d) + (1) * B (a, B, 6 * d )........... the coefficient before B (a, B, I) is the value of the Mobius function. Then the formula can be written as: A (a, B, 1) = u (1) * B (a, B, 1) + u (2) * B (a, B, 2) + u (3) * B (a, B, 3) + u (4) * B (a, B, 4) + u (5) * B (, b, 5) + u (6) * B (a, B, 6 )........... A (a, B, 2) = u (1) * B (a, B, 2) + u (2) * B (a, B, 4) + u (3) * B (a, B, 6) + u (4) * B (a, B, 8) + u (5) * B (a, B, 10) + u (6) * B (a, B, 12 )........... A (a, B, 3) = u (1) * B (a, B, 3) + u (2) * B (a, B, 6) + u (3) * B (a, B, 9) + u (4) * B (a, B, 12) + u (5) * B (a, B, 15) + u (6) * B (a, B, 18 )...........
A (a, B, 4) = u (1) * B (a, B, 4) + u (2) * B (a, B, 8) + u (3) * B (a, B, 12) + u (4) * B (a, B, 16) + u (5) * B (a, B, 20) + u (6) * B (a, B, 24 )...........
The answer is A (a, B, 1) + A (a, B, 2) + A (a, B, 3) + ...... A (a, B, d) = u (1) * B (a, B, 1) + (u (1) + u (2) * B (, b, 2) + ....... (u (1) + u (2) + u (3) + u (6) * B (a, B, 6 )........ it can be seen that the coefficient before A (a, B, d) is sigma (u (I) (I is the approximate number of d) = C (a, B, d)
Then, there is another constraint for this question: to make the number of prime factors less than or equal to p, then we define this function D (a, B, d, p) it indicates the coefficient before B (a, B, d), then we only need to select some coefficients that meet the conditions from C (a, B, d. Use an array F [d] [cnt] (cnt is the number of prime factors) to record. The array represents the size of the influencing factor where the number of prime factors of d is cnt. Calculate a single prefix and (useful next ). Next, we find that for a d, B (a, B, d) = (B, a, B, d + x) is satisfied ), in addition, if x = min (a/d), B/(B/d), the computation of the entire formula is segmented, therefore, you can use the prefix and when calculating this interval.

#include 
 
  #include 
  
   #include 
   
    #include 
    
     #include 
     
      #include #include 
      
       using namespace std;typedef long long ll;const int maxn = 5e5+10;int mobi[maxn];int preMobi[maxn][20];int pricnt[maxn];bool isPrime[maxn];vector
       
         prime;void getMobi(){ memset(isPrime,1,sizeof isPrime); memset(pricnt,0,sizeof pricnt); mobi[1] = 1; for(int i = 2; i < maxn; i++){ if(isPrime[i]){ mobi[i] = -1; pricnt[i] = 1; prime.push_back(i); } for(int j = 0; j < prime.size() && i*prime[j] < maxn; j++){ pricnt[i*prime[j]] = pricnt[i]+1; isPrime[i*prime[j]] = false; if(i%prime[j]==0){ mobi[i*prime[j]] = 0; break; }else{ mobi[i*prime[j]] = -mobi[i]; } } }}void getpreMobi(){ memset(preMobi,0,sizeof preMobi); for(int i = 1; i < maxn; i++){ for(int j = i; j < maxn; j += i){ preMobi[j][pricnt[i]] += mobi[j/i]; } } for(int i = 1; i < maxn; i++){ for(int j = 0; j < 20; j++){ preMobi[i][j] += preMobi[i-1][j]; } } for(int i = 0; i < maxn; i++){ for(int j = 1; j < 20; j++){ preMobi[i][j] += preMobi[i][j-1]; } }}int n,m,p;ll ans;void solve(){ ans = 0; for(int i = 1; i <= n; i++){ int ed = min(n/(n/i),m/(m/i)); ans += ll(preMobi[ed][p]-preMobi[i-1][p])*(n/i)*(m/i); i = ed; } cout<> ncase; while(ncase--){ scanf("%d%d%d",&n,&m,&p); if(p>=19){ cout<<(ll)n*m<
        
          m) swap(n,m); solve(); } return 0;}