【BZOJ3960】DZY Loves Math V（数论）

题目：

BZOJ3560

分析：

orz跳瓜。

欧拉函数的公式：

[phi(n)=n(prod frac{p_i-1}{p_i}) ]

其中 (p_i) 取遍 (n) 的所有质因子。

考虑原式，把欧拉函数展开，得到：

[sum_{b_1|a_1}sum_{b_2|a_2}cdotssum_{b_n|a_n}prod b_i prod frac{(p_j-1)}{p_j}= sum_{b_1|a_1}sum_{b_2|a_2}cdotssum_{b_n|a_n}prod p_j^{a_j}frac{(p_j-1)}{p_j}]

其中 (p_j) 取遍 (prod b_i) 的所有质因子， (prod p_j^{a_j}=prod b_i)。

可以看出，对于每个质数 (p) ，它的贡献是独立的。设 (a_j) 中 (p_i) 的次数为 (c_{ij}) ，则答案为：

[prod_i left(1+frac{p_i-1}{p_i}prod_{j=1}^n sum_{k_j=0}^{c_{ij}}left[sum k_j>0 ight]p_i^{sum k_j} ight) ]

其中 (k_j) 表示在第 (j) 个数中取了 (p_i) 的 (k_j) 次幂。当 (sum k_j=0) （即 (prod b_i) 不含质因子 (p_i) ）时不乘 (frac{p_i-1}{p_i}) ，这种情况中 (p_i) 对答案的贡献是乘 (1) （即最前面加的 (1) ）。

上面的式子相当于（括号比较鬼畜，凑合看吧……）：

[prod_i left(1+frac{p_i-1}{p_i}left(left(prod_{j=1}^n sum_{k_j=0}^{c_{ij}}p_i^{sum k_j} ight)-1 ight) ight) ]

（即直接减去不合法的 (p_i^{sum k_j}=p_i^0=1) ）

然后可以暴力算。对于每一个质数预处理出它的幂的前缀和，并可以只枚举有质因子 (p_i) 的 (a_j) 。由于每个数的质因子种类数最多 (10) 个左右，所以总复杂度大约 (O(10n))，可以过。具体实现详见代码。

代码：

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cctype>
using namespace std;
 
namespace zyt
{
    template<typename T>
    inline bool read(T &x)
    {
        char c;
        bool f = false;
        x = 0;
        do
            c = getchar();
        while (c != EOF && c != '-' && !isdigit(c));
        if (c == EOF)
            return false;
        if (c == '-')
            f = true, c = getchar();
        do
            x = x * 10 + c - '0', c = getchar();
        while (isdigit(c));
        if (f)
            x = -x;
        return true;
    }
    template<typename T>
    inline void write(T x)
    {
        static char buf[20];
        char *pos = buf;
        do
            *pos++ = x % 10 + '0';
        while (x /= 10);
        while (pos > buf)
            putchar(*--pos);
    }
    typedef long long ll;
    const int N = 1e5 + 10, MX = 1e6 + 10, B = 30, mod = 1e9 + 7;
    int n, cnt, pow[B];
    pair<int, int> prime[MX];
    void get_prime(int a)
    {
        for (int i = 2; (ll)i * i <= a; i++)
            if (a % i == 0)
            {
                prime[cnt++] = make_pair(i, 0);
                while (a % i == 0)
                    ++prime[cnt - 1].second, a /= i;
            }
        if (a > 1)
            prime[cnt++] = make_pair(a, 1);
    }
    inline int power(int a, int b)
    {
        int ans = 1;
        while (b)
        {
            if (b & 1)
                ans = (ll)ans * a % mod;
            a = (ll)a * a % mod;
            b >>= 1;
        }
        return ans;
    }
    inline int inv(const int a)
    {
        return power(a, mod - 2);
    }
    int cal(const int l, const int r)
    {
        int p = prime[l].first;
        int ans = 1, mx = 0;
        for (int i = l; i <= r; i++)
            mx = max(mx, prime[i].second);
        pow[0] = 1;
        for (int i = 1; i <= mx; i++)
            pow[i] = (ll)pow[i - 1] * p % mod;
        for (int i = 1; i <= mx; i++)
            pow[i] = (pow[i] + pow[i - 1]) % mod;
        for (int i = l; i <= r; i++)
            ans = (ll)ans * pow[prime[i].second] % mod;
        ans = (ans - 1 + mod) % mod;
        ans = (ll)ans * (p - 1) % mod * inv(p) % mod;
        return (ans + 1) % mod;
    }
    int work()
    {
        int n;
        read(n);
        for (int i = 0; i < n; i++)
        {
            int a;
            read(a);
            get_prime(a);
        }
        sort(prime, prime + cnt);
        int pre = 0, ans = 1;
        for (int i = 0; i < cnt; i++)
            if (prime[i].first != prime[i + 1].first)
                ans = (ll)ans * cal(pre, i) % mod, pre = i + 1;
        write(ans);
        return 0;
    }
}
int main()
{
    return zyt::work();
}