HDU/HDOJ 3875 Euclidean Algorithm 多校联合第四场
来源:互联网 发布:python 最大回撤 编辑:程序博客网 时间:2024/05/17 05:08
题目连接:http://acm.hdu.edu.cn/showproblem.php?pid=3875
Every scenario consists of a single line containing two integers n and c separated by a space.
if the difference is the multiple of c print "yes",otherwise print "no"
26 1004 2
noyes
这个题应该是数论中比较综合的一道题目。
做这个题需要知道大数的质因数分解方法(米勒罗宾测试),欧拉函数的积性函数性质,除法取模的方法等等
这个题由于数据量很大,所以要分解质因数必须要用到米勒罗宾的随机算法,快速分解。
然后后面就是公式推导了。
首先看gcd sum的快速求解方法。
这里我们假设f(x)=sigma(gcd(i,n))
由于这是一个积性函数那么必然有:
f(x)=f(p1^k1)*f(p2^k2)*f(p3^k3)*....*f(pn^kn)
单单这样,似乎规律并不明显,不过第一步我们已经明确。就是将n进行质因数分解
然后接下来看形如f(p^k)有没有什么规律
注意到p^k的约数只有1,p,p^2,p^3.....p^k
那么必然的,f(p^k)=1*eular(p^k)+p*eular(p^(k-1))+....+p^k*eular(1)//这个题的基础是那个电子科大的数论中的一道题目演变而来的。这个公式不在附加说明。
注意到:eular(p^k)=p^k-p^(k-1)这个公式来自于数论概论,不过利用欧拉函数性质也很好推导。
把这个式子代入f(p^k)得到
f(p^k)=k*(p^k-p^(k-1))+p^k
于是对于gcd的和我们可以再分解完因式之后快速计算出来。
接下来是比较困难的lcm求和
有一个较为朴素的公式,相信大家都能推导出来:
ans=sigma(n/d*eular(n/d)*n/2)
n/2是每一个与n的最大公约数为d的数的平均数。(理论依据是所有小于n的且与n互质的数的和为(1+2+3+..+n)/2,详细内容可以查阅相关资料)
然后说说这个式子怎么转化吧。
其实,n/d*eular(n/d)=d*eular(d)。原因就是d会出现的数字,n/d也一定会出现。
举个例子:n=6,那么d=1,2,3,6 n/d=6,3,2,1
完全一样的说~
然后于是ans=sigma(d*eular(d)*n/2)
这里说说当d=1的时候,这个时候是一个特殊情况。
因为我们这个时候平均数事实上是n,而非n/2
所以式子需要修正:ans=sigma(d*eular(d)*n/2+n/2)
ans=n*sigma(d*eular(d)+1)
注意到上面两个式子在数学上面并不相等,但是由于这里,除号是一个整除法,所以这里,其实两个式子是相等的(至少计算机算出来是一样的)
然后这里之后,我们还可以继续化简
假设f(x)=x*eular(x)
这个也是一个积性函数
那么f(n)=f(p1^k1)*f(p2*k2)*f(p3^k3)....*f(pn*kn)
然后同样的,我们单独考虑
f(p^k)=1*1+p*eular(p)+p^2*eular(p^2)+....+p^k*eular(p^k)
继续化简:(等比数列求和)
得到:f(p^k)=(p^2k+1)/(p+1)
这样我们就也可以在分解完素因数之后快速求的lcm sum
最后注意一点由于lcm sum的公式里面有一个1/2
所以我们必须把原来的数扩大两倍,然后再把c扩大两倍,取余的结果最后除以2就是原来的答案
详细见代码:
#include<iostream>#include<cstring>#include<cstdlib>#include<cstdio>using namespace std;typedef unsigned __int64 u64;const int MAX=100;u64 f0[100],f1[100],ff,n,tmp,ret,ret1,p,pp;u64 myrandom(){ u64 a; a=rand(); a*=rand(); a*=rand(); a*=rand(); return a;}u64 mulmod(u64 a,u64 b,u64 c){ u64 ret=0; while(b){ if(b&1){ ret+=a; if(ret>c) ret-=c; } a<<=1; if(a>c) a-=c; b>>=1; } return ret;}u64 powmod(u64 a,u64 b,u64 c){ u64 ret=1; while(b){ if(b&1) ret=mulmod(ret,a,c); a=mulmod(a,a,c); b>>=1; } return ret;}int miller(u64 base,u64 n){ u64 m=n-1,k=0; while(m%2==0){ m>>=1; k++; } if(powmod(base,m,n)==1) return 1; for(int i=0;i<k;i++){ if(powmod(base,m<<i,n)==n-1) return 1; } return 0;}int Miller_Rabin(u64 n){ for(int i=2;i<100;++i) if(n%i==0) if(n==i) return 1; else return 0; for( i=0;i<MAX;++i){ u64 base=myrandom()%(n-1)+1; if(!miller(base,n)) return 0; } return 1;}u64 gcd(u64 a, u64 b){ if(b==0) return a; return gcd(b,a%b);}u64 f(u64 a,u64 b){ return (mulmod(a,a,b)+1)%b;}u64 Pollard_Rho(u64 n){ if(n<=2) return 0; for(u64 i=2;i<100;++i) if(n%i==0) return i; u64 p,x,xx; for( i=1;i<MAX;i ++){ x=myrandom()%n; xx=f(x,n); p=gcd((xx+n-x)%n,n); while(p==1){ x=f(x,n); xx=f(f(xx,n),n); p=gcd((xx+n-x)%n,n)%n; } if(p) return p; } return 0;}u64 Prime(u64 a){ if(Miller_Rabin(a)) return 0; u64 t=Pollard_Rho(a); u64 p=Prime(t); if(p) return p; return t;}int main(){ int tt; scanf("%d",&tt); while(tt--){ scanf("%I64d %I64d",&n,&p);pp=p*2; u64 old=n%pp; ret=1,ret1=1,ff=0; while(n>1){ if(Miller_Rabin(n)) break; tmp=Prime(n); f0[ff]=tmp; f1[ff]=1; n/=tmp; while(n%tmp==0){ n/=tmp; f1[ff]++; } ff++; } if(n>0){ f0[ff]=n; f1[ff++]=1; } for(int i=0;i<ff;++i){ u64 tmp=1; for(int j=0;j<f1[i];++j) tmp=tmp*f0[i]; ret1=ret1*(f1[i]*(tmp-tmp/f0[i])+tmp); ret1=ret1%p; tmp=1; for(j=0;j<2*f1[i]+1;++j) tmp=tmp*f0[i];ret=ret*((tmp+1)/(1+f0[i]));ret=ret%pp; } ret=((old*(ret+1))%pp)/2; if(ret==ret1) printf("yes\n"); else printf("no\n"); } return 0;}
- HDU/HDOJ 3875 Euclidean Algorithm 多校联合第四场
- HDOJ多校联合第四场
- hdu/hoj 3875 Euclidean Algorithm
- hdu 4642 fliping game 多校联合训练第四场
- hdu 4639 hehe 多校联合训练第四场
- HDU 6069 Counting Divisors 多校联合第四场
- HDOJ 4632 - Palindrome subsequence/2013多校联合第四场A 区间DP
- 2013多校联合训练第四场
- 2016多校联合第四场 HDU5768
- HDOJ多校联合第六场
- HDOJ多校联合第五场
- hdu 4632 Palindrome subsequence 区间dp 多校联合训练第四场
- 2016多校联合训练赛 第四场1001 Another Meaning hdu 5763
- 2016多校联合训练赛 第四场1012 Bubble Sort hdu 5775
- 2017多校联合第四场/HDU 6068 Classic Quotation(kmp+dp)
- HDU/HDOJ 3519 矩阵二分幂 2010多校联合第九场
- HDU/HDOJ 3609 Up-up 2010多校联合17场ZSTU
- HDOJ 4635 - Strongly connected/2013多校联合第四场D Tarjan求强联通图.缩点.
- 生产者-消费者问题实现 (linux下C语言)
- c语言的一些有用的库,记录一下。
- 二、转换成3NF的保持函数依赖的分解
- 生产者-消费者问题实现 (windows)
- linux c常用字符串处理函数
- HDU/HDOJ 3875 Euclidean Algorithm 多校联合第四场
- C# foreach或许你所不知道的
- sizeof的疑问
- Windows 2000/XP 注册表终极修改(转载) .
- 三、转换成3NF的保持无损连接和函数依赖的分解
- highcharts document
- 【VC++】如何解决unexpected end of file while looking for precompiled header directive
- HashMap遍历的两种方式
- 四、转换成BCNF的保持无损连接的分解