Description

Input

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

Output

Sample Input

样例1

3 3 2

样例2

5000000 5000000 4000000

Sample Output

样例1

样例2

913111630

Data Constraint

Solution

直接上莫比乌斯反演：

Ans=∑i=1n∑j=1mgcd(i,j)kAns=∑i=1n∑j=1mgcd(i,j)k

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

Ans=∑ddk∑i=1n∑j=1m[gcd(i,j)=d]Ans=∑ddk∑i=1n∑j=1m[gcd(i,j)=d]

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

Ans=∑ddk∑i=1⌊nd⌋∑j=1⌊md⌋[gcd(i,j)=1]Ans=∑ddk∑i=1⌊nd⌋∑j=1⌊md⌋[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}[gcd(i,j)=1]

Ans=∑ddk∑i=1⌊nd⌋∑j=1⌊md⌋∑u|gcd(i,j)μ(u)Ans=∑ddk∑i=1⌊nd⌋∑j=1⌊md⌋∑u|gcd(i,j)μ(u)

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)

这一步是因为：

∑d|nμ(d)=[n=1]∑d|nμ(d)=[n=1]

\sum_{d|n}μ(d)=[n=1]
于是我们继续变换：

Ans=∑ddk∑uμ(u)⌊ndu⌋⌊mdu⌋Ans=∑ddk∑uμ(u)⌊ndu⌋⌊mdu⌋

Ans=\sum_{d}d^k\sum_{u}μ(u)\lfloor\frac{n}{du}\rfloor\lfloor\frac{m}{du}\rfloor
令 D=duD=duD=du ，再将 DDD 提到前面（这是为了凑形式）：

Ans=∑D⌊nD⌋⌊mD⌋∑d|Ddk·μ(Dd)" role="presentation">Ans=∑D⌊nD⌋⌊mD⌋∑d|Ddk⋅μ(Dd)Ans=∑D⌊nD⌋⌊mD⌋∑d|Ddk·μ(Dd)

Ans=\sum_{D}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{d|D}d^k·μ(\frac{D}{d})
令单位幂函数 idk(d)=dkidk(d)=dkid^k(d)=d^k ，那么有：

Ans=∑D⌊nD⌋⌊mD⌋f(D)Ans=∑D⌊nD⌋⌊mD⌋f(D)

Ans=\sum_{D}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor f(D)
其中：

f(D)=∑d|Didk(d)⋅μ(Dd)f(D)=∑d|Didk(d)·μ(Dd)

f(D)=\sum_{d|D}id^k(d)·μ(\frac{D}{d})
是狄利克雷卷积！即：f=idk∗μf=idk∗μf=id^k*μ
对于 fff 这个积性函数，我们可以用线筛求出（跟筛 φ" role="presentation" style="position: relative;">φφφ 差不多）。
之后做出前缀和就可以分块求答案了。
这种算法明显优于其它什么两次分块的算法，可以做多组询问。
就是提前给出 kkk ，每次询问给出一组 n,m" role="presentation" style="position: relative;">n,mn,mn,m ，范围还是 n,m,k≤5∗106n,m,k≤5∗106n,m,k\le5*10^6 ，询问数 T≤2000T≤2000T\le2000 。
那先预处理出 fff ，每次 O(n)" role="presentation" style="position: relative;">O(n−−√)O(n)O(\sqrt n) 回答询问，时间复杂度 O(n+Tn−−√)O(n+Tn)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;
}

JZOJ 4161. 于神之怒相关推荐

GMOJ 4161 / Luogu P4449 于神之怒 (加强版) 题解
于神之怒 (加强版) 题解 Description 求 F ( n , m ) = ∑ i = 1 n ∑ j = 1 m gcd ⁡ ( i , j ) k F(n,m)=\sum_{i=1}^n\ ...
JZOJ 5461 购物 —— 贪心
题目:https://jzoj.net/senior/#main/show/5461 贪心,原来想了个思路,优先选优惠价最小的 K 个,然后其他按原价排序遍历: 如果当前物品没选过,原价选上,如果选过 ...
BZOJ 4407: 于神之怒加强版
4407: 于神之怒加强版 Time Limit: 80 Sec Memory Limit: 512 MB Submit: 560 Solved: 271 [Submit][Status][Dis ...
JZOJ 5372. 【NOIP2017提高A组模拟9.17】猫
Description 信息组最近猫成灾了!隔壁物理组也拿猫没办法.信息组组长只好去请神刀手来帮他们消灭猫.信息组现在共有n 只猫(n 为正整数),编号为1 到n,站成了一个环,第i 只猫的左边是第i ...
uwsgi部署到nginx出现invalid request block size: 4161 (max 4096)...skip问题(亲测)
http://50vip.com/419.html 使用Flask制作一个网页平台之后,登陆使用openid登陆,然后使用uwsgi服务部署到nginx上,运行起来没有什么问题,但是偶尔在登陆的时候出 ...
[jzoj NOIP2018模拟 11.01]
很庆幸打了这场模拟赛,因为这一场爆零好像上次纪中的某场比赛我也出现了同样的问题,光是计算时间复杂度而忘记了空间的限制.想必是比上次惨的,考场上就写了两题而这两题都因为MLE爆零了.而且我T2还码了7 ...
BZOJ 4407 于神之怒加强版
题目链接:于神之怒加强版这个式子还是很妙的,只是我已经思维僵化了 \begin{aligned} &\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k \\ =& ...
P4449 于神之怒加强版
P4449 于神之怒加强版推式子 ∑i=1n∑j=1ngcd(i,j)h\sum_{i = 1} ^{n} \sum_{j = 1} ^{n} gcd(i, j) ^ h i=1∑nj=1∑ng ...
【DP】小学生语文题（jzoj 5102）
正题 jzoj 5102 题目大意给你两个串A,B,字母个数相等,可以把B的一个字符移到前面某个位置,问你最少移多少次可以使A,B相等解题思路设fi,jf_{i,j}fi,j为A匹配了i-n, ...

JZOJ 4161. 于神之怒