JZOJ 4161. 于神之怒

Description

Input

一行三个数，m, n, k。

Output

Sample Input

样例1

3 3 2

样例2

5000000 5000000 4000000

Sample Output

样例1

20

样例2

​913111630

Data Constraint

Solution

• 直接上莫比乌斯反演：

$$Ans=\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)^k$$

$$Ans=\sum_{d}d^k\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]$$

$$Ans=\sum_{d}d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]$$

$$Ans=\sum_{d}d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{u|gcd(i,j)}μ(u)$$

• 这一步是因为：

  $$\sum_{d|n}μ(d)=[n=1]$$

• 于是我们继续变换：
  

  $$Ans=\sum_{d}d^k\sum_{u}μ(u)\lfloor\frac{n}{du}\rfloor\lfloor\frac{m}{du}\rfloor$$

• 令 $D=du$ ，再将 $D$ 提到前面（这是为了凑形式）：

  $$Ans=\sum_{D}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{d|D}d^k·μ(\frac{D}{d})$$

• 令单位幂函数 $id^k(d)=d^k$ ，那么有：

  $$Ans=\sum_{D}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor f(D)$$

• 其中：

  $$f(D)=\sum_{d|D}id^k(d)·μ(\frac{D}{d})$$

• 是狄利克雷卷积！即：$f=id^k*μ$

• 对于 $f$ 这个积性函数，我们可以用线筛求出（跟筛 $φ$ 差不多）。

• 之后做出前缀和就可以分块求答案了。

• 这种算法明显优于其它什么两次分块的算法，可以做多组询问。

• 就是提前给出 $k$ ，每次询问给出一组 $n,m$ ，范围还是 $n,m,k\le5*10^6$ ，询问数 $T\le2000$ 。

• 那先预处理出 $f$ ，每次 $O(\sqrt n)$ 回答询问，时间复杂度 $O(n+T\sqrt n)$ 。

Code

#include<cstdio>
#include<algorithm>
#include<cctype>
using namespace std;
typedef long long LL;
const int N=5e6+5,mo=1e9+7;
int f[N],g[N];
bool bz[N];
inline int read()
{
    int X=0,w=0; char ch=0;
    while(!isdigit(ch)) w|=ch=='-',ch=getchar();
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    return w?-X:X;
}
inline int ksm(int x,int y)
{
    int s=1;
    while(y)
    {
        if(y&1) s=(LL)s*x%mo;
        x=(LL)x*x%mo;
        y>>=1;
    }
    return s;
}
inline int min(int x,int y)
{
    return x<y?x:y;
}
int main()
{
    int n=read(),m=read(),k=read(),ans=0;
    g[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!bz[i])
        {
            f[++f[0]]=i;
            g[i]=ksm(i,k)-1;
        }
        for(int j=1;j<=f[0] && i*f[j]<N;j++)
        {
            bz[i*f[j]]=true;
            if(i%f[j]==0)
            {
                g[i*f[j]]=(LL)g[i]*(g[f[j]]+1)%mo;
                break;
            }
            g[i*f[j]]=(LL)g[i]*g[f[j]]%mo;
        }
    }
    for(int i=1;i<N;i++) g[i]=(g[i-1]+g[i])%mo;
    if(n>m) swap(n,m);
    for(int i=1,p;i<=n;i=p+1)
    {
        p=min(n/(n/i),m/(m/i));
        ans=(ans+(LL)(n/p)*(m/p)%mo*(g[p]-g[i-1]+mo)%mo)%mo;
    }
    printf("%d\n",ans);
    return 0;
}