【UOJ448】【集訓隊做業2018】人類的本質 min_25篩

題目大意

　　給你 \(n,m\)，求
\[ \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)) \]
　　對 \({10}^9+7\) 取模。

　　\(nm\leq {10}^9\)

題解

　　先推一下式子：
\[ 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)}\)。對於一個 \(y_j\)，有多少個知足要求的 \(x_j\) 呢？
\[ \frac{i}{\gcd(i,x_j)}=y_j\\ \gcd(i,x_j)=\frac{i}{y_j}\\ \gcd(y_j,\frac{x_j}{y_j})=1\\ \]
　　因此有 \(\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)} \]
　　容易發現 \(f(n)\) 是積性函數。
\[ f(p^k)=\sum_{i=0}^kp^{k-i}\sum_{j=i}^k\mu(\frac{j}{i})g(j)\\ =\sum_{j=0}^kg(j)\sum_{i=0}^jp^{k-i}\mu(\frac{j}{i})\\ g(j)=\begin{cases} p^{km}&,j=1\\ (p^k-p^{j-1})^m&,2\leq j\leq k \end{cases}\\ \]

\[ f(p^k)=\sum_{x_1,x_2,\ldots,x_m=1}^{p^k}lcm(\gcd(p^k,x_1),\gcd(p^k,x_2),\ldots,\gcd(p^k,x_m))\\ =\sum_{x_1,x_2,\ldots,x_m=0}^{k}p^{\max(x_1,x_2,\ldots,x_m)}\prod_{j=1}^m\varphi(p^{k-x_j})\\ =\sum_{i=0}^kp^i(g(i)-g(i-1))\\ g(i)=\sum_{x_1,x_2,\dots,x_m=0}^i\prod_{j=1}^m\varphi(p^{k-x_j})\\ =\prod_{j=1}^m(\sum_{l=0}^i\varphi(p^{k-l}))\\= \begin{cases} p^{km}&,i=k\\ (p^k-p^{k-1-i})^m&,0\leq i<k \end{cases}\\ f(p^{k+1})=p^mf(p^k)-p^k(p^{(k+1)m}-(p^{k+1}-1)^m)+p^{k+1}(p^{(k+1)m}-(p^{k+1}-1)^m)\\ =p^mf(p^k)+p^k(p-1)(p^{(k+1)m}-(p^{k+1}-1)^m) \]

　　當 \(n\leq {10}^7\) 的時候能夠線性篩算出 \(f(1)\sim f(n)\)。

　　時間複雜度：\(O(n(1+\frac{\log m}{\log n}))\)

　　當 \(n\) 大的時候：
\[ f(p)=g(0)+pg(1)-pg(0)\\ =(1-p)(p-1)^m+p^{m+1}\\ =p^{m+1}-(p-1)^{m+1}\\ =-\sum_{i=0}^m{(-1)}^{m+1-i}p^i\binom{m+1}{i}\\ =\sum_{i=0}^m{(-1)}^{m-i}p^i\binom{m+1}{i} \]
　　就能夠上 min_25 篩了。

　　時間複雜度：\(O(\frac{n^{\frac{3}{4}}}{\log n}m+\sqrt nm^2)\)

代碼

