题目地址：BZOJ：Problem 3560. — DZY Loves Math V

题目描述

给定n个正整数a1,a2,…,an，求 \displaystyle \sum_{i_1|a_1} \sum_{i_2|a_2} \cdots \sum_{i_n|a_n} \varphi(i_1 i_2 \cdots i_n)的值（答案模10^9+7）。

输入输出格式

输入格式：
第一行一个正整数n。
接下来n行，每行一个正整数，分别为a1,a2,…,an。

输出格式：
仅一行答案。

输入输出样例

输入样例#1：

3
6
10
15

输出样例#1：

1595

说明

1<=n<=10^5，1<=ai<=10^7。共3组数据。

题解

首先考虑到\varphi(x)是积性函数，因此我们可以对ai中的每一个质因数讨论，如果某个因数p在ai中的个数是bi，则计算出每个因数的
\sum_{i_1=0}^{b_1} \sum_{i_2=0}^{b_2} \cdots \sum_{i_n=0}^{b_n} \varphi(p^{i_1+i_2+\cdots+i_n})
并乘起来就得到了最终结果。然而强行计算这个值复杂度不可接受，我们需要给出一种优化算法，要对这个式子化简。
我们知道\varphi(x)有一种求法，如果x的所有质因数分别为pi，则
\varphi(x) = x(1 - \frac{1}{p_1})(1 - \frac{1}{p_2}) \cdots (1 - \frac{1}{p_n})
那么我们观察每个因数的结果式，显然1-分数只有一项即 (1-\frac{1}{p}) ，前面的x可以求和，即化简成下面这种形式
\left[ \left( \sum_{i_1=0}^{b_1} \sum_{i_2=0}^{b_2} \cdots \sum_{i_n=0}^{b_n} p^{i_1+i_2+\cdots+i_n} \right) - 1 \right] (1-\frac{1}{p}) + 1
括号里的-1是减去p_0的情况，这种情况我们在外面计算，即\varphi(1) = 1。前面的一坨东西仍然会使复杂度变大，但是这种形式已经好优化了
\left[ \left( \prod_{i=1}^n \sum_{j=0}^{b_i} p^j \right) - 1 \right] (1-\frac{1}{p}) + 1
这样我们只需要对p的乘方求前缀和，就可以实现O(n)计算每个因数了。总复杂度O(n \sqrt{n})。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>

#include <algorithm>
#include <vector>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for (; !isdigit(c); c = fgc()) if (c == '-') neg = -1;
    for (; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

const int MAXN = 1000005, MO = 1e9 + 7;

inline LL fpow(LL x, LL k) {
    LL res = 1;
    for(; k; k >>= 1) {
        if(k & 1) res = res * x % MO;
        x = x * x % MO;
    }
    return res;
}

int n, a[100005];
typedef std::pair<int, int> PII;
std::vector<PII> fact;

inline void calfact(int x) {
    for(int i = 2; i * i <= x; i++) {
        if(x % i) continue;
        int cnt = 0;
        while(x % i == 0) {
            x /= i; cnt++;
        }
        fact.push_back(PII(i, cnt));
    }
    if(x > 1) fact.push_back(PII(x, 1));
}

LL sum[50];

inline LL calc(int l, int r) {
    sum[0] = 1; int p = fact[l].first;
    for(int i = 1; i <= fact[r].second; i++) {
        sum[i] = sum[i - 1] * p % MO;
    }
    for(int i = 1; i <= fact[r].second; i++) {
        sum[i] = (sum[i - 1] + sum[i]) % MO;
    }
    LL res = 1;
    for(int i = l; i <= r; i++) {
        res = res * sum[fact[i].second] % MO;
    }
    res--; 
    res = res * (p - 1) % MO; res = res * fpow(p, MO - 2) % MO;
    res++; res %= MO;
    return res;
}

int main() {
    n = readint();
    for(int i = 1; i <= n; i++) {
        a[i] = readint(); calfact(a[i]);
    }
    std::sort(fact.begin(), fact.end());
    LL ans = 1; int lst = -1;
    for(int i = 0; i < fact.size(); i++) {
        if(i == fact.size() - 1 || fact[i].first != fact[i + 1].first) {
            ans = ans * calc(lst + 1, i) % MO; lst = i;
        }
    }
    printf("%lld", ans);
    return 0;
}