JZOJ 5261 求和

发表于 | 分类于 |

神仙数论。
真难写……

随手推一下，原式变为 \(\sum\limits_{x=1}^n \sum\limits_{d=1}^k f_d(x) (2\sum\limits_{i=1}^{\lfloor\frac nx\rfloor}\varphi(i) - 1)\)。
对 \(\varphi\) 杜教筛，然后考虑设 \(F_d(n) = \sum\limits_{i=1}^n f_d(i)\)。

考虑容斥，可以强行钦定几个质因子的指数大于 \(d\)，然后删去贡献（\(f_{\infty}\)），然后再把重复算的贡献补回。
即 \[ \begin{align*} F_d(n) =&\sum\limits_{i=1}^n f_d(i) \\ =&\sum\limits_{i=1}^n \mu(i) \sum\limits_{j=1}^{\left\lfloor\frac n{i^{d+1}}\right\rfloor} f_{\infty}(i^{d+1}j) \\ =&\sum\limits_{i=1}^n f_{\infty}^{d+1}(i) \mu(i) F_{\infty}\left(\left\lfloor\frac n{i^{d+1}}\right\rfloor\right) \\ =&\sum\limits_{i=1}^{\lfloor n^{1/(d+1)}\rfloor} f_{\infty}^{d+1}(i) \mu(i) F_{\infty}\left(\left\lfloor\frac n{i^{d+1}}\right\rfloor\right) \end{align*} \]

这是因为 \(f_{\infty}\) 是完全积性函数。
于是如果能求得 \(F_{\infty}\)，\(F_d\) 就可以在 \(O(n^{2/3}k^{2/3})\) 内求得，即预处理 + 暴力算（实际上由于空间问题达不到这个最优值）。

考虑杜教筛，注意到若 \(n = \prod\limits_{i=1}^m p_i^{c_i}\)，则 \((f_{\infty} * 1)(n) = \prod\limits_{i=1}^m [2\mathop{|}c_i]\)，即 \(n\) 为完全平方数。
于是就做完了。
（当然直接 Min_25 筛或者洲阁筛也不错）

但是有一点，即 \(F_d\) 的预处理。
实际上是可以直接把 \(\sum\limits_{d=1}^k f_d(n)\) 的前缀和线性筛出来的。
具体的说，由于每个 \(f_d\) 的贡献非 \(0\) 即 \(f_{\infty}\)，故只需要筛出每个数的质因子最大指数，然后算出有多少个 \(d\) 使得 \(f_d\) 有贡献，乘上 \(f_{\infty}\) 的值即可。

上文中有提到求 \(\sum\limits_{d=1}^k f_d\) 的复杂度，注意到若预处理 \((nk)^{2/3} = n^{2/3}k^{2/3}\) 内的值，可以达到最优复杂度 \(O(n^{2/3}k^{2/3})\)，但是由于空间问题无法达到。
所以就大概少预处理个 \(2 \times 10^7\) 的样子吧（逃

