Loading

BZOJ3518 点组计数

莫比乌斯反演做法。

将横线竖线斜线分开考虑。这题有个很好的特点,就是以四个角落为线段端点所涵盖的的线的种类是最多的。那么我们先将枚举的线的一个端点固定在左上角,其坐标为 (0,0),然后再枚举另外一个端点 (i,j),那么这种线在本题中的贡献为 (gcd(i,j)1)×(ni+1)×(mj+1),为什么呢?因为这种线段在图像中有 (ni+1)(mj+1) 个位置可以存放,并且中间的点有 (gcd(i,j)1) 个位置可以选择,所以我们将本题转化为计算 i=1nj=1m(gcd(i,j)1)(ni+1)(mj+1)

我们将式子拆开,分别计算出

a=i=1nj=1m(gcd(i,j)1)(n+1)(m+1)

b=i=1nj=1mi(gcd(i,j)1)(m+1)

c=i=1nj=1mj(gcd(i,j)1)(n+1)

d=i=1nj=1m(gcd(i,j)1)ij

ans=abc+d

即可。

计算可以利用莫比乌斯反演,时间复杂度 O(n)

注意,因为我左上角坐标为 (0,0),故开始 n,m 都要减一。最后再加回去来计算横竖线产生的贡献。

代码:

#include <bits/stdc++.h>
#define rep(i, l, r) for (int i = l; i <= r; ++ i)
#define rrp(i, l, r) for (int i = r; i >= l; -- i)
#define eb emplace_back
#define int long long
using namespace std;
constexpr int N = 5e4 + 5, P = 1e9 + 7;
inline int rd ()
{
  int x = 0, f = 1;
  char ch = getchar ();
  while (! isdigit (ch))
  {
    if (ch == '-') f = - 1;
    ch = getchar ();
  }
  while (isdigit (ch))
  {
    x = (x << 1) + (x << 3) + (ch ^ 48);
    ch = getchar ();
  }
  return x * f;
}
int qpow (int x, int y)
{
  int ret = 1;
  for (; y; y >>= 1, x = x * x % P) if (y & 1) ret = ret * x % P;
  return ret;
}
int T;
int mu[N], pri[N], cnt;
int sum[N], s1[N], s2[N], s3[N], s4[N], s5[N];
bool ip[N];
int func (int a, int b, int t)
{
  int n = min (a, b);
  int ret = 0;
  int l1, l2, r1, r2;
  l1 = l2 = 1;
  while (l1 <= n && l2 <= n)
  {
    int k1 = a / l1;
    r1 = min (a / k1, n);
    int k2 = b / l2;
    r2 = min (b / k2, n);
    if (t == 1)
    (ret += (sum[min (r1, r2)] - sum[max (l1, l2) - 1]) * k1 * k2) %= P;
    if (t == 2)
    (ret += (s4[min (r1, r2)] - s4[max (l1, l2) - 1]) * (k1 * (k1 + 1) / 2 % P) % P * k2) %= P;
    if (t == 3)
    (ret += (s4[min (r1, r2)] - s4[max (l1, l2) - 1]) * (k2 * (k2 + 1) / 2 % P) % P * k1) %= P;
    if (t == 4)
    (ret += (s5[min (r1, r2)] - s5[max (l1, l2) - 1]) * (k1 * (k1 + 1) / 2 % P) % P * (k2 * (k2 + 1) / 2 % P)) %= P;
    if (r1 < r2) l1 = r1 + 1; else l2 = r2 + 1;
  }
  return ret;
}
signed main ()
{
  mu[1] = 1;
  rep (i, 2, N - 1)
  {
    if (! ip[i])
    {
      pri[++ cnt] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= cnt && i * pri[j] <= N - 1; ++ j)
    {
      ip[i * pri[j]] = 1;
      if (i % pri[j] == 0)
      {
        mu[i * pri[j]] = 0;
        break;
      }
      mu[i * pri[j]] = - mu[i];
    }
  }
  rep (i, 0, N - 1) sum[i] = sum[i - 1] + mu[i];
  rep (i, 1, N - 1)
  	s1[i] = (s1[i - 1] + max (0ll, i - 1)) % P,
  	s2[i] = (s2[i - 1] + i * max (0ll, i - 1)) % P,
	s3[i] = (s3[i - 1] + i * i * max (0ll, i - 1)) % P,
	s4[i] = (s4[i - 1] + mu[i] * i) % P,
	s5[i] = (s5[i - 1] + mu[i] * i * i) % P;
  int n = rd () - 1, m = rd () - 1;
  if (n > m) swap (n, m);
  int l1 = 1, l2 = 1, ret = 0;
  while (l1 <= n && l2 <= n)
  {
    int r1 = n / (n / l1);
    int r2 = m / (m / l2);
    int l = max (l1, l2), r = min (r1, r2);
    (ret +=
	func (n / l, m / r, 1) * (m + 1) % P * (n + 1) % P * (s1[r] - s1[l - 1]) -
	func (n / l, m / r, 2) * (m + 1) % P * (s2[r] - s2[l - 1]) - 
	func (n / l, m / r, 3) * (n + 1) % P * (s2[r] - s2[l - 1]) + 
	func (n / l, m / r, 4) * (s3[r] - s3[l - 1])
	);
    if (r1 < r2) l1 = r1 + 1; else l2 = r2 + 1;
  }
  ret <<= 1;
  ++ n, ++ m;
  rep (i, 3, n) (ret += (i - 2) * m * (n - i + 1)) %= P;
  rep (i, 3, m) (ret += (i - 2) * n * (m - i + 1)) %= P;
  printf ("%lld\n", (ret + P) % P);
}
 

作者:lalaouye

出处:https://www.cnblogs.com/lalaouyehome/p/18264277

版权:本作品采用「114514」许可协议进行许可。

posted @   lalaouye  阅读(4)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· Manus重磅发布:全球首款通用AI代理技术深度解析与实战指南
· 被坑几百块钱后,我竟然真的恢复了删除的微信聊天记录!
· 没有Manus邀请码?试试免邀请码的MGX或者开源的OpenManus吧
· 园子的第一款AI主题卫衣上架——"HELLO! HOW CAN I ASSIST YOU TODAY
· 【自荐】一款简洁、开源的在线白板工具 Drawnix
more_horiz
keyboard_arrow_up dark_mode palette
选择主题
点击右上角即可分享
微信分享提示