spoj 4491. Primes in GCD Table 莫比乌斯反演

http://www.spoj.com/problems/PGCD/

题意： 给出 a ，b  （1<=a,b <=10^7） ,求 gcd（x，y） =prime 的 方案数，   其中 1=<x<=a ,  1<=x<=b 。

如果gcd(x ,y) = 1; 这个问题比较好求。

g（n） 为  gcd(x ,y) = n 的方案数。

f( n )  为  gcd(x ,y) = b ,其中 n|b 的所有方案数

f( a ) =  sigma( g( b ))   条件  a|b；

g( n ) = sigma( u(b)* f(b) )   n|b

u(b) =1   ,   b  是 Square-free Number 且 b是偶数个素因子的积

u(b) = 0      b  是 Square-free Number 且 b是寄数个素因子的积

 u(b) = -1   b  是 不是  Square-free Number

f( n ) = a/n * b/n

g( n ) =  sigma(  u(t/n )*a/t* b/t);

ans  =  sigma( g( d )) ,d 是素数

ans = sigma{  u[ j ]* a/( d* j)* b/( d *j ) }   d 是素数。

对于  x=p1^k1 * p2^k2 * p3^k3 ... 其中pi为素数。

当所有 ki全为1。那么从里面拿出一个素数d，会剩下一个Square-Free-Number。

1设x有t个素因子。取出一个后（其实就是   对于（a，b） 中， gcd（n）=  取值范围(a / n, b/ n) 中的gcd( 1)  ）而取法有t种于是ans+= ( ((t-1)&1)?-1:1 ) * t * (m/x) * (n/x)

2 当ki有一个是2，其它全是1 时， 所取的值只能是 ki等于2 的那个数，其余的都不会构成 square-free Number 。 所以取法仅一种。

那么显然 ans+=( (t&1)?-1:1 )*(m/x)*(n/x)

3 对于有多个2 的 x  不用考虑， 无论取那个都不是 square-free Number。

解法 对 u 进行打表， 求解过程中 分块求求值。

打表过程比较慢  程序 45s  。。  不解，本地跑得挺快的。。

#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <stack>
#include <cstring>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
 
#include <cstdlib>
#include <ctime>
#include <assert.h>
#include <queue>
#define REP(i,n) for(int i=0;i<n;i++)
#define TR(i,x) for(typeof(x.begin()) i=x.begin();i!=x.end();i++)
#define ALLL(x) x.begin(),x.end()
#define SORT(x) sort(ALLL(x))
#define CLEAR(x) memset(x,0,sizeof(x))
#define FILLL(x,c) memset(x,c,sizeof(x))
using namespace std;
 
const double EPS = 1e-8;
 
#define LL long long 
#define pb push_back
const int maxn = 10000001;
int f[maxn];
int sum[maxn];
int two[maxn];
void init(){
    for(int i =1 ;i<maxn;i++){
        f[i] = 0;
        two[i]=0;
    }
    for(int i=2;i<maxn;i++){
        if(!f[i]){
             for(int j = i;j<maxn;j+=i){
                    f[j] = i;
             }            
        }
    }
    for(int i=2;i<maxn;i++){
            if(f[i] ==0){
                f[i] =1;
            }else{
                int j = i/f[i];
                if(j%f[i] ==0){
                    two[i] = two[j]+1;
                }else{
                    two[i] = two[j];
                }
                f[i] = f[j]+1;
            }    
    }
    for(int i=2;i<=maxn;i++){
        if(two[i] == 1){
                if(f[i]&1){
                    f[i] =1;
                }else{
                    f[i] = -1;
                }
        }else if(two[i] == 0){
            
            if((f[i])&1){
            
            }else{
                f[i] = -f[i];    
            }
        }else{
            f[i] = 0;
        }
    }
    sum[0]= 0;
    sum[1] = 0;
    for(int i=2;i<maxn;i++){
        sum[i] = sum[i-1] + f[i];
    }
}
int a,b, c;
LL ans;
void gcd(int a,int b){
    int k = min(a,b);
    for(int i = k , j ; i > 1 ; i = j ){
        LL ta = a/i, tb = b/i;
        j =max(a/ (ta+1) , b /(tb+1));
        ans += (LL)ta*tb*(sum[i] - sum[j]);
        // cout << ta<<" "<<tb<< " "<<sum[i] - sum[j]<<endl;
  }
}
 
int main(){
    init();
    int t ;
    cin>>t;
    while(t--){
    
        scanf("%d%d",&a,&b);
        ans = 0;
        gcd(a,b);
        printf("%lld\n",ans);
    }
    return 0;
}