博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
【BZOJ】3309: DZY Loves Math 莫比乌斯反演优化
阅读量:4552 次
发布时间:2019-06-08

本文共 2913 字,大约阅读时间需要 9 分钟。

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

参考博文:

盗的

我主要讲讲为什么只有在当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最小质因子的幂指数乘积即p1^a1,这是为了在之后得到g[i]的时候,判断a[i] ?= a[j/val[i]]即最小的质因子的幂指数是否和倒数第二小的幂指数相等。注意这里是递推关系,即并不是只比较了最小的两个幂指数,就相当于马尔科夫链一样~~继承下来了~~(只是看别人的code的理解)巧妙 orz

#include
#include
#include
#include
#include
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 0x7ffffffftypedef long long ll;template
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
void read2(T &a,T &b){read1(a);read1(b);}template
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;}

 

转载于:https://www.cnblogs.com/hxer/p/5290059.html

你可能感兴趣的文章
proxy和proxy-no的策略取值区别
查看>>
Silverlight代码编写对控件的PlaneProjection.RotationY属性控制动画
查看>>
AFNetworking
查看>>
unity3d Start执行不同时问题
查看>>
session
查看>>
JS只能输入数字
查看>>
Laravel 数据库连接, 数据库名,配置文件修改
查看>>
屌丝接盘侠们,孩子可能不是你们亲生的!
查看>>
BZOJ 1854 【SCOI2010】 游戏
查看>>
JavaScript - 匿名函数和闭包
查看>>
负载均衡下的资源文件配置/多站点下的资源文件夹共享(Windows IIS)
查看>>
MySQL firstmatch strategy
查看>>
MS SQL server 2014 创建用户及权限
查看>>
office很抱歉遇到一些临时服务器问题
查看>>
禁止键盘上的刷新键F5等
查看>>
SAP中对于获取订单的状态
查看>>
oracle PL/SQL块
查看>>
sklearn.preprocessing.LabelBinarizer
查看>>
C teaching
查看>>
分隔指定内容,提取章节数
查看>>