HDU 4675 (莫比乌斯反演)

HDU 4675 (莫比乌斯反演)

题目链接:点击这里

题意:给定n个数,每个数不超过m,和k.问从$d$从$1$到$m$能够构造多少个$b$数组使得$b$数组的$gcd$为$d$并且$b$数组每个数都不超过$m$,并且恰好有$k$个数和下表对应的a数组的数不相等.

用$F(x)$表示gcd为x的倍数的方案数,$f(x)$表示gcd为x的方案数.先考虑$F(d)$怎么计算.可以把a数组中的数分成两类,第一类是必须对应下标不等的($a[i]$不是$d$的倍数),其他的就是第二类,假设第二类的数量是$p$,第一类的数量就是$n-p$,因为要选择$k$个不同的,第一类必须不同,所以需要在第二类中选择$k-n+p$个,而$b$数组中每一个数都有$[m/p]$种选择.所以最终的结果就是 $$F(d)=C_p^{k-n+p}*([m/d]-1)^{k-n+p}*([m/d]^{n-p})$$ .

然后暴力计算每一个$f(d)$就好了.总的复杂度是$n(lgn)$.

    #include <cstdio>
    #include <cstring>
    #include <cmath>
    #include <algorithm>
    using namespace std;
    #define maxn 311111
    #define mod 1000000007

    int n, m, k;

    int prime[maxn], prime_cnt, mu[maxn];
    bool vis[maxn];
    long long fac[maxn], rev[maxn];

    long long extend_gcd(long long a,long long b,long long &x,long long &y)
    {
        if(a==0&&b==0) return -1;
        if(b==0){x=1;y=0;return a;}
        long long d=extend_gcd(b,a%b,y,x);
        y-=a/b*x;
        return d;
    }
    long long mod_rev(long long a,long long n)
    {
        long long x,y;
        long long d=extend_gcd(a,n,x,y);
        if(d==1) return (x%n+n)%n;
        else return -1;
    }

    void init () {
        prime_cnt = 0;
        mu[1] = 1;
        for (int i = 2; i < maxn; i++) {
            if (!vis[i]) {
                prime[prime_cnt++] = i;
                mu[i] = -1;
            }
            for (int j = 0; j < prime_cnt; j++) {
                if (i*prime[j] >= maxn)
                    break;
                vis[i*prime[j]] = 1;
                if (i%prime[j] == 0) {
                    mu[i*prime[j]] = 0;
                    break;
                }
                else {
                    mu[i*prime[j]] = -mu[i];
                }
            }
        }
        fac[0] = 1;
        rev[0] = 1;
        for (int i = 1; i < maxn; i++) {
            fac[i] = fac[i-1]*i%mod;
            rev[i] = mod_rev (fac[i], mod);
        }
    }

    long long quick_mod(long long a, long long b)
    {
        long long ans = 1;
        a %= mod;
        while(b)
        {
            if(b & 1)
            {
                ans = ans * a % mod;
                b--;
            }
            b >>= 1;
            a = a * a % mod;
        }
        return ans;
    }

    long long C (long long m, long long n) {
        if (n == 0)
            return 1;
        long long ans = fac[m] * rev[n] % mod * rev[m-n] % mod;
        return ans;
    }

    long long F[maxn], f[maxn];
    int a[maxn];
    int cnt[maxn];//对于i,有多少a[j]是i的倍数

    int scan () {
        int num = 0;
        char ch;
        while (1) {
            ch = getchar ();
            while (ch <= '9' && ch >= '0') {
                num = num*10 + ch - '0';
                ch = getchar ();
                if (ch > '9' || ch < '0')
                    return num;
            }
        }
        return num;
    }

    int main () {
        //freopen ("in.txt", "r", stdin);
        //freopen ("out.txt", "w", stdout);
        init ();
        while (scanf ("%d%d%d", &n, &m, &k) == 3) {
            memset (cnt, 0, sizeof cnt);
            for (int i = 1; i <= n; i++) {
                a[i] = scan ();
                cnt[a[i]]++;
            }
            for (int i = 1; i <= m; i++) {
                for (int j = i+i; j <= m; j += i) {
                    cnt[i] += cnt[j];
                }
            }
            for (int i = 1; i <= m; i++) {
                long long p = cnt[i];
                if (k-n+p < 0) {
                    F[i] = 0;
                    continue;
                }
                F[i] = C (p, k-n+p) * quick_mod (m/i-1, k-n+p) % mod * quick_mod (m/i, n-p) % mod;
            }
            memset (f, 0, sizeof f);
            for (int i = 1; i <= m; i++) {
                if (F[i] == 0) {
                    f[i] = 0;
                }
                else for (int j = i; j <= m; j += i) {
                    f[i] += mu[j/i]*F[j];
                    f[i] %= mod;
                }
                printf ("%lld%c", (f[i] + mod) % mod, (i == m ? '\n' : ' '));
            }
        }
        return 0;
    }