【学术篇】分析矿洞--省选数论从入门到放弃-优快云博客

本文详细解析了一道数论题目的解题思路，重点介绍了如何使用杜教筛算法解决大范围求和问题，并给出了具体实现代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

数论什么的都去死吧！

看着题解我都能化式子用完4页草纸。。。
另外吐槽一句出题人的拼音学的是真好, 不知道是不是故意的.
其实题解已经写得挺详细的了.
我就是提一些出题人觉得太easy没必要提但是做题还是需要的一些东西….(因为这些东西我基本都是现学的)

然而之前刚刚学完mobius反演就暂时性脱坑的我啥也不会啊..
看到前排dp和曲神在水luogu的欢(bao/du)乐(ling/liu)赛, 就想去看看.
然后就点了报名但是发现自己什么都不会.

去看了看T1. 就是这道题.
然后成功的化出了第一步的式子.
这样就可以水30pts了.
一眼看出应该是反演类型的题目, 但是真的tmd不会啊,, 80pts都水不到.
(部分分给的也是有点迷, 80pts和100pts完全不是一样东西好么= =)

30pts:
通过一眼看出法可以得到激光扫到的第一个点的坐标是

(x g c d ( x , y ), y g c d ( x , y ))

$(\frac x{gcd(x,y)},\frac y{gcd(x,y)})$
所以

(x + y x g c d ( x , y ) + y g c d ( x , y )) 2 = g c d 2 (x, y)

$(\frac{x+y}{\frac x{gcd(x,y)}+\frac y{gcd(x,y)}})^2=gcd^2(x,y)$
于是很显然地就是要求出

\sum i = 1 N \sum j = 1 N φ (g c d 2 (i, j)) \dots ①

$\sum_{i=1}^N\sum_{j=1}^N\varphi(gcd^2(i,j))\cdots①$
这个东西, 于是就变成了一道纯数论题. (本来就是一道纯数论题不是?!)

然后就 $O(n^2logn)$ 乱搞就30pts到手了.

80pts:
继续化式子.
对于 $\varphi$ 我们有这么一种操作:

φ (n) = \prod i p k i - 1 i (p i - 1)

$\varphi(n)=\prod_ip_i^{k_i-1}(p_i-1)$
所以可以得到

φ (n m) = \prod i p m k i - 1 i (p i - 1) = \prod i p k i - 1 i (p i - 1) * p (m - 1) k i i = φ (n) * n m - 1

$\varphi(n^m)=\prod_ip_i^{mk_i-1}(p_i-1)=\prod_ip_i^{k_i-1}(p_i-1)*p_i^{(m-1)k_i}=\varphi(n)*n^{m-1}$
我们就可以把①中的gcd拆出来, 把平方去掉, 变成

① = \sum i = 1 N \sum j = 1 N φ (g c d (i, j)) * g c d (i, j)

$①=\sum_{i=1}^{N}\sum_{j=1}^{N}\varphi(gcd(i,j))*gcd(i,j)$
然后反演的常见套路之 枚举公因数

= \sum i = 1 N \sum j = 1 N \sum d | i, d | j φ (d) * d * b [g c d (i, j) = d]

$=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{d|i,d|j}\varphi(d)*d*b[gcd(i,j)=d]$
其中

b[x] $b[x]$ 表示

x $x$ 的真假性,

x $x$ 真则

b[x]=1 $b[x]=1$ , 否则

b[x]=0 $b[x]=0$
也就等价于在一个边长为

⌊nd⌋ $\left \lfloor\frac nd\right \rfloor$ 的方阵中找互质的

(i,j) $(i,j)$ , 然后对

d∗φ(d) $d*\varphi(d)$ 求和.

= \sum d = 1 N φ (d) * d \sum i = 1 ⌊ N d ⌋ \sum j = 1 ⌊ N d ⌋ b [g c d (i, j) = 1] \dots ②

$=\sum_{d=1}^N\varphi(d)*d\sum_{i=1}^{\left \lfloor\frac Nd\right \rfloor}\sum_{j=1}^{\left \lfloor\frac Nd\right \rfloor} b[gcd(i,j)=1]\cdots②$
后半边式子有点眼熟?? 我们好像在哪里见过这种形式.

