传送门:http://acm.hdu.edu.cn/showproblem.php?pid=6833
借鉴大佬: https://blog.csdn.net/weixin_44282912/article/details/107844614 https://blog.csdn.net/oampamp1/article/details/107865300
题意
求解n∑a1=1n∑a2=1…n∑ax=1(x∏j=1akj)f(gcd(a1,a2…ax))gcd(a1,a2…ax)
思路
令d=gcd(a1,a2…ax)得:
n∑a1=1n∑a2=1…n∑ax=1(x∏j=1akj)f(d)d
枚举d得:
n∑d=1df(d)n∑a1=1n∑a2=1…n∑ax=1(x∏j=1akj)[gcd(a1,a2…ax)=d]
n∑d=1dkx+1f(d)⌊nd⌋∑a1=1⌊nd⌋∑a2=1…⌊nd⌋∑ax=1(x∏j=1akj)[gcd(a1,a2…ax)=1]
n∑d=1dkx+1f(d)⌊nd⌋∑a1=1⌊nd⌋∑a2=1…⌊nd⌋∑ax=1(x∏j=1akj)∑ta1ta2…taxμ(t)
n∑d=1dkx+1f(d)⎛⎜
⎜⎝⌊nd⌋∑i=1ik⎞⎟
⎟⎠x∑ta1ta2…taxμ(t)
枚举t得:
n∑d=1dkx+1f(d)⌊nd⌋∑t=1μ(t)⎛⎜
⎜⎝⌊ndt⌋∑i=1(it)k⎞⎟
⎟⎠x
n∑d=1dkx+1f(d)⌊nd⌋∑t=1μ(t)tkx⎛⎜
⎜⎝⌊ndt⌋∑i=1ik⎞⎟
⎟⎠x
令T=dt并枚举T得:
n∑T=1⎛⎜
⎜⎝⌊ndt⌋∑i=1ik⎞⎟
⎟⎠xTkx∑dTdf(d)μ(Td)
令g(T)=∑dTdf(d)μ(Td),则
n∑T=1⎛⎜
⎜⎝⌊nT⌋∑i=1ik⎞⎟
⎟⎠xTkxg(T)
O(nlogn)筛f和g,并对Tkxg(T)和∑⌊nT⌋i=1ik做前缀和,最后√n分块计算,复杂度为O(nlogn+T√n)
Code(2121MS)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
| ##include <bits/stdc++.h>
using namespace std;
typedef long long ll; typedef long double ld; typedef pair<int, int> pdd;
##define INF 0x3f3f3f3f ##define lowbit(x) x & (-x) ##define mem(a, b) memset(a , b , sizeof(a)) ##define FOR(i, x, n) for(int i = x;i <= n; i++)
const ll mod = 1e9 + 7;
const int N = 2e5 + 10;
ll T, k, x;
int mu[N]; bool is_prime[N]; int prime[N]; int cnt; ll g[N], f[N]; ll sumG[N], sumi[N];
ll quick_pow(ll a, ll b) { ll ans = 1; while(b) { if(b & 1) ans = ans * a % mod; a = a * a % mod; b >>= 1; } return ans % mod; }
void Mobi() { f[1] = mu[1] = 1; is_prime[0] = is_prime[1] = true; for(int i = 2;i < N; i++) { f[i] = 1; if (!is_prime[i]) { mu[i] = -1; prime[++cnt] = i; } for (int j = 1; j <= cnt && i * prime[j] < N; j++) { is_prime[i * prime[j]] = true; if (i % prime[j] == 0) { mu[i * prime[j]] = 0; break; } mu[i * prime[j]] = -mu[i]; } } for(int d = 2;d * d < N; d++) { for(int i = d * d;i < N; i += d * d) { f[i] = 0; } } for(int d = 1;d < N; d++) { for(int i = d;i < N; i += d) { g[i] = (g[i] + 1ll * d * f[d] % mod * mu[i / d] % mod + mod) % mod; } } for(int i = 1;i < N; i++) { ll t = quick_pow(i, k); sumi[i] = (sumi[i - 1] + t) % mod; sumG[i] = (sumG[i - 1] + quick_pow(t, x) * g[i] % mod) % mod; } }
void solve() { cin >> T >> k >> x; Mobi(); while(T--) { int n; cin >> n; ll ans = 0; for(int l = 1, r;l <= n; l = r + 1) { r = min(n, n / (n / l)); ans = (ans + (sumG[r] - sumG[l - 1] + mod) % mod * quick_pow(sumi[n / l], x) % mod) % mod; }
cout << ans << endl; } }
signed main() { ios_base::sync_with_stdio(false); ##ifdef FZT_ACM_LOCAL freopen("in.txt", "r", stdin); freopen("out.txt", "w", stdout); signed test_index_for_debug = 1; char acm_local_for_debug = 0; do { if (acm_local_for_debug == '$') exit(0); if (test_index_for_debug > 20) throw runtime_error("Check the stdin!!!"); auto start_clock_for_debug = clock(); solve(); auto end_clock_for_debug = clock(); cout << "Test " << test_index_for_debug << " successful" << endl; cerr << "Test " << test_index_for_debug++ << " Run Time: " << double(end_clock_for_debug - start_clock_for_debug) / CLOCKS_PER_SEC << "s" << endl; cout << "--------------------------------------------------" << endl; } while (cin >> acm_local_for_debug && cin.putback(acm_local_for_debug)); ##else solve(); ##endif return 0; }
|
本文作者:jujimeizuo
本文地址: https://blog.jujimeizuo.cn/2020/10/07/hdu-6833/
本博客所有文章除特别声明外,均采用 CC BY-SA 3.0 协议。转载请注明出处!