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

最新推荐文章于 2019-11-06 20:09:08 发布

Here_jiaxinwei

最新推荐文章于 2019-11-06 20:09:08 发布

阅读量169

点赞数

分类专栏：【数论-------------------------】文章标签： hdu

本文链接：https://blog.youkuaiyun.com/baidu_35643793/article/details/76704601

版权

【数论-------------------------】专栏收录该内容

18 篇文章

订阅专栏

本文介绍了 TrickGCD 问题的两种解法，利用莫比乌斯反演求解 gcd 相关组合计数问题。通过对给定数组 A 构造满足特定条件的数组 B，探讨了如何计算满足 gcd 大于等于 2 的序列数量。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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

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

方法二：

$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 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.youkuaiyun.com/qq_32570675/article/details/76212197、http://blog.youkuaiyun.com/wing_wuchen/article/details/76298196*******************/