BZOJ3561 DZY Loves Math VI 数论 快速幂 莫比乌斯反演
原文链接http://www.cnblogs.com/zhouzhendong/p/8116330.html
UPD(2018-03-26):回来重新学数论啦。之前的博客版面放在更新之后的后面。
题目传送门 - BZOJ3561
题意概括
给出$n,m$,求$\Large\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)^{\gcd(i, j)}$。
$1\leq n,m\leq 500000$
题解
先推式子：（假设$n\leq m$）
$$\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)^{\gcd(i, j)}\\=\sum_{d=1}^{n}\sum_{i=1}^{\left\lfloor\frac nd\right\rfloor}\sum_{j=1}^{\left\lfloor\frac md\right\rfloor}(ijd)^d\cdot[\gcd(i,j)=1]\\=\sum_{d=1}^{n}\sum_{i=1}^{\left\lfloor\frac nd\right\rfloor}\sum_{j=1}^{\left\lfloor\frac md\right\rfloor}(ijd)^d\cdot\sum_{p|i,p|j}\mu(p)\\=\sum_{d=1}^{n}d^{d}\sum_{p=1}^{\left\lfloor\frac nd\right\rfloor}\mu(p)\sum_{i=1}^{\left\lfloor\frac{n}{pd}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{pd}\right\rfloor}(ijp^2)^d\\=\sum_{d=1}^{n}d^{d}\sum_{p=1}^{\left\lfloor\frac nd\right\rfloor}\mu(p)p^{2d}\sum_{i=1}^{\left\lfloor\frac {n}{pd}\right\rfloor}i^d\sum_{j=1}^{\left\lfloor\frac{m}{pd}\right\rfloor}j^d$$
然后发现对于$d^d$可以直接快速幂。对于某一个$d$，要枚举的$p$有$O(\frac nd)$个，对于后面的一堆数的幂和，只要前缀和预处理，要处理的个数也是$O(\frac md)$的。所以总复杂度为$O(n \log n)$。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=500005;
const LL mod=1e9+7;
int n,m,prime[N],u[N],pcnt=0;
bool f[N];
void init(int n){
memset(f,true,sizeof f);
f[0]=f[1]=0,u[1]=1;
for (int i=2;i<=n;i++){
if (f[i])
prime[++pcnt]=i,u[i]=-1;
for (int j=1;j<=pcnt&&i*prime[j]<=n;j++){
f[i*prime[j]]=0;
if (i%prime[j])
u[i*prime[j]]=-u[i];
else {
u[i*prime[j]]=0;
break;
}
}
}
}
LL Pow(LL x,LL y){
if (!y)
return 1LL;
LL xx=Pow(x,y/2);
xx=xx*xx%mod;
if (y&1LL)
xx=xx*x%mod;
return xx;
}
LL pows[N],sum[N];
int main(){
scanf("%d%d",&n,&m);
if (n>m)
swap(n,m);
init(n);
for (int i=1;i<=m;i++)
pows[i]=1;
LL ans=0;
for (int d=1;d<=n;d++){
LL now=0;
sum[0]=0;
for (int i=1;i<=m/d;i++)
pows[i]=pows[i]*i%mod,sum[i]=(sum[i-1]+pows[i])%mod;
for (int p=1;p<=n/d;p++)
now=(now+u[p]*pows[p]*pows[p]%mod*sum[n/p/d]%mod*sum[m/p/d])%mod;
now=(now%mod+mod)%mod;
ans=(ans+Pow(d,d)*now)%mod;
}
printf("%lld",ans);
return 0;
}
——————old——————
题意概括
题解
博主越来越懒了。
http://blog.csdn.net/lych_cys/article/details/50721642?locationNum=1&fps=1
代码
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long LL;
const int N=500005;
const LL mod=1e9+7;
LL n,m,u[N],prime[N],pcnt,v[N],sum[N];
bool isprime[N];
LL Pow(LL x,LL y){
if (y==0)
return 1LL;
LL xx=Pow(x,y/2);
xx=xx*xx%mod;
if (y&1LL)
xx=xx*x%mod;
return xx;
}
void Get_Mobius(){
memset(isprime,true,sizeof isprime);
isprime[0]=isprime[1]=pcnt=0;
u[1]=1;
for (LL i=2;i<=n;i++){
if (isprime[i])
u[i]=-1,prime[++pcnt]=i;
for (LL j=1;j<=pcnt&&i*prime[j]<=n;j++){
isprime[i*prime[j]]=0;
if (i%prime[j])
u[i*prime[j]]=-u[i];
else {
u[i*prime[j]]=0;
break;
}
}
}
}
int main(){
scanf("%lld%lld",&n,&m);
if (n<m)
swap(n,m);
Get_Mobius();
for (int i=1;i<=n;i++)
v[i]=1;
LL ans=0;
for (LL d=1;d<=m;d++){
sum[0]=0;
for (LL i=1;i<=(LL)(n/d);i++)
v[i]=v[i]*i%mod,sum[i]=(v[i]+sum[i-1])%mod;
LL res=0;
for (LL p=1;p<=(LL)(m/d);p++)
res=(res+v[p]*v[p]%mod*u[p]*sum[n/d/p]%mod*sum[m/d/p]%mod+mod)%mod;
ans=(ans+res*Pow(d,d))%mod;
}
printf("%lld",ans);
return 0;
}
