1. 程式人生 > >BZOJ3561 - DZY Loves Math VI

BZOJ3561 - DZY Loves Math VI

ble php .com 根據 1=1 script source else http

Portal

Description

\(T(T\leq3)\)組測試數據。給出\(n,m(n\leq5\times10^5)\),求
\[ \sum_{i=1}^n \sum_{j=1}^m lcm(i,j)^{gcd(i,j)}\]

Solution

喜聞樂見推推推。
\[\begin{align*} ans &= \sum_{i=1}^n \sum_{j=1}^m lcm(i,j)^{gcd(i,j)} \&= \sum_{d=1}^{+∞} \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=d](\frac{ij}{d})^d \&= \sum_{d=1}^{+∞} \sum_{i=1}^{?\frac{n}{d}?} \sum_{j=1}^{?\frac{m}{d}?} [gcd(i,j)=1](ijd)^d \&= \sum_{d=1}^{+∞} d^d \sum_{d_1=1}^{+∞}\mu(d_1) \sum_{d_1|i}^{?\frac{n}{d}?} \sum_{d_1|j}^{?\frac{m}{d}?} (ij)^d \&= \sum_{d=1}^{+∞} d^d \sum_{d_1=1}^{+∞}\mu(d_1) \sum_{i=1}^{?\frac{n}{dd_1}?} \sum_{j=1}^{?\frac{m}{dd_1}?} (id_1\cdot jd_1)^d \&= \sum_{d=1}^{+∞} d^d \sum_{d_1=1}^{+∞}\mu(d_1)d_1^{2d} \sum_{i=1}^{?\frac{n}{dd_1}?} i^d\sum_{j=1}^{?\frac{m}{dd_1}?} j^d \\end{align*}\]

於是我們枚舉\(d\),預處理出\(f(x)=\sum_{i=1}^x i^d\),然後枚舉\(d_1\)

根據調和級數總復雜度為\(O(nlogn)\)

Code

//DZY Loves Math VI
#include <algorithm>
#include <cstdio>
using std::swap;
typedef long long lint;
const int N=5e5+10;
const int P=1e9+7;
int prCnt,pr[N]; bool prNot[N];
int mu[N];
void init(int n)
{
    mu[1
]=1; for(int i=2;i<=n;i++) { if(!prNot[i]) pr[++prCnt]=i,mu[i]=-1; for(int j=1;j<=prCnt;j++) { int x=i*pr[j]; if(x>n) break; prNot[x]=true; if(i%pr[j]) mu[x]=-mu[i]; else break; } } } int powD[N],S[N]; int pow(int
x,int y) { lint r=1,t=x; for(int i=y;i;i>>=1,t=(t*t)%P) if(i&1) r=(r*t)%P; return r; } int main() { int task=1; init(5e5); while(task--) { int n,m; scanf("%d%d",&n,&m); if(n>m) swap(n,m); int ans=0; for(int i=1;i<=m;i++) powD[i]=1; for(int d=1;d<=n;d++) { int k1=n/d,k2=m/d,ansD=0; for(int i=1;i<=k2;i++) powD[i]=(1LL*powD[i]*i)%P; for(int i=1;i<=k2;i++) S[i]=(S[i-1]+powD[i])%P; for(int d1=1;d1<=k1;d1++) { lint t=1LL*(mu[d1]+P)%P*powD[d1]%P*powD[d1]%P; t=t*S[k1/d1]%P*S[k2/d1]%P; ansD=(ansD+t)%P; } ans=(ans+1LL*ansD*pow(d,d))%P; } printf("%d\n",ans); } return 0; }

BZOJ3561 - DZY Loves Math VI