您的位置:首页 > 其它

hdu6053 TrickGCD

2017-07-29 21:19 225 查看
莫比乌斯反演

#include<bits/stdc++.h>
using namespace std;
#define N 100010
#define LL long long int
const LL MOD=1e9+7;
int vis
, mu
, cnt, prime
;
LL CNT
, f
;
//CNT[i]表示<=i的数有几个,f[i]表示GCD为i的倍数的组数有多少个
void Init()
{
cnt=0;
memset(vis,0,sizeof(vis));
mu[1] = 1;
cnt = 0;
for(int i=2; i<N; i++)
{
if(!vis[i])
{
prime[cnt++] = i;
mu[i] = -1;
}
for(int j=0; j<cnt&&i*prime[j]<N; j++)
{
vis[i*prime[j]] = 1;
if(i%prime[j]) mu[i*prime[j]] = -mu[i];
else
{
mu[i*prime[j]] = 0;
break;
}
}
}
}
LL quick_pow(LL  a, LL b)
{
LL ans=1;
LL base =a;
while(b!=0)
{
if(b&1)
ans=(ans*base)%MOD;
base=(base*base)%MOD;
b>>=1;
}
return ans%MOD;
}
void cal(LL m)
{
/*i相当于GCD,GCD最小为2(题目条件),最大为输入里的最小值(显然),枚举i*/
for(int i=2; i<=m; i++)
{
f[i]=1;
/* 对于某一个i来说,不同的数可能是i不同的倍数。比如i是2,那么对4来说,4/2=2, 对5来说,5/2=2, 对6来说,6/2=3。因为对2来说,被4除和被5除是一样。*/
/*所以可以简化,将一段与i相除后值相同的乘起来*/
for(int j=2*i-1, k=1; j<=100000; j+=i, k++)
{
f[i]=(f[i]*quick_pow(k, CNT[j]-CNT[j-i])%MOD)%MOD;
/*下面的步骤是因为防止j<1000000但是j+k又>100000时,有一部分数据没有被算到。*/
if(j+i>100000&&j!=100000)
f[i]=(f[i]*quick_pow(k+1, CNT[100000]-CNT[j])%MOD)%MOD;
}
//  printf("f[%d]=%d\n", i, f[i]);
}
}
int main()
{
Init();
int T, n, a, cas=1;
LL Max, Min;
scanf("%d", &T);
while(T--)
{
Max=-1;
Min=10000000;
memset(CNT, 0, sizeof(CNT));
memset(f, 0, sizeof(f));
scanf("%d", &n);
for(int i=0; i<n; i++)
{
scanf("%d", &a);
CNT[a]++;
if(a>Max)
Max=a;
if(a<Min)
Min=a;
}
for(int i=1; i<=1e5; i++)
{
CNT[i]+=CNT[i-1];
//printf("CNT[%d]=%d\n", i, CNT[i]);
}
//   printf("Max=%d, Mi=%d\n", Max, Min);
cal(Min);
LL ans=0, temp;
/*对于每个数i,算出GCD为i时的组数(而不是为i的倍数!),然后相加,j表示的其实是倍数*/
for(int i=2; i<=100000; i++)
{
temp=0;
for(int j=1; i*j<=100000 ;j++)
{
temp+=mu[j]*f[i*j];
temp%=MOD;
temp=(temp+MOD)%MOD;
}
ans+=temp;
ans%=MOD;
ans=(ans+MOD)%MOD;
}
printf("Case #%d: %I64d\n",cas++, ans);
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: