对于正整数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)。
标签：
对于正整数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)。
第一行一个数T，表示询问数。
接下来T行，每行两个数a,b，表示一个询问。
对于每一个询问，输出一行一个非负整数作为回答。
【数据规模】
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 的取法就为2k(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最小质因子的幂指数乘积即p[1]^a[i]，这是为了在之后得到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; }
【BZOJ】3309: DZY Loves Math 莫比乌斯反演优化
标签：
原文地址：http://www.cnblogs.com/hxer/p/5290059.html
