BZOJ3518 点组计数

莫比乌斯反演做法。

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

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

$$a=\sum _{i=1}^n\sum _{j=1}^m(gcd(i,j)-1)(n+1)(m+1)$$

$$b=\sum _{i=1}^n\sum _{j=1}^{m}i(gcd(i,j)-1)(m+1)$$

$$c=\sum _{i=1}^n\sum _{j=1}^{m}j(gcd(i,j)-1)(n+1)$$

$$d=\sum _{i=1}^n\sum _{j=1}^m(gcd(i,j)-1)ij$$

$$ans=a-b-c+d$$

即可。

计算可以利用莫比乌斯反演，时间复杂度 $\mathcal{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);
}