Project Euler 465

发表于 分类于 阅读次数： 本文字数： 3.1k 阅读时长 ≈ 3 分钟

Project Euler 465

题目

Polar polygons

The kernel of a polygon is defined by the set of points from which the entire polygon’s boundary is visible. We define a polar polygon as a polygon for which the origin is strictly contained inside its kernel.

For this problem, a polygon can have collinear consecutive vertices. However, a polygon still cannot have self-intersection and cannot have zero area.

For example, only the first of the following is a polar polygon (the kernels of the second, third, and fourth do not strictly contain the origin, and the fifth does not have a kernel at all):

Notice that the first polygon has three consecutive collinear vertices.

Let \(P(n)\) be the number of polar polygons such that the vertices \((x, y)\) have integer coordinates whose absolute values are not greater than \(n\).

Note that polygons should be counted as different if they have different set of edges, even if they enclose the same area. For example, the polygon with vertices \([(0,0),(0,3),(1,1),(3,0)]\) is distinct from the polygon with vertices \([(0,0),(0,3),(1,1),(3,0),(1,0)]\).

For example, \(P(1) = 131, P(2) = 1648531, P(3) = 1099461296175\) and \(P(343) \bmod 1000000007 = 937293740\).

Find \(P(7^{13}) \bmod 1000000007\).

解决方案

设点集来自方格\(\mathcal G_n=\{(x,y)\in\mathbb Z^2:|x|\le n,|y|\le n\}\setminus\{(0,0)\}.\)把任意候选多边形的顶点按极角升序连接。原点严格在核心内可等价转化为两个条件：

  1. 同一条过原点射线上至多取一个点。否则远点会被近点在同方向遮挡，边界不可能从原点全部可见。
  2. 所选方向不被任何过原点的闭半平面完全包含。等价说法是：按极角排序后，相邻方向夹角都严格小于 \(\pi\)。

在这两个条件下，按极角连接得到的就是一个无自交、面积非零且原点在核心内的极多边形；反过来极多边形也必满足这两个条件。因此问题变成：在每条方向射线上选 0 或 1 个点，并且选出的方向集合不落在任何闭半平面内。

把方向按互质向量 \((a,b)\)（\(\gcd(|a|,|b|)=1\)）分类。若该方向满足 \(\max(|a|,|b|)=d\)，则该方向上可选点数\(c(d)=\left\lfloor \dfrac n d\right\rfloor .\) 并且满足 \(\max(|a|,|b|)=d\) 的方向数是\(8\varphi(d).\)

定义\(\displaystyle A=\prod_{u\in\mathcal U}(c_u+1),N=\sum_{u\in\mathcal U}c_u,Q=\sum_{u\in\mathcal U}c_u^2.\)其中 \(\mathcal U\) 表示全部方向的索引集合。则有 \(\displaystyle A=\prod_{d=1}^n\left(\left\lfloor\frac nd\right\rfloor+1\right)^{8\varphi(d)}=M^8,\)其中\(\displaystyle M=\prod_{d=1}^n\left(\left\lfloor\frac nd\right\rfloor+1\right)^{\varphi(d)}.\)

同时可得\(\displaystyle N=(2n+1)^2-1=4n(n+1),Q=\sum_{d=1}^n 8\varphi(d)\left\lfloor\frac nd\right\rfloor^2=8S,S=\sum_{d=1}^n\varphi(d)\left\lfloor\frac nd\right\rfloor^2.\)

先数每个方向至多选一个点的非空集合，可得数量为\(A-1\)。接着减去落在某个闭半平面内的坏集合。用锚点计数：选一个具体点 \(P\) 作为半平面边界上的锚点（共 \(N\) 种）；其余点只能从与 \(P\) 对应的闭半平面内各方向选或不选，这一乘积恒为 \(\sqrt A\)（这是因为对向方向成对、半圈取一侧）。所以锚点计数总数是 \(N\sqrt A\)。

但坏集合中，含一对关于原点对向点的集合会被双计一次。双计总数为\(\dfrac Q2\)（每对对向方向贡献 \(c_u c_{-u}=c_u^2\)，再除以 \(2\) 去重）。于是\(P(n)=(A-1)-\left(N\sqrt A-\dfrac Q2\right)= A-N\sqrt A+\dfrac Q2-1.\)代入 \(A=M^8,\sqrt A=M^4,N=4n(n+1),Q=8S\)，得到

\[ P(n)=M^8-4n(n+1)M^4+4S-1. \]

记欧拉函数前缀和为\(\displaystyle\Phi(x)=\sum_{k=1}^x\varphi(k).\)那么我们考虑使用数论分块来计算。我们按 \(\left\lfloor \dfrac n l\right\rfloor\) 分块：若\(q=\left\lfloor\dfrac nl\right\rfloor,r=\left\lfloor\dfrac n q\right\rfloor,\)则在区间 \([l,r]\) 上商都等于 \(q\)，且权重\(\displaystyle \sum_{k=l}^r\varphi(k)=\Phi(r)-\Phi(l-1).\)于是