仪仗队!!!!!!

我们可以得出

\sum i = 1 N \sum j = 1 N b [g c d (i, j) = 1]

$\sum_{i=1}^N\sum_{j=1}^Nb[gcd(i,j)=1]$
这个式子的结果是

2∗∑Ni=1φ(N)−1 $2*\sum_{i=1}^N\varphi(N)-1$ .
所以

② = \sum d = 1 N φ (d) * d * (2 * \sum i = 1 ⌊ N d ⌋ φ (i) - 1) \dots ③

$②=\sum_{d=1}^N\varphi(d)*d*(2*\sum_{i=1}^{\left \lfloor\frac Nd\right \rfloor}\varphi(i)-1)\cdots③$
这样这个式子就化到头了.
而此时我们枚举

d $d$ 就可以做到

O(nn√) $O(n\sqrt n)$ 求解了..
这样就能水到80pts. (个人感觉80pts部分分给得略高了)

100pts
满分做法就要用到一种高端科技了..

杜教筛!

顾名思义是一种筛法. 但是要比线筛快一些.
举个栗子, 我们来求一下

\sum i = 1 N φ (i)

$\sum_{i=1}^N\varphi(i)$ (其实跟上面的式子是有联系的OvO)
那么我们看数据范围想算法:
-

N<=1000 $N<=1000$ ?
- 枚举, 每次从头扫求一遍欧拉函数都能过.
-

N<=10000000 $N<=10000000$ ?
-

φ(x) $\varphi(x)$ 是个积性函数, 直接线筛一下就好咯.
-

N<=10000000000 $N<=10000000000$ ?
- 这个……

O(n) $O(n)$ 过不了啊..
这就是说我们必须要想别的方法了.
比如 杜教筛..
我们先来化一波式子, 尽量把N变小到能做的范围.
对于

φ $\varphi$ 函数有一条:

\sum d | n φ (d) = n

$\sum_{d|n}\varphi(d)=n$
那么

\sum d | n, d < n φ (d) + φ (n) = n φ (n) = n - \sum d | n, d < n φ (d)

$\sum_{d|n,d<n}\varphi(d)+\varphi(n)=n\\ \varphi(n)=n-\sum_{d|n,d<n}\varphi(d)$
我们

ϕ (i) = \sum i = 1 N φ (i) = \sum i = 1 N (i - \sum d | i, d < i φ (d)) = \sum i = 1 N i - \sum i = 2 N \sum d | i, d < i φ (d) = \sum i = 1 N i - \sum i d = 2 N \sum d = 1 ⎢ ⎣ ⎢ ⎢ n i d ⎥ ⎦ ⎥ ⎥ φ (d) 令 j = i d, 则 ϕ (i) \sum i = 1 N φ (i) = \sum i = 1 N i - \sum j = 2 N \sum d = 1 ⌊ n j ⌋ φ (d) = \sum i = 1 N - \sum j = 2 N ϕ (⌊ n j ⌋)

$\phi(i)= \sum_{i=1}^N\varphi(i)=\sum_{i=1}^N(i-\sum_{d|i,d<i}\varphi(d))=\sum_{i=1}^Ni-\sum_{i=2}^N\sum_{d|i,d<i}\varphi(d) =\sum_{i=1}^Ni-\sum_{\frac id=2}^N\sum_{d=1}^{\left\lfloor\frac n{\frac id}\right\rfloor}\varphi(d)\\ 令j=\frac id,则\phi(i)\sum_{i=1}^N\varphi(i)=\sum_{i=1}^Ni-\sum_{j=2}^N\sum_{d=1}^{\left\lfloor\frac nj\right\rfloor}\varphi(d)=\sum_{i=1}^N-\sum_{j=2}^N\phi(\left\lfloor\frac nj\right\rfloor)$
其中减号前面的显然是可以

O(1) $O(1)$ 计算的(别说你不会), 后面的值是不会超过

