HDU 5321 Beautiful Set 美丽集合

题意：给定一个集合，含有n个数。浙理工先生和杭电先生各自有计算这个集合美丽值的方法。

浙理工先生的计算方法是：对于这个n个数的某个排列，此排列的美丽值为这个排列所有的区间最大公约数之和。然后这个集合的美丽值为n个数的所有排列的美丽值之和。

杭电先生的计算方法是：在这个n个数中选出k(1 ≤ k ≤ n)个数，对于某种选取方案，这种方案的美丽值为k个数的最大公约数乘上k。然后这个集合的美丽值为所有选数方案的美丽值之和。

然后他们想比比谁得到的美丽值更大。因为数很大，所以他们只比较各自的结果对258280327取余的大小。

做此题的过程真是艰难。。。。。我来分享下。。。。。

看起来完全没法做。不过数论题么，总应该是可以怎么乱搞一下或者推公式的。。。。。推出公式来要很半天，然后敲代码10分钟基本上就这样。。。。。

当然，按照题意去写公式很好写，然而公式是(n^2)的复杂度，需要利用各种数论定理或者化简技巧去化简。

我们先看浙理工先生的计算方法。对任意一个排列，计算所有区间GCD的和。所以直接枚举排列再枚举区间不就完了嘛，然而这显然是不可能的。于是我们转化一下方法，既然归根结底是求所有的区间GCD，我们不妨看看有哪些的区间可以用来计算。显然，在n个数中任取k个数，这k个数连续的放在一起，可以构成一个长度为k区间。因为k个数任意排列都不影响整个区间的GCD，所以我们把这整个区间看成一个物品，和剩下n - k个数一共组成n - k + 1个物品，又有(n - k + 1)!种排列方法（也就是捆绑法思想）。所以这个区间在所有的排列中，一共会被计算k!*(n - k + 1)!次。而每次计算都会加上这个区间的GCD值，不妨设为x，于是GCD为x，且长度为k的这样一个区间对整个美丽值的贡献即为x*k!*(n - k + 1)!。我们不妨设n个数中最大的数为MAXN。我们定义一个函数G(k, x)表示从n个数中选k个数且这个k个数最大公约数为x的方案数。于是我们得到了计算浙理工先生的计算公式：

ZSTU = sigma(1 ≤ k ≤ n，sigma(1 ≤ x ≤ MAXN，G(k，x)*x*k!*(n - k + 1)! ))。

然而，这样暴力枚举k和x是n^2复杂度，必然需要将公式化简。直观看来，两个sigma求和似乎都找不到什么好的公式。于是我们对G(k，x)下手。设法计算G(k，x)的值，看能不能有什么突破。考虑到G(k，x)是表示k个数最大公约数恰为x的方案数，恰为x不好计算，我们可以定义一个辅助函数——F(k，x)，表示n个数中选k个数的最大公约数为x的倍数的方案数。为什么用这个辅助函数呢。因为这个函数很好计算——既然gcd是x的倍数，那么就在所有x的倍数中选k个不就完了嘛。于是我们定一个cnt[x]表示n个数中x的倍数的个数，就得到F(k，x) = C(cnt[x]，k)（C是组合数，下同）。我们再来看看F函数和G函数的关系——很简单，其实就是F(k，x) = sigma(x | d，G(k，d))看到这种形式的等式，不由得就想到了莫比乌斯反演！于是G函数的计算就被转化为F了。

然而，貌似好像前进了一大步，不过，剩下来的问题更加困难——怎么化简公式。打字太麻烦。。。。。我就直接上图。。。。。计算中的μ函数是莫比乌斯函数。

注意这里第二次变形的时候把x和d的地位互换了，为什么呢。因为枚举约数x计算倍数d有个限制，那就是d不能超过MAXN，这样就比较麻烦，而转化为枚举倍数计算约数时，每个数的约数是固定的，于是计算起来就比较方便。所以对上式最后一个sigma，我们可以对任一的d直接预处理。我们记最后一个sigma所表示的值为H(d)。这个H(d)是可以离线打表计算的。继续化简——