代码（杜教筛）： 

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
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
const long long N = 1e10;
const int MX = 2e7;
const int K = 40;
long long n;
int k;
int vis[MX + 5],cnt,prime[MX + 5];
int mn[MX + 5],mx[MX + 5],f[MX + 5],phi[MX + 5],mu[MX + 5],lambda[MX + 5];
int w[3][MX + 5];
inline int &mem_slambda(long long x)
{
    return w[0][n / x];
}
inline int &mem_sphi(long long x)
{
    return w[1][n / x];
}
inline int &mem_sf(long long x)
{
    return w[2][n / x];
}
int slambda(long long n)
{
    if(n <= MX)
        return lambda[n];
    if(~mem_slambda(n))
        return mem_slambda(n);
    int ret = sqrt(n);
    for(register long long l = 2,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret -= (r - l + 1) * slambda(n / l);
    }
    return mem_slambda(n) = ret;
}
int sphi(long long n)
{
    if(n <= MX)
        return phi[n];
    if(~mem_sphi(n))
        return mem_sphi(n);
    int ret = (n & 1) ? (n + 1) / 2 * n : n / 2 * (n + 1);
    for(register long long l = 2,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret -= (r - l + 1) * sphi(n / l);
    }
    return mem_sphi(n) = ret;
}
int sf(long long n)
{
    if(n <= MX)
        return f[n];
    if(~mem_sf(n))
        return mem_sf(n);
    int ret = 0;
    for(register int i = 1;(long long)i * i <= n;++i)
        if(mu[i])
        {
            long long j = (long long)i * i;
            int s = (lambda[i] - lambda[i - 1]) * (lambda[i] - lambda[i - 1]);
            for(register int d = 2;d <= k + 1 && j <= n;++d,j *= i,s *= (lambda[i] - lambda[i - 1]))
                ret += s * mu[i] * slambda(n / j);
        }
    return mem_sf(n) = ret;
}
int ans;
int main()
{
    scanf("%lld%d",&n,&k);
    memset(w,-1,sizeof w),phi[1] = mu[1] = lambda[1] = 1,f[1] = k;
    for(register int i = 2;i <= MX;++i)
    {
        if(!vis[i])
            mn[prime[++cnt] = i] = mx[i] = 1,phi[i] = i - 1,lambda[i] = mu[i] = -1;
        for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
        {
            vis[i * prime[j]] = 1,lambda[i * prime[j]] = -lambda[i];
            if(!(i % prime[j]))
            {
                mn[i * prime[j]] = mn[i] + 1,mx[i * prime[j]] = max(mx[i],mn[i] + 1),phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            mn[i * prime[j]] = 1,mx[i * prime[j]] = mx[i],phi[i * prime[j]] = phi[i] * (prime[j] - 1),mu[i * prime[j]] = -mu[i];
        }
        f[i] = f[i - 1] + lambda[i] * (mx[i] <= k ? k - mx[i] + 1 : 0),phi[i] += phi[i - 1],lambda[i] += lambda[i - 1];
    }
    for(register long long l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ans += (sf(r) - sf(l - 1)) * (2 * sphi(n / l) - 1);
    }
    printf("%d\n",ans & (1 << 30) - 1);
}

代码（Min_25 筛）： 

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
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
const long long N = 1e10;
const int MX = 2e7;
const int K = 40;
long long n;
int k;
int vis[MX + 5],cnt,prime[MX + 5];
int mn[MX + 5],mx[MX + 5],f[MX + 5],mu[MX + 5],lambda[MX + 5];
int w[MX + 5];
namespace Min_25
{
    const long long N = 1e10;
    const int MX = 1e5;
    long long n;
    int lim;
    int vis[MX + 5],prime[MX + 5],cnt,Gprime[MX + 5];
    int tot,le[MX + 5],ge[MX + 5];
    long long lis[2 * MX + 5];
    int G[2 * MX + 5][2],F[2 * MX + 5][2];
    inline int &id(long long x)
    {
        return x <= lim ? le[x] : ge[n / x];
    }
    void init(long long m)
    {
        lim = sqrt(n = m);
        for(register int i = 2;i <= lim;++i)
        {
            if(!vis[i])
                prime[++cnt] = i,Gprime[cnt] = Gprime[cnt - 1] + i;
            for(register int j = 1;j <= cnt && i * prime[j] <= lim;++j)
            {
                vis[i * prime[j]] = 1;
                if(!(i % prime[j]))
                    break;
            }
        }
        for(register long long l = 1,r;l <= n;l = r + 1)
        {
            r = n / (n / l);
            lis[id(n / l) = ++tot] = n / l;
            G[tot][0] = n / l - 1,G[tot][1] = (n / l & 1) ? (n / l - 1) / 2 * (n / l + 2) : (n / l + 2) / 2 * (n / l - 1);
        }
        for(register int k = 1;k <= cnt;++k)
        {
            int p = prime[k];
            long long s = (long long)p * p;
            for(register int i = 1;lis[i] >= s;++i)
                G[i][0] -= G[id(lis[i] / p)][0] - (k - 1),
                G[i][1] -= (long long)p * (G[id(lis[i] / p)][1] - Gprime[k - 1]);
        }
        for(register int i = 1;i <= tot;++i)
            F[i][0] = -G[i][0],F[i][1] = G[i][1] - G[i][0];
        for(register int k = cnt;k;--k)
        {
            int p = prime[k];
            long long s = (long long)p * p;
            for(register int i = 1;lis[i] >= s;++i)
            {
                long long pw = p,pw2 = s,pw0 = 1;
                for(register int c = 1;pw2 <= lis[i];++c,pw0 = pw,pw = pw2,pw2 *= p)
                    F[i][0] += ((c & 1) ? -1 : 1) * (F[id(lis[i] / pw)][0] + k) + ((c & 1) ? 1 : -1),
                    F[i][1] += (pw - pw0) * (F[id(lis[i] / pw)][1] - (Gprime[k] - k)) + (pw2 - pw);
            }
        }
    }
    inline int slambda(long long n)
    {
        return n ? F[id(n)][0] + 1 : 0;
    }
    inline int sphi(long long n)
    {
        return n ? F[id(n)][1] + 1 : 0;
    }
}
inline long long fpow(long long a,int b)
{
    long long ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret *= a),a *= a;
    return ret;
}
inline int &mem_sf(long long x)
{
    return w[n / x];
}
int sf(long long n)
{
    if(n <= MX)
        return f[n];
    if(~mem_sf(n))
        return mem_sf(n);
    int ret = 0;
    for(register int i = 1;(long long)i * i <= n;++i)
        if(mu[i])
        {
            long long j = (long long)i * i;
            int s = lambda[i] * lambda[i];
            for(register int d = 2;d <= k + 1 && j <= n;++d,j *= i,s *= lambda[i])
                ret += s * mu[i] * Min_25::slambda(n / j);
        }
    return mem_sf(n) = ret;
}
int ans;
int main()
{
    scanf("%lld%d",&n,&k),Min_25::init(n);
    memset(w,-1,sizeof w),mu[1] = lambda[1] = 1,f[1] = k;
    for(register int i = 2;i <= MX;++i)
    {
        if(!vis[i])
            mn[prime[++cnt] = i] = mx[i] = 1,lambda[i] = mu[i] = -1;
        for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
        {
            vis[i * prime[j]] = 1,lambda[i * prime[j]] = -lambda[i];
            if(!(i % prime[j]))
            {
                mn[i * prime[j]] = mn[i] + 1,mx[i * prime[j]] = max(mx[i],mn[i] + 1);
                break;
            }
            mn[i * prime[j]] = 1,mx[i * prime[j]] = mx[i],mu[i * prime[j]] = -mu[i];
        }
        f[i] = f[i - 1] + lambda[i] * (mx[i] <= k ? k - mx[i] + 1 : 0);
    }
    for(register long long l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ans += (sf(r) - sf(l - 1)) * (2 * Min_25::sphi(n / l) - 1);
    }
    printf("%d\n",ans & (1 << 30) - 1);
}