noi.ac 716 答案是整数 题解

发表于 | 阅读次数：
本文字数： 1.5k 字 | 阅读时长 ≈ 7 分钟

题意简述

求
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\dfrac{ij}{\gcd (i,j)}$。

对 $998244353$ 取模。$n<=10^7,m<=10^{14}$。

注：空间限制只有64MB，只够开 $10^7$ 的int型数组。

思路

基础操作

设 $n\le m$，因为 $n,m$ 对称。
设 $g=\gcd(i,j)$ （简写）
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} \dfrac{ij}{g^2}$
$\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}[g=1]ij$ （枚举$\gcd$）
$\sum\limits_{q=1}^{n}\mu(q)\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}[q|g]ij$
$\sum\limits_{q=1}^{n}\mu(q)\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/dq}\sum\limits_{j=1}^{m/dq}iq\times jq$
$\sum\limits_{q=1}^{n}\mu(q)q^2\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/dq}\sum\limits_{j=1}^{m/dq}ij$
设 $s(x)=\sum\limits_{i=1}^x i$
原式化为：
$\sum\limits_{q=1}^n\mu(q)q^2\sum\limits_{d=1}^{n}s(n/dq)s(m/dq)$
我们发现， $dq>n$ 时，$s(n/dq)=0$，也就不会有贡献了。那么我们珂以缩小 $d$ 的范围：
$\sum\limits_{q=1}^{n}\mu(q)q^2\sum\limits_{d=1}^{n/q}s(n/dq)s(m/dq)$

整除分块两次，是 $O(n)$ 的（？）。

但是似乎这个复杂度还和 $m$ 有点关系，我不是很清楚…总之，它TLE了，只有 $80$ 分。（亲测）

那么我们考虑枚举 $q\times d$，设为 $T$。它的贡献就是 $s(n/T)s(m/T)$ 。而它被算到的此时，显然， $T=dq$ 是 $q$ 的倍数，所以 $q$ 是 $T=dq$ 的因数。那 $dq=T$ 被算到的次数就是 $\sum\limits_{q|T} \mu(q)q^2$。将这个值设为 $f(T)$。

式子化一下，变成：$\sum\limits_{T=1}^{n}f(T)\times s(n/T)\times s(m/T)$。

那么问题就在于我们怎么筛这个 $f(T)$ ，还要严格线性，因为 $n<=10^7$。

如何筛f函数

显然，$f(1)=1$，而且当 $p$ 和 $q$ 互质时，$f(p\times q)=f(p)\times f(q)$。这说明它也许能线性筛。

考虑质数情况。当 $p$ 为质数时，$f(p)=1-p^2$。

接着，$f(p^k)=\sum\limits_{d|p^k}\mu(d)d^2$。显然，$p^k$ 的因数只有 $p^0,p^1\cdots p^k$。而对于 $p^j$ （$2\le j \le k$），显然，$\mu(p^j)=0$，也就不会有贡献了。因此，有贡献的只有 $p^0$ 和 $p^1$ 。

那也就是说，$f(p^k)=f(p)$！！！（换句话说，同一个质因子乘多少遍都不会改变 $f$ 的值）

那么我们在线性筛的时候有这样一步：枚举 $i$，找一个质数 $u$。如果 $u$ 不是 $i$ 的因数，那么显然 $i$ 和 $u$ 互质，$f(i\times u)=f(i)\times f(u)$。否则 $u$ 是 $i$ 的质因子，然后我们要计算 $f(i\times u)$ 的值，并 break。

由上面那个性质得，$u$ 是 $i$ 的质因子，而我们又把它乘了一遍，并不会改变 $f$ 的值。也就是 $f(i\times u)=f(i)\times f(u)$。

知道了这些，就能线性筛这个 $f$ 函数了！！！

线性筛出来了 $f$ 之后，整除分块都不用，直接 $O(n)$ 暴力就能求出原式的值了。而且是很稳的 $O(n)$ 哦！

极 限 卡 常

上面说了，空间限制只有 64MB ，只能开的下一个 $10^7$ 个int数组，两个就不行了。

说的具体点，能开的下 $16 777 216$个 int。

但是我们常见的线性筛，需要一个int primes[1e7]保存质数，还有一个bool notp[1e7]标记是否不是素数，然后才是我们的f数组。这可怎么办呢？

首先，$[1,n]$ 中的质数个数和 $n$ 绝对不是同一个级别。具体的说，大概是 $\dfrac{n}{\ln(n)}$ 级别的。$[1,1e7]$ 中的质数，经过测试，还不足 $10^6$ 个。那么我们的 primes 数组就珂以开的小一点了。

接着，bool notp[1e7] 珂以用一个 bitset 存储，空间直接少掉 $32$ 倍，你说爽不爽。

然后我们就能卡进空间常数了！nice！

还有，记得不要全开long long 啊！

代码

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
#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
    #define N 10000007    
    #define mod 998244353ll
    #define i6  166374059ll
    #define i2  499122177ll
    // 模数，2的逆元，6的逆元 
    // 只是打个板子，并不是所有的数都有用上
    #define ll long long
    #define F(i,l,r) for(int i=l;i<=r;++i)
    #define D(i,r,l) for(int i=r;i>=l;--i)
    #define Fs(i,l,r,c) for(int i=l;i<=r;c)
    #define Ds(i,r,l,c) for(int i=r;i>=l;c)
    #define MEM(x,a) memset(x,a,sizeof(x))
    #define FK(x) MEM(x,0)
    #define Tra(i,u) for(int i=G.Start(u),__v=G.To(i);~i;i=G.Next(i),__v=G.To(i))
    #define p_b push_back
    #define sz(a) ((int)a.size())
    #define iter(a,p) (a.begin()+p)
    //并不想用快读，数据小

    ll n,m;
    void Input()
    {
        scanf("%lld%lld",&n,&m);
    }

    int f[N];
    int primes[N/10]; //N/10足够了
    bitset<N> notp; //STL万岁
    void Init()
    {
        int &cnt=primes[0];
        int n=1e7;
        f[1]=notp[1]=1;

        F(i,2,n)
        {
            if (!notp[i])
            {
                f[i]=(mod+1-1ll*i*i%mod)%mod;
                // 上面说了，i为质数的时候，f[i]=1-i*i
                // 别忘了取模啊
                primes[++cnt]=i;
            }

            for(int j=1;j<=cnt and i*primes[j]<=n;++j)
            {
                int u=primes[j];
                notp[i*u]=1;
                if (i%u!=0)
                {
                    f[i*u]=1ll*f[i]*f[u]%mod;
                    //利用积性性质，f[i*u]=f[i]*f[u]
                }
                else
                {
                    f[i*u]=f[i];
                    //上面特意讲了，此时f[i*u]=f[i]
                    break;
                }
            }
        }
    }
    ll s(ll x){x%=mod;return x*(x+1)%mod*i2%mod;}
    //计算1+2+3...+x的和，也就是x*(x+1)/2
    void Soviet()
    {
        if (n>m) swap(n,m);
        //令n<=m
        ll ans=0;
        F(i,1,n)
        {
            ans+=1ll*f[i]*s(n/i)%mod*s(m/i)%mod;
            //每次ans+=f[i]*s(n/i)*s(m/i)
            //别忘了取模啊！
            ans%=mod;
        }
        printf("%lld\n",ans);
    }

    #define Flan void
    Flan IsMyWife()
    {
        Init();
        Input();
        Soviet();
    }
}
int main()
{
    Flandre_Scarlet::IsMyWife();
    getchar();getchar();
    return 0;
}
w