　　這個代碼有些東西沒有預處理，時間複雜度可能比上面寫的高一些。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<functional>
#include<cmath>
#include<vector>
#include<assert.h>
//using namespace std;
using std::min;
using std::max;
using std::swap;
using std::sort;
using std::reverse;
using std::random_shuffle;
using std::lower_bound;
using std::upper_bound;
using std::unique;
using std::vector;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef std::pair<int,int> pii;
typedef std::pair<ll,ll> pll;
void open(const char *s){
#ifndef ONLINE_JUDGE
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
void open2(const char *s){
#ifdef DEBUG
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
int rd(){int s=0,c,b=0;while(((c=getchar())<'0'||c>'9')&&c!='-');if(c=='-'){c=getchar();b=1;}do{s=s*10+c-'0';}while((c=getchar())>='0'&&c<='9');return b?-s:s;}
void put(int x){if(!x){putchar('0');return;}static int c[20];int t=0;while(x){c[++t]=x%10;x/=10;}while(t)putchar(c[t--]+'0');}
int upmin(int &a,int b){if(b<a){a=b;return 1;}return 0;}
int upmax(int &a,int b){if(b>a){a=b;return 1;}return 0;}
const ll p=1000000007;
ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}
int n,m;
namespace gao1
{
    const int N=10000010;
    int f[N];
    int pw[N];
    bool b[N];
    int pri[N];
    int c[N];
    int cnt;
    void gao()
    {
        f[1]=1;
        pw[1]=1;
        for(int i=2;i<=n;i++)
        {
            if(!b[i])
            {
                pri[++cnt]=i;
                pw[i]=fp(i,m);
            }
            for(int j=1;j<=cnt&&i*pri[j]<=n;j++)
            {
                b[i*pri[j]]=1;
                pw[i*pri[j]]=(ll)pw[i]*pw[pri[j]]%p;
                if(i%pri[j]==0)
                    break;
            }
        }
        for(int i=2;i<=n;i++)
        {
            if(!b[i])
            {
                c[i]=i;
                f[i]=fp(i,m+1)-fp(i-1,m+1);
            }
            for(int j=1;j<=cnt&&pri[j]<=i&&i*pri[j]<=n;j++)
            {
                int v=i*pri[j];
                if(i%pri[j]==0)
                {
                    c[v]=c[i]*pri[j];
                    if(c[v]==v)
                        f[v]=((ll)pw[pri[j]]*f[i]+(ll)i*(pri[j]-1)*(pw[v]-pw[v-1]))%p;
                    else
                        f[v]=(ll)f[c[v]]*f[v/c[v]]%p;
                    break;
                }
                c[v]=pri[j];
                f[v]=(ll)f[i]*f[pri[j]]%p;
            }
        }
        ll ans=0;
        for(int i=1;i<=n;i++)
            ans+=f[i];
        ans=(ans%p+p)%p;
        printf("%lld\n",ans);
    }
}
namespace gao2
{
    const int N=100010;
    ll inv[N],fac[N],ifac[N];
    ll binom(int x,int y)
    {
        return x>=y&&y>=0?fac[x]*ifac[y]%p*ifac[x-y]%p:0;
    }
    int _;
    int b[N],pri[N],cnt;
    void init()
    {
        inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
        for(int i=2;i<=100000;i++)
        {
            inv[i]=-p/i*inv[p%i]%p;
            fac[i]=fac[i-1]*i%p;
            ifac[i]=ifac[i-1]*inv[i]%p;
        }
        for(int i=2;i<=_;i++)
        {
            if(!b[i])
                pri[++cnt]=i;
            for(int j=1;j<=cnt&&i*pri[j]<=_;j++)
            {
                b[i*pri[j]]=1;
                if(i%pri[j]==0)
                    break;
            }
        }
    }
    ll pw[N],s[N];
    void init(int x)
    {
        pw[1]=1;
        for(int i=2;i<=x+2;i++)
        {
            if(!b[i])
                pw[i]=fp(i,x);
            for(int j=1;j<=cnt&&i*pri[j]<=n&&pri[j]<=i;j++)
            {
                pw[i*pri[j]]=pw[i]*pw[pri[j]]%p;
                if(i%pri[j]==0)
                    break;
            }
        }
        for(int i=1;i<=x+2;i++)
            s[i]=(s[i-1]+pw[i])%p;
    }
    ll pre[N],suf[N];
    ll 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]=pre[i-1]*(x-i)%p;
        suf[y+3]=1;
        for(int i=y+2;i>=1;i--)
            suf[i]=suf[i+1]*(x-i)%p;
        ll res=0;
        for(int i=1;i<=y+2;i++)
            res=(res+s[i]*pre[i-1]%p*suf[i+1]%p*ifac[i-1]%p*ifac[y+2-i]%p*((y+2-i)&1?-1:1))%p;
        return res;
    }
    ll s1[N],s2[N],ans1[N],ans2[N];
    void gao(int x)
    {
        init(x);
        calc(3,x);
        for(int i=1;i<=_;i++)
            s1[i]=(fp(i,x)+s1[i-1])%p;
        for(int i=1;n/i>_;i++)
            s2[i]=calc(n/i,x);
        for(int i=1;i<=cnt;i++)
        {
            ll v=fp(pri[i],x);
            int j;
            for(j=1;n/j/pri[i]>_&&n/j>=(ll)pri[i]*pri[i];j++)
                s2[j]=(s2[j]-v*(s2[j*pri[i]]-s1[pri[i]-1]))%p;
            for(;n/j>_&&n/j>=(ll)pri[i]*pri[i];j++)
                s2[j]=(s2[j]-v*(s1[n/j/pri[i]]-s1[pri[i]-1]))%p;
            for(j=_;j>=(ll)pri[i]*pri[i];j--)
                s1[j]=(s1[j]-v*(s1[j/pri[i]]-s1[pri[i]-1]))%p;
        }
    }
    ll g(int z,int x,int y)
    {
        if(z<0)
            return 0;
        return z==y?fp(x,(ll)y*m):fp(fp(x,y)-fp(x,y-z-1),m);
    }
    ll get(int x,int y)
    {
        ll res=0;
        for(int i=0;i<=y;i++)
            res=(res+fp(x,i)*(g(i,x,y)-g(i-1,x,y)))%p;
        return res;
    }
    ll solve(ll x,int y)
    {
        if(x<=1||x<pri[y])
            return 0;
        if(y>cnt)
            return (x<=_?ans1[x]:ans2[n/x])-ans1[_];
        ll res=(x<=_?ans1[x]:ans2[n/x])-ans1[pri[y]-1];
        for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
        {
            ll x1=x/pri[i];
            for(int j=1;x1>=pri[i];j++,x1/=pri[i])
                res=(res+solve(x1,i+1)*get(pri[i],j)+get(pri[i],j+1))%p;
        }
        return res;
    }
    void gao()
    {
        _=sqrt(n)+0.5;
        init();
        for(int i=0;i<=m;i++)
        {
            gao(i);
            ll v=binom(m+1,i)*((m-i)&1?-1:1);
            for(int i=1;i<=_;i++)
                ans1[i]=(ans1[i]+v*s1[i])%p;
            for(int i=1;n/i>_;i++)
                ans2[i]=(ans2[i]+v*s2[i])%p;
        }
        ll ans=solve(n,1);
        ans++;
        ans=(ans%p+p)%p;
        printf("%lld\n",ans);
    }
}
int main()
{
    open("c");
    scanf("%d%d",&n,&m);
    if(n<=1e7)
        gao1::gao();
    else
        gao2::gao();
    return 0;
}