起因：单纯的不想写代码，做不动题，状态非常差……于是就写点清新莫反。
记 (f(x)) 表示 (x) 的约数个数。
[sum_{i=1}^{n}sum_{j=1}^{m}f(i)f(j)f(gcd(i,j))\
sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=d]f(d)f(i)f(j)\
sum_{d=1}^{n}sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{m}{d}}[gcd(i,j)=1]f(d)f(id)f(jd)\
=sum_{d=1}^{n}f(d)sum_{i=1}^{frac{n}{d}}f(id)sum_{j=1}^{frac{m}{d}}f(jd)sum_{x|i,x|j}mu(x)\
=sum_{d=1}^{n}f(d)sum_{x=1}^{frac{n}{d}}mu(x)sum_{i=1}^{frac{n}{dx}}f(idx)sum_{j=1}^{frac{m}{dx}}f(jdx)\
=sum_{d=1}^{n}f(d)sum_{x=1}^{frac{n}{d}}mu(x)g(dx)h(dx)\
]
其中 (g(x)=sumlimits_{i=1}^{n}[x|i]f(i),h(x)=sumlimits_{i=1}^{m}[x|i]f(i))。
于是就可以 (O(nlog n)) 暴力求了，预处理 (f,g,h)，最后统计答案都是调和级数。不过跑的非常的慢，大力卡常之后才贴着时限跑。
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N = 2000005;
int f[N], g[N], h[N], n, m, mod, ans;
int pri[N], pct, mu[N];
bool vis[N];
inline void Sieve(const int &lim = N - 1) {
for(int i = 1; i <= lim; ++i)
for(int j = 1; i * j <= lim; ++j)
++f[i * j];
for(int i = 1; i <= n; ++i) {
LL tmp = 0;
for(int j = i; j <= n; j += i)
tmp += f[j];
g[i] = tmp % mod;
}
for(int i = 1; i <= m; ++i) {
LL tmp = 0;
for(int j = i; j <= m; j += i)
tmp += f[j];
h[i] = tmp % mod;
}
mu[1] = 1;
for(int i = 2; i <= lim; ++i) {
if(!vis[i]) mu[i] = -1, pri[++pct] = i;
for(int j = 1; j <= pct && i * pri[j] <= lim; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
mu[i * pri[j]] = 0;
break;
} else mu[i * pri[j]] = -mu[i];
}
}
}
signed main() {
cin >> n >> m >> mod;
//n = 2e6, m = 2e6, mod = 998244353;//264180408
if(n > m) n ^= m ^= n ^= m;
Sieve();
rep(d, 1, n) {
int res = 0;
rep(x, 1, n / d) {
if(mu[x] == 1)
res = (res + 1ll * g[d * x] * h[d * x]) % mod;
else if(mu[x] == -1)
res = (res + mod - 1ll * g[d * x] * h[d * x] % mod) % mod;
}
ans = (ans + 1ll * f[d] * res) % mod;
}
cout << ans << '
';
}
然后发现提交记录有人不到 100ms，就去看了下。好 nb 啊！接下去才是重点。
这题是可以 (O(nloglog n)) 的，一部分一部分看。
首先约数个数和线性筛，多维护每个数最小质因子的次数，降为 (O(n))。
(g,h) 可以狄利克雷后缀和 (O(nloglog n))。
统计答案的式子可以再化，可以看出有个狄利克雷卷积：
[=sum_{d=1}^{n}f(d)sum_{x=1}^{frac{n}{d}}mu(x)g(dx)h(dx)\
=sum_{i=1}^{n}g(i)h(i)sum_{j|i}f(j)mu(dfrac{i}{j})
]
看一看 (d*mu) 是什么：(d=I*I,d*mu=I*I*mu=I*epsilon=I)。
于是后面那个恒为 (1)，答案就是 (sumlimits_{i=1}^{n}g(i)h(i))。
甚至不用筛 (mu)……
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N = 2000005;
int f[N], g[N], h[N], n, m, mod, ans, sf[N];
int pri[N], pct, low[N];
bool vis[N];
inline int add(int x, int y) { return (x += y) >= mod ? x - mod : x; }
inline void Sieve(const int &lim = N - 1) {
f[1] = 1;
for(int i = 2; i <= lim; ++i) {
if(!vis[i]) low[i] = 1, f[i] = 2, pri[++pct] = i;
for(int j = 1; j <= pct && i * pri[j] <= lim; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
f[i * pri[j]] = f[i] / (low[i] + 1) * (low[i] + 2),
low[i * pri[j]] = low[i] + 1;
break;
} else {
f[i * pri[j]] = f[i] * 2,
low[i * pri[j]] = 1;
}
}
}
rep(i, 1, n) g[i] = f[i];
rep(i, 1, m) h[i] = f[i];
for(int i = 1; i <= lim; ++i) sf[i] = add(sf[i], f[i - 1]);
for(int i = 1; i <= pct; ++i)
for(int j = n / pri[i]; j >= 1; --j)
g[j] = add(g[j], g[j * pri[i]]);
for(int i =1; i <= pct; ++i)
for(int j = m / pri[i]; j >= 1; --j)
h[j] = add(h[j], h[j * pri[i]]);
}
signed main() {
cin >> n >> m >> mod;
// n = 2e6, m = 2e6, mod = 998244353;//264180408
if(n > m) n ^= m ^= n ^= m;
Sieve();
rep(i, 1, n) ans = add(ans, 1ll * g[i] * h[i] % mod);
cout << ans << '
';
}
[sum_{i_1=L}^{H}sum_{i_2=L}^{H}cdotssum_{i_+n=L}^{H}[gcd(i_1,i_2,cdots,i_n)=K]\
]
记 (A=lfloordfrac{L-1}{K} floor +1,B=lfloordfrac{R}{K} floor)，上式可以化简：
[sum_{i_1=A}^{B}sum_{i_2=A}^{B}cdotssum_{i_n=A}^{B}[gcd(i_1,i_2,cdots,i_n)=1]\=sum_{i_1=A}^{B}sum_{i_2=A}^{B}cdotssum_{i_n=A}^{B}sum_{x|i_1,x|i_2,cdots,x|i_n}mu(x)\=sum_{x=1}^{B}mu(x)sum_{i_1=A}^{B}[x|i_1]sum_{i_2=A}^{B}[x|i_2]cdotssum_{i_n=A}^{B}[x|i_n]\=sum_{x=1}^{B}mu(x)(lfloordfrac{B}{x}
floor-lfloordfrac{A-1}{x}
floor)^n
]
后面整除分块前面杜教筛，不懂 (H-Lle 10^5) 有什么用……
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
#define mod 1000000007
const int N = 1000005;
int n, K, L, H, ans;
int pri[N], pct, smu[N];
bool vis[N];
inline void Sieve(const int&n = N - 1) {
smu[1] = 1;
for(int i = 2; i <= n; ++i) {
if(!vis[i]) pri[++pct] = i, smu[i] = mod - 1;
for(int j = 1; i * pri[j] <= n; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
smu[i * pri[j]] = 0;
break;
}
smu[i * pri[j]] = mod - smu[i];
}
}
for(int i = 1; i <= n; ++i) smu[i] = (smu[i] + smu[i - 1]) % mod;
}
unordered_map<int, int> _mu;
int sum_mu(int n) {
if(n < N) return smu[n];
if(_mu.find(n) != _mu.end()) return _mu[n];
int res = 1;
for(int l = 2, r; l <= n; l = r + 1)
r = n / (n / l), res = (res + mod - 1ll * (r - l + 1) * sum_mu(n / l) % mod) % mod;
return _mu[n] = res;
}
inline int qpow(int n, int k) {
int res = 1;
for(; k; k >>= 1, n = 1ll * n * n % mod)
if(k & 1) res = 1ll * n * res % mod;
return res;
}
signed main() {
Sieve();
cin >> n >> K >> L >> H, --L;
int A = L / K, B = H / K;
for(int l = 1, r; l <= B; l = r + 1) {
r = B / (B / l);
if(A >= l) ckmin(r, A / (A / l));
ans = (ans + 1ll * qpow(B / l - A / l, n) * (sum_mu(r) - sum_mu(l - 1) + mod)) % mod;
}
cout << ans << '
';
}
挺不错的一道题，然而我没做出来/kk，体现了我只会做板子题的事实。
[dfrac{1}{a}+dfrac{1}{b}=dfrac{1}{c}\
bc+ac=ab\
ab-bc-ac+c^2=c^2\
(a-c)(b-c)=c^2
]
因为 (gcd(a,b,c)=1)，所以 (gcd(a-c,b-c,c)=1)。
然后就是非常神奇的一步，我就卡在这里了：
[gcd(a-c,b-c,c)=1Rightarrow gcd(a-c,b-c)=1
]
原因：如果存在一个数 (x) 满足 (x|a-c,x|b-c)，那么根据 ((a-c)(b-c)=c^2)，左边必然有因子 (x^2)，那么就有 (x|c)。
由于 (gcd(a-c,b-c)=1)，那么 (c^2) 的每一种质因子只能全给 (a-c) 或 (b-c)。所以 (a-c,b-c) 都是完全平方数。
设 (x^2=a-c,y^2=b-c)，那么 (c=xy)。
因为 (gcd(a-c,b-c)=1)，所以 (gcd(x,y)=1)。
所以我们可以枚举 (x)，这样枚举的范围只有 (sqrt{n})。记 (m=sqrt{n})
式子是：
[sum_{x=1}^{m}sum_{yge 1}[gcd(x,y)=1][1le ale n][1le ble n][1le cle n]\=sum_{x=1}^{m}sum_{yge 1}[gcd(x,y)=1][1le x^2+xyle n][1le y^2+xyle n][1le xyle n]
]
对 (c) 的限制显然弱于对 (a,b) 的限制，直接去掉：
[=sum_{x=1}^{m}sum_{yge 1}[gcd(x,y)=1][1le x^2+xyle n][1le y^2+xyle n]
]
发现后两个限制只是限制了 (y) 的范围，不妨直接解出来。
[x^2+xyle nRightarrow yledfrac{n-x^2}{x}\y^2+xyle nLeftrightarrow y^2+xy-nle 0Leftrightarrow \dfrac{-x-sqrt{x^2+4n}}{2}le yle dfrac{-x+sqrt{x^2+4n}}{2}
]
令 (a_x=min(dfrac{n-x^2}{x},dfrac{-x+sqrt{x^2+4n}}{2}))，那么答案就是
[sum_{i=1}^{m}sum_{j=1}^{a_i}[gcd(i,j)=1]
]
一个人口普查形式的莫反。不过注意不能根号。算了还是写到底吧。
[sum_{i=1}^{m}sum_{j=1}^{a_i}[gcd(i,j)=1]\=sum_{i=1}^{m}sum_{x|i}mu(x)sum_{j=1}^{a_i}[x|j]\=sum_{x=1}^{m}mu(x)sum_{x|i}dfrac{a_i}{x}
]
可以调和级数计算。
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N = 1000005;
int mu[N], pct, pri[N];
bool vis[N];
LL n, ans, a[N];
int m;
inline void Sieve(const int &n = N - 1) {
mu[1] = 1;
for(int i = 2; i <= n; ++i) {
if(!vis[i]) pri[++pct] = i, mu[i] = -1;
for(int j = 1; j <= pct && i * pri[j] <= n; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
mu[i * pri[j]] = 0;
break;
} else mu[i * pri[j]] = -mu[i];
}
}
}
signed main() {
cin >> n, m = sqrt(n);
Sieve(m);
rep(i, 1, m) a[i] = min((n - 1ll * i * i) / i, (LL)(sqrt(1ll * i * i + 4ll * n) - i) / 2);
rep(i, 1, m) {
LL tmp = 0;
for(int j = 1; i * j <= m; ++j) tmp += 1ll * a[i * j] / i;
ans += tmp * mu[i];
}
cout << ans << '
';
}
好像是曾经的月赛题，场上没做出来，现在还是不会做……
考虑怎么化简这个 (p(a,b))。
首先注意到，(gcd(a,b)>1) 就不行。又注意到 (amod 2 = bmod 2) 也不行。剩下的凭直觉就感觉都可以啊，打个表感觉很对啊，往上一交就 20 了啊。
也就是这个东西：
inline int p(int x, int y) {
return ((x & 1) ^ (y & 1)) && (__gcd(x, y) == 1);
}
所以答案就是
[sum_{i=1}^{n}sum_{j=1}^{n}[gcd(i,j)=1][imod 2
ot= jmod 2]
]
不妨钦定 (i) 为偶数 (j) 为奇数，问题转化成 (sum_{i=1}^{frac{n}{2}}sum_{j=1}^{n}[gcd(i,j)=1])，后面那个东西可以根号，人口普查形式的莫反不再多说，复杂度 (O(nsqrt{n})) 得到 35。
感觉这个思路没得救，换个思路。
考虑计算 (f(x)=sum_{i=1}^{x}[imod 2 ot= x mod 2][gcd(x,i)=1])。
对于 (xmod 2 = 0)，直接 (f(x)=varphi(x)) 即可，因为所有偶数都不会被统计在欧拉函数内。
对于 (xmod 2 = 1)，(f(x)=dfrac{varphi(x)}{2})，因为欧拉函数每统计到一个与 (x) 互质的奇数数 (y)，都会再统计到一个与 (x) 互质的偶数 (x-y)，所以除以二。
特判 (x=1)，最后答案乘 (2) 即可，线性筛 (varphi)，可以得到 50 分。
看看现在的式子（把 (2) 先乘进去）：
[sum_{i=1}^{n}left([imod 2 = 0]2varphi(i)+[imod 2 = 1]varphi(i)
ight)-1
]
设 (f(n)=sum_{i=1}^{n}[imod 2 = 0]2varphi(i)+[imod 2 = 1]varphi(i))
注意到可以变形成 (f(n)=sum_{i=1}^{n}varphi(i)+sum_{i=1}^{n}[imod 2 = 0]varphi(i))。
而 (sum_{i=1}^{n}[imod 2 = 0]varphi(i))，可以发现它就等于 (f(lfloordfrac{n}{2} floor))！
具体分析就是，如果 (imod 2 = 0)，那么 (varphi(2i)=2varphi(i))，否则 (varphi(2i)=varphi(i))。分讨一下 (le lfloordfrac{n}{2} floor) 的奇偶性，发现刚好就是。
所以 (f(n)=sum_{i=1}^{n}varphi(i)+f(lfloordfrac{n}{2} floor))，直接递归，加个杜教筛，就好了！
无限 orz 出题人 QuantAsk！！！
Code
#include <cstdio>
#include <tr1/unordered_map>
using namespace std::tr1;
const int N = 10000005;
typedef unsigned long long ull;
int pri[N], pct;
ull phi[N], g[N];
bool vis[N];
inline void Sieve(const int &n = N - 1) {
phi[1] = 1;
for(int i = 2; i <= n; ++i) {
if(!vis[i]) pri[++pct] = i, phi[i] = i - 1;
for(int j = 1; j <= pct && i * pri[j] <= n; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
}
}
for(int i = 1; i <= n; ++i)
g[i] = g[i - 1] + (i & 1 ? phi[i] : phi[i] + phi[i]), phi[i] += phi[i - 1];
}
unordered_map<ull, ull> _phi;
ull sum_phi(ull n) {
if(n < N) return phi[n];
if(_phi.find(n) != _phi.end()) return _phi[n];
ull res = n & 1 ? (n + 1) / 2 * n : n / 2 * (n + 1);
for(ull l = 2, r; l <= n; l = r + 1)
r = n / (n / l), res -= (r - l + 1) * sum_phi(n / l);
return _phi[n] = res;
}
ull f(ull n) {
if(n < N) return g[n];
return f(n / 2) + sum_phi(n);
}
signed main() {
Sieve();
int T; ull n;
for(scanf("%d", &T); T--;)
scanf("%lld", &n), printf("%llu
", f(n) - 1);
}