n√ $\sqrt n$ 个的, 我们枚举因数递归计算即可.
代码太长而且基本是这题代码的子集就不糊在这里了 ..留个传送门自己去看吧.

然后我们回到这道题. 我们化出了③式, 为了防止忘掉, 我们再贴一遍.

\sum d = 1 N φ (d) * d * (2 * \sum i = 1 ⌊ N d ⌋ φ (i) - 1)

$\sum_{d=1}^N\varphi(d)*d*(2*\sum_{i=1}^{\left \lfloor\frac Nd\right \rfloor}\varphi(i)-1)$
这样括号里面的刚刚学习了怎么筛(所以说是子集嘛), 所以问题就集中在了前面的

\sum d = 1 N φ (d) * d

$\sum_{d=1}^N\varphi(d)*d$
怎么快速的筛出来. 而这个题解已经说的挺清楚了的 ~~(说你是不是懒得继续化了←_←)~~
我们令

f(i)=∑Nd=1φ(d)∗d $f(i)=\sum_{d=1}^N\varphi(d)*d$
这个

φ(d)∗d=φ(d2) $\varphi(d)*d=\varphi(d^2)$ , 我们就猜测和

∑Ni=1i2 $\sum_{i=1}^Ni^2$ (这个式子可以用平方和公式

O(1) $O(1)$ 求哟~)有什么联♂系.
(这个理由是蒙的, 比赛的时候怎么凑知道该怎么凑出来不是很急但是会在线等= =)
那我们就化一下

\sum i = 1 N i 2 = \sum i = 1 N (i * \sum d | i φ (d)) = \sum i = 1 N (i d * \sum d | n φ (d) * d)

$\sum_{i=1}^Ni^2=\sum_{i=1}^N(i*\sum_{d|i}\varphi(d))=\sum_{i=1}^N(\frac id*\sum_{d|n}\varphi(d)*d)$
这样凑出了

φ(d)∗d $\varphi(d)*d$ 的形式. 但是还没有

∑ni=1 $\sum_{i=1}^n$ 的形式, 考虑枚举.
我们令

t=id $t=\frac id$ , 然后枚举

t $t$ . 因为要求1..N的和, 所以如果有

t∗d>N $t*d>N$ 的

d $d$ 显然不能对答案做出贡献, 所以我们枚举

d $d$ 的时候枚举到

⌊Nt⌋ $\left\lfloor\frac Nt\right\rfloor$ 即可. 也就是说

\sum i = 1 N i 2 = \sum t = 1 N t * \sum d = 1 ⌊ N t ⌋ φ (d) * d = \sum d = 1 N φ (d) * d + \sum t = 2 N t * \sum d = 1 ⌊ N t ⌋ φ (d) * d

$\sum_{i=1}^Ni^2=\sum_{t=1}^Nt*\sum_{d=1}^{\left\lfloor\frac Nt\right\rfloor}\varphi(d)*d=\sum_{d=1}^N\varphi(d)*d+\sum_{t=2}^Nt*\sum_{d=1}^{\left\lfloor\frac Nt\right\rfloor}\varphi(d)*d$
那我们就可以得到

f (i) = \sum i = 1 N i 2 - \sum t = 2 N t * \sum d = 1 ⌊ N t ⌋ φ (d) * d = \sum i = 1 N i 2 - \sum t = 2 N f (⌊ N t ⌋)

$f(i)=\sum_{i=1}^Ni^2-\sum_{t=2}^Nt*\sum_{d=1}^{\left\lfloor\frac Nt\right\rfloor}\varphi(d)*d=\sum_{i=1}^Ni^2-\sum_{t=2}^Nf(\left\lfloor\frac Nt\right\rfloor)$
就这么得到了一个杜教筛的形式, 就可以仿照上面做咯~

md化式子化到吐系列……

