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;}