HDU/HDOJ 3875 Euclidean Algorithm 多校联合第四场

题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=3875

Problem Description

wr recently learned the Euclidean algorithm to solve the greatest common divisor,then he realized that the algorithm can also be able to solve the least common multiple.In order to train the skill of code,he find a number p(1<=p<=2^25)and a number q(1<=q<=2^25) then he got n=p*q,he calculated a ₁=∑gcd(i, n) 1<=i <=n and a ₂= ∑lcm(i, n) 1<=i <=n,he wants to know if a2-a1 is the multiple of c(1<=c<=10^9)

Input

The first line contains the number of scenarios.
Every scenario consists of a single line containing two integers n and c separated by a space.

Output

for every scenario print a single line contain "yes" or "no"
if the difference is the multiple of c print "yes",otherwise print "no"

Sample Input

2 6 100 4 2

Sample Output

no yes

Source

2011 Multi-University Training Contest 4 - Host by SDU

Recommend

lcy

这个题应该是数论中比较综合的一道题目。

做这个题需要知道大数的质因数分解方法（米勒罗宾测试），欧拉函数的积性函数性质，除法取模的方法等等

这个题由于数据量很大，所以要分解质因数必须要用到米勒罗宾的随机算法，快速分解。

然后后面就是公式推导了。

首先看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; }