\[ M=\prod_{b} (q+1)^{\Phi(r)-\Phi(l-1)}, S=\sum_{b} q^2\bigl(\Phi(r)-\Phi(l-1)\bigr). \]

块数是 \(O(\sqrt n)\)。

利用欧拉函数的和恒等函数的Dirichlet 卷积：\(\displaystyle\sum_{d=1}^m \varphi(d)\left\lfloor\frac md\right\rfloor=\frac{m(m+1)}2,\)可得递推\(\displaystyle\Phi(m)=\frac{m(m+1)}2-\sum_{l=2}^m \Phi\left(\left\lfloor\frac ml\right\rfloor\right).\) 对相同商分组（\(l\to r=\lfloor m/(m/l)\rfloor\)）可在 \(O(\sqrt m)\) 时间求一个 \(\Phi(m)\)。

代码

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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

ll MOD = 1000000007;
ll MODM1 = MOD - 1;
typedef long long ll;

struct PairMod
{
    int a, b;
};

int mulmod(ll x, ll y, int m) { return (int)(x * y % m); }

int LIM;
vector<int> phiSmall, prefMod, prefModM1;
unordered_map<ll, PairMod> memo;

PairMod tri(ll n)
{
    ll a = n, b = n + 1;
    if (a & 1)
        b >>= 1;
    else
        a >>= 1;
    return {
        (int)((a % MOD) * (b % MOD) % MOD),
        (int)((a % MODM1) * (b % MODM1) % MODM1)};
}

PairMod Phi(ll n)
{
    if (n <= LIM)
        return {prefMod[n], prefModM1[n]};
    auto it = memo.find(n);
    if (it != memo.end())
        return it->second;
    PairMod res = tri(n);
    for (ll l = 2, r; l <= n; l = r + 1)
    {
        ll v = n / l;
        r = n / v;
        PairMod t = Phi(v);
        ll cnt = r - l + 1;
        res.a = (res.a - mulmod(cnt % MOD, t.a, MOD) + MOD) % MOD;
        res.b = (res.b - mulmod(cnt % MODM1, t.b, MODM1) + MODM1) % MODM1;
    }
    memo[n] = res;
    return res;
}

int solve(ll n)
{
    LIM = (int)pow((long double)n, 2.0L / 3.0L);
    LIM = max(LIM, 1000000);

    phiSmall.assign(LIM + 1, 0);
    prefMod.assign(LIM + 1, 0);
    prefModM1.assign(LIM + 1, 0);

    vector<int> lp(LIM + 1, 0), primes;
    phiSmall[1] = 1;
    for (int i = 2; i <= LIM; i++)
    {
        if (lp[i] == 0)
        {
            lp[i] = i;
            primes.push_back(i);
            phiSmall[i] = i - 1;
        }
        for (int p : primes)
        {
            ll v = 1LL * p * i;
            if (v > LIM)
                break;
            lp[(int)v] = p;
            if (i % p == 0)
            {
                phiSmall[(int)v] = phiSmall[i] * p;
                break;
            }
            else
            {
                phiSmall[(int)v] = phiSmall[i] * (p - 1);
            }
        }
    }
    if (LIM >= 1)
        phiSmall[1] = 1;
    for (int i = 1; i <= LIM; i++)
    {
        prefMod[i] = (prefMod[i - 1] + phiSmall[i]) % MOD;
        prefModM1[i] = (prefModM1[i - 1] + phiSmall[i]) % MODM1;
    }

    int M = 1;
    ll S = 0;

    for (ll l = 1, r; l <= n; l = r + 1)
    {
        ll q = n / l;
        r = n / q;

        PairMod pir = Phi(r);
        PairMod pil = Phi(l - 1);

        int cntMod = (pir.a - pil.a + MOD) % MOD;
        int cntExp = (pir.b - pil.b + MODM1) % MODM1;

        int base = (int)((q + 1) % MOD);
        if (base == 0)
        {
            if (cntMod != 0)
                M = 0;
        }
        else
        {
            M = mulmod(M, qpow(base, cntExp, MOD), MOD);
        }

        int qmod = (int)(q % MOD);
        int q2 = mulmod(qmod, qmod, MOD);
        S = (S + mulmod(q2, cntMod, MOD) + MOD) % MOD;
    }

    ll m4 = qpow(M, 4, MOD);
    ll m8 = mulmod(m4, m4, MOD);

    int nmod = n % MOD;
    int termN = mulmod(mulmod(4, nmod, MOD), (nmod + 1) % MOD, MOD);

    ll ans = (m8 - termN * m4 % MOD + 4ll * S % MOD - 1 + MOD * 2) % MOD;
    return ans;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    ll n = 1;
    for (int i = 0; i < 13; i++)
        n *= 7;
    cout << solve(n) << '\n';
    return 0;
}