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

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

  
  

   
   1
4
4 4 4 4
  
  

Sample Output

  
  

   
   Case #1: 17
  
  

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

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

                               :gcd(i,j)=d的倍数的数对的个数；

                             ：gcd(i,j)=d的数对的个数；

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

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

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

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

那么这样就好办了：

方法一：

发现系数满足莫比乌斯函数；

其中把转化成了;

我来解释一下这个式子：就是把枚举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 long

const 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;
}

方法二：

其中，

           

             

关于m的解释同方法一；

代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;

#define mod 1000000007
#define ll long long

const 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*******************/