LibreOJ 6686 Stupid GCD

发表于 | 分类于 | 12

真 TM 毒瘤……
典型的式子好推不可实现题……

一个显然的想法是先枚举 r=⌊3√i⌋，然后再考虑 i 的贡献。
即 ⌊3√i⌋∑i=1gcd(⌊3√i⌋,i)=⌊3√i⌋∑r=1n∑i=1[⌊3√i⌋=r]gcd(r,i)。

又因为 ⌊3√i⌋=r 等价于 r3≤i<(r+1)3，相当于把这些数分成了 ⌊3√i⌋ 块。
得 ⌊3√i⌋∑i=1gcd(⌊3√i⌋,i)=n∑i=⌊3√n⌋3gcd(⌊3√i⌋,i)+⌊3√i⌋−1∑r=1(r+1)3−1∑i=r3gcd(r,i)。

根据基本反演套路易得 n∑i=1gcd(r,i)=∑d|r⌊nd⌋φ(d)。
于是两个和式分别算即可。

第一个和式需要 ⌊3√n⌋ 的约数的欧拉函数值，可以考虑类似杜教筛的思想，预处理 n2/9 之内的值，然后别的值分解质因数硬算（只需要 n1/6 以下的质数）。
复杂度比较难算，姑且认为是 O(n2/9)（实际上跑不满）。
（实际上可以通过 DFS 做到 O(n1/6)，因为都是 ⌊3√n⌋ 的约数，但是我懒得写）

第二个和式，根据上式得 ⌊3√i⌋−1∑r=1(r+1)3−1∑i=r3gcd(r,i)=⌊3√i⌋−1∑r=1∑d|r(⌊(r+1)3−1d⌋−⌊r3−1d⌋)φ(d)=⌊3√i⌋−1∑d=1φ(d)⌊⌊3√i⌋−1d⌋∑r=1(⌊(dr+1)3−1d⌋−⌊(dr)3−1d⌋)

注意到 ⌊(dr+1)3−1d⌋=⌊d3r3+3d2r2+3drd⌋=d2r3+3dr2+3r，且 ⌊(dr)3−1d⌋=⌊d3r3−1d⌋=d2r3−1。
故 ⌊3√i⌋−1∑d=1φ(d)⌊⌊3√i⌋−1d⌋∑r=1(⌊(dr+1)3−1d⌋−⌊(dr)3−1d⌋)=⌊3√i⌋−1∑d=1φ(d)⌊⌊3√i⌋−1d⌋∑r=1(3dr2+3r+1)

然后稍微杜教筛一下就行了。

代码： 

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
#include <cstdio>
#include <cstring>
using namespace std;
template<class T>
inline void read(T &x)
{
    char ch = 0,w = 0;
    for(;ch < '0' || ch > '9';w |= ch == '-',ch = getchar());
    for(x = 0;ch >= '0' && ch <= '9';x = (x << 3) + (x << 1) + (ch ^ '0'),ch = getchar());
    w && (x = -x);
}
const int MX = 5e6;
const int mod = 998244353;
const int inv2 = 499122177;
const int inv6 = 166374059;
__int128 n;
long long rt;
int vis[MX + 5],cnt,prime[MX + 5],phi[MX + 5],sum[MX + 5];
inline long long cubert(__int128 n)
{
    __int128 l = 1,r = n,mid,ret;
    while(l <= r)
        mid = l + r >> 1,mid <= n / mid / mid ? (l = mid + 1,ret = mid) : (r = mid - 1);
    return ret;
}
int w[3][MX + 5];
inline int &mem_calc(long long x)
{
    return w[0][rt / x];
}
inline int &mem_sphi(long long x)
{
    return w[1][rt / x];
}
inline int &mem_sidphi(long long x)
{
    return w[2][rt / x];
}
int calc(long long n)
{
    if(n <= MX)
        return phi[n] - phi[n - 1];
    if(~mem_calc(n))
        return mem_calc(n);
    long long ret = n,t = n;
    for(register int i = 1;(long long)prime[i] * prime[i] <= n;++i)
    {
        !(t % prime[i]) && (ret -= ret / prime[i]);
        for(;!(t % prime[i]);t /= prime[i]);
    }
    t > 1 && (ret -= ret / t);
    return mem_calc(n) = ret % mod;
}
int sphi(long long n)
{
    if(n <= MX)
        return phi[n];
    if(~mem_sphi(n))
        return mem_sphi(n);
    long long ret = (n % mod) * ((n + 1) % mod) % mod * inv2 % mod;
    for(register long long l = 2,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret = (ret - (r - l + 1) % mod * sphi(n / l) % mod + mod) % mod;
    }
    return mem_sphi(n) = ret;
}
int sidphi(long long n)
{
    if(n <= MX)
        return sum[n];
    if(~mem_sidphi(n))
        return mem_sidphi(n);
    long long ret = (n % mod) * ((n + 1) % mod) % mod * ((2 * n + 1) % mod) % mod * inv6 % mod;
    for(register long long l = 2,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret = (ret - ((l + r) % mod) * ((r - l + 1) % mod) % mod * inv2 % mod * sidphi(n / l) % mod + mod) % mod;
    }
    return mem_sidphi(n) = ret;
}
int ans;
int main()
{
    memset(w,-1,sizeof w),sum[1] = phi[1] = 1;
    for(register int i = 2;i <= MX;++i)
    {
        if(!vis[i])
            phi[prime[++cnt] = i] = i - 1;
        for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
            {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
        sum[i] = (sum[i - 1] + (long long)i * phi[i] % mod) % mod,phi[i] = (phi[i] + phi[i - 1]) % mod;
    }
    read(n),rt = cubert(n);
    for(register int i = 1;(long long)i * i <= rt;++i)
        if(!(rt % i))
        {
            ans = (ans + (long long)((int)(n / i % mod) - (int)(((__int128)rt * rt * rt - 1) / i % mod) + mod) * calc(i) % mod) % mod;
            if(i ^ rt / i)
                ans = (ans + (long long)((int)(n / (rt / i) % mod) - (int)(((__int128)rt * rt * rt - 1) / (rt / i) % mod) + mod) * calc(rt / i) % mod) % mod;
        }
    --rt;
    for(register long long l = 1,r;l <= rt;l = r + 1)
    {
        r = rt / (rt / l);
        ans = (ans + 3LL * (sidphi(r) - sidphi(l - 1) + mod) % mod * (rt / l % mod) % mod * ((rt / l + 1) % mod) % mod * ((rt / l * 2 + 1) % mod) % mod * inv6 % mod) % mod;
        ans = (ans + (long long)(sphi(r) - sphi(l - 1) + mod) * (3LL * (rt / l % mod) % mod * ((rt / l + 1) % mod) % mod * inv2 % mod + (rt / l) % mod) % mod) % mod;
    }
    printf("%d\n",ans);
}

Related Issues not found

Please contact @Alpha1022 to initialize the comment