集训队作业2018人类的本质

题目链接

UOJ

代码参考了yww的

sol

\[ans=\sum_{i=1}^n\sum_{x_1,x_2,\ldots,x_m=1}^i\operatorname{lcm}(\gcd(i,x_1),\gcd(i,x_2),\ldots,\gcd(i,x_m))\\ =\sum_{i=1}^n\sum_{x_1,x_2,\ldots,x_m=1}^i\frac{i}{\gcd(\frac{i}{\gcd(i,x_1)},\frac{i}{\gcd(i,x_2)},\ldots,\frac{i}{\gcd(i,x_m)})}\\\]

然后考虑\(y_j=\frac i{gcd(i,x_j)}\)，容易发现它出现的次数等于\(\varphi(y_j)\)。（左式简单变换即可

即求
\[ans=\sum_{i=1}^n\sum_{y_1,y_2,\ldots,y_m\mid i}\varphi(y_1)\varphi(y_2)\cdots\varphi(y_m)\frac{i}{\gcd(y_1,y_2,\ldots,y_m)}\]

令\(f(n)=\sum_{y_1,y_2,\ldots,y_m\mid n}\varphi(y_1)\varphi(y_2)\cdots\varphi(y_m)\frac{n}{\gcd(y_1,y_2,\ldots,y_m)}\)

然后这好像就是个积性函数？讲真我看不出来有谁会证评论教教我好不好呀\(QwQ\)。

下面我们假装它是积性函数。

实际上我们只要会算\(f(p^k)\)就行了，和他们一样，我们暂时用\(d\)来代替\(k\)。

\[f(p^k)=\sum p^{max(x1,x2,...,xd)}\prod\varphi(p^{k-x_i})\]

设\(mx\)为最大指数
\[g(mx)=\sum_{i=0}^{mx}(\varphi(x^{k-i}))^d\\ =\begin{cases} p^{kd}&,mx=1\\ (p^k-p^{mx-1})^d&,2\leq mx\leq k \end{cases}\]

则
\[f(p^k)=\sum_{mx=0}^kp^{mx}(g(mx)-g(mx-1))\\ f(p^{k+1})=p^df(p^k)-p^k(p^{(k+1)d}-(p^{k+1}-1)^d)+p^{k+1}(p^{(k+1)d}-(p^{k+1}-1)^d)\\ =p^df(p^k)+p^k(p-1)(p^{(k+1)d}-(p^{k+1}-1)^d)\]
就是考虑\(g(mx)\)中的\(k\)的变化，暴力推。

现在我们可以做\(n\le10^7\)了

关于\(n>10^7\)，有\(k<100\)于是我们想到\(Min\_25\)筛。

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

\(f(p^k)\)暴力算。

然后就可以\(Min\_25\)筛了。

关于这份代码，我一开始写斯特林数处理自然数幂和，莫名\(WA\)了，于是蒯了yww的插值做法，恩。

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define gt getchar()
#define ll long long
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
inline int in()
{
    int k=0;char ch=gt;
    while(ch<'-')ch=gt;
    while(ch>'-')k=k*10+ch-'0',ch=gt;
    return k;
}
const int YL=1e9+7;
inline int MO(const int &x){return x>=YL?x-YL:x;}
inline int OD(const int &x){return x<0?x+YL:x;}
inline int ksm(int a,ll k){int r=1;while(k){if(k&1)r=1ll*r*a%YL;a=1ll*a*a%YL,k>>=1;}return r;}
int n,k,d;
namespace w1
{
    const int N=1e7+5;bool np[N];
    int pr[N],tot,pw[N],f[N],c[N];
    void work()
    {
        f[1]=pw[1]=1;
        for(int i=2;i<=n;++i)
        {
            if(!np[i])pr[++tot]=i,pw[i]=ksm(i,k);
            for(int j=1;i*pr[j]<=n;++j)
            {
                np[i*pr[j]]=1;pw[i*pr[j]]=1ll*pw[i]*pw[pr[j]]%YL;
                if(i%pr[j]==0)break;
            }
        }
        for(int i=2;i<=n;++i)
        {
            if(!np[i])c[i]=i,f[i]=OD(ksm(i,k+1)-ksm(i-1,k+1));
            for(int j=1;i*pr[j]<=n;++j)
            {
                int v=i*pr[j];
                if(i%pr[j]==0)
                {
                    c[v]=c[i]*pr[j];
                    if(c[v]==v)
                        f[v]=(1ll*pw[pr[j]]*f[i]%YL+1ll*i*(pr[j]-1)%YL*OD(pw[v]-pw[v-1])%YL)%YL;
                    else f[v]=1ll*f[c[v]]*f[v/c[v]]%YL;break;
                }
                c[v]=pr[j],f[v]=1ll*f[i]*f[pr[j]]%YL;
            }
        }
        ll ans=0;for(int i=1;i<=n;++i)ans+=f[i];
        printf("%lld\n",ans%YL);
    }
}

