0%

ECNU3618 数三角形

题目链接

题意:求在 $n\cdot n$ (包括边界)整点三角形有多少个。由于答案可能很大,对于 $10^9+7$ 取模。($1 \leq n \leq 2 \cdot 10^6$)

前置知识

  • 莫比乌斯反演公式:
    $F(n)=\sum\limits_{d\mid n}f(d)\Rightarrow f(n)=\sum\limits_{d\mid n}\mu(d)F(\frac{n}{d})$
    $F(n)=\sum\limits_{n\mid d}f(d)\Rightarrow f(n)=\sum\limits_{n\mid d}\mu(\frac{d}{n})F(d)$

  • 一些预处理技巧(莫比乌斯函数的预处理等)

  • 一系列的推公式技巧

推导出原始公式

一看到这样的题,我最先想到的是暴力打个表再说。
打表的程序如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from itertools import combinations

def isTri(a, b, c):
return (a[1] - b[1]) * (c[0] - b[0]) != (c[1] - b[1]) * (a[0] - b[0])

def calc(n):
v, ret = [], 0
for i in range(n + 1):
v += [(i, j) for j in range(n + 1)]
for tri in combinations(v, 3):
ret += 1 if isTri(tri[0], tri[1], tri[2]) else 0
return ret

for i in range(1, 10):
print(calc(i))

输出:

1
2
3
4
5
6
7
8
9
4
76
516
2148
6768
17600
40120
82608
157252

然后我就在oeis找到了他的原始公式
(我不知道是怎么推出来的)

原始公式: $A(n)=\frac{n^{2}(n+1)^{2}(n+2)^{2}}{6}-2\cdot \sum\limits_{l=1}^{n}\sum\limits_{k=1}^{n}(n-k+1)(n-l+1)\gcd(k,l)$

简化公式

我们可以发现前面部分是比较好求的,所以我们令 $$B(n)=\sum\limits_{l=1}^{n-1}\sum\limits_{k=1}^{n-1}(n-k)(n-l)\gcd(k,l)$$
那么 $$A(n)=\frac{n^{2}(n+1)^{2}(n+2)^{2}}{6}-2\cdot B(n+1)$$
枚举所有的$\gcd=d$的情况,那么 $$B(n)=\sum\limits_{d=1}^{n-1}d\sum\limits_{l=1}^{n-1}\sum\limits_{k=1}^{n-1}(n-k)(n-l)[\gcd(k,l)==d]$$ ($[\gcd(k,l)==d]$ 表示当gcd(k,l)==d时为1,否则为0)
这时令$$f(d)=\sum\limits_{l=1}^{n-1}\sum\limits_{k=1}^{n-1}(n-k)(n-l)[\gcd(k,l)==d]$$
$$=\sum\limits_{l=1}^{\lfloor\frac{n-1}{d}\rfloor}\sum\limits_{k=1}^{\lfloor\frac{n-1}{d}\rfloor}(n-ld)(n-kd)[\gcd(k,l)==1]$$ (假设n是个常量)

再令$$F(d)=\sum\limits_{l=1}^{n-1}\sum\limits_{k=1}^{n-1}(n-k)(n-l)[\gcd(k,l)为d的倍数]$$

所以$$F(d)=\sum\limits_{d\mid n}f(n)$$
易知$$F(d)=\sum\limits_{l=1}^{\lfloor\frac{n-1}{d}\rfloor}\sum\limits_{k=1}^{\lfloor\frac{n-1}{d}\rfloor}(n-ld)(n-kd)$$
$$=n^2\sum\limits_{l=1}^{\lfloor\frac{n-1}{d}\rfloor}\sum\limits_{k=1}^{\lfloor\frac{n-1}{d}\rfloor}1-dn\sum\limits_{l=1}^{\lfloor\frac{n-1}{d}\rfloor}\sum\limits_{k=1}^{\lfloor\frac{n-1}{d}\rfloor}(l+r)+d^2\sum\limits_{l=1}^{\lfloor\frac{n-1}{d}\rfloor}l\sum\limits_{k=1}^{\lfloor\frac{n-1}{d}\rfloor}k$$
$$=\lfloor\frac{n-1}{d}\rfloor^2n^2-(\lfloor\frac{n-1}{d}\rfloor^3+\lfloor\frac{n-1}{d}\rfloor^2)dn+\frac{\lfloor\frac{n-1}{d}\rfloor^2(\lfloor\frac{n-1}{d}\rfloor+1)^2}{4} d^2$$
由莫比乌斯反演公式得$$f(d)=\sum\limits_{d\mid m, m\lt n}\mu(\frac{m}{d})F(m)$$
$$=\sum\limits_{i=1}^{\lfloor\frac{n-1}{d}\rfloor}\mu(i)F(id)$$
则$$B(n)=\sum\limits_{d=1}^{n-1}d\sum\limits_{i=1}^{\lfloor\frac{n-1}{d}\rfloor}\mu(i)F(id)$$

在可以预处理前n个$\mu(n)$和前n个$F(n)$的情况下,复杂度为$O(n+\sum\limits_{i=1}^{n}\lfloor\frac{n}{i}\rfloor)\leq O(33333333)$在合理的复杂度范围内。

小技巧

在处理$(a\cdot b\cdot …)%m$时粗心的小伙伴总会因为太复杂而出错。
这样写基本上就不会出错了。
虽然可能会比较慢,但一般不会卡这个常吧。

1
2
3
4
5
6
7
8
9
10
11
long long mod(long long x, long long p) {
return ((x % p) + p) % p;
}

long long mmod() {
return 1ll;
}
template<typename T, typename ...U>
long long mmod(const T & head, const U &... tail) {
return mod((head % MOD) * mmod(tail...), MOD);
}

使用方法long long res = mmod(a, b, c, d);

最终代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include<bits/stdc++.h>
using namespace std;

const long long MOD = 1e9 + 7;
const int MAXN = 2e6 + 5;
bool isnpri[MAXN];
vector<int>pri;
int mu[MAXN];
void mobius() {
mu[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!isnpri[i]) {
pri.push_back(i);
mu[i] = -1;
}
for (int j = 0; j < int(pri.size()); j++) {
const int cur = i * pri[j];
if (cur >= MAXN) {
break;
}
isnpri[cur] = true;
if (i % pri[j] == 0) {
mu[cur] = 0;
break;
} else {
mu[cur] = -mu[i];
}
}
}
}

long long mod(long long x, long long p) {
return ((x % p) + p) % p;
}

long long mmod() {
return 1ll;
}
template<typename T, typename ...U>
long long mmod(const T & head, const U &... tail) {
return mod((head % MOD) * mmod(tail...), MOD);
}

long long F[MAXN];
long long n;
void init() {
for(int i = 1; i <= n; i++) {
long long cur = (n - 1) / i;
F[i] = mod(mmod(cur, cur, n, n) - mmod(cur * cur * cur + cur * cur, i, n) + mmod(cur, cur, cur + 1, cur + 1, 250000002ll, i, i), MOD);
}
}

signed main() {
mobius();
cin >> n;
long long sub = mmod(n, n, n + 1, n + 1, n + 2, n + 2, 166666668ll);

n++;
init();
long long ans = 0;
for (int d = 1; d <= n; d++) {
long long tmp = 0;
for (int j = 1; j <= n / d; j++) {
tmp = (tmp + F[j * d] * mu[j]) % MOD;
}
ans = (ans + tmp * d) % MOD;
}
cout << mod(sub - ans * 2, MOD) << endl;
}

回到开头