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_backconst 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;}