博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
luogu4240 毒瘤之神的考验(毒瘤乌斯反演)
阅读量:4355 次
发布时间:2019-06-07

本文共 4806 字,大约阅读时间需要 16 分钟。

题意:求出\(\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)\),对998244353取模

多组数据,\(T\le 10^4,n,m\le 10^5\)

前置知识:\(\varphi(ij)=\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))}\)

证明:我是口胡呢还是好好证呢还是口胡吧

按照欧拉函数的计算式展开,会发现,左边是\(ij\prod_{p|i \mathrm{\color{green}{or}}p|j}\frac{p-1}p\)

右边是\(\frac{i\prod_{p|i}\frac{p-1}pj\prod_{p|j}\frac{p-1}p\gcd(i,j)}{\gcd(i,j)\prod_{p|i\mathrm{\color{green}{and}}p|j}\frac{p-1}p}\)

显然,根据容斥原理,两边是相等的

然后推式子

\(\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)\)

\(=\sum_{i=1}^n\sum_{j=1}^m\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(ij))}\)

\(=\sum_{p=1}^n\frac p{\varphi(p)}\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)[\gcd(i,j)=p]\)

\(=\sum_{p=1}^n\frac p{\varphi(p)}\sum_{i=1}^{n/p}\sum_{j=1}^{m/p}\varphi(ip)\varphi(jp)[\gcd(i,j)=1]\)

\(=\sum_{p=1}^n\frac p{\varphi(p)}\sum_{i=1}^{n/p}\sum_{j=1}^{m/p}\varphi(ip)\varphi(jp)\sum_{d|i,d|j}\mu(d)\)

\(=\sum_{p=1}^n\frac p{\varphi(p)}\sum_{d=1}^n\mu(d)\sum_{i=1}^{n/dp}\sum_{j=1}^{m/dp}\varphi(idp)\varphi(jdp)\)

\(=\sum_{q=1}^n\sum_{p|q}\frac{p\mu(\frac qp)}{\varphi(p)}\sum_{i=1}^{n/q}\sum_{j=1}^{m/q}\varphi(iq)\varphi(jq)\)

\(=\sum_{q=1}^n\left(\sum_{p|q}\frac{p\mu(\frac qp)}{\varphi(p)}\right)\left(\sum_{i=1}^{n/q}\varphi(iq)\right)\left(\sum_{i=1}^{m/q}\varphi(iq)\right)\)

前面这一部分好处理--\(O(n\log n)\)枚举倍数。后面?按照套路?数论分块?怎么分????

观察了你谷题解后,终于懂了

\(sum(q)=\sum_{p|q}\frac{p\mu(\frac qp)}{\varphi(p)}\),显然可以在\(O(n\log n)\)的时间复杂度内处理出来。

\(g(x,y)=\sum_{i=1}^x\varphi(iy)\),显然有递推式\(g(x,y)=g(x-1,y)+\varphi(xy)\)

由于\(xy<=n\),对于每个\(x\),有\(\frac nx\)的数值,我们可以通过动态申请内存,在\(O(n\log n)\)的时间复杂度和空间复杂度内求出\(g\)数组。

\(T(n,a,b)=\sum_{q=1}^n\left(\sum_{p|q}\frac{p\mu(\frac qp)}{\varphi(p)}\right)\left(\sum_{i=1}^{a}\varphi(iq)\right)\left(\sum_{i=1}^{b}\varphi(iq)\right)=\sum_{q=1}^nsum(q)g(a,q)g(b,q)\)

显然T的递推式为\(T(n,a,b)=T(n-1,a,b)+sum(n)g(a,n)g(b,n)\)

根据数论分块那套理论,对于一个\(n/q\)\(m/q\)相同的\(q\)的区间,当\(n/q=a,m/q=b\)时,这一区间的\(ans=T(r,a,b)-T(l-1,a,b)\),r和l是这一区间内的最大值和最小值

