P1587 [NOI2016]循环之美

发表于 | 阅读数

$\sum_{x=1}^n\sum_{y=1}^m \frac{x}{y}$在$k$进制下能表示成循环节从第一位小数开始的无限循环小数或整数的最简分数个数

设循环长度为$p$

$gcd(k,y)=1$

先老老实实化简

边界$F(a,b,1)=\sum_{i=1}^{min(a,b})\mu(d)\frac{a}{d}\frac{b}{d}$ 杜教筛也只需要$\sqrt{n}$的空间

$O(\log n \log m \sqrt{k} + n^{\frac{2}{3}})$

代码 
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

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define pii pair<int, int>
#define mk make_pair
const int N = 2e6 + 10;
const int mod = 1e9 + 7;

int pcnt, pr[N], npr[N];
int mu[N], usum[N];

void Prime_init(int X)
{
    npr[1] = 1;
    mu[1] = 1;
    for (int i = 1; i <= X; i++)
    {
        if (!npr[i])
            pr[++pcnt] = i, mu[i] = -1;
        for (int j = 1; j <= pcnt && pr[j] * i <= X; j++)
        {
            npr[pr[j] * i] = 1;

            if (i % pr[j] == 0)
            {
                mu[pr[j] * i] = 0;
                break;
            }
            else
            {
                mu[pr[j] * i] = mu[i] * (-1);
            }
        }
    }
    for (int i = 1; i <= X; i++)
        usum[i] = usum[i - 1] + mu[i];
}

map<int, ll> su;
ll Smu(int x)
{
    if (x <= 1e6 + 10)
        return usum[x];
    if (su[x])
        return su[x];
    ll res = 1;
    for (int i = 2, j; i <= x; i = j + 1)
    {
        j = x / (x / i);
        res -= 1ll * (j - i + 1) * Smu(x / i);
    }
    return su[x] = res;
}
struct node
{
    int n, m, k;
    inline bool operator<(const node &b) const { return n == b.n ? (m == b.m ? k < b.k : m < b.m) : n < b.n; }
};

map<node, ll> mp;
ll S(int n, int m, int k)
{
    node now = {n, m, k};
    if (mp[now])
        return mp[now];
    if (!n || !m)
        return 0;
    ll res = 0;
    if (k == 1)
    {

        int d = min(n, m);

        for (int i = 1, j; i <= d; i = j + 1)
        {
            j = min(m / (m / i), n / (n / i));
            res += 1ll * (n / i) * (m / i) * (Smu(j) - Smu(i - 1));
        }
    }
    else
    {
        vector<int> g;
        //g.clear();
        int P = sqrt(k);
        for (int i = 1; i <= P; i++)
        {
            if (k % i == 0)
            {
                g.push_back(i);
                if (k / i != i)
                    g.push_back(k / i);
            }
        }

        for (int x : g)
            if (mu[x] != 0)
            {
                res += 1ll * mu[x] * S(m / x, n, x);
            }
    }
    return mp[now] = res;
}
int main()
{
    int n, m, k;
    cin >> n >> m >> k;
    Prime_init(1e6 + 10);
    Smu(m), Smu(n);
    printf("%lld", S(n, m, k));
}