【数论-莫比乌斯】hdu 6053 TrickGCD

来源：互联网发布：jsp中引用js文件编辑：程序博客网时间：2024/06/05 06:43

TrickGCD

Time Limit: 5000/2500 MS (Java/Others) Memory Limit: 262144/262144 K (Java/Others)
Total Submission(s): 2846 Accepted Submission(s): 1071

Problem Description
You are given an array A , and Zhu wants to know there are how many different array B satisfy the following conditions?
* 1≤Bi≤Ai
* For each pair( l , r ) (1≤l≤r≤n) , gcd(bl,bl+1...br)≥2

Input
The first line is an integer T(1≤T≤10) describe the number of test cases.
Each test case begins with an integer number n describe the size of array A.
Then a line contains n numbers describe each element of A
You can assume that 1≤n,Ai≤105

Output
For the kth test case , first output "Case #k: " , then output an integer as answer in a single line . because the answer may be large , so you are only need to output answermod109+7

Sample Input

144 4 4 4

Sample Output

Case #1: 17

题意：已知序列A，请你有多少种方式构造一个序列B，并且使得gcd(B)>=2，其中Bi<=Ai；

思路：利用莫比乌斯反演变形求gcd：

$F(d)=\sum_{x=1}^{\frac{n}{d}}f(d\cdot x)$ :gcd(i,j)=d的倍数的数对的个数；

$f(d)=\sum_{x=1}^{\frac{n}{d}}\mu (x)F(d\cdot x)$ ：gcd(i,j)=d的数对的个数；

之前的问题是给出两个区间，问gcd=d的数对有多少个；

推广之，现在的问题是给出n个区间，问gcd=d的序列有多少个。

本题是给出 Ai，你要选择的 Bi 在（1~Ai）之间，

那不就相当于给出了n个区间：(1~A1)，(1~A2)，(1~A3)...(1~An)，然后找出gcd>=2的序列有多少个；

那么这样就好办了：

方法一：

$ANS=F(2)+F(3)+F(5)-F(6)+...$ 发现系数满足莫比乌斯函数；

$ANS=\sum_{i=2}^{min(A)}-\mu (i)\prod_{j=1}^{n}\frac{Aj}{i}$ $=\sum_{i=2}^{min(A)}-\mu (i)\prod_{j=1}^{\frac{max(A)}{i}}j^{m}$

其中把 $\prod_{j=1}^{n}\frac{Aj}{i}$ 转化成了 $\prod_{j=1}^{\frac{max(A)}{i}}j^{m}$ ;

我来解释一下这个式子：就是把枚举n个 ai 转化成了枚举 ai/gcd 的值，这个点值一定在1和maxa/gcd之间，这个化简大概优化到了logn。然后就是m的含义，m代表n个 ai 中除以gcd为i的数量，这个m我们可以通过求一遍 ai 的前缀和得到（就是开个数组，把 ai 当作下标，记录 ai出现了多少次，然后求这个数组的前缀和，这样就可以O（1）的求出一个范围内 ai 出现了多少次了，也就得到了m的值）;

代码：

