HDU 6390 GuGuFishtion 【gcd+欧拉函数+mobius函数】

题目链接

题意:已知

Gu(a,b)=ϕ(ab)ϕ(a)ϕ(b) G u ( a , b ) = ϕ ( a b ) ϕ ( a ) ϕ ( b )

(a=1mb=1nGu(a,b))(modp) ( ∑ a = 1 m ∑ b = 1 n G u ( a , b ) ) ( mod p )

先推公式……必不可少的一步

由欧拉函数的定义,

ϕ(x)=xi=1n(11pi) ϕ ( x ) = x ∏ i = 1 n ( 1 − 1 p i )

Gu(a,b) G u ( a , b ) 展开,发现分子是a、b和ab的质因数,分母是a、b、a的质因数和b的质因数。而ab的质因数在a的质因数中有出现,在b的质因数中也有出现,也就是说 gcd(a,b) g c d ( a , b ) 这部分在分母多出现了一次。约去其他部分后可以得到:

Gu(a,b)=gcd(a,b)ϕ(gcd(a,b)) G u ( a , b ) = g c d ( a , b ) ϕ ( g c d ( a , b ) )

a[i]=iϕ(i) a [ i ] = i ϕ ( i ) ,所求式

(a=1mb=1nGu(a,b))(modp) ( ∑ a = 1 m ∑ b = 1 n G u ( a , b ) ) ( mod p )

=i=1mj=1na[gcd(i,j)] = ∑ i = 1 m ∑ j = 1 n a [ g c d ( i , j ) ]

=k=1min(m,n)i=1mj=1na[k][gcd(i,j)==k] = ∑ k = 1 m i n ( m , n ) ∑ i = 1 m ∑ j = 1 n a [ k ] ∗ [ g c d ( i , j ) == k ]

=k=1min(m,n)a[k]i=1mj=1n[gcd(i,j)==k] = ∑ k = 1 m i n ( m , n ) a [ k ] ∑ i = 1 m ∑ j = 1 n [ g c d ( i , j ) == k ]

后面这部分可以用mobius函数求出……不会的话可以看看HDU 1695 题解

然后这题还卡常数……能不开long long的尽量用int,并且线性预处理出逆元就可以过了

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1000000+10;

bool check[maxn];
int prime[maxn],mu[maxn];

void Mobius(){
    memset(check,false,sizeof(check));
    mu[1]=1;
    int tot=0;
    for (int i=2;i<=1000000;i++){
        if (!check[i]){
            prime[tot++]=i;
            mu[i]=-1;
        }
        for (int j=0;j<tot;j++){
            if (i*prime[j]>1000000) break;
            check[i*prime[j]]=true;
            if (i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            else{
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
}

int euler[maxn];

void getEuler(){
    memset(euler,0,sizeof(euler));
    euler[1]=1;
    for (int i=2;i<=1000000;i++){
        if (!euler[i]){
            for (int j=i;j<=1000000;j+=i){
                if (!euler[j])  euler[j]=j;
                euler[j]=euler[j]/i*(i-1);
            }
        }
    }
}

ll mod;

ll solve(int a,int b,ll k){
    ll ans=0;
    if (a>b)    swap(a,b);
    a/=k,b/=k;
    for (int i=1;i<=a;i++){
        ans+=1ll*mu[i]*(a/i)*(b/i);
        ans%=mod;
    }
    ll tmp=0;
//  for (int i=1;i<=a;i++){
//      tmp+=1ll*mu[i]*(a/i)*(a/i);
//  }
    return ans-tmp/2;
}

ll qpow(ll x,ll n,ll mod){
    ll ans=1;
    while (n){
        if (n&1)    ans=ans*x%mod;
        n>>=1;
        x=x*x%mod;
    }
    return ans%mod;
}

int T;
int m,n;
ll inv[maxn];
ll a[maxn];

void init(){
    inv[0]=inv[1]=1;
    for (int i=2;i<=min(m,n);i++){
        inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
    }
    for (int i=1;i<=min(m,n);i++){
        a[i]=1ll*i*inv[euler[i]]%mod;
    }
}

int main(){
#ifdef __APPLE__
    freopen("1.in","r",stdin);
#endif
    Mobius();
    getEuler();
    scanf("%d",&T);
    while (T--){
        ll ans=0;
        scanf("%d %d %lld",&m,&n,&mod);
        init();
        for (int k=1;k<=min(m,n);k++){
            ans+=1ll*a[k]*solve(m,n,k);
            ans%=mod;
        }
        printf("%lld\n",ans%mod);
    }
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值