HDU 6053 TrickGCD

原题连接：
http://acm.hdu.edu.cn/showproblem.php?pid=6053
（对于不了解，莫比乌斯反演的同学不需要纠结。直接用结论就好了）

莫比乌斯反演公式

若对于 定义域 为正整数，值域为复数子集 的函数f：

定义它的和函数

F(n)=∑d|nf(d)=∑d|nf(nd)

则对于任意正整数n：

f(n)=∑d|nμ(d)F(nd)=∑d|nμ(nd)F(d)

同时 。也有一种无穷形式的反演公式：

F(n)=∑n|df(d)

这里，d取n的倍数。 所以这是一个无穷和式子。有：

f(n)=∑n|dμ(dn)F(d)

反演的形式很多。都是很巧妙的（可以参阅组合数学，有更一般的反演与容斥）。

上面的的莫比乌斯反演的定义是很强的。- -！一切算术函数

————————————————————————————————————————————————

对于这一题。gcd(Bl,Bl+1,...Br)>=2,对于任意的l<=r,l,r∈[1,n]

设：gcd(B)=n的B数列个数为f(n)

设

F(n)=∑n|dnf(d)

由 第二个莫比乌斯反演得：

f(1)=∑d>=1μ(d)F(d)

这里F(n)就是gcd=n的倍数的所有方案数量。

虽然 得到的式子是一个无穷和式子。

但。当d>min(A1,A2,A3,...An)时，F(d)=0

我么可以通过预处理与快速幂快速计算出F函数

对于F(k)那么任意B都应该含有k这个因子。

那么Bi可选择的方案数就是 ⌊Aik⌋

那么乘法原理得：

F(k)=∏i=1n⌊Aik⌋

令m=min(A1,A2,..An)

令h=max(A1,A2,...An)

暴力算所有的F从1到m的时间就是O(nm)

设cnt[u]是A中A[i]<=u的元素个数。

对于A中满足⌊Aik⌋=t的所有的元素个数为：

cnt[(t+1)k−1]−cnt[tk−1]

那么 A中与k的商等于t 的元素对F(k)的贡献为：

tcnt[(t+1)k−1]−cnt[tk−1]

那么计算F(k)的复杂度为O(hklog(cnt[(t+1)k−1]−cnt[tk−1]))

计算出1到m的F(k)所用时间为：

O(h∑mk=1log(cnt)k)=O(h∗logm∗logh)<O(105log2(105))

利用反演公式：

f(1)=∑d>=1μ(d)F(d)

计算出f(1)

answer=F(1)−f(1)

#include <algorithm>#include <string.h>#include <stdio.h>#include <set>#include <math.h>#define MAXN 100005using namespace std;typedef long long LL;const LL mod=1e9+7;int cnt[MAXN];LL f[MAXN];const int N=MAXN;int mu[N];void getMu(){    for(int i=1;i<=N;i++)    {        int target=i==1?1:0;        int delta=target-mu[i];        mu[i]=delta;        for(int j=i+i;j<=N;j+=i)        {            mu[j]+=delta;        }    }}LL pow(LL a,int b){    LL tmp=1;    while(b)    {        if(b&1)            tmp=tmp*a%mod;        a=a*a%mod;        b>>=1;    }    return tmp;}int main (){    int T,ca=1;    getMu();    scanf("%d",&T);    while(T--)    {        int A;        LL ans=0;        memset(cnt,0,sizeof cnt);        int n,m=MAXN,ma=-1;        scanf("%d",&n);        for(int i=1;i<=n;i++)        {            scanf("%d",&A);            m=min(A,m);            ma=max(A,ma);            cnt[A]++;        }        for(int i=1;i<=ma;i++) cnt[i]+=cnt[i-1];        for(int k=1;k<=m;k++)        {            f[k]=1;            for(int l=-1,j=k-1,t=1; ;t++)            {                j+=k;                l+=k;                if(j>ma)j=ma;                f[k]*=pow((LL)t,cnt[j]-cnt[l]);                f[k]%=mod;                if(j==ma)break;            }        }        for(int i=1;i<=m;i++)            ans=(ans+f[i]*mu[i]+mod)%mod;        printf("Case #%d: %lld\n",ca++,(f[1]-ans+mod)%mod);    }    return 0;}