听说这题卡常数? 但也没怎么卡嘛感觉随便一跑就轻松第一页了?
(好像出题人改了题, 所以只算改题之后的话应该就rank1了= =之前q1~q4数组开大了memset废掉好多时间= =
用了一些小trick 比如全程能开int绝对不开long long… (所以因为少%了一次第一次交80pts嘤嘤嘤…
比如把结果记忆化一下. 看到std这一步是用map做的.
但是因为都是 $N$ 的因数, 我们就可以把这些因数分比一个数大的和比这个数小的两类分
这个数大约取 $\sqrt N$ 即可. (大点也能存下, 但是小点好像可能会出现冲突)
开int的话一定要记得:步步取模, 强转long long!
~~Emmmm还有一种非常无良的针对数据的压常trick就是n<=1e5的时候线筛可以少筛一点...(这样就稳稳rank1了)~~

const SQ=1e6/7;     //142857
int p1[SQ],p2[SQ];
int& getaddr(LL x){
    if(x<SQ) return p1[x];
    return p2[n/x];
} //返回储存地址

这样就能少个log, 就会非常快…

代码~~(为什么不去看std呢)~~

#include <cstdio>
#include <cstring>
const int N=5e6+5;
const int P=1e9+7;
const int SQ=1e6/7;
typedef long long LL;
LL n;
int p1[SQ],p2[SQ],p3[SQ],p4[SQ];
int prime[N],euler[N],mu[N],eusum[N],eumul[N],tot;
bool notp[N];
void shai(){
    notp[1]=euler[1]=mu[1]=eusum[1]=eumul[1]=1;
    for(int i=2;i<=5e6;++i){
        if(!notp[i]){
            prime[++tot]=i;
            euler[i]=i-1;
            mu[i]=-1;
        }
        for(int j=1;j<=tot&&i*prime[j]<5e6;++j){
            int k=i*prime[j];
            notp[k]=1;
            if(i%prime[j]==0){
                mu[k]=0;
                euler[k]=euler[i]*prime[j];
                break;
            }
            else{
                mu[k]=-mu[i];
                euler[k]=euler[i]*(prime[j]-1);
            }
        }
        eusum[i]=(eusum[i-1]+euler[i])%P;
        eumul[i]=(eumul[i-1]+1LL*euler[i]*i%P)%P;
    }
}
inline int qpow(int a,int b,int s=1){
    for(;b;b>>=1,a=1LL*a*a%P)
        if(b&1) s=1LL*s*a%P;
    return s;
}
int inv2=qpow(2,P-2),inv6=qpow(6,P-2);
inline int& getaddr(LL x,bool flag){
    if(flag){
        if(x<SQ) return p1[x];
        return p2[n/x];
    }
    if(x<SQ) return p3[x];
    return p4[n/x]; 
}
int eulersum(LL x){
    if(x<=5e6) return eusum[x];
    int& addr=getaddr(x,1);
    if(addr!=-1) return addr;
    int ans; LL last;
    ans=1LL*(x%P)*(x%P+1)%P*inv2%P;
    for(LL i=2;i<=x;i=last+1){
        last=x/(x/i);
        ans=(ans-1LL*(last-i+1)*eulersum(x/i)%P)%P;
    }
    return addr=(ans%P+P)%P;
}
int eumulsum(LL x){
    if(x<=5e6) return eumul[x];
    int& addr=getaddr(x,0);
    if(addr!=-1) return addr;
    int ans; LL last;
    ans=1LL*(x%P)*(x%P+1)%P*((x+x+1)%P)%P*inv6%P;
    for(LL i=2;i<=x;i=last+1){
        last=x/(x/i);
        ans=(ans-1LL*(i+last)%P*(last-i+1)%P*inv2%P*eumulsum(x/i)%P)%P;
    }
    return addr=(ans%P+P)%P;
}
int main(){
    memset(p1,-1,sizeof(p1));
    memset(p2,-1,sizeof(p2));
    memset(p3,-1,sizeof(p3));
    memset(p4,-1,sizeof(p4));
    shai(); scanf("%lld",&n);
    int s=0; LL last;
    for(LL i=1;i<=n;i=last+1){
        last=n/(n/i);
        s=(s+1LL*(eumulsum(last)-eumulsum(i-1)+P)*((2*eulersum(n/i)%P-1+P)%P)%P)%P;
    }
    printf("%d",s);
}

真是要吐了完结撒花~