1. 程式人生 > >HDU 6053 TrickGCD 【容斥定理】【莫比烏斯函式】

HDU 6053 TrickGCD 【容斥定理】【莫比烏斯函式】

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 answer mod 109+7

Sample Input

1
4
4 4 4 4

Sample Output

Case #1: 17

題意:

給定數列a,求b數列的不同組成種數,對於b數列的任意子集的gcd≥2。

想法:

從2遍歷到amin的數作為gcd,看每個ai中能拿出幾個gcd的倍數來。
舉個栗子,4 5 6 8
x=2時,4和5中都能拿出2,4;6中能拿出2,4,6;8中能拿出2,4,6,8,所以sum1=2*2*3*4=48;
x=3時,4和5中都能拿出3;6和8中都能拿出3,6,所以sum2=2*2=4;
x=4時,4,5,6中都能拿出4;8中能拿出4,8,所以sum3=1*2=2;
x>=4時, 4中已經拿不出數了。
所以sum=sum1+sum2+sum3=54
但會有重複的情況。
比如x=2時,可以取出4 4 4 4,x=4時,也可以取出4 4 4 4。

由此我們發現,我們得出公式sum=aminx=2ajj=1
為了處理重複的情況,我們採用專門用來處理重複情況的容斥定理。

我們知道任何≥2的數都可以分解成質因數相乘的情況。
而每一個倍數都只能統計存在一次。
假設我們討論n, n=p*a,p是最小素因子
這裡寫圖片描述
①代表此倍數在之前p的次數中已經出現過了,比如n=4=2*2;
②代表每個素數都只會本身出現一次,之前未出現過,比如n=6=2*3;
③比如n=12=2*6,再將6分解成2*3帶回最初討論。

即根據容斥思想,考慮一個數n,假設他能分解成三個質因數相乘:a*b*c,那麼在算gcd為a、b、c時分別算了一遍,所以我們要減去gcd為a*b、a*c、b*c的,再加上gcd為a*b*c的。
推廣:我們需要加上質因數為奇數個的,減去質因數為偶數個的。當某個質因數的指數大於等於2,那麼這個數我們之前肯定算過,所以這個不計入最後答案。
而取反後的莫比烏斯函式正好符合我們的需要,
所以最終公式為sum=aminx=2-μ(x)aii=1aix

然鵝數很大,這樣直接做會超時,所以下面我們來想一想優化的方式。
之前我們是在求出每個數能取到此時的gcd倍數有幾個,現在我們把aii相同的ai放到一個組來看,求遍歷i時每一個區域內能取到的個數,然後用快速冪降低時間複雜度,再用莫比烏斯函式去重。
優化後的式子為:sum=amini=2-μ(i)ji100000i=1jnum[ji+i1]num[ji1]

比如前面栗子中的4 5 6 8
x=2時,取出一個2的範圍內(2,3),不存在;取出兩個2的範圍內(4,5),有4 5 存在;取出三個2的範圍內(6,7),6存在;取出四個2的範圍(8,9)內,有8存在。所以sum1=102231*41=48。
接著再累加求和去重~和優化前一個思路辣~

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
using namespace std;
const int maxn=2e6+10;
const int mod=1e9+7;
typedef long long ll;
bool vis[maxn];
ll prime[maxn],mu[maxn],num[maxn],a[maxn];
void mobius()//mu[i]表示莫比烏斯函式
{
    memset(vis,false,sizeof vis);
    mu[1]=1;
    int tot=0;
    for(ll i=2;i<maxn;++i)
    {
        if(!vis[i])
        {
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot;++j)
        {
            if(i*prime[j]>maxn)
                break;
            vis[i*prime[j]]=true;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else
            {
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
}
ll qmod(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b%2==1)
            ans=ans*a%mod;
        a=(a*a)%mod;
        b/=2;
    }
    ans=ans%mod;
    return ans;
}
int main()
{
    mobius();
    int t,n;
    scanf("%d",&t);
    for(int cnt=1;cnt<=t;++cnt)
    {
        scanf("%d",&n);
        ll minn=maxn;
        memset(num,0,sizeof num);
        for(int i=1;i<=n;++i)
        {
            scanf("%lld",&a[i]);
            num[a[i]]++;
            minn=min(a[i],minn);
        }
        for(int i=1;i<maxn;++i)
        {//num[i]表示≤i的在a數列中的個數
            num[i]+=num[i-1];
        }
        ll ans=0;
        for(ll i=2;i<=minn;++i)
        {
            ll sum=1;
            for(ll j=1;j*i<=100000;++j)
            {//j表示倍數
                sum=(sum*qmod(j,num[(j+1)*i-1]-num[j*i-1])%mod)%mod;
            }
            ans=(ans-sum*mu[i]%mod+mod)%mod;
        }
        printf("Case #%d: %lld\n",cnt,ans);
    }
    return 0;
}

ps:沃德天 這個題 說多了都是淚/(ㄒoㄒ)/~~