#include<iostream>#include<cstdio>#include<cstring>#include<algorithm>using namespace std;#define mod 1000000007#define ll long longconst int N= 202000;bool check[N+10];int prime[N+10];int mu[N+10];int A[N+10],sum[N+10];int minn,maxx;void Moblus(){    memset(check,false,sizeof(check));    mu[1]=1;    int tot=0;    for(int i=2; i<=N; i++)    {        if(!check[i])        {            prime[tot++]=i;            mu[i]=-1;        }        for(int j=0; j<tot; j++)        {            if(i*prime[j]>N) break;            check[i*prime[j]]=true;            if(i%prime[j]==0)            {                mu[i*prime[j]]=0;                break;            }            else            {                mu[i*prime[j]]=-mu[i];            }        }    }}ll mod_pow(ll a,ll b){    ll res=1;    while(b)    {        if(b&1)            res=res*a%mod;        a=a*a%mod;        b>>=1;    }    return res;}ll solve()     //公式{    ll ans=0;    for(int i=2;i<=minn;i++)    {        ll an=1;        for(int j=1;j*i<=maxx;j++)            an=an*mod_pow(j,sum[j*i+i-1]-sum[j*i-1])%mod;        ans=ans-(an*mu[i]);        ans%=mod;    }    return ans;}int main(){    Moblus();    int t,cas=0;    scanf("%d",&t);    while(t--)    {        memset(A,0,sizeof(A));        memset(sum,0,sizeof(sum));        minn=1000000,maxx=-1;        int n,val;        scanf("%d",&n);        for(int i=1;i<=n;i++) //预处理        {            scanf("%d",&val);            A[val]++;            maxx=max(maxx,val);            minn=min(minn,val);        }        for(int i=1;i<=N;i++)            sum[i]=sum[i-1]+A[i];  //                    ll ans=solve();          printf("Case #%d: %lld\n",++cas,(ans+mod)%mod);    }    return 0;}

方法二：

$ANS=F(1)-f(1)$

其中，

$F(1)=\prod_{i=1}^{n}Ai$

$f(1)=\sum_{x=1}^{min(A)}\mu (x)F(1*x)$

$=\sum_{x=1}^{min(A)}\mu (x)\prod_{i=1}^{n}\frac{Ai}{x}$

$=\sum_{x=1}^{min(A)}\mu (x)\prod_{j=1}^\frac{max(A)}{x}j^{m}$

关于m的解释同方法一；

代码：

#include<cstdio>#include<cstring>#include<algorithm>using namespace std;#define mod 1000000007#define ll long longconst int N=300100;int A[N+10],sum[N+10];int maxx=-1,minn=1000000;bool check[N+10];int prime[N+10];int mu[N+10];void Moblus(){    memset(check,false,sizeof(check));    mu[1]=1;    int tot=0;    for(int i=2; i<=N; i++)    {        if(!check[i])        {            prime[tot++]=i;            mu[i]=-1;        }        for(int j=0; j<tot; j++)        {            if(i*prime[j]>N) break;            check[i*prime[j]]=true;            if(i%prime[j]==0)            {                mu[i*prime[j]]=0;                break;            }            else            {                mu[i*prime[j]]=-mu[i];            }        }    }}inline ll mod_pow(ll x,ll n){    ll res=1;    while(n>0)    {        if(n&1)            res=res*x%mod;        x=x*x%mod;        n>>=1;    }    return res;}ll solve()       //f(1){    ll ans=0;    for(int i=1; i<=minn; i++)    {        ll an=1;        for(int j=1; j*i<=maxx; j++)            an=an*mod_pow(j,(sum[j*i+i-1]-sum[j*i-1]))%mod;        ans=(ans+(mu[i]*an))%mod;    }    return ans;}int main(){    Moblus();    int t,cas=0;    scanf("%d",&t);    while(t--)    {        memset(A,0,sizeof(A));        memset(sum,0,sizeof(sum));        maxx=-1,minn=1000000;        int n;        scanf("%d",&n);        for(int i=1; i<=n; i++)        {            int val;            scanf("%d",&val);            A[val]++;            maxx=max(maxx,val);            minn=min(minn,val);        }        for(int i=1; i<=N; i++)            sum[i]=sum[i-1]+A[i];        ll F1=1,f1=0;         for(int i=minn; i<=maxx; i++)    //F1        {            F1=F1*mod_pow(i,A[i]);            F1%=mod;        }         f1=solve();                     //f1        f1%=mod;        ll ans=((F1-f1)+mod)%mod;        printf("Case #%d: %lld\n",++cas,ans);    }    return 0;}

/****************参考blog：http://blog.csdn.net/qq_32570675/article/details/76212197、http://blog.csdn.net/wing_wuchen/article/details/76298196*******************/

阅读全文

0 0