【BZOJ】3309: DZY Loves Math 莫比乌斯反演优化

3309: DZY Loves Math

Description

对于正整数n，定义f(n)为n所含质因子的最大幂指数。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。
给定正整数a,b，求sigma(sigma(f(gcd(i,j)))) (i=1..a, j=1..b)。

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<=10^7

参考博文:PoPoQQQ

盗的

我主要讲讲为什么只有在当T素因数分解所有的幂指数都相等时才会有(-1)^k+1的贡献;

为了能将内层循环优化掉，就需要将a,b提出来，预处理；这时就变成第三步了，对T质因数分解为 T = p1^a1*p2^a2*...*pk^ak;

有莫比乌斯反演公式易知，d|T要使得mu[T/d] != 0只有每一个素因子的系数均 <= 1才行，这时总的d 的取法就为2^k(k为T质因子的个数)；

<1>当a1 = ... = ak = a时，对于f[d]易知只有当全部取的都是a-1时，才为a-1,这时只有一种情况；贡献为(a-1)*(-1)^k = a*(-1)^k + (-1)^k;

f[d] = a需要用组合数求解；(组合数为C(k,i)，其中i表示取a-1的个数，即i = 0,1,...,k-1),这里i != k;(i = k时就是f[d] = a-1了)，

这时的贡献为a*((-1)^0*C(k,0)+(-1)^1*C(k,1)+...+(-1)^(k-1)*C(k,k-1)) = ?

直接看(1-1)^K就会发现少了a*(-1^k*C(K,K)) = a*(-1)^k,所以加上第一种情况之后，需要-(-1)^k = (-1)^k+1

<2> 当存在ai != aj时，最大幂无论去a还是a-1都是f[d]的取值，所以不可能出现分情况讨论，这时其他的去aj | aj-1都会造成奇偶的变化而使得f[d]变为0；

所以结论就是只有a1 = a2 = ... = ak = a时，Σf[d]*μ(T/d)才不为0，且f[d]是不需要计算的，要得到的只是k的奇偶 = > (-1)^(k+1)

预处理时间复杂度与筛法相同O(n);之后每组数据处理时间为O(sqrt(n)),即分块加速;

编码的细节:里面使用了递推来得到g[]:得到当前值为0还是1；即(-1)^(k+1)，之后使用前缀和是为了分块加速；

同样递推还用在了a[i] : i的最小质因子的幂指数；val[i]：i最小质因子的幂指数乘积即p1^a1，这是为了在之后得到g[i]的时候，判断a[i] ?= a[j/val[i]]即最小的质因子的幂指数是否和倒数第二小的幂指数相等。注意这里是递推关系，即并不是只比较了最小的两个幂指数，就相当于马尔科夫链一样~~继承下来了~~(只是看别人的code的理解)巧妙 orz

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<stdlib.h>
using namespace std;
#define rep0(i,l,r) for(int i = (l);i < (r);i++)
#define rep1(i,l,r) for(int i = (l);i <= (r);i++)
#define rep_0(i,r,l) for(int i = (r);i > (l);i--)
#define rep_1(i,r,l) for(int i = (r);i >= (l);i--)
#define inf 0x7fffffff
typedef long long ll;
template<typename T>
void read1(T &m)
{
    T x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    m = x*f;
}
template<typename T>
void read2(T &a,T &b){read1(a);read1(b);}
template<typename T>
void read3(T &a,T &b,T &c){read1(a);read1(b);read1(c);}
const int N = 1e7+7;
int vis[N],g[N],p[N],a[N],val[N];
void init()
{
    rep0(i,2,N){
        if(vis[i] == 0){
            p[++p[0]] = i;
            g[i] = 1;
            a[i] = 1;
            val[i] = i;
        }
        for(int j = 1;p[j]*i < N;j++){
            vis[i*p[j]] = 1;
            if(i%p[j] == 0){
                a[i*p[j]] = a[i]+1;//最小质因数的幂指数;
                val[i*p[j]] = val[i]*p[j];
                int tmp = i/val[i];
                if(tmp == 1) g[i*p[j]] = 1;//只有一个质因数
                else
                    g[i*p[j]] = (a[tmp] == a[i*p[j]]?-g[tmp]:0);                
                break;
            }
            a[p[j]*i] = 1;//表示p[j]的次数为1;
            val[i*p[j]] = p[j];
            g[p[j]*i] = (a[i] == 1?-g[i]:0);
        }
    }
    rep0(i,1,N) g[i] += g[i-1];
}
int main()
{
    init();
    int n,m,T,i,j;
    read1(T);
    while(T--){
        read2(n,m);
        ll ans = 0;
        if(n > m) swap(n,m);
        for(i = 1;i <= n;i = j + 1){
            j = min(n/(n/i),m/(m/i));
            ans += 1LL*(g[j] - g[i-1])*(n/i)*(m/i);
        }
        printf("%lld\n",ans);// I64d是会WA的，运行系统不同的结果..
    }
    return 0;
}