题意:

计算:

\[\sum\limits_{a = 1}^{m}\sum\limits_{b = 1}^{n} \frac{\varphi(ab)}{\varphi(a)\varphi(b)} (\bmod p)
\]

思路:

考虑算术基本定理和\(\varphi(x)\)函数积性将式子化简:

令\(a = p_1^{t_1}p_2^{t_2} \cdots p_n^{t_n}\),\(b = p_1^{q_1}p_2^{q_2} \cdots p_n^{q_n}\)。

那么原式有:

\[\begin{eqnarray*}
\frac{\varphi(ab)}{\varphi(a)\varphi(b)} (\bmod p) = \frac{\varphi(p_1^{t_1 + q_1} \cdots p_n^{t_n + q_n})}{\varphi(p_1^{t_1} \cdots \varphi(p_n^{t_n})) \cdot \varphi(p_1^{q_1} \cdots p_n^{q_n})} (\bmod p)
\end{eqnarray*}
\]

我们单独考虑一下\(p_1\),那么有:

\[\begin{eqnarray*}
\frac{\varphi(p_1^{t_1 + q_1})}{\varphi(p_1^{t_1}) \cdot \varphi(p_1^{q_1})} = \frac{p_1^{t_1 + q_1} \cdot (1 - \frac{1}{p_1})} {p_1^{t_1} (1 - \frac{1}{p_1})\cdot p_1^{q_1}(1 - \frac{1}{p_1})}
\end{eqnarray*}
\]

我们令\(t_1 < p_1\),即\(p_1^{t_1}是gcd(a, b)\)的一部分,那么约分之后有:

\[\begin{eqnarray*}
\frac{p_1^{t_1}}{p_1^{t_1} (1 - \frac{1}{p_1})}
\end{eqnarray*}
\]

我们再同理考虑\(p_1 \cdots p_n\),我们发现分子刚好是\(gcd(a, b)\), 而分母是\(\varphi(gcd(a, b))\),即:

\[\begin{eqnarray*}
\frac{\varphi(ab)}{\varphi(a)\varphi(b)} (\bmod p) &=& \frac{\varphi(p_1^{t_1 + q_1} \cdots p_n^{t_n + q_n})}{\varphi(p_1^{t_1} \cdots \varphi(p_n^{t_n})) \cdot \varphi(p_1^{q_1} \cdots p_n^{q_n})} (\bmod p) \\
&=& \frac{gcd(a, b)}{\varphi(gcd(a, b))}
\end{eqnarray*}
\]

所以现在我们的问题转化成了求解:

\[\begin{eqnarray*}
\sum\limits_{a = 1}^{m}\sum\limits_{b = 1}^{n} \frac{gcd(a, b)}{\varphi(gcd(a, b))} (\bmod p)
\end{eqnarray*}
\]

令\(gcd(a, b) = d\),并且令\(n <= m\),有:

\[\begin{eqnarray*}
\sum\limits_{a = 1}^{m} \sum\limits_{b = 1}^{n} \frac{d}{\varphi(d)} = \sum\limits_{d = 1}^{n} d \cdot \varphi(d)^{-1} \sum\limits_{a = 1}^{n} \sum\limits_{b = 1}^{m} [gcd(a, b) == d] \\
\end{eqnarray*}
\]

我们令:

\[\begin{eqnarray*}
f(d) &=& \sum\limits_{a = 1}^{n} \sum\limits_{b = 1}^{m} [gcd(a, b) == d] \\
g(d) &=& \sum\limits_{d|x}f(x) \\
&=& \sum\limits_{a = 1}^{n} \sum\limits_{b = 1}^{m} [d | gcd(a, b)] \\
&=& \sum\limits_{a = 1}^{n/d}\sum\limits_{b = 1}^{m/d} [1 | gcd(a, b)] \\
&=& \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor
\end{eqnarray*}
\]

进行莫比乌斯反演,有:

