【BZOJ-4407】于神之怒加强版 莫比乌斯反演 + 线性筛

4407: 于神之怒加强版

Time Limit: 80 Sec  Memory Limit: 512 MB
Submit: 241  Solved: 119
[Submit][Status][Discuss]

Description

给下N,M,K.求

Input

输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示，下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。

Output

如题

Sample Input

1 2
3 3

Sample Output

20

HINT

1<=N,M,K<=5000000,1<=T<=2000

题解：JudgeOnline/upload/201603/4407.rar

Source

命题人：成都七中张耀楠，鸣谢excited上传。

 

Solution

首先变换一下式子：

$$sum_{d=1}^{n}d^{k}sum_{i=1}^{n}sum_{j=1}^{m}left lfloor gcdleft ( i,j ight )= d ight floor$$

那么我们设$fleft ( d ight )$表示$gcdleft ( i,j ight )= d$的点对的数目，那么可以莫比乌斯反演得到：

$$fleft ( d ight )= sum_{x=1}^{left lfloor frac{n}{d} ight floor}mu left ( x ight )left lfloor frac{n}{dx} ight floorleft lfloor frac{m}{dx} ight floor$$

那么就有：

 $$Ans= sum_{d=1}^{n}d^{k} imes f(d)$$

但这还不够求解，那么令$y= dx$代换一下可以得到：

$$Ans= sum_{y}^{n}left lfloor frac{n}{y} ight floorleft lfloor frac{m}{y} ight floorsum_{d|y}d^{k}mu left ( frac{y}{d} ight )$$

到这一步就已经可以求解了：

令$gleft ( y ight )= sum_{d|y}d^{k}mu left ( frac{y}{d} ight )$，发现是积性函数，那么线性筛处理出来即可

然后分块求解即可。

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
int read()
{
    int x=0,f=1; char ch=getchar();
    while (ch<'0' || ch>'9') {if (ch=='-') f=-1; ch=getchar();}
    while (ch>='0' && ch<='9') {x=x*10+ch-'0'; ch=getchar();}
    return x*f;
}
#define maxn 5000010
#define p 1000000007
int T,K,N,M;
long long quick_pow(long long x,int y)
{
    long long re=1; x=x%p; y%=p;
    for (int i=y; i; i>>=1,x=x*x%p)
        if (i&1) re=re*x%p;
    return re;
}
bool flag[maxn];long long F[maxn],prime[maxn],cnt,sum[maxn];
void prework()
{
    flag[1]=1; F[1]=1; sum[1]=1;
    for (int i=2; i<maxn; i++)
        {
            if (!flag[i]) prime[++cnt]=i,F[i]=quick_pow(i,K)-1;
            for (int j=1; j<=cnt && i*prime[j]<maxn; j++)
                {
                    flag[i*prime[j]]=1;
                    if (!(i%prime[j]))
                        {F[i*prime[j]]=F[i]*quick_pow(prime[j],K)%p;break;}
                    else F[i*prime[j]]=F[i]*F[prime[j]]%p;
                }
            sum[i]=sum[i-1]+F[i]%p;
        }
}
void work(int n,int m)
{
    if (n>m) swap(n,m); 
    long long ans=0;
    for (int j,i=1; i<=n; i=j+1)
        j=min(m/(m/i),n/(n/i)),
        ans+=(sum[j]-sum[i-1]+p)%p*(n/i)%p*(m/i)%p,ans%=p;
    printf("%lld
",ans);
}
int main()
{
    T=read(),K=read();
    prework();
    while (T--)
        {
            N=read(),M=read();
            work(N,M);
        }
    return 0;
}

数论题做的巨心累，推了半天，毫无头绪，最后默默看题解....zky学长说这是裸题...