DZY Loves Math [Bzoj 3309]

题目地址请点击——

DZY Loves Math

Description

对于正整数 n，定义 f(n) 为 n 所含质因子的最大幂指数。

例如 f(1960)=f(23∗51∗72)=3 , f(10007)=1 , f(1)=0 。

给定正整数 a , b，求 ∑ai=1∑bj=1f(gcd(i,j))。

Input

第一行一个数 T，表示询问数。

接下来 T 行，每行两个数 a,b，表示一个询问。

Output

对于每一个询问，输出一行一个非负整数作为回答。

Sample Input

4

7558588 9653114

6514903 4451211

7425644 1189442

6335198 4957

Sample Output

35793453939901

14225956593420

4332838845846

15400094813

Hint

T<=10000

1<=a,b<=107

Solution

令 n≤m 。

设 g(d) 表示 ∑ni=1∑mj=1[gcd(i,j)=1]⋅f(d)，则

g(d)=∑d|p⌊np⌋⌊mp⌋μ(pd)f(d)=∑i=1⌊nd⌋⌊nid⌋⌊mid⌋μ(i)f(d)

ans=∑i=1dg(d)=∑i=1n⌊ni⌋⌊mi⌋∑t|iμ(t)f(it)

令 F(i)=∑t|iμ(t)f(it)，设 i 的唯一分解式为 pq11×pq22×…×pqss 。

要使 μ(t) 不为零，则

t=pm11×pm22×…×pmss , it=pq1−m11×pq2−m22×…×pqs−mss

0≤mj≤1(1≤j≤s)

设 Max(q1,q2,q3,…,qs)=a，number(qj=a)=b 。

所以 f(it)=a 或 a−1。

(1)f(it)=a，

∴sum=a⋅∑i=1b[Cib⋅∑j=0s−bCjs−b⋅(−1)b−i+j]=a⋅∑i=1b[Cib⋅(−1)b−i⋅∑j=0s−bCjs−b⋅(−1)j]

当且仅当 s=b 时，

sum=a⋅∑i=1s[Cis⋅(−1)s−i]=a⋅(∑i=0sCis⋅(−1)s−i−C0s⋅(−1)s)=−a(−1)s

否则，∑s−bj=0Cjs−b⋅(−1)j=0，则 sum=0。

(2)f(it)=a−1，

∴sum=(a−1)⋅∑j=0s−bCjs−b⋅(−1)b+j=(a−1)⋅(−1)b⋅∑j=0s−bCjs−b⋅(−1)j

当且仅当 s=b 时，

sum=(a−1)(−1)s

否则，∑s−bj=0Cjs−b⋅(−1)j=0，则 sum=0。

综上，当 s=b 时，F(i)=(a−1)(−1)s−a(−1)s=(−1)s+1，否则，F(i)=0。

求出 F 的前缀和，然后分块优化即可。

Code

#include <iostream>
#include <cstdio>

#define LL long long
#define MAXN 10000000
#define Min(x,y) ((x)<(y)?(x):(y))

using namespace std;

LL T,n,m;
bool no_prime[MAXN+10];
LL prime[MAXN+10];
LL sum[MAXN+10];
LL g[MAXN+10];
LL ci[MAXN+10];
LL rci[MAXN+10];

int main(){

sum[1]=0;

for(LL i=2;i<=MAXN;i++){
if(!no_prime[i]){prime[++prime[0]]=i;ci[i]=i;g[i]=1;rci[i]=1;}
for(LL j=1;prime[j]*i<=MAXN;j++){
no_prime[prime[j]*i]=true;
if(i%prime[j]==0){
rci[prime[j]*i]=rci[i]+1;
ci[prime[j]*i]=ci[i]*prime[j];
if(prime[j]*i/ci[prime[j]*i]==1)g[prime[j]*i]=1;
else if(rci[prime[j]*i]==rci[prime[j]*i/ci[prime[j]*i]])g[prime[j]*i]=g[prime[j]*i/ci[prime[j]*i]]*(-1);
break;
}
rci[prime[j]*i]=1;
ci[prime[j]*i]=prime[j];
if(prime[j]*i/ci[prime[j]*i]==1)g[prime[j]*i]=1;
else if(rci[prime[j]*i]==rci[prime[j]*i/ci[prime[j]*i]])g[prime[j]*i]=g[prime[j]*i/ci[prime[j]*i]]*(-1);
}
sum[i]=sum[i-1]+g[i];
}

scanf("%lld",&T);

while(T--){
LL ans=0;
scanf("%lld%lld",&n,&m);
if(n==0||m==0){
printf("0\n");
continue;
}
if(n>m)swap(n,m);
for(LL i=1,it;i<=n;i=it+1){
it=Min(n/(n/i),m/(m/i));
ans+=(sum[it]-sum[i-1])*(n/i)*(m/i);
}
printf("%lld\n",ans);
}

return 0;
}