题意

给你一个\(n*m\)方格图,统计上面有多少个格点三角形,除了三个顶点,不覆盖其他的格点(包括边和内部).

答案对于\(998244353\)取模... (\(n,m \le 5 * 10^9\))

题解

这个题十分的巧妙... 集训时是大佬ztzshiwo出的..

据他所说,是不那么杜教筛的杜教筛QAQ

考试时候提示了一个皮克定理...

皮克定理:

\[S=a+\frac{b}{2}-1
\]

\(S\)为格点多边形面积,\(a\)为多边形内部点数,\(b\)为多边形边上点数.

然而我还是只会暴力,正解是真的太神了啊QAQ.

我们考虑一个\(a*b\)的矩形,以它对顶点为端点的三角形,只当\(a \bot b\)时存在四个解.

这个我是听wearry证明的qwq 十分巧妙

简单的证明:

代入皮克定理,可知

当且仅当格点三角形面积为\(S=0+\frac{3}{2}-1=\frac{1}{2}\)时才计入答案

我们可以用叉积的形式表示这个面积. \(S=\frac{1}{2} |AB| |BC| \sin \theta =\frac{1}{2} |\vec {AB} \times \vec{AC}|\).(高中必修内容QAQ)

![image](http://images.cnblogs.com/cnblogs_com/zjp-shadow/1056673/t_sol.png)

我们令之前那个矩形的一条对角线为三角形的一条边,

令左下角为向量出发点,也就是其中一个顶点,然后这条边的向量坐标表达就为\((a,b)\).

我们令另外一条边为\((i,j)\),然后三角形面积就是\(\frac{1}{2}|(a,b) \times (i,j)| = \frac{1}{2}|aj - bi|\).

这个为\(\frac{1}{2}\)所以\(|aj-bi|=1\). \(\therefore aj-bi= \pm 1\)

我们是要求解\((i,j)\) 所以不难发现这是一个扩欧的形式,当且仅当\(a \bot b\)时有整数解.

又\(\because 0 < i < a, 0 < j < b\).. \(\therefore\)可以通过扩欧的相邻解确定在这个区域仅一解.

所以\(\pm 1\)各有一解,换个对角线又有对称的一组解.所以最后总共\(2*2=4\)组解.

所以我们要求的就是原图中每个矩形的贡献就行了...

此处\(n\),\(m\)都是要减一的... (至于为啥...手推就知道了QAQ)

\[\displaystyle \mathrm {ans}=\sum_{i=1}^{n} \sum _{j=1}^{m} 4 \cdot (n-i)(m-j) [i \bot j]
\]

\[\displaystyle =4\sum_{i=1}^{n} \sum_{j=1}^{m} (n-i)(m-j) \sum _ {x|\gcd (i,j)} \mu(x)
\]

\[=\displaystyle 4 \sum_{x=1}^{min(n,m)} \mu(x) \sum_{i=1}^{\lfloor \frac{n}{x} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{x} \rfloor} nm - mix - njx + ijx^2
\]

令$$\displaystyle S(i)=\sum_{i=1}^{n} i = \frac{n(n+1)}{2}$$.

\[\displaystyle \mathrm{sum1}=nm \sum_{x=1}^{min(n,m)} \mu(x) \lfloor \frac{nm}{x^2} \rfloor
\]

\[\displaystyle \mathrm{sum2}=m \sum_{x=1}^{min(n,m)} (\mu(x) \cdot x) S(\lfloor \frac{n}{x} \rfloor) \lfloor \frac{m}{x} \rfloor
\]

\[\displaystyle \mathrm{sum3}=n \sum_{x=1}^{min(n,m)} (\mu(x) \cdot x) \lfloor \frac{n}{x} \rfloor S(\lfloor \frac{m}{x} \rfloor)
\]

\[\displaystyle \mathrm{sum4}=\sum_{x=1}^{min(n,m)} (\mu(x) \cdot x^2) S(\lfloor \frac{n}{x} \rfloor) S(\lfloor \frac{m}{x} \rfloor)
\]

\[\mathrm{ans}=\mathrm{sum1}-\mathrm{sum2}-\mathrm{sum3}+\mathrm{sum4}
\]

这个用根号分块就能做到\(\Theta (n+\sqrt {n})\)复杂度啦... 具体推导证明看我的一篇博客线性筛与莫比乌斯反演.

然而这并不能满分...fuck

所以就有杜教筛卡了30分.

\(\displaystyle \sum _{i=1}^{n} \mu(i)\)之前那篇博客杜教筛小结中有介绍.

然后就介绍另外两个套路求的东西吧..

令\(\displaystyle id(x)=x, mx(x)= \mu(i) i\). (然后之后默认把第一个字母大写记作前缀和比如\(\displaystyle Id(x)=\sum_{i=1}^{x} id(i) = \frac{x(x+1)}{2}\))

所以就有

\[mx * id (n)
\]

\[\displaystyle = \sum _{d|n} \mu(d) \cdot d \cdot \frac{n}{d} =\displaystyle \sum _{d|n} \mu(d) \cdot n
\]

\[=\displaystyle n \sum_{d|n} \mu(d) = n \cdot [n=1] = \epsilon
\]

代入之前的套路式子就有

\[\displaystyle 1 - Mx(n) = \sum _{i=2}^{n} i \cdot Mx(\lfloor \frac{n}{i} \rfloor)
\]

然后就可以尝试推出\(\displaystyle \sum _{i=1}^{n} \mu(i) \cdot i \cdot i\).

这个也不麻烦QAQ...

然后本人比较懒 就直接用c++11 中的unordered_map了(这个基于哈希算法)

有些地方有点细节\(5*10^9 * 5*10^9 = 2.5 * 10^{19}\)会爆long long所以很多地方都要记得先取模!!!

代码

#include <bits/stdc++.h>
#define For(i, l, r) for(register ll i = (l), _end_ = (ll)(r); i <= _end_; ++i)
#define Fordown(i, r, l) for(register ll i = (r), _end_ = (ll)(l); i >= _end_; --i)
#define Set(a, v) memset(a, v, sizeof(a))
using namespace std; typedef long long ll;
inline ll read() {
ll x = 0, fh = 1; char ch = getchar();
for (; !isdigit(ch); ch = getchar() ) if (ch == '-') fh = -1;
for (; isdigit(ch); ch = getchar() ) x = (x<<1) + (x<<3) + (ch ^ '0');
return x * fh;
} void File() {
#ifdef zjp_shadow
freopen ("1456.in", "r", stdin);
freopen ("1456.out", "w", stdout);
#endif
} const ll Mod = 998244353; ll n, m; const int N = 1e7 + 1e3;
int prime[N], cnt = 0;
int Limit = N - 1e3; ll mux[N], muxx[N], mu[N];
bitset<N> is_prime; void Init(int maxn) {
int res;
mu[1] = 1;
is_prime.set(); is_prime[0] = is_prime[1] = false;
For (i, 2, maxn) {
if (is_prime[i]) { prime[++cnt] = i; mu[i] = -1; }
For (j, 1, cnt) {
res = prime[j] * i;
if (res > maxn) break;
is_prime[res] = false;
if (i % prime[j]) mu[res] = -mu[i];
else { mu[res] = 0; break ; }
}
}
For (i, 1, maxn) {
mux[i] = mux[i - 1] + 1ll * mu[i] * i % Mod;
mux[i] = (mux[i] % Mod + Mod) % Mod; muxx[i] = muxx[i - 1] + 1ll* mu[i] * i % Mod * i % Mod;
muxx[i] = (muxx[i] % Mod + Mod) % Mod; mu[i] += mu[i - 1];
mu[i] = (mu[i] % Mod + Mod) % Mod;
}
} ll fpm(ll x, ll power) { ll res = 1; x %= Mod; for (; power; power >>= 1, (x *= x) %= Mod) if (power & 1) (res *= x) %= Mod; return res; }
const ll inv2 = fpm(2, Mod - 2), inv6 = fpm(6, Mod - 2);
ll Sum(ll x) { x %= Mod; return (x + 1) * x % Mod * inv2 % Mod; } ll sum1, sum2, sum3, sum4;
ll Nextx; unordered_map<ll, ll> MU, MUX, MUXX;
ll mu_(ll x) {
if (x <= Limit) return mu[x];
if (MU.count(x)) return MU[x];
ll res = 1, Nextx;
For (i, 2, x) { Nextx = x / (x / i); (res += Mod - (Nextx - i + 1) * mu_(x / i) % Mod) %= Mod; i = Nextx; }
return (MU[x] = res);
} inline ll Sum1(ll x) { x %= Mod; return x * (x + 1) % Mod * inv2 % Mod; }
ll mux_(ll x) {
if (x <= Limit) return mux[x];
if (MUX.count(x)) return MUX[x];
ll res = 1, Nextx;
For (i, 2, x) { Nextx = x / (x / i); (res += Mod - (Sum1(Nextx) - Sum1(i - 1) + Mod) * mux_(x / i) % Mod) %= Mod; i = Nextx; }
return (MUX[x] = res);
} inline ll Sum2(ll x) { x %= Mod ; return x * (x + 1) % Mod * (2 * x + 1) % Mod * inv6 % Mod; }
ll muxx_(ll x) {
if (x <= Limit) return muxx[x];
if (MUXX.count(x)) return MUXX[x];
ll res = 1, Nextx;
For (i, 2, x) { Nextx = x / (x / i); (res += Mod - (Sum2(Nextx) - Sum2(i - 1) + Mod) * muxx_(x / i) % Mod) %= Mod; i = Nextx; }
return (MUXX[x] = res);
} int main () {
File();
n = read(); m = read();
if (n > m) swap(n, m);
Init(Limit); For (x, 1, n) {
Nextx = min(n / (n / x), m / (m / x));
(sum1 += (mu_(Nextx) - mu_(x - 1)) * (n / x) % Mod * (m / x) % Mod * n % Mod * m % Mod) %= Mod;
(sum2 += (mux_(Nextx) - mux_(x - 1)) * Sum(n / x) % Mod * (m / x) % Mod) %= Mod;
(sum3 += (mux_(Nextx) - mux_(x - 1)) * Sum(m / x) % Mod * (n / x) % Mod) %= Mod;
(sum4 += (muxx_(Nextx) - muxx_(x - 1)) * Sum(n / x) % Mod * Sum(m / x) % Mod) %= Mod;
x = Nextx;
}
ll ans = sum1 - 1ll * m * sum2 % Mod - 1ll * n * sum3 % Mod + sum4;
ans = (ans % Mod + Mod) % Mod;
ans = ans * 4ll % Mod;
printf ("%lld\n", ans);
return 0;
}

最新文章

  1. JavaScript闭包(Closure)学习笔记
  2. spring功能总结
  3. 工作需求----表单select多选交互
  4. 【转载】Linux启动过程
  5. [NOIP2008] 提高组 洛谷P1155 双栈排序
  6. Linux配置SSH公钥认证与Jenkins远程登录进行自动发布
  7. C#基础篇03
  8. iOS 里面如何使用第三方应用程序打开自己的文件,调用wps其他应用打开当前应用里面的的ppt doc xls
  9. img 中的src的应用
  10. putty 使用 注意事项
  11. nginx的https环境如何配置
  12. 学习web的第二天
  13. Dynamics CRM 日常使用JS整理(二)
  14. JPA(Hibernate)
  15. Linux驱动之USB总线驱动程序框架简析
  16. 0.2:Game and Art
  17. elasticsearch6.7 05. Document APIs(8)Multi Get API
  18. 《大道至简》第一章读后感Java伪代码
  19. -bash: /etc/init.d/nginx: /bin/bash^M: bad interpreter: No such file or directory
  20. composer 安装yii2 The package is not available in a stable-enough version解决办法

热门文章

  1. PhpStorm的破解 汉化
  2. python学习:字典排序
  3. 洛谷P1171 售货员的难题【状压DP】
  4. angular-file-upload封装为指令+图片尺寸限制
  5. 加快compser install 和update的方法
  6. Qt 5.9.4 如何静态编译和部署?
  7. java网络编程(3)——UDP
  8. 在SpringBoot中存放session到Redis
  9. 基于Python的Flask的开发实战(第二节程序的基本结构)
  10. c#获取文件MD5算法