化简到这里，貌似还是个(n^2)公式。不过内层的sigma已经是一个非常有规律的公式了——我们把它写出来看看：答案记为P(cnt[d])。然后我们来求这个P函数。为了简单点写。。。。。就暂时用x代替cnt[d]。

P(x)

= sigma(1 ≤ k ≤ x，(n - k + 1)! / (x - k)!

= n! / (x - 1)! + (n - 1)! / (x - 2)! + ...... + (n - x + 1)! / 0!

= n*(n - 1)*(n - 2)* ...... *x +

   (n - 1)*(n - 2)* ...... *(x - 1) +

   (n - 2)*(n - 3)* ...... *(x - 2) +

   (n - x + 1)*(n - x)* ...... *1

可以发现此和式的每一项都是连续(n - x + 1)个数相乘的结果，且每一项的连续数的开头都等于上一项的开头 + 1。对于这种类型的和式，有一种非常巧妙的求和办法——我举个简单例子：求 1*2*3 + 2*3*4 + ..... + n*(n + 1)*(n + 2)。我们记答案为ans。

对于每个i*(i + 1)*(i + 2)这一项，我们都让它乘上4（4的由来是每一项中连乘数的个数 + 1）然后，对于每一项中的4，都拆成( (i + 3) - (i - 1) )。为什么要这么拆呢？拆完带进去你就发现答案了——

4*ans

= 1*2*3*(4 - 0) + 2*3*4*(5 - 1) + ...... +(n*(n + 1)*(n + 2)*( (n + 3) - (n - 1) )

= 1*2*3*4 - 0*1*2*3 + 2*3*4*5 - 1*2*3*4 + 3*4*5*6 - 2*3*4*5 + ...... + n*(n + 1)*(n + 2)*(n + 3) - (n - 1)*n*(n + 1)*(n + 2)

中间项全都抵消了！于是答案为ans = n*(n + 1)*(n + 2)*(n + 3)/4。

这样，我们得到了P(x)的公式——P(x) = (n + 1)*n*(n - 1)* ...... *x / (n - x + 1) = (n + 1)! / (x - 1)! / (n - x + 1) —— 在预处理出逆元、阶乘取模、阶乘逆元后，这是个O(1)的公式，除法可以用逆元计算。至于逆元的求法不要用扩欧，因为会重复用到很多个数的逆元。这里我们用O(n)递推打表的公式求每个数的逆元inv[i]——公式为inv[i] = (MOD - MOD/i)*inv[MOD%i]%MOD，约定inv[0] = inv[1] = 1。这个公式我就不证了。。。。。这样，我们先预处理出莫比乌斯函数、H函数、逆元、阶乘取模和阶乘逆元等，最后的公式就完全成了个O(n)的公式。至此浙理工先生的答案已成功解出。

那么杭电先生的公式呢，比较简单。我就直接写了——其中内层G函数还是转化为F函数去求，方法和计算浙理工先生时一模一样。

HDU

= sigma(1 ≤ k ≤ n，sigma(1 ≤ x ≤ MAXN，G(k，x)*k))

= sigma(1 ≤ d ≤ n，H(d)*sigma(1 ≤ k ≤ cnt[d]，k*C(cnt[d]，k)))

因为k*C(cnt[d]，k) = cnt[d]*C(cnt[d] - 1, k - 1)，所以内层sigma

= sigma(1 ≤ k ≤ cnt[d]，k*C(cnt[d]，k)))

= sigma(1 ≤ k ≤ cnt[d]，cnt[d]*C(cnt[d] - 1，k - 1))

= cnt[d]*(sigma(0 ≤ k ≤ cnt[d] - 1，C(cnt[d] - 1，k))

= cnt[d]*2^(cnt[d] - 1)。

带回去后得到HDU

= sigma(1 ≤ d ≤ n，H(d)*cnt[d]*2^(cnt[d] - 1)。

OK，两位先生的答案都算出来了，再比较一下输出就可以了！

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <string>
#include <vector>
#include <map>
#include <stack>
#include <queue>
using namespace std;

const int MAX = 1e5 + 5;
const int MOD = 258280327;
bool u[MAX]; //素数筛
int prime[MAX], num; //素数表
int mobius[MAX]; //莫比乌斯函数表
long long inv[MAX]; //逆元表
long long factorial[MAX]; //阶乘表
long long fac_inv[MAX]; //阶乘逆元表
long long pow2[MAX]; //2的次方表
long long sum[MAX]; //H函数表

int n, max_value;
int cnt[MAX], cnt_temp[MAX]; //cnt_temp存储每个数出现了多少次，方便cnt数组的计算
long long P[MAX];

void initial() //预处理各种表
{
    memset(u, true, sizeof(u));
    num = 0;
    mobius[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        if(u[i])
        {
            prime[num++] = i;
            mobius[i] = -1;
        }
        for(int j = 0; j < num; j++)
        {
            if(i*prime[j] >= MAX)
                break;
            u[i*prime[j]] = false;
            if(i%prime[j] == 0)
            {
                mobius[i*prime[j]] = 0;
                break;
            }
            else
                mobius[i*prime[j]] = -mobius[i];
        }
    }

    inv[0] = inv[1] = 1;
    factorial[0] = factorial[1] = 1;
    fac_inv[0] = fac_inv[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        inv[i] = (MOD - MOD/i)*inv[MOD%i]%MOD;
        factorial[i] = factorial[i - 1]*i%MOD;
        fac_inv[i] = fac_inv[i - 1]*inv[i]%MOD;
    }

    memset(sum, 0, sizeof(sum));
    sum[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        for(int j = 1; j*j <= i; j++)
        {
            if(i%j == 0)
            {
                sum[i] += j*mobius[i/j];
                if(j*j != i)
                    sum[i] += i/j*mobius[j];
            }
        }
    }

    pow2[0] = 1;
    for(int i = 1; i < MAX; i++)
        pow2[i] = pow2[i - 1]*2%MOD;
}

void input()
{
    int x;
    memset(cnt_temp, 0, sizeof(cnt_temp));
    memset(cnt, 0, sizeof(cnt));
    max_value = 0;
    for(int i = 0; i < n; i++)
    {
        scanf("%d", &x);
        cnt_temp[x]++;
        max_value = max(max_value, x);
    }
}

void solve()
{
    for(int i = 1; i <= max_value; i++) //计算cnt数组
    {
        for(int j = i; j <= max_value; j += i) //枚举i的倍数
            cnt[i] += cnt_temp[j];
    }
    P[0] = 0; //不要忘了特殊值。
    for(int i = 1; i <= n; i++)
        P[i] = factorial[n + 1]*fac_inv[i - 1]%MOD*inv[n - i + 2]%MOD;
    long long ZSTU = 0, HDU = 0;
    for(int i = 1; i <= max_value; i++)
    {
        HDU = (HDU + sum[i]*cnt[i]%MOD*(pow2[cnt[i] - 1]))%MOD;
        ZSTU = (ZSTU + sum[i]*factorial[cnt[i]]%MOD*P[cnt[i]]%MOD)%MOD;
    }
    if(ZSTU > HDU)
        printf("Mr. Zstu %I64d\n", ZSTU);
    if(HDU > ZSTU)
        printf("Mr. Hdu %I64d\n", HDU);
    if(ZSTU == HDU)
        printf("Equal %I64d\n", HDU);
}

int main()
{
    //for(int i = 2; i*i <= MOD; i++)if(MOD%i == 0){printf("NO\n");return 0;}printf("YES\n"); 好吧它是个质数 ^_^
    initial();
    while(scanf("%d", &n) != EOF)
    {
        input();
        solve();
    }
    return 0;
}