P3768 简单的数学题
题解
题意
给定n,p 求:
(i=1∑nj=1∑nijgcd(i,j))modp
其中,n≤1010,5.8×108≤p≤1.1×109
解法
首先观察到n 的数据范围,肯定是要以一个亚线性时间复杂度跑过去的,须有用到杜教筛或Min_25筛,这里用的杜教筛。
首先化简式子:
i=1∑nj=1∑nijgcd(i,j)=d=1∑ndi=1∑nj=1∑nij[gcd(i,j)==d]=d=1∑nd3i=1∑⌊dn⌋j=1∑⌊dn⌋ij[gcd(i,j)==1]
后面是个经典的式子,设为S(n),接着化简:
S(n)=i=1∑nj=1∑nij[gcd(i,j)==1]=i=1∑nj=1∑nijd∣gcd(i,j)∑μ(d)=d=1∑nμ(d)d2i=1∑⌊dn⌋j=1∑⌊dn⌋1
后面的式子可以O(1)处理,再套上一个数论分块,然后就可以得到60分了,然而这并不能解决这道题,其实可以在其中的一步换一种方式去卷mu,但是有更简单的做法,重新考虑式子,对其使用欧拉反演:
有:
φ∗1=Id
即:
d∣n∑φ(d)=n
原来的式子可以化为:
i=1∑nj=1∑nijgcd(i,j)=i=1∑nj=1∑nijd∣i,d∣j∑φ(d)=d=1∑nφ(d)d2i=1∑⌊dn⌋j=1∑⌊dn⌋ij
显然后面的式子可以O(1)计算,
4n2(n+1)2
现在考虑如何计算前面的
d=1∑nφ(d)d2
既然φ∗1=Id,将φ(d)d2设为f,设g=Id2显然有
(f∗g)n=d∣n∑φ(d)d2(dn)2=n3
看来问题已经解决了,S(n)为我们想要求的和,
i=1∑n(f∗g)i=i=1∑ng(i)S(⌊in⌋)=i=1∑ni3
然后使用技巧,
g(1)S(n)=i=1∑n(f∗g)i−i=2∑ng(i)S(⌊in⌋)
又因为g(1)=1,得到递归式,
S(n)=i=1∑n(f∗g)i−i=2∑ng(i)S(⌊in⌋)
需要用到:
i=1∑ni2=6n(n+1)(2n+1)i=1∑ni3=4n2(n+1)2
用杜教筛,时间复杂度O(n32)
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
| #include <bits/stdc++.h> #include <unordered_map> using namespace std;
#define int long long
const int N = 5e6 + 10;
int mod;
int inv4, inv6;
inline int qpow(int a, int b) { int t = 1; while(b != 0) { if(b & 1)t = t * a%mod; a = a * a%mod; b >>= 1; } return t%mod; }
int cnt, prime[N], phi[N];
int S_phi[N];
bool vis[N];
inline void init(int n) { phi[1] = 1; for(int i = 2; i <= n; i++) { if(!vis[i]) { prime[++cnt] = i; phi[i] = i - 1; } for(int j = 1; j <= cnt && prime[j] * i <= n; j++) { vis[i * prime[j]] = true; if(i % prime[j] == 0) { phi[i * prime[j]] = prime[j] * phi[i]; break; } phi[i * prime[j]] = phi[i] * phi[prime[j]]; } } for(int i = 1; i <= n; i++) S_phi[i] = (S_phi[i-1] + phi[i]%mod * i%mod * i%mod)%mod; }
inline int f(int n) { return n %mod* n%mod * (n + 1ll)%mod * (n + 1)%mod *inv4 %mod; }
inline int g(int n) { return n%mod * (n + 1)%mod * (n * 2ll%mod + 1ll)%mod *inv6 %mod; }
unordered_map<int, int> Phi;
int S(int n) { if(n <= 5e6)return S_phi[n]%mod; if(Phi.count(n))return Phi[n]; int sum = f(n)%mod; for(int l = 2, r; l <= n; l = r + 1) { r = n/(n/l); sum = (sum - S(n/l)%mod * (g(r) - g(l-1)) %mod + mod)%mod; } return Phi[n] = sum; }
signed main() { int n; scanf("%lld%lld", &mod, &n); inv4 = qpow(4, mod-2); inv6 = qpow(6, mod-2); init(5e6); int ans = 0; for(int l = 1, r; l <= n; l = r + 1) { r = n/(n/l); ans = (ans + f(n/l)%mod * (S(r) - S(l-1) + mod)%mod)%mod; } printf("%lld", ans); return 0; }
|
tips: 欧拉反演的证明
d∣n∑φ(d)=n
设f(n)=∑d∣nφ(d), m⊥n,根据欧拉函数的积性函数的性质有:
f(n)×f(m)=i∣n∑φ(i)j∣m∑φ(j)=i∣n∑j∣m∑φ(i)φ(j)=i∣n∑j∣m∑φ(ij)=d∣nm∑φ(d)=f(nm)
将n质因数分解p1c1×p2c2×p3c3×⋯×pkck。
所以有,
f(n)=f(p1c1)×f(p2c2)×⋯×f(pkck)
其中的每一项有:
f(pc)=i=0∑cφ(pi)=pc
所以得证:
f(n)=f(p1c1)×f(p2c2)×⋯×f(pkck)=Πi=1kpici=n=d∣n∑φ(d)