Jackiesteed

www.github.com/jackiesteed

  博客园  :: 首页  :: 新随笔  :: 联系 :: 订阅 订阅  :: 管理

PollardRho算法的应用.

输入a,b,处理为b/=a.

先使用Pollard算法求出b的所有的质因子,求出的信息中包含每个质因子在b中的个数.

设结果为x,y,

那么b=(x1^p1)*(x2^p2)*(x3^p3)*...,化简为b=X1*X2*X3*...*Xn,其中Xi=xi^pi.(因为x和y是互质的,所以x和y会独占Xi,这是化简的理由.)

通过打表可以推出n最大为15.

那么,那么位运算一下即可,(O(1<<n))

for(i=0;i<=(1<<n);i++)

{

  LL tmp=1;

  for(int j=0;j<n;j++)

  {if(i&(1<<j)) tmp*=Xj;}

  //下面是判断处理...求出最有结果即可...

}

最后x*=a,y*=a.

附上代码:

View Code
#include <map>
#include
<iostream>
#include
<fstream>
#include
<algorithm>
#include
<cmath>
#include
<cstring>
#include
<climits>
#include
<cstdio>
#include
<ctime>




using namespace std;

typedef
long long LL;

LL factor[
1000];
LL cnt;


inline LL gcd(LL a,LL b)
{
LL tmp;
while(b)
{
tmp
=a%b;
a
=b;
b
=tmp;
}
return a;
}

//求a*b%n,(a*b可能大于LLONG_MAX)
inline LL mulMod(LL a,LL b,LL n)
{
LL res
=0,e=a%n;
while(b)
{
if(b&1)
{
res
+=e;
if(res>n)
res
-=n;
}
b
>>=1;
e
<<=1;
if(e>n)
e
-=n;
}
return res;
}

//幂取模
LL expMod(LL a,LL b,LL n)
{
LL res
=1,e=a%n;
while(b)
{
if(b&1)
{
res
=mulMod(res,e,n);
}
b
>>=1;
e
=mulMod(e,e,n);
}

return res;
}


bool millerRabin(LL n,LL times)
{
if(n<2)
return false;
if(n==2)
return true;
if(!(n&1))
return false;

long long l=0,a,x,y;
long long m=n-1;
while(!(m&1))
l
++,m>>=1;
srand(time(NULL));
while(times--)
{
a
=rand()%(n-1)+1;
x
=expMod(a,m,n);
if(x==1)
continue;
for(int i=0;i<l;i++)
{
y
=mulMod(x,x,n);
if(y==1 && 1!=x && (n-1)!=x)
return false;
x
=y;
}
if(y!=1)
return false;
}
return true;
}

//质因子分解
LL pollardRho(long long n,long long c)
{
LL x,y,d,i
=1,k=2;
srand(time(NULL));

x
=rand()%(n-1)+1;
y
=x;
while(true)
{
i
++;
x
=(mulMod(x,x,n)+c)%n;
d
=gcd(y-x,n);
if(d>1 && d<n)
return d;
if(x==y)
return n;
if(i==k)
{
y
=x;
k
<<=1;
}
}
return n;
}

//寻找所有的质因子
void findFactor(LL n,LL k)
{
if(n==1)
return;

if(millerRabin(n,10))
{
factor[cnt
++]=n;
return;
}
LL p
=n;
while(p>=n)
{
p
=pollardRho(p,k--);
}
findFactor(p,k);
findFactor(n
/p,k);
}

int main()
{

long long a,b;
long long nMin=LLONG_MAX;
LL x,y;
map
<LL,LL> hash;
while(scanf("%I64d %I64d",&a,&b)!=EOF)
{
if(a==b)
{
printf(
"%I64d %I64d\n",a,b);
continue;
}
hash.clear();
b
/=a;
factor[
0]=1;
cnt
=0;
findFactor(b,
200);
if(0==cnt)
{
printf(
"%I64d %I64d\n",a,a*b);
continue;
}
for(int i=0;i<cnt;i++)
{
if(hash.find(factor[i])==hash.end())
{
hash[factor[i]]
=factor[i];
}
else
{
hash[factor[i]]
*=factor[i];
}
}
cnt
=0;
for(map<LL,LL>::iterator itr=hash.begin();itr!=hash.end();itr++)
{
factor[cnt
++]=itr->second;
}
LL tmp;
nMin
=LLONG_MAX;
for(int i=0;i<(1<<cnt);i++)
{
tmp
=1;
for(int j=0;j<cnt;j++)
{
if(i&(1<<j))
tmp
*=factor[j];
}
//printf("tmp %I64d\n",tmp);
if(nMin>tmp+b/tmp)
{
nMin
=tmp+b/tmp;
x
=min(tmp,b/tmp);
y
=max(tmp,b/tmp);
}
}
printf(
"%I64d %I64d\n",x*a,y*a);





}


return 0;
}
posted on 2011-04-18 15:05  Jackiesteed  阅读(442)  评论(0编辑  收藏  举报