BZOJ3560: DZY Loves Math V(欧拉函数)

Time Limit: 10 Sec  Memory Limit: 256 MB
Submit: 557  Solved: 318
[Submit][Status][Discuss]

Description

给定n个正整数a1,a2,…,an，求

的值（答案模10^9+7）。

Input

第一行一个正整数n。

接下来n行，每行一个正整数，分别为a1,a2,…,an。

Output

仅一行答案。

Sample Input

3
6
10
15

Sample Output

1595

HINT

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

Source

By Jcvb

将$a_i$分解质因数后$p$的出现次数设为$b_i$

那么我们要求的就是

$sum_{i_1 = 0}^{b_1} sum_{i_2 = 0}^{b_2} dots sum_{i_n = 0}^{b_n} phi( p^{sum_{j = 1}^n i_j})$

考虑到$phi(p^k) = p^k - p^{k - 1}$

同时我们需要特判一下$1$的情况！

上式可以化为

$sum_{i_1 = 0}^{b_1} sum_{i_2 = 0}^{b_2} dots sum_{i_n = 0}^{b_n}( p^{sum_{j = 1}^n i_j} - 1) * frac{p - 1}{p} + 1$

再重新考虑每一个$p$

$[(prod_{i = 1}^{n} sum_{i = 0}^{b_i} p^j) - 1] * frac{p - 1}{p} + 1$

#include<cstdio>
#include<algorithm>
#define LL long long 
#define int long long 
using namespace std;
const int MAXN = 1e6 + 10, mod = 1e9 + 7;
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int T, N, M;
struct Node {
    int p, a;
    bool operator < (const Node &rhs) const {
        return p == rhs.p ? a < rhs.a : p < rhs.p;
    }
}s[MAXN];
int tot = 0;
void insert(int x) {
    for(int i = 2; i * i <= x; i++) {
        if(!(x % i)) {
            s[++tot].p = i;
            while(!(x % i)) 
                x /= i, s[tot].a++;
        }
    }
    if(x > 1) s[++tot] = (Node) {x, 1};
}
int fastpow(int a, int p, int mod) {
    int base = 1;
    while(p) {
        if(p & 1) base = (base * a) % mod;
        a = (a * a) % mod; p >>= 1;
    }
    return base % mod;
}
int calc(int l, int r) {
    static int sum[31] = {}; sum[0] = 1;
    for(int i = 1; i <= s[r].a; i++) sum[i] = sum[i - 1] * s[l].p % mod;
    for(int i = 1; i <= s[r].a; i++) sum[i] = (sum[i - 1] + sum[i]) % mod;
    int rt = 1;
    for(int i = l; i <= r; i++) rt = (rt * sum[s[i].a]) % mod;
    rt--; 
    rt = (rt * fastpow(s[l].p, mod - 2, mod) + mod) % mod;
    rt = (rt * (s[l].p - 1) + mod) % mod;
    return (rt + 1) % mod;
}
main() {
//    freopen("a.in", "r", stdin);
    N = read();
    for(int i = 1; i <= N; i++) 
        insert(read());
    sort(s + 1, s + tot + 1);
    int pre = 0, ans = 1;
    for(int i = 1; i <= tot; i++) 
        if(s[i].p != s[i + 1].p || i == tot)
            ans = (ans * calc(pre + 1, i)) % mod, pre = i;
    printf("%lld", (ans + mod) % mod);
}
/*
3
600
1010
15010

*/