题面

题解

肝了整整一天……膜拜和巨巨……（虽然它们的题解里我就没看懂几个字）

请备好草稿纸和笔，这种题目就是需要耐心推倒

题目所求是这么一个东西

\[ \begin{aligned} ans &=\sum_{i=1}^n\sum_{x_1=1}^i\sum_{x_2=1}^i...\sum_{x_k=1}^ilcm(\gcd(i,x_1),\gcd(i,x_2),...,\gcd(i,x_k))\\ &=\sum_{i=1}^n\sum_{x_1=1}^i\sum_{x_2=1}^i...\sum_{x_k=1}^i\frac{i}{\gcd(\frac{i}{\gcd(i,x_1)},\frac{i}{\gcd(i,x_2)},\ldots,\frac{i}{\gcd(i,x_k)})}\\ \end{aligned} \]

对于每一个\(y_j={i\over\gcd(i,x_j)}\)来说，我们考虑它会出现多少次

\[ \begin{aligned} {i\over\gcd(i,x_j)}=y_j\\ \gcd(i,x_j)={i\over y_j}\\ \gcd(y_j,{x_jy_j\over i})=1\\ \end{aligned} \]

然后我们发现，满足条件的\(x_j\)恰好有\(\varphi(y_j)\)个

那么答案就可以改成

\[ ans =\sum_{i=1}^n\sum_{y_1,y_2,...,y_k|i}\varphi(y_1)\varphi(y_2)...\varphi(y_k)\frac{i}{\gcd(y_1,y_2,\ldots,y_k)} \]

我们记

\[f_i=\sum_{y_1,y_2,...,y_k|i}\varphi(y_1)\varphi(y_2)...\varphi(y_k)\frac{i}{\gcd(y_1,y_2,\ldots,y_k)}\]

那么这道题就变成要我们对\(f_i\)求和了

首先我们可以发现，\(f_i\)是一个积性函数

考虑如何证明。设\(a,b\)互质，那么\(a\)的任何一个因子和\(b\)的任何一个因子也互质，于是有

\[ \begin{aligned} f_af_b &=\sum_{y_1,y_2,...,y_k|a}\sum_{z_1,z_2,...,z_k|b}\varphi(y_1)\varphi(y_2)...\varphi(y_k)\varphi(z_1)\varphi(z_2)...\varphi(z_k)\frac{a}{\gcd(y_1,y_2,\ldots,y_k)}\frac{b}{\gcd(z_1,z_2,\ldots,z_k)}\\ &=\sum_{y_1,y_2,...,y_k|a}\sum_{z_1,z_2,...,z_k|b}\varphi(y_1z_1)\varphi(y_2z_2)...\varphi(y_kz_k){ab\over \gcd(y_1z_1,y_2z_2,...,y_kz_k)}\\ &=\sum_{y_1z_1,y_2z_2,...,y_kz_k|ab}\varphi(y_1z_1)\varphi(y_2z_2)...\varphi(y_kz_k){ab\over \gcd(y_1z_1,y_2z_2,...,y_kz_k)}\\ &=f_{ab}\\ \end{aligned} \]

可以这么理解，就是每一个\(ab\)的因子分解成\(a\)的因子和\(b\)的因子的乘积的方法是唯一的，所以可以不重不漏地枚举所有\(ab\)的因子，而且\(y\)和\(z\)全部互质，所以两者积的\(\gcd\)就等于两者\(\gcd\)的积

那么我们需要能够快速求出\(f(p^c)\)

\[ \begin{aligned} f(p^c) &=\sum_{x_1,x_2,\ldots,x_k=0}^{c}p^{\max(x_1,x_2,\ldots,x_k)}\prod_{j=1}^k\varphi(p^{c-x_j})\\ &=\sum_{i=0}^cp^i(g(i)-g(i-1))\\ \end{aligned} \]

其中\(g(i)\)表示\(\max(x_1,x_2,...,x_k)\)小于等于\(i\)的总方案数

显然有

\[ g(i)= \begin{cases} p^{ck}&,i=c\\ (p^c-p^{c-1-i})^k&,0\leq i< c \end{cases} \]

然后再来考虑一下该怎么递推

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

这个直接根据公式就可以得出转移了

那么当\(n\)范围很小的时候我们就可以直接递推了

然而如果\(n\)很大怎么办？

这就是个积性函数求和啦！

显然，除非脑子被咱踢了，否则这种奇怪的式子不可能用杜教筛搞的，我们只能考虑\(Min\_25\)筛了

然而\(Min\_25\)筛需要\(f(p)\)是一个完全积性函数，或者可以拆成若干个完全积性函数的和，但现在的\(f(p)\)是个啥呢？