我们考虑预处理\(n*B*B\)范围的答案,B是我们钦定的一个数字,T数组开的空间复杂度为\(O(nB^2)\)(实际上由于\(a*n,b*n\le 10^5\)的限制,应该开不到\(O(nB^2)\)

对于每次询问,我们只能在\(n/q\le B\)时候进行数论分块操作通过\(T\)数组计算答案,复杂度根据数论分块那套理论为\(O(\sqrt n)\)

对于\(n/q>B\)的部分,有\(q<n/B\),暴力枚举\(q\),通过\(g\)数组计算答案,这一部分单次计算的复杂度为\(O(n/B)\)

总复杂度为\(O(n\log n+nB^2+T(\sqrt n+n/B))\)。实测,B开到50左右跑的快一点,且内存占用超小。

下面是乱七八糟的代码= =

注意讲文明,new来的内存要主动回收垃圾

注意取模(这题如果写的复杂度没错的话不卡常,开#define int long long也是没问题的

#include 
#include
using namespace std;const int p = 998244353;const int b = 50;bool vis[100010];int prime[100010], tot, fuck = 100000;int mu[100010], phi[100010], invphi[100010];int sum[100010];int *g[100010], *t[100][100]; //注意这里t数组下标是[2][3][1]int qpow(int x, int y){ int res = 1; for (x %= p; y > 0; x = x * (long long)x % p, y >>= 1) if (y & 1) res = res * (long long)x % p; return res;}int main(){ //线性筛phi,mu,预处理前面的部分 phi[1] = mu[1] = invphi[1] = 1; for (int i = 2; i <= fuck; i++) { if (vis[i] == false) prime[++tot] = i, phi[i] = i + (mu[i] = -1); for (int j = 1; j <= tot && i * prime[j] <= fuck; j++) { vis[i * prime[j]] = true; if (i % prime[j] == 0) { phi[i * prime[j]] = phi[i] * prime[j]; break; } phi[i * prime[j]] = phi[i] * (prime[j] - 1); mu[i * prime[j]] = -mu[i]; } invphi[i] = qpow(phi[i], p - 2); if (phi[i] * (long long)invphi[i] % p != 1) { return -233; printf("cnm\n"); } } for (int pp = 1; pp <= fuck; pp++) for (int q = pp, d = 1; q <= fuck; q += pp, d++) sum[q] = (sum[q] + pp * (long long)invphi[pp] % p * mu[d]) % p, sum[q] += (sum[q] < 0 ? p : 0); //处理g数组 for (int i = 1; i <= fuck; i++) { g[i] = new int[(fuck / i) + 1], g[i][0] = 0; for (int j = 1, sb = fuck / i; j <= sb; j++) g[i][j] = (g[i][j - 1] + phi[i * j]) % p; } //处理t数组 注意有第一维<=第二维,因为下面我们强制n<=m了 for (int j = 1; j <= b; j++) for (int k = j; k <= b; k++) { int len = fuck / max(j, k); t[j][k] = new int[len + 1], t[j][k][0] = 0; for (int i = 1; i <= len; i++) t[j][k][i] = (t[j][k][i - 1] + sum[i] * (long long)g[i][j] % p * g[i][k] % p) % p; } //处理询问 int tat; scanf("%d", &tat); while (tat --> 0) { int n, m, res = 0; scanf("%d%d", &n, &m); if (n > m) swap(n, m); //对于n/q>b的部分,暴力,通过g数组和sum数组计算计算 for (int i = 1, sb = m / b; i <= sb; i++) res = (res + sum[i] * (long long)g[i][n / i] % p * g[i][m / i] % p) % p; //对于n/q
< 0 ? p : 0); } printf("%d\n", res); } //垃圾回收 for (int i = 1; i <= fuck; i++) delete []g[i], g[i] = 0; for (int i = 1; i <= b; i++) for (int j = i; j <= b; j++) delete[] t[i][j], t[i][j] = 0; return 0;}

转载于:https://www.cnblogs.com/oier/p/10307241.html

你可能感兴趣的文章
Bootstrap相关优质项目学习清单
查看>>
GS给客户单发包以及m_queGcWait(所有GC共享)
查看>>
Ulipad运行源码出错
查看>>
Windows Live Writer 网易博客配置
查看>>
逻辑运算符
查看>>
WPA2无线网络加密协议,部分受到影响产品以及对应补丁
查看>>
jvm
查看>>
安卓学习-界面-ui-ImageView
查看>>
利用键值对来找对应值的信息
查看>>
JNI全局对象,及原生线程JNIENV传递
查看>>
POJ 1159 回文LCS滚动数组优化
查看>>
ASP.Net MVC3 - The easier to run Unit Tests by moq #Reprinted#
查看>>
《Python从入门基础到实践》
查看>>
for循环
查看>>
BZOJ2132: 圈地计划
查看>>
PHP数据库连接mysql与mysqli的区别与用法
查看>>
char * 与char []探究理解
查看>>
QT窗体显示在屏幕中间位置
查看>>
emmet使用技巧
查看>>
RPC-Thrift(二)
查看>>