min25筛学习笔记

目录

一、素数部分
二、合数部分
三、相关例题
例1、\(\texttt{P5325 【模板】Min\_25筛}\)
例2、\(\texttt{LOJ6235 区间素数个数}\)
例3、\(\texttt{LOJ6053 简单的函数}\)
例4、\(\texttt{UOJ188 【UR \#13】Sanrd}\)
例5、\(\texttt{LOJ572 「LibreOJ Round \#11」Misaka Network 与求和}\)
例6、\(\texttt{UOJ448 【集训队作业2018】人类的本质}\)

\(\min25\) 筛可以在亚线性时间内求积性函数 \(f\) 的前缀和。

一、素数部分

一些约定：

• \(p_i\) 表示从小到大第 \(i\) 个素数。
• \(\text{minp}(i)\) 表示 \(i\) 的最小素因子。
• \(v_p(i)\) 表示 \(i\) 含素数 \(p\) 的幂次。
• 若未特别说明， \(\frac xy\) 表示 \(\lfloor\frac xy\rfloor\) 。

目标对 \(n\) 的 \(2\sqrt n\) 个整除分块节点 \(m\) ，我们希望求出 \(f\) 在 \([1,m]\) 中素数位置点值之和。

即 \(\forall x\in\{\frac n1,\cdots,\frac nn\}\) ，求 \(G(m)=\sum_{i=1}^m[i\in\text{prime}]f(i)\) 。

我们需要找一个函数 \(f'\) 满足以下条件：

• \(f'\) 和 \(f\) 在素数处的点值全部相同。

• \(f'\) 是完全积性函数。

• \(f'\) 可以快速求前缀和。

目标转化为求 \(G(m)=\sum_{i=1}^m[i\in\text{prime}]f'(i)\) 。

构造 \(f'\) 的目的就是快速计算素数点值的前缀和，对于合数 \(x\) ， \(f'(x)\) 和 \(f(x)\) 可以不相等。当然，在合数部分也不会用到 \(f'\) 。

有时直接构造 \(f'\) 会有困难，可以将 \(f\) 拆成若干个完全积性函数的和，最后将所得结果相加，具体请看模板题。

---

考虑动态规划，定义：

\[g(n,j)=\sum_{i=1}^n[i\in\text{prime}\vee\text{minp}(i)\gt p_j]f'(i)\\ \]

即经过 \(j\) 轮埃氏筛后剩下的所有数之和。

初始值 \(g(n,0)=\sum_{i=2}^nf'(i)\) ，别忘了删掉 \(f'(1)\) 的贡献。

再来考虑转移。

• 如果 \(p_j^2\gt n\) ，那么第 \(j\) 轮不会筛掉任何数。

• 如果 \(p_j^2\le n\) ，那么第 \(j\) 轮被筛掉的是最小素因子为 \(p_j\) 的所有数， \(p_j\) 自身除外。

  \[\sum_{i=2\cdot p_j}^n[p_j\mid i\wedge\text{minp}(i)\ge p_j]f'(i)=f'(p_j)\sum_{i=2}^\frac n{p_j}[\text{minp}(i)\ge p_j]f'(i)\\ \]

  这里用到了 \(f'\) 是完全积性函数的条件，提出素因子 \(p_j\) 时，不需要考虑剩余部分是否与 \(p_j\) 互素。

  预处理 \([1,\sqrt n]\) 中， \(f'\) 在素数处点值的前缀和，记 \(s_k=\sum_{i=1}^kf'(p_i)\) 。

  继续化简：

  \[\sum_{i=1}^\frac n{p_j}[\text{minp}(i)\ge p_j]f'(i)=g(\frac n{p_j},j-1)-s_{j-1}\\ \]

完整的转移方程如下：

\[g(n,j)=\begin{cases} g(n,j-1)&p_j^2\gt n\\ g(n,j-1)-f'(p_j)\big(g(\frac n{p_j},j-1)-s_{j-1}\big)&p_j^2\le n\\ \end{cases} \]

代码实现时，可以用滚动数组优化掉第二维，倒序枚举 \(n\) 进行转移。

第一维虽然只有 \(2\sqrt n\) 个数，但还涉及到下标存储的问题。

开两个长为 \(B=\sqrt n\) 的数组 id1[],id2[] ，令 \(x\) 的编号为 x<=B?id1[x]:id2[n/x] 。

记 \([1,\sqrt n]\) 中素数个数为 \(cnt\) ，则 \(G(m)=g(m,cnt)\) 。

int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
for(int l=1,r=0;l<=n;l=r+1)
{
    int x=n/l;
    r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
    /*预处理g[x]=\sum_{i=2}^xf'(x),记得减掉f'(1)的贡献*/
}
for(int j=1;j<=cnt;j++)
    for(int i=1;p[j]*p[j]<=val[i];i++)
    {///val[]倒序存储所有整除分块节点,所以顺序枚举i
        int k=id(val[i]/p[j]);
        /*g[x]-=f'(p[j])*(g[k]-s[j-1]);*/
    }
/*G(x)=g[id(x)]*/

时间复杂度分析：

对每个 \(x=\lfloor\frac ni\rfloor\) ，它会被 \(\sqrt x\) 以内的每个素数枚举到，代价为 \(\mathcal O(\frac{\sqrt x}{\log\sqrt x})\) 。

因此这一部分的代价为：

\[T_1(n)=\sum_{i=1}^\sqrt n\mathcal O(\frac{\sqrt i}{\log\sqrt i})+\sum_{i=1}^\sqrt n\mathcal O(\frac{\sqrt\frac ni}{\log\sqrt\frac ni})\\ \]

显然瓶颈在后半部分，将分母视为 \(\log n\) ，对分子积分：

\[T_1(n)=\frac{\sqrt n}{\log n}\mathcal O(\sum_{i=1}^\sqrt ni^{-\frac 12}) =\frac{\sqrt n}{\log n}\mathcal O(\sqrt x\mid_0^\sqrt n) =\mathcal O(\frac{n^{0.75}}{\log n})\\ \]

因此素数部分的时间复杂度为 \(\mathcal O(\frac{n^{0.75}}{\log n})\) 。

二、合数部分

定义：

\[S(n,j)=\sum_{i=2}^n[\text{minp}(i)\gt p_j]f(i)\\ \]

即最小素因子严格大于 \(p_j\) 的所有数的贡献之和。

显然最终答案为 \(f(1)+S(n,0)\) 。

转移分别计算素数和合数的贡献：

• 显然素数的贡献为 \(G(n)-s_j\) 。

• 计算合数的贡献需要枚举 \(i\) 的最小素因子 \(p\) ，及其幂次 \(e\) ：

\[\sum_{i=1}^n[\text{minp}(i)=p\and v_p(i)=e]f(i)\\ =f(p^e)\cdot\bigg([e\neq 1]+\sum_{i=1}^\frac n{p^e}[\text{minp}(i)\gt p]f(i)\bigg)\\ =f(p^e)\cdot\bigg([e\neq 1]+S(\frac n{p^e},id_p)\bigg)\\ \]

这里用到了\(f\) 为积性函数的条件，注意 \(p^e\) 与后面一大堆东西互素。

由于仅考虑合数，所以枚举上界为\(p\le\sqrt n\)。

完整转移方程如下：

\[S(n,j)=G(n)-s_j+\sum_{p_j<p_i\le\sqrt n}\sum_{e=1}^{\log_{p_i}n}f(p_i^e)\cdot\bigg(S(\frac n{p_i^e},i)+[e\neq1]\bigg) \]

注意代码实现时不需要记忆化。

时间复杂度分析：

当 \(n\to\infty\) 时，可以证明合数部分时间复杂度为 \(\mathcal O(n^{1-\epsilon})\) ，其中 \(\epsilon\) 为无穷小量。

但是当 \(n\le 10^{13}\) 时，可以证明时间复杂度为 \(\mathcal O(\frac{n^{0.75}}{\log n})\) ，证明参考 2018 年集训队论文 《一些特殊的数论函数求和问题》 。

三、相关例题

例1、\(\texttt{P5325 【模板】Min\_25筛}\)

题目描述

积性函数 \(f\) 满足 \(f(p^k)=p^k(p^k-1)\) ，求 \(\sum_{i=1}^nf(i)\) ，对 \(10^9+7\) 取模。

数据范围

• \(1\le n\le 10^{10}\) 。

时间限制 \(\texttt{2s}\) ，空间限制 \(\texttt{250MB}\) 。

分析

素数处的点值 \(f(p)=p(p-1)=p^2-p\) ，拆成完全积性函数 \(f_1(p)=p^2\) 和 \(f_2(p)=p\) 之差。

低次多项式都可以用这个做法，需要拆的原因是多项式不满足完全积性，但单项式满足。

时间复杂度 \(\mathcal O(\frac{n^{0.75}}{\log n})\) 。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=2e5+5,mod=1e9+7,inv2=(mod+1)/2,inv6=(mod+1)/6;///maxn=2\sqrt n
int B,m,n,cnt;
int b[maxn],p[maxn];
int id1[maxn],id2[maxn],val[maxn];
int g1[maxn],g2[maxn],s1[maxn],s2[maxn];
void init(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
    for(int i=1;i<=cnt;i++)
    {
        s1[i]=(s1[i-1]+p[i])%mod;
        s2[i]=(s2[i-1]+p[i]*p[i])%mod;
    }
}
int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
int f(int x)
{
    x%=mod;
    return x*(x-1)%mod;
}
int S(int n,int j)
{
    if(p[j]>=n) return 0;
    int k=id(n),res=((g2[k]-g1[k])-(s2[j]-s1[j]))%mod;
    for(int i=j+1;i<=cnt&&p[i]*p[i]<=n;i++)
        for(int e=1,x=p[i];x<=n;e++,x*=p[i])
            res=(res+1ll*f(x)*(S(n/x,i)+(e!=1)))%mod;
    return res;
}
signed main()
{
    scanf("%lld",&n),init(B=sqrt(n));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        x%=mod;
        g1[m]=(x*(x+1)%mod*inv2-1)%mod;
        g2[m]=(x*(x+1)%mod*(2*x+1)%mod*inv6-1)%mod;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)///val[]从大到小存储整除分块的O(\sqrt n)个值
        {
            int k=id(val[i]/p[j]);
            g1[i]=(g1[i]-p[j]*(g1[k]-s1[j-1]))%mod;
            g2[i]=(g2[i]-p[j]*p[j]%mod*(g2[k]-s2[j-1]))%mod;
        }
    printf("%lld\n",(1+S(n,0)+mod)%mod);
    return 0;
}

例2、\(\texttt{LOJ6235 区间素数个数}\)

题目描述

给定 \(n\) ，求 \(1\sim n\) 中素数个数。

数据范围

• \(2\le n\le 10^{11}\) 。

时间限制 \(\texttt{1s}\) ，空间限制 \(\texttt{512MB}\) 。

分析

本题其实比模板题还要简单，令 \(f'(p)=1\) 即可，甚至都不需要合数部分的操作。

时间复杂度 \(\mathcal O(\frac{n^{0.75}}{\log n})\) 。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=7e5+5;
int B,m,n,cnt;
int b[maxn],p[maxn];
int id1[maxn],id2[maxn],val[maxn];
int g[maxn];
void init(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
}
int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
signed main()
{
    scanf("%lld",&n),init(B=sqrt(n));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        g[m]=x-1;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
            g[i]-=g[id(val[i]/p[j])]-(j-1);
    printf("%lld\n",g[id(n)]);
    return 0;
}

例3、\(\texttt{LOJ6053 简单的函数}\)

题目描述

积性函数 \(f\) 满足 \(f(p^c)=p\oplus c\) ，求 \(\sum_{i=1}^nf(i)\) ，对 \(10^9+7\) 取模。

数据范围

• \(1\le n\le 10^{10}\) 。

时间限制 \(\texttt{2s}\) ，空间限制 \(\texttt{256MB}\) 。

分析

在素数部分把 \(f(2)\) 当成 \(1\) （注意合数部分无影响）求前缀和，如果 \(n\ge 2\) 再补上 \(n=2\) 的贡献。

因此 \(f(p)=p-1\) ，拆成 \(f_1(p)=p\) 和 \(f_2(p)=1\) 之差即可。

时间复杂度 \(\mathcal O(\frac{n^{0.75}}{\log n})\) 。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=2e5+5,mod=1e9+7,inv2=(mod+1)/2;
int B,m,n,cnt;
int b[maxn],p[maxn];
int id1[maxn],id2[maxn],val[maxn];
int g1[maxn],g2[maxn],s1[maxn],s2[maxn];
void init(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
    for(int i=1;i<=cnt;i++) s1[i]=(s1[i-1]+p[i])%mod,s2[i]=i;
}
int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
int S(int n,int j)
{
    if(p[j]>=n) return 0;
    int k=id(n),res=((g1[k]-g2[k])-(s1[j]-s2[j]))%mod;
    for(int i=j+1;i<=cnt&&p[i]*p[i]<=n;i++)
        for(int e=1,x=p[i];x<=n;e++,x*=p[i])
            res=(res+(p[i]^e)*(S(n/x,i)+(e!=1)))%mod;
    return res;
}
signed main()
{
    scanf("%lld",&n),init(B=sqrt(n));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        x%=mod;
        g1[m]=(x*(x+1)%mod*inv2-1)%mod,g2[m]=x-1;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
        {
            int k=id(val[i]/p[j]);
            g1[i]=(g1[i]-p[j]*(g1[k]-s1[j-1]))%mod;
            g2[i]=(g2[i]-(g2[k]-s2[j-1]))%mod;
        }
    printf("%lld\n",(1+S(n,0)+2*(n>=2)+mod)%mod);
    return 0;
}

例4、\(\texttt{UOJ188 【UR \#13】Sanrd}\)

题目描述

定义 \(f(x)\) 为 \(x\) 的非严格次大质因子，规定 \(f(1)=f(p)=0\) 。

给定 \(l,r\) ，求 \(\sum_{i=l}^rf(i)\) 。

数据范围

• \(1\le l\le r\le 10^{11},l+r\le 10^{11}\) 。

时间限制 \(\texttt{5s}\) ，空间限制 \(\texttt{512MB}\) 。

分析

模仿 \(\min25\) 筛合数部分的 \(dfs\) ，定义 \(S(n,j)=\sum_{i=1}^n[\text{minp(i)}\gt p_j]f(i)\) 。

枚举最小素因子 \(p_i\) 及其次数 \(e\) ，有两种情况：

• 如果除掉 \(p_i^e\) 后仍为合数，可以递归算。

• 如果除掉 \(p_i^e\) 后为素数，贡献为 \(p_i\) 。

  这种情况的出现次数为 \(|[p_i,\frac n{p_i^e}]\cap\text{prime}|=G(\frac n{p_i^e})-(i-1)\) ，其中 \(G(x)\) 表示\(1\sim x\)中素数个数，可以用 \(\min25\) 筛预处理。

  注意剩下的素数可以为 \(p_i\) ， \(f(p_i^e)\) 的贡献会在枚举 \(p_i^{e-1}\) 时被计算。

转移方程如下：

\[S(n,j)=\sum_{p_j\lt p_i\le\sqrt n}\sum_{e=1}^{\log_{p_i}n}\bigg(S(\frac n{p_i^e},i)+p_i\cdot\max\big(G(\frac n{p_i^e})-i+1,0\big)\bigg) \]

时间复杂度 \(O(\frac{r^{0.75}}{\log r})\) 。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=7e5+5;
int B,l,m,n,r,cnt;
int b[maxn],p[maxn];
int id1[maxn],id2[maxn],val[maxn];
int g[maxn];
void init(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
}
int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
int S(int n,int j)
{
    if(p[j]>=n) return 0;
    int res=0;
    for(int i=j+1;i<=cnt&&p[i]*p[i]<=n;i++)
        for(int e=1,x=p[i];x<=n;e++,x*=p[i])
            res+=S(n/x,i)+p[i]*max(g[id(n/x)]-i+1,0ll);
    return res;
}
int calc(int _n)
{
    n=_n,m=cnt=0,init(B=sqrt(n));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        g[m]=x-1;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
            g[i]-=g[id(val[i]/p[j])]-(j-1);
    return S(n,0);
}
signed main()
{
    scanf("%lld%lld",&l,&r);
    printf("%lld\n",calc(r)-calc(l-1));
    return 0;
}

例5、\(\texttt{LOJ572 「LibreOJ Round \#11」Misaka Network 与求和}\)

题目描述

定义 \(f(x)\) 为 \(x\) 的非严格次大质因子，规定 \(f(1)=0,f(p)=1\) ，求 \(\sum_{i=1}^n\sum_{j=1}^nf^k(\gcd(i,j))\) ，对 \(2^{32}\) 取模。

数据范围

• \(1\le k,n\le 2\cdot 10^9\) 。

时间限制 \(\texttt{4s}\) ，空间限制 \(\texttt{512MB}\) 。

分析

先套路莫反：

\[\begin{aligned} &\sum_{i=1}^n\sum_{j=1}^nf^k(\gcd(i,j))\\ =&\sum_{d=1}^nf^k(d)\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]\\ =&\sum_{d=1}^nf^k(d)\sum_{d|x}\mu(\frac xd)\lfloor\frac nx\rfloor^2\\ =&\sum_{x=1}^n\lfloor\frac nx\rfloor^2\sum_{d|x}f^k(d)\mu(\frac xd)\\ \end{aligned} \]

记 \(F(x)=\sum_{d|x}f(d)\mu(\frac xd)\) ，外层整除分块后内层要求 \(F\) 的前缀和。

杜教筛， \(F*1=f*\mu*1=f*\epsilon=f\) 。

\(1\) 的前缀和可以直接求， \(f\) 的前缀和需要用 \(\min25\) 筛，具体做法和上题类似。

计算 \(S(n,j)\) 时先不考虑素数的贡献（否则递推式会变的复杂很多），最后加上 \([1,n]\) 中素数个数即可。

瓶颈在于杜教筛没法预处理，时间复杂度 \(\mathcal O(n^{0.75})\) 。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=9e4+5;
int B,k,m,n,cnt,res;
int b[maxn],p[maxn];
int id1[maxn],id2[maxn],val[maxn];
int g[maxn],mp[maxn],pw[maxn];
bool vis[maxn];
inline int qpow(int a,int k)
{
    int res=1;
    while(k)
    {
        if(k&1) res=res*a;
        a=a*a,k>>=1;
    }
    return res;
}
inline void init(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!b[i]) p[++cnt]=i;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
    for(int i=1;i<=n;i++) pw[i]=qpow(i,k);
}
inline int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
inline int S(int n,int j)
{
    if(p[j]>=n) return 0;
    int res=0;
    for(int i=j+1;i<=cnt&&p[i]*p[i]<=n;i++)
        for(int e=1,x=p[i];x<=n;e++,x*=p[i])
            res+=S(n/x,i)+pw[p[i]]*max(g[id(n/x)]-i+1,0ll);
    return res;
}
inline int sum(int n)
{
    if(vis[id(n)]) return mp[id(n)];
    vis[id(n)]=true;
    int res=S(n,0)+g[id(n)];
    for(int l=2,r=0;l<=n;l=r+1) r=n/(n/l),res-=(r-l+1)*sum(n/l);
    return mp[id(n)]=res;
}
signed main()
{
    scanf("%lld%lld",&n,&k),init(B=sqrt(n));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        g[m]=x-1;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
            g[i]-=g[id(val[i]/p[j])]-(j-1);
    for(int l=1,r=0;l<=n;l=r+1) r=n/(n/l),res+=(n/l)*(n/l)*(sum(r)-sum(l-1));
    printf("%lld\n",res&(-1u));
    return 0;
}

例6、\(\texttt{UOJ448 【集训队作业2018】人类的本质}\)

题目描述

给定 \(n,k\) ，求：

\[\sum_{i=1}^n\sum_{1\le x_j\le i}\text{lcm}(\gcd(i,x_1),\cdots,\gcd(i,x_k)) \]

对 \(10^9+7\) 取模。

数据范围

• \(n,k\ge 1,n\cdot k\le 10^9\) 。

时间限制 \(\texttt{2s}\) ，空间限制 \(\texttt{512MB}\) 。

分析

固定 \(i\) ，先化简内层式子：

\[\sum_{1\le x_j\le i}\text{lcm}(\gcd(i,x_1),\cdots,\gcd(i,x_k))=\sum_{1\le x_j\le i}\frac i{\gcd(\frac i{\gcd(i,x_1)},\cdots,\frac i{\gcd(i,x_k)})}\\ \]

枚举 \(y_j=\frac i{\gcd(i,x_j)}\) ，对应的 \(x_j\) 有多少个？

\[y_j=\frac i{\gcd(i,x_j)}\iff\gcd(i,x_j)=\frac i{y_j}\iff\gcd(y_j,x_j\cdot\frac{y_j}i)=1\\ \]

由 \(1\le x_i\le i\) ，知 \(1\le x_j\cdot\frac{y_j}i\le y_j\) ，因此满足条件的 \(x_j\) 个数为 \(\varphi(y_j)\) 。

于是原式可以化简为：

\[\sum_{i=1}^n\sum_{y_j|i}\bigg(\varphi(y_1)\cdots\varphi(y_k)\frac i{\gcd(y_1,\cdots,y_k)} \bigg) \]

---

记\(f(x)=\sum\limits_{y_j|x}\bigg(\varphi(y_1)\cdots\varphi(y_k)\frac x{\gcd(y_1,\cdots,y_k)}\bigg)\)，目标变为求\(f\)的前缀和。

本题最有迷惑性的地方在于，看到\(\gcd\)很容易往反演上面想，但如果真的去反演就意味着走进死胡同了。

关键性质：\(f(x)\)是积性函数！

证明：

\(\forall\gcd(a,b)=1,\forall y_j|a,z_j|b\)，对\(f(a)\cdot f(b)\)的贡献为：

\[\varphi(y_1)\cdots\varphi(y_k)\frac a{\gcd(y_1,\cdots,y_k)} \cdot \varphi(z_1)\cdots\varphi(z_k)\frac b{\gcd(z_1,\cdots,z_k)} \]

对\(f(a\cdot b)\)的贡献为：

\[\varphi(y_1\cdot z_1)\cdots\varphi(y_k\cdot z_k)\frac{ab}{\gcd(y_1z_1,\cdots,y_kz_k)} \]

由\(\gcd(a,b)=1\)可知上面两式相同。

于是本题等价于求积性函数\(f(x)\)的前缀和！

---

注意到 \(f(x)\) 是积性函数，所以只需关注如何求 \(f(p^c)\) 。

枚举 \(y_j=p^{x_j}\) ：

\[\begin{aligned} f(p^c)&=\sum_{x_j=0}^c\varphi(p^{x_1})\cdots\varphi(p^{x_k})\cdot p^{c-\min(x_1,\cdots,x_k)}\\ &=\sum_{x_j=0}^c\varphi(p^{c-x_1})\cdots\varphi(p^{c-x_k})\cdot p^{\max(x_1,\cdots,x_k)}\\ \end{aligned} \]

定义 \(g(i)\) 为 \(\max x_j\le i\) 时的求和结果：

\[\begin{aligned} g(i)&=\sum_{x_j=0}^i\varphi(p^{c-x_1})\cdots\varphi(p^{c-x_k})\\ &=\big(\sum_{x=0}^i\varphi(p^{c-x})\big)^k\\ &=\begin{cases} (p^c-p^{c-i-1})^k&i\neq c\\ (p^c)^k&i=c\\ \end{cases} \end{aligned} \]

然后可以化简 \(f(p^c)\) 了：

\[\begin{aligned} f(p^c)&=\sum_{i=0}^cp^i(g(i)-g(i-1))\\ &\color{red}{=p^c(p^{ck}-(p^c-1)^k)+\sum_{i=0}^{c-1}p^i\bigg((p^c-p^{c-i-1})^k-(p^c-p^{c-i})^k\bigg)}\\ \end{aligned} \]

注意 \(g(-1)=0\) ，所以不必单独考虑。

考虑用 \(f(p^c)\) 递推 \(f(p^{c+1})\) ：

\[\begin{aligned} f(p^{c+1})&=p^{c+1}(p^{(c+1)k}-(p^{c+1}-1)^k)+p^k\bigg(f(p^c)-p^c(p^{ck}-(p^c-1)^k)\bigg) +p^c\bigg((p^{c+1}-p^{(c+1)-c-1})^k-(p^{c+1}-p^{(c+1)-c})^k\bigg)\\ &=p^kf(p^c)+p^{(c+1)(k+1)}-p^{c+1}(p^{c+1}-1)^k-p^{ck+c+k}+p^{c+k}(p^c-1)^k +p^c(p^{c+1}-1)^k-p^c(p^{c+1}-p)^k\\ &=p^kf(p^c)+p^{ck+c+k}(p-1)-p^c(p-1)(p^{c+1}-1)^k\\ &=p^kf(p^c)+p^c(p-1)\cdot\bigg(p^{(c+1)k}-(p^{c+1}-1)^k\bigg)\\ \end{aligned} \]

初始值可以通过在上式中代入 \(c=0\) 得到：

\[\begin{aligned} f(p)&=p^k+(p-1)\cdot(p^k-(p-1)^k)\\ &=p^{k+1}-(p-1)^{k+1}\\ &=\sum_{i=0}^k\binom{k+1}i(-1)^{k-i}p^i\\ \end{aligned} \]

---

至此 \(n\le B\) 时已经可以线性筛了，时间复杂度 \(\mathcal O(n)\) 。

在线性筛的同时维护 \(x^k\) ，可以省掉快速幂的一只 \(\log\) 。

如果 \(n\gt B\) ，由于 \(k\cdot n\le 10^9\) ，所以 \(k\) 很小。

考虑 \(\min25\) 筛，素数部分 \(\texttt{dp}\) 初始值需要用到自然数幂和，可以用斯特林数或拉格朗日插值预处理。

时间复杂度 \(\mathcal O(k\cdot\frac{n^{0.75}}{\log n})\) 。

实测取 \(B=1.5\cdot 10^7\) 时效率最高。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=1.5e7+5,mod=1e9+7;
int k,n;
inline int qpow(int a,int k)
{
    int res=1;
    while(k)
    {
        if(k&1) res=res*a%mod;
        a=a*a%mod,k>>=1;
    }
    return res;
}
namespace task1
{
    int cnt,res;
    int b[maxn],p[maxn],pw[maxn];
    int f[maxn],g[maxn];///g[i]表示i的最小素因子及其幂次
    void init(int n)
    {
        pw[1]=1;
        for(int i=2;i<=n;i++)
        {
            if(!b[i]) p[++cnt]=i,pw[i]=qpow(i,k);
            for(int j=1;j<=cnt&&i*p[j]<=n;j++)
            {
                b[i*p[j]]=1,pw[i*p[j]]=pw[i]*pw[p[j]]%mod;
                if(i%p[j]==0) break;
            }
        }
    }
    void solve()
    {
        init(n);
        for(int i=2;i<=n;i++) b[i]=0;
        f[1]=g[1]=1;
        for(int i=2;i<=n;i++)
        {
            if(!b[i]) f[i]=(i*pw[i]-(i-1)*pw[i-1])%mod,g[i]=i;
            for(int j=1;j<=cnt&&i*p[j]<=n;j++)
            {
                int x=i*p[j];
                b[x]=1;
                if(i%p[j]==0)
                {
                    g[x]=g[i]*p[j];
                    if(x==g[x]) f[x]=(pw[p[j]]*f[i]+i*(p[j]-1)%mod*(pw[x]-pw[x-1]))%mod;
                    else f[x]=f[x/g[x]]*f[g[x]]%mod;
                    break;
                }
                f[x]=f[i]*f[p[j]]%mod,g[x]=p[j];
            }
        }
        for(int i=1;i<=n;i++) res+=f[i];
        printf("%lld\n",(res%mod+mod)%mod);
    }
}
namespace task2
{
    const int maxn=7e4+5;
    int B,m,cnt;
    int b[maxn],p[maxn];
    int fac[maxn],inv[maxn];
    int id1[maxn],id2[maxn],val[maxn];
    int f[maxn][30];
    int g[maxn][70],s[maxn][70],pw[maxn][70];
    int pre[maxn],suf[maxn];
    int G[maxn],sum[maxn];///G[i]表示1~val[i]素数贡献,sum[i]表示素数点值前缀和
    void init(int n)
    {
        for(int i=2;i<=n;i++)
        {
            if(!b[i]) p[++cnt]=i;
            for(int j=1;j<=cnt&&i*p[j]<=n;j++)
            {
                b[i*p[j]]=1;
                if(i%p[j]==0) break;
            }
        }
        fac[0]=1;
        for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
        inv[n]=qpow(fac[n],mod-2);
        for(int i=n-1;i>=0;i--) inv[i]=inv[i+1]*(i+1)%mod;
        for(int i=1;i<=n;i++)
        {
            pw[i][0]=1;
            for(int j=1;j<=k;j++) pw[i][j]=pw[i][j-1]*i%mod;
        }
        for(int i=1;i<=cnt;i++)
            for(int j=0;j<=k;j++)
                s[i][j]=(s[i-1][j]+pw[p[i]][j])%mod;
    }
    int c(int n,int m)
    {
        if(m<0||n<m) return 0;
        return 1ll*fac[n]*inv[m]%mod*inv[n-m]%mod;
    }
    int id(int x)
    {
        return x<=B?id1[x]:id2[n/x];
    }
    int calc(int n,int k)
    {
        static int sum[70];
        int m=k+2,res=0;
        pre[0]=suf[m+1]=1;
        for(int i=1;i<=m;i++) sum[i]=(sum[i-1]+pw[i][k])%mod;
        for(int i=1;i<=m;i++) pre[i]=pre[i-1]*(n-i)%mod;
        for(int i=m;i>=1;i--) suf[i]=suf[i+1]*(n-i)%mod;
        for(int i=1;i<=m;i++) res=(res+(m-i&1?mod-1:1)*sum[i]%mod
            *pre[i-1]%mod*suf[i+1]%mod*inv[i-1]%mod*inv[m-i])%mod;
        return res;
    }
    int S(int n,int j)
    {
        if(p[j]>=n) return 0;
        int res=G[id(n)]-sum[j];
        for(int i=j+1;i<=cnt&&p[i]*p[i]<=n;i++)
            for(int e=1,x=p[i];x<=n;e++,x*=p[i])
                res=(res+f[i][e]*(S(n/x,i)+(e!=1)))%mod;
        return res;
    }
    void solve()
    {
        init(B=sqrt(n));
        for(int i=1;i<=cnt;i++)///预处理f(p^c)
        {
            f[i][0]=1;
            for(int e=1,x=p[i];x<=n;e++,x*=p[i])
                f[i][e]=(pw[p[i]][k]*f[i][e-1]+x/p[i]*(p[i]-1)*(qpow(x,k)-qpow(x-1,k)))%mod;
            sum[i]=(sum[i-1]+f[i][1])%mod;
        }
        for(int l=1,r=0;l<=n;l=r+1)
        {
            int x=n/l;
            r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
            for(int i=0;i<=k;i++) g[m][i]=calc(x,i)-1;
        }
        for(int j=1;j<=cnt;j++)
            for(int i=1;p[j]*p[j]<=val[i];i++)
                for(int l=0;l<=k;l++)
                    g[i][l]=(g[i][l]-pw[p[j]][l]*(g[id(val[i]/p[j])][l]-s[j-1][l]))%mod;
        for(int i=1;i<=m;i++)
            for(int l=0;l<=k;l++)
                G[i]=(G[i]+(k-l&1?-1:1)*c(k+1,l)*g[i][l])%mod;
        printf("%lld\n",(1+S(n,0)+mod)%mod);
    }
}
signed main()
{
    scanf("%lld%lld",&n,&k);
    if(n<=maxn-5) task1::solve();
    else task2::solve();
    return 0;
}