namespace w2
{
    const int N=2e5+5,M=N-5;bool np[N];
    int pr[N],tot,pw[N],ans[N],g[N],id1[N],id2[N],sp[N];
    int fnv[N],fac[N],s[N],w[N],m,S[105][105],sum[N];
    inline int C(int x,int y){return y>x?0:1ll*fac[x]*fnv[y]%YL*fnv[x-y]%YL;}
    void init()
    {
        fac[0]=1;
        for(int i=1;i<=M;++i)fac[i]=1ll*fac[i-1]*i%YL;fnv[M]=ksm(fac[M],YL-2);
        for(int i=M;i>=1;--i)fnv[i-1]=1ll*fnv[i]*i%YL;
        for(int i=2;i<=M;++i)
        {
            if(!np[i])pr[++tot]=i,sum[tot]=MO(sum[tot-1]+OD(ksm(i,k+1)-ksm(i-1,k+1)));
            for(int j=1;i*pr[j]<=M;++j)
            {np[i*pr[j]]=1;if(i%pr[j]==0)break;}
        }   
    }
    void pre_(int x)
    {
        int cnt=0;s[1]=1;
        for(int i=2;i<=M;++i)
        {
            if(!np[i])pw[i]=ksm(i,x),++cnt,sp[cnt]=MO(sp[cnt-1]+pw[i]);
            for(int j=1;i*pr[j]<=M;++j)
            {pw[i*pr[j]]=1ll*pw[i]*pw[pr[j]]%YL;if(i%pr[j]==0)break;}
            s[i]=MO(s[i-1]+pw[i]);
        }
    }
    int pre[N],suf[N];
    int calc(int x,int y)
    {
        if(x<=y+2)return s[x];
        pre[0]=1;
        for(int i=1;i<=y+2;i++)pre[i]=1ll*pre[i-1]*(x-i)%YL;
        suf[y+3]=1;
        for(int i=y+2;i>=1;i--)suf[i]=1ll*suf[i+1]*(x-i)%YL;
        ll res=0;
        for(int i=1;i<=y+2;i++)res=(res+1ll*s[i]*pre[i-1]%YL*suf[i+1]%YL*fnv[i-1]%YL*fnv[y+2-i]%YL*((y+2-i)&1?YL-1:1))%YL;
        return (res%YL+YL)%YL;
    }
#define ID(x) ((x)<=d?id1[x]:id2[n/(x)])
    void deal(int x)
    {
        pre_(x);
        for(int i=1;i<=m;++i)g[i]=((w[i]<=d)?s[w[i]]:calc(w[i],x))-1;
        for(int i=1;pr[i]<=n/pr[i];++i)
            for(int j=1;j<=m&&pr[i]<=w[j]/pr[i];++j)
                g[j]=OD(g[j]-1ll*pw[pr[i]]*OD(g[ID(w[j]/pr[i])]-sp[i-1])%YL);
    }
    inline int G(int z,int x,int y)
    {
        if(z<0)return 0;if(z==y)return ksm(x,1ll*y*k);
        return ksm(OD(ksm(x,y)-ksm(x,y-z-1)),k);
    }
    inline int Get(int x,int y)
    {
        int res=0;
        for(int i=0;i<=y;++i)
            res=(1ll*ksm(x,i)*OD(G(i,x,y)-G(i-1,x,y))+res)%YL;
        return res;
    }
    int solve(int x,int y)
    {
        if(x<=1||pr[y]>x)return 0;
        int res=OD(ans[ID(x)]-sum[y-1]);
        for(int i=y;i<=tot&&1ll*pr[i]*pr[i]<=x;++i)
            for(int e=1,p=pr[i];1ll*p*pr[i]<=x;++e,p*=pr[i])
                res=MO(res+(1ll*solve(x/p,i+1)*Get(pr[i],e)+Get(pr[i],e+1))%YL);
        return res;
    }
    void work()
    {
        d=sqrt(n);init();pw[1]=1;
        for(int i=1;i<=n;i=n/(n/i)+1)w[++m]=n/i,ID(w[m])=m;
        for(int i=0;i<=k;++i)
        {
            deal(i);int v=1ll*C(k+1,i)*((k-i)&1?YL-1:1)%YL;
            for(int j=1;j<=m;++j)ans[j]=(1ll*v*g[j]+ans[j])%YL;
        }
        int ans=solve(n,1)+1;printf("%d\n",ans);
    }
}
int main()
{
    n=in(),k=in();
    if(n<=1e7)w1::work();
    else w2::work();
    return 0;
}

转载于:https://www.cnblogs.com/cx233666/p/10194741.html