[BZOJ3309]DZY Loves Math
[BZOJ3309]DZY Loves Math
试题描述
对于正整数 n,定义 f(n) 为 n 所含质因子的最大幂指数。例如 f(1960)=f(23×51×72)=3,f(10007)=1,f(1)=0。
给定正整数 a,b,求 ∑ni=1∑mj=1f(gcd(i,j))。
输入
第一行一个数 T,表示询问数。
接下来 T 行,每行两个数 a,b,表示一个询问。
输出
对于每一个询问,输出一行一个非负整数作为回答。
输入示例
4
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957
输出示例
35793453939901
14225956593420
4332838845846
15400094813
数据规模及约定
T≤10000
1≤a,b≤107
题解
按照套路来,枚举最大公约数(以下有些过程的求和符号的上界应该是 ⌊min{n,m}x⌋,不失正确性,以下都写成 ⌊nx⌋)
n∑i=1m∑j=1f(gcd(i,j))=n∑g=1f(g)n∑i=1m∑j=1[gcd(i,j)=g]=n∑g=1f(g)⌊ng⌋∑i=1⌊mg⌋∑j=1∑d|i,d|jμ(d)=n∑g=1f(g)⌊ng⌋∑d=1μ(d)⌊ngd⌋⌊mgd⌋
继续按套路来,令 T=gd
=n∑T=1⌊nT⌋⌊mT⌋∑d|Tμ(d)f(Td)
令 g(T)=∑d|Tμ(d)f(Td),那么现在的问题就是快速求出 g(T) 的前缀和。
n∑T=1g(T)=n∑T=1∑d|Tμ(d)f(Td)=n∑d=1μ(d)⌊nd⌋∑T=1
f(T),μ(d) 的前缀和可以线性筛得到,那么我们只需要预处理 g(T) 的前 n23 项就可以杜教筛了。这个做法总时间复杂度 O(Tn23),在某些 OJ 上能过……
然而 BZOJ 上是不行的,所以需要想办法线性筛出 g(T)。
接下来是不知道怎么想出来的部分。
令 T=∏ki=1paii,由于 μ(d) 不为 0 的 d 的每个质因子都只有 1 次幂,所以有如下结论:
- 若 ∃i≠j,ai≠aj,那么 g(T)=0,这个可以用 (1−1)k 的二项式展开证明:把贡献分为 d 是否含有所有的指数最大的质因子,两类贡献都是形如 k′ 个数中取偶数个数的方案减去 k′ 个数中取奇数个数的方案,它们都是 0;
- 否则 g(T)=(−1)k+1,如果有偶数个不同的质因子,会在“全选”的时候少加 1,因为 f 值变小了 1(注意这里是所有 ai 都相等的情况),同理奇数个不同质因子会在“全选”时少减 1。
那么就可以线性筛出 g(T) 了,复杂度降为 O(T√n)。
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)
const int BufferSize = 1 << 16;
char buffer[BufferSize], *Head, *Tail;
inline char Getchar() {
if(Head == Tail) {
int l = fread(buffer, 1, BufferSize, stdin);
Tail = (Head = buffer) + l;
}
return *Head++;
}
int read() {
int x = 0, f = 1; char c = Getchar();
while(!isdigit(c)){ if(c == '-') f = -1; c = Getchar(); }
while(isdigit(c)){ x = x * 10 + c - '0'; c = Getchar(); }
return x * f;
}
#define maxn 10000010
#define LL long long
bool vis[maxn];
int prime[maxn], cp, mu[maxn], f[maxn], g[maxn], smu[maxn], nos[maxn];
LL sf[maxn], sg[maxn];
void init() {
int n = maxn - 1;
mu[1] = smu[1] = 1; f[1] = sf[1] = g[1] = sg[1] = 0;
rep(i, 2, n) {
if(!vis[i]) prime[++cp] = i, mu[i] = -1, f[i] = 1, nos[i] = 1, g[i] = 1;
for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
vis[i*prime[j]] = 1;
if(i % prime[j] == 0) {
mu[i*prime[j]] = 0;
f[i*prime[j]] = max(f[nos[i]], f[i/nos[i]] + 1);
if(nos[i] > 1) g[i*prime[j]] = f[nos[i]] == f[i/nos[i]] + 1 ? -g[nos[i]] : 0;
else g[i*prime[j]] = 1;
nos[i*prime[j]] = nos[i];
break;
}
mu[i*prime[j]] = -mu[i];
f[i*prime[j]] = max(f[i], 1);
g[i*prime[j]] = f[i] == 1 ? -g[i] : 0;
nos[i*prime[j]] = i;
}
smu[i] = smu[i-1] + mu[i];
sf[i] = sf[i-1] + f[i];
sg[i] = sg[i-1] + g[i];
}
return ;
}
int num[50], cntn;
void putLL(LL x) {
if(!x) return (void)puts("0");
cntn = 0;
while(x) num[cntn++] = x % 10, x /= 10;
dwn(i, cntn - 1, 0) putchar(num[i] + '0'); putchar('\n');
return ;
}
void work() {
int n = read(), m = read();
LL ans = 0;
for(int T = 1; T <= min(n, m); ) {
int r = min(min(n / (n / T), m / (m / T)), min(n, m));
ans += (LL)(n / T) * (m / T) * (sg[r] - sg[T-1]);
T = r + 1;
}
putLL(ans);
return ;
}
int main() {
init();
int T = read();
while(T--) work();
return 0;
}
再贴一下杜教筛的暴力
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)
const int BufferSize = 1 << 16;
char buffer[BufferSize], *Head, *Tail;
inline char Getchar() {
if(Head == Tail) {
int l = fread(buffer, 1, BufferSize, stdin);
Tail = (Head = buffer) + l;
}
return *Head++;
}
int read() {
int x = 0, f = 1; char c = Getchar();
while(!isdigit(c)){ if(c == '-') f = -1; c = Getchar(); }
while(isdigit(c)){ x = x * 10 + c - '0'; c = Getchar(); }
return x * f;
}
#define maxn 10000010
#define maxm 4000010
#define LL long long
bool vis[maxn];
int prime[maxn], cp, mu[maxn], f[maxn], smu[maxn], nos[maxn];
LL sf[maxn], g[maxm], sg[maxm];
void init() {
int n = maxn - 1;
mu[1] = smu[1] = 1; f[1] = sf[1] = 0;
rep(i, 2, n) {
if(!vis[i]) prime[++cp] = i, mu[i] = -1, f[i] = 1, nos[i] = 1;
for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
vis[i*prime[j]] = 1;
if(i % prime[j] == 0) {
mu[i*prime[j]] = 0;
f[i*prime[j]] = max(f[nos[i]], f[i/nos[i]] + 1);
nos[i*prime[j]] = nos[i];
break;
}
mu[i*prime[j]] = -mu[i];
f[i*prime[j]] = max(f[i], 1);
nos[i*prime[j]] = i;
}
smu[i] = smu[i-1] + mu[i];
sf[i] = sf[i-1] + f[i];
}
n = maxm - 1;
rep(i, 1, n) for(int j = 1; i * j <= n; j++) g[i*j] += mu[i] * f[j];
rep(i, 1, n) sg[i] = sg[i-1] + g[i];
return ;
}
const int HMOD = 1000037;
struct Hash {
int ToT, head[HMOD], nxt[maxn], key[maxn];
LL val[maxn];
Hash() { ToT = 0; memset(head, 0, sizeof(head)); }
LL Find(int x) {
int u = x % HMOD;
for(int e = head[u]; e; e = nxt[e]) if(key[e] == x) return val[e];
return -1;
}
void Insert(int x, LL v) {
int u = x % HMOD;
key[++ToT] = x; val[ToT] = v; nxt[ToT] = head[u]; head[u] = ToT;
return ;
}
} hG;
LL G(int n) {
if(n < maxm) return sg[n];
LL getg = hG.Find(n);
if(getg >= 0) return getg;
LL ans = 0;
for(int d = 1; d <= n; ) {
int r = min(n / (n / d), n);
ans += sf[n/d] * (smu[r] - smu[d-1]);
d = r + 1;
}
hG.Insert(n, ans);
return ans;
}
int num[50], cntn;
void putLL(LL x) {
if(!x) return (void)puts("0");
cntn = 0;
while(x) num[cntn++] = x % 10, x /= 10;
dwn(i, cntn - 1, 0) putchar(num[i] + '0'); putchar('\n');
return ;
}
void work() {
int n = read(), m = read();
LL ans = 0;
for(int T = 1; T <= min(n, m); ) {
int r = min(min(n / (n / T), m / (m / T)), min(n, m));
ans += (LL)(n / T) * (m / T) * (G(r) - G(T - 1));
T = r + 1;
}
putLL(ans);
return ;
}
int main() {
init();
int T = read();
while(T--) work();
return 0;
}
可以看出卡常的痕迹。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 你所不知道的 C/C++ 宏知识
· 聊一聊 操作系统蓝屏 c0000102 的故障分析
· SQL Server 内存占用高分析
· .NET Core GC计划阶段(plan_phase)底层原理浅谈
· .NET开发智能桌面机器人:用.NET IoT库编写驱动控制两个屏幕
· 我干了两个月的大项目,开源了!
· 推荐一款非常好用的在线 SSH 管理工具
· 千万级的大表,如何做性能调优?
· 聊一聊 操作系统蓝屏 c0000102 的故障分析
· .NET周刊【1月第1期 2025-01-05】