\[ \begin{aligned} f(p) &=g(0)+pg(1)-pg(0)\\ &=(1-p)(p-1)^{k+1}-p^{k+1}\\ &=p^{k+1}-(p-1)^{k+1} \end{aligned} \]

这式子怎么看都不像个完全积性函数啊……

我们考虑把\((p-1)^{k+1}\)用二项式定理展开一下，等于\(\sum_{i=0}^{k+1}(-1)^{k+1-i}p^i{k+1\choose i}\)，然后用\(p^{k+1}\)减去它，发现

\[f(p)=-\sum_{i=0}^k(-1)^{k+1-i}p^i{k+1\choose i}=\sum_{i=0}^k(-1)^{k-i}p^i{k+1\choose i}\]

那么\(f(p)\)就可以拆成\(k+1\)个完全积性函数的和了。我们每一次枚举\(i\)，把记\(f(p)\)的值为\((-1)p^i\)，拉格朗日插值求一下前缀和作为初始值，算出答案之后再乘上后面的组合数加到\(g\)数组里就行了

然而\(Min\_25\)筛还需要能够快速求出\(f(p^c)\)的值，这里没有好的办法，只能暴力了

实测当\(n\leq 5000000\)的时候用方法\(1\)，其他情况用方法\(2\)求和复杂度最优

明明是个完全抄题解的却因为调了一下啥时候改用方法\(1\)啥时候该用方法\(2\)跑到了\(rank\ 1\)我有点慌

//minamoto#include
    
     #define R register#define ll long long#define pi pair
     
      #define IT map
      
       ::iterator#define fp(i,a,b) for(R int i=a,I=b+1;i
       
        I;--i)#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)using namespace std;const int P=1e9+7;inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}int ksm(R int x,R ll y){    R int res=1;    for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);    return res;}int n,k;namespace chino{    const int N=1e7+5;    bitset
        
         vis;int p[N/3+5],pw[N],f[N],c[N],tot,v,res;    void MAIN(){        f[1]=pw[1]=1;        fp(i,2,n){            if(!vis[i])p[++tot]=i,pw[i]=ksm(i,k);            for(R int j=1;j<=tot&&1ll*i*p[j]<=n;++j){                vis[v=i*p[j]]=1,pw[v]=mul(pw[i],pw[p[j]]);                if(i%p[j]==0)break;            }        }        fp(i,2,n){            if(!vis[i])c[i]=i,f[i]=dec(ksm(i,k+1),ksm(i-1,k+1));            for(R int j=1;j<=tot&&1ll*i*p[j]<=n;++j){                v=i*p[j];                if(i%p[j]==0){                    c[v]=c[i]*p[j];                    f[v]=(c[v]==v)?add(mul(pw[p[j]],f[i]),1ll*i*(p[j]-1)%P*dec(pw[v],pw[v-1])%P):mul(f[c[v]],f[v/c[v]]);                    break;                }                c[v]=p[j],f[v]=mul(f[i],f[p[j]]);            }        }        fp(i,1,n)res=add(res,f[i]);        printf("%d\n",res);    }}namespace yosino{    #define ID(x) (x<=sqr?id1[x]:id2[n/(x)])    const int N=1e5+5;    bitset
         
          vis;int p[N],ans[N],g[N],id1[N],id2[N],sp[N]; int fac[N],ifac[N],w[N],sum[N],Pre[N],suf[N],s[N],pw[N]; int tot,cnt,sqr,m;map
          
           mp; inline int C(R int n,R int m){return n
           
            second; int res=0,las=0,now; fp(i,0,k)now=G(i,p,k),res=add(res,mul(ksm(p,i),dec(now,las))),las=now; return mp[pi(p,k)]=res; } int S(int x,int y){ if(x<=1||p[y]>x)return 0; int res=dec(ans[ID(x)],sum[y-1]); for(R int i=y;i<=tot&&1ll*p[i]*p[i]<=x;++i) for(R int e=1,pr=p[i],las=F(p[i],1),now;1ll*pr*p[i]<=x;++e,pr*=p[i],las=now) now=F(p[i],e+1),res=add(res,add(mul(S(x/pr,i+1),las),now)); return res; } void MAIN(){ sqr=sqrt(n),init(),pw[1]=1; for(R int i=1,j;i<=n;i=j+1)j=n/(n/i),w[++m]=n/i,ID(w[m])=m; fp(i,0,k){ solve(i);int v=((k-i)&1)?P-C(k+1,i):C(k+1,i); fp(j,1,m)ans[j]=add(ans[j],mul(v,g[j])); } printf("%d\n",S(n,1)+1); }}int main(){// freopen("testdata.in","r",stdin); scanf("%d%d",&n,&k); n<=5e6?chino::MAIN():yosino::MAIN(); return 0;}