HDU 1695

看见别人的用的莫比乌斯来做，我看了好久也没明白，实在佩服，看到是组合数学的内容，只好先留着，待我学了组合数学后再用莫比乌斯来写。

 

求GCD(X,Y)=K.其实即是在[1,X/K]和[1,Y/K]的区间内求GCD(X,Y)=1的对数。这样，假设X/K<=Y/K，在[1,X/K]的区间内可以用欧拉函数求出对数，这是显然的。

而在[X/K+1,Y/K]的区间内，则可以用容斥原理求出不与区间内的数互质的对数，再后来减去即可。

怎么样去求区间内不互质的数呢，设一个P在[X/K+1,Y/K]，分解质因数为P1*P2*P3.....*PS，然后用对[1,X/K]应用容斥原理即可。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define LL __int64
using namespace std;

const int Max=1000005;

bool isprime[Max];
int prime[Max],np;
int phi[Max];
int num[Max],no;
LL all,ans;

void exchange(int &a,int &b){
    if(a>b){
        int tmp=a;
        a=b; b=tmp;
    }
}

void do_prime(){
    np=0;
    memset(isprime,true, sizeof(isprime));
    isprime[1]=false;
    for(LL i=2;i<Max;i++){
        if(isprime[i]){
            prime[np++]=i;
            for(LL j=i*i;j<Max;j+=i){
                isprime[j]=false;
            }
        }
    }
}

void do_phi(){
    for(int i=1;i<Max;i++) phi[i]=i;
    for(int i=2;i<Max;i+=2) phi[i]/=2;
    for(int i=3;i<Max;i+=2){
        if(phi[i]==i){
            for(int j=i;j<Max;j+=i)
            phi[j]=phi[j]/i*(i-1);
        }
    }
}

void initial(){
    do_prime();
    do_phi();
}


void dfs(int i,int nu,int x,int mu,int b){
    if(nu==x){
        all+=(LL)(b/mu); return;
    }
    if(i==no) return ;
    dfs(i+1,nu+1,x,mu*num[i],b);
    dfs(i+1,nu,x,mu,b);
}

int Rong(int b,int d){
    int s=0;
    for(int i=1;i<=no;i++){
        all=0;
        dfs(0,0,i,1,b);
        if(i&1) s+=(int)all;
        else s-=(int)all;
    }
    return b-s;
}

int main(){
    int T,a,b,c,d,k;
    scanf("%d",&T);
    initial();
    int kase=0;
    while(T--){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        if(k>b||k>d||k==0){
            printf("Case %d: 0\n",++kase); 
            continue;
        }
        exchange(b,d);
        b/=k; d/=k;
        all=ans=0;
        for(int i=1;i<=b;i++)
        ans+=(LL)phi[i];
        for(int i=b+1;i<=d;i++){
            no=0;
            int tmp=i;
            for(int j=0;(LL)prime[j]*(LL)prime[j]<=(LL)tmp&&j<np;j++){
                if(tmp%prime[j]==0){
                    num[no++]=prime[j];
                    while(tmp%prime[j]==0){
                        tmp/=prime[j];
                    }
                }
            }
            if(tmp>1)
            num[no++]=tmp;
            ans+=(LL)Rong(b,d);
        }
        printf("Case %d: %I64d\n",++kase,ans);
    }
    return 0;
}