BZOJ 4407: 于神之怒加强版

简单反演题，都是套路，直接写式子了……

\[\sum_{i=1}^n\sum_{j=1}^m \gcd(i,j)^k\\=\sum_{g=1}^{\min(n,m)} g^k\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor} [\gcd(i,j)=1]\\=\sum_{g=1}^{\min(n,m)} g^k\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor} \sum_{d|\gcd(i,j)} \mu(d)\\=\sum_{g=1}^{\min(n,m)} g^k\sum_{d=1}^{\min(\lfloor\frac{n}{g}\rfloor,\lfloor\frac{m}{g}\rfloor)} \mu(d)\times\lfloor\frac{n}{dg}\rfloor\lfloor\frac{m}{dg}\rfloor)\\=\sum_{D=1}^{\min(n,m)} \lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor \sum_{g|D} \mu(\frac{D}{g})\times g^k\\ \]

推导这里我们发现只要设\(f(D)=\sum_{g|D} \mu(\frac{D}{g})\times g^k\)，预处理出\(f\)就可以\(O(\sqrt n+\sqrt m)\)除法分块得到答案了

那么\(f\)怎么搞呢，首先这是个狄利克雷卷积的形式，同时里面的\(\mu,g^k\)都是积性函数，因此\(f\)也是积性函数

既然是积性函数就考虑当\(D=p_i^{a_i}\)时的情况，显然只有两种取值方法（注意\(\mu\)有平方质因子时值为\(0\)），我们容易算出：

\[f(n)=\begin{cases}1&n=1\\n^k-1&n\in Prime\end{cases} \]

然后容易推出线性筛的时候：

\[f(i\times P_j)=\begin{cases}f(i)\times f(P_j)&P_j\not| \ i\\f(i)\times P_j^k&P_j|i\end{cases} \]

然后就做完了，预处理是\(O(n)\)的

#include<cstdio>
#include<iostream>
#define RI register int
#define CI const int&
using namespace std;
const int N=5000000,mod=1e9+7;
int n,m,k,t,prime[N+5],cnt,f[N+5],ans; bool vis[N+5];
inline int quick_pow(int x,int p,int mul=1)
{
    for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline int sum(CI x,CI y)
{
    int t=x+y; return t>=mod?t-mod:t;
}
inline int sub(CI x,CI y)
{
    int t=x-y; return t<0?t+mod:t;
}
#define P(x) prime[x]
inline void init(CI n)
{
    RI i,j; for (vis[1]=f[1]=1,i=2;i<=n;++i)
    {
        if (!vis[i]) prime[++cnt]=i,f[i]=sub(quick_pow(i,k),1);
        int tp; for (j=1;j<=cnt&&(tp=i*P(j))<=n;++j)
        {
            vis[tp]=1; if (i%P(j)) f[tp]=1LL*f[i]*f[P(j)]%mod;
            else f[tp]=1LL*f[i]*quick_pow(P(j),k)%mod;
        }
    }
    for (i=2;i<=n;++i) f[i]=sum(f[i],f[i-1]);
}
#undef P
int main()
{
    for (scanf("%d%d",&t,&k),init(N);t;--t)
    {
        scanf("%d%d",&n,&m); ans=0; if (n>m) swap(n,m);
        for (RI l=1,r;l<=n;l=r+1) r=min(n/(n/l),m/(m/l)),
        ans=sum(ans,1LL*(n/l)*(m/l)%mod*sub(f[r],f[l-1])%mod);
        printf("%d\n",ans);
    }
    return 0;
}