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

时间：2024-01-02 13:57:44

4407: 于神之怒加强版

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

Description

给下N,M,K.求
【BZOJ-4407】于神之怒加强版 莫比乌斯反演 + 线性筛

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 gcd\left ( i,j \right )= d \right \rfloor$$

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

$$f\left ( d \right )= \sum_{x=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu \left ( x \right )\left \lfloor \frac{n}{dx} \right \rfloor\left \lfloor \frac{m}{dx} \right \rfloor$$

那么就有：

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

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

$$Ans= \sum_{y}^{n}\left \lfloor \frac{n}{y} \right \rfloor\left \lfloor \frac{m}{y} \right \rfloor\sum_{d|y}d^{k}\mu \left ( \frac{y}{d} \right )$$

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

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

然后分块求解即可。

Code

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cmath>

using namespace std;

int read()

{

    int x=,f=; char ch=getchar();

    while (ch<'' || ch>'') {if (ch=='-') f=-; ch=getchar();}

    while (ch>='' && ch<='') {x=x*+ch-''; 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=; x=x%p; y%=p;

    for (int i=y; i; i>>=,x=x*x%p)

        if (i&) re=re*x%p;

    return re;

}

bool flag[maxn];long long F[maxn],prime[maxn],cnt,sum[maxn];

void prework()

{

    flag[]=; F[]=; sum[]=;

    for (int i=; i<maxn; i++)

        {

            if (!flag[i]) prime[++cnt]=i,F[i]=quick_pow(i,K)-;

            for (int j=; j<=cnt && i*prime[j]<maxn; j++)

                {

                    flag[i*prime[j]]=;

                    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-]+F[i]%p;

        }

}

void work(int n,int m)

{

    if (n>m) swap(n,m);

    long long ans=;

    for (int j,i=; i<=n; i=j+)

        j=min(m/(m/i),n/(n/i)),

        ans+=(sum[j]-sum[i-]+p)%p*(n/i)%p*(m/i)%p,ans%=p;

    printf("%lld\n",ans);

}

int main()

{

    T=read(),K=read();

    prework();

    while (T--)

        {

            N=read(),M=read();

            work(N,M);

        }

    return ;

}

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

相关文章