\[\begin{eqnarray*}
f(d) &=& \sum\limits_{d|x} \mu(\frac{x}{d}) g(d) \\
&=& \sum\limits_{d|x} \mu(\frac{x}{d}) \cdot \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \\
&=& \sum\limits_{x = 1}^{n/d} \mu(x) \cdot \lfloor \frac{n}{xd} \rfloor \lfloor \frac{m}{xd} \rfloor \\
\end{eqnarray*}
\]

所以,原式为:

\[\begin{eqnarray*}
\sum\limits_{i = 1}^{n} i \cdot \varphi(i)^{-1} \sum\limits_{d = 1}^{n|i} \mu(d) \lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor
\end{eqnarray*}
\]

预处理逆元,\(\varphi()\)函数,\(\mu()\)函数,然后直接算即可。

复杂度为\(\sum\limits_{i = 1}^{n} \sqrt{(i)}\)

#include <bits/stdc++.h>
using namespace std; #define ll long long
#define N 1000010
ll p;
int n, m;
int prime[N], mu[N];
int phi[N], inv[N], g[N];
bool check[N]; void init()
{
memset(check, 0, sizeof check);
prime[0] = 0;
phi[1] = 1;
mu[1] = 1;
for (int i = 2; i < N; ++i)
{
if (!check[i])
{
prime[++prime[0]] = i;
phi[i] = i - 1;
mu[i] = -1;
}
for (int j = 1; j <= prime[0]; ++j)
{
if (1ll * i * prime[j] >= N)
break;
check[i * prime[j]] = 1;
if (i % prime[j] == 0)
{
phi[i * prime[j]] = phi[i] * prime[j];
mu[i * prime[j]] = 0;
break;
}
else
{
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
mu[i * prime[j]] = -mu[i];
}
}
}
} void work()
{
inv[1] = 1;
for (int i = 2; i <= n; ++i)
inv[i] = 1ll * inv[p % i] * (p - p / i) % p;
for (int i = 1; i <= n; ++i)
g[i] = 1ll * i * inv[phi[i]] % p;
} ll get(int n, int m)
{
ll res = 0;
for (int i = 1; i <= n; ++i)
res = (res + 1ll * mu[i] * (n / i) * (m / i)) % p;
return res;
} int main()
{
init();
int T; cin >> T;
while (T--)
{
scanf("%d %d %lld\n", &n, &m, &p);
if (n > m) swap(n, m);
work();
ll res = 0;
for (int i = 1; i <= n; ++i)
res = (res + g[i] * get(n / i, m / i)) % p;
printf("%lld\n", res);
}
return 0;
}

最新文章

  1. 2016年8月ios面试问题总结
  2. 十一个行为模式之观察者模式(Observer Pattern)
  3. Maven 常用的命令
  4. js中遍历出查询后的listmodel(下拉框系列)
  5. Vector和Stack源码分析/List集合的总结
  6. 关于webapi post 使用多个参数的间接用法
  7. ImageLoader实现图片异步加载
  8. (转)深入探讨在集群环境中使用 EhCache 缓存系统
  9. 【转】Android与JNI(二) -- 不错
  10. Eclipse vim插件安装使用
  11. Codeforces Round #277.5 (Div. 2)A——SwapSort
  12. C#:判断当前程序是否通过管理员运行
  13. LeetCode题目----求中位数---标签:Array
  14. 多线程之Java线程阻塞与唤醒
  15. DOM4j 修改和删除
  16. Response.AddHeader小结
  17. [转]客户端js判断文件类型和文件大小即限制上传大小
  18. spring boot配置druid数据源和监控配置
  19. 黑帆第四季/全集Black Sails迅雷下载
  20. opengl redbook Fatal Error in XXXXXX pixel format with necessary capabilities not found

热门文章

  1. node爬虫(转)
  2. tomcat启动报错:serializer.jar (系统找不到指定的文件。)
  3. Hexo - 把word转成markdown
  4. Jmeter与Jenkins结合进行Web接口测试
  5. codeforces 586B/C
  6. 【紫书】Play on Words UVA - 10129 欧拉回路
  7. 安装XP时BIOS的设置(ahci ide)
  8. Buffer和cahce的区别
  9. 最大生成树——LCA
  10. 洛谷P3966 单词 [TJOI2013] AC自动机