hihocoder #1456 : Rikka with Lattice(杜教筛)

题意

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

答案对于998244353取模... (n,m5109)

题解

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

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

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

皮克定理:

S=a+b21

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

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

我们考虑一个ab的矩形,以它对顶点为端点的三角形,只当ab时存在四个解.

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

简单的证明:

代入皮克定理,可知

当且仅当格点三角形面积为S=0+321=12时才计入答案

我们可以用叉积的形式表示这个面积. S=12|AB||BC|sinθ=12|AB×AC|.(高中必修内容QAQ)

![image](http://images.cnblogs.com/cnblogs_com/zjp-shadow/1056673/t_sol.png)
我们令之前那个矩形的一条对角线为三角形的一条边,

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

我们令另外一条边为(i,j),然后三角形面积就是12|(a,b)×(i,j)|=12|ajbi|.

这个为12所以|ajbi|=1. ajbi=±1

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

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

所以±1各有一解,换个对角线又有对称的一组解.所以最后总共22=4组解.

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

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

ans=i=1nj=1m4(ni)(mj)[ij]

=4i=1nj=1m(ni)(mj)x|gcd(i,j)μ(x)

=4x=1min(n,m)μ(x)i=1nxj=1mxnmmixnjx+ijx2

S(i)=i=1ni=n(n+1)2.

sum1=nmx=1min(n,m)μ(x)nmx2

sum2=mx=1min(n,m)(μ(x)x)S(nx)mx

sum3=nx=1min(n,m)(μ(x)x)nxS(mx)

sum4=x=1min(n,m)(μ(x)x2)S(nx)S(mx)

ans=sum1sum2sum3+sum4

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

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

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

i=1nμ(i)之前那篇博客杜教筛小结中有介绍.

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

id(x)=x,mx(x)=μ(i)i. (然后之后默认把第一个字母大写记作前缀和比如Id(x)=i=1xid(i)=x(x+1)2)

所以就有

mxid(n)

=d|nμ(d)dnd=d|nμ(d)n

=nd|nμ(d)=n[n=1]=ϵ

代入之前的套路式子就有

1Mx(n)=i=2niMx(ni)

然后就可以尝试推出i=1nμ(i)ii.

这个也不麻烦QAQ...

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

有些地方有点细节51095109=2.51019会爆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; }

__EOF__

本文作者zjp_shadow
本文链接https://www.cnblogs.com/zjp-shadow/p/8492876.html
关于博主:评论和私信会在第一时间回复。或者直接私信我。
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   zjp_shadow  阅读(594)  评论(0编辑  收藏  举报
编辑推荐:
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
阅读排行:
· 地球OL攻略 —— 某应届生求职总结
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】
点击右上角即可分享
微信分享提示