Pollard Rho prime factorization
in [General Programming]
|
From: mukesh tiwari on 4 Feb 2010 15:36 Hello Everyone. I am trying to solve a problem on SPOJ (https:// www.spoj.pl/problems/FACT0/)but i am getting time limit exceed.I implemented pollard rho algorithm given on wikipedia.(http:// en.wikipedia.org/wiki/Pollard%27s_rho_algorithm) Is there any other faster version of this algorithm.Here is the full code but only major part is Pollard and factor function. Thank you #include <vector> #include <list> #include <map> #include <set> #include <deque> #include <queue> #include <stack> #include <bitset> #include <algorithm> #include <functional> #include <numeric> #include <utility> #include <sstream> #include <iostream> #include <iomanip> #include <cstdio> #include <cmath> #include <cstdlib> #include <cctype> #include <string> #include <cstring> #include <cstdio> #include <cmath> #include <cstdlib> #include <ctime> using namespace std; typedef long long ll; bool B[1000000]; vector<long long>v,p; void isprime() { B[0]=B[1]=1; for(int i=2;i*i<1000000;i++) if(!B[i]) for(int j=i*i;j<1000000;j+=i) B[j]=1; for(int i=0;i<1000000;i++) if(!B[i]) p.push_back((long long)i); //for(int i=0;i<10;i++) cout<<p[i]<<" "; //cout<<endl; } long long mulmod(long long a,long long b,long long c){ long long x = 0,y=a%c; while(b > 0){ if(b%2 == 1){ x = (x+y)%c; } y = (y*2)%c; b /= 2; } return x%c; } ll modulo(ll a,ll b,ll c){ ll x=1,y=a; // ll is taken to avoid overflow of intermediate results while(b > 0){ if(b%2 == 1){ x=mulmod(x,y,c); } y = mulmod(y,y,c); // squaring the base b /= 2; } return x%c; } bool Miller(long long p,int iteration){ if(p<2){ return false; } if(p!=2 && p%2==0){ return false; } long long s=p-1; while(s%2==0){ s/=2; } for(int i=0;i<iteration;i++){ long long a=rand()%(p-1)+1,temp=s; long long mod=modulo(a,temp,p); while(temp!=p-1 && mod!=1 && mod!=p-1){ mod=mulmod(mod,mod,p); temp *= 2; } if(mod!=p-1 && temp%2==0){ return false; } } return true; } ll PollardRho(ll n) { //cout<<n<<endl; if(n%2==0) return 2LL; ll x=(ll) random(); ll y=x; ll c=(ll)random(); ll d; do { x=((x*x)%n+c)%n; //if(x<0) cout<<x<<endl; y=((y*y)%n+c)%n; y=((y*y)%n+c)%n; d=__gcd(abs(x-y),n); //cout<<"d= "<<d<<endl; }while(d==1); return d; } void factor(ll n) { if(n==1) return ; if(Miller(n,1)) { v.push_back(n); return ;} ll div=PollardRho(n); factor(div); factor(n/div); } int main() { ll n; isprime(); while(scanf("%lld",&n) && n!=0) { v.clear(); for(int i=0;i<p.size() && p[i]*p[i]<=n;) { if(n%p[i]==0) { v.push_back(p[i]); n/=p[i]; } else i++; } //cout<<n<<endl; factor(n); sort(v.begin(),v.end()); v.push_back(0); int p=1; for(int i=0;i+1<v.size();i++) { if(v[i]==v[i+1]) p++; else { cout<<v[i]<<"^"<<p; if(i+2<v.size()) printf(" "); else printf("\n"); p=1; } } } }
From: Richard Heathfield on 4 Feb 2010 18:52 mukesh tiwari wrote: > Hello Everyone. I am trying to solve a problem on SPOJ (https:// > www.spoj.pl/problems/FACT0/)but i am getting time limit exceed.I > implemented pollard rho algorithm given on wikipedia.(http:// > en.wikipedia.org/wiki/Pollard%27s_rho_algorithm) Is there any other > faster version of this algorithm.Here is the full code but only major > part is Pollard and factor function. > Thank you <snip> > void isprime() > { > B[0]=B[1]=1; > for(int i=2;i*i<1000000;i++) > if(!B[i]) for(int j=i*i;j<1000000;j+=i) B[j]=1; Once you've dealt with 2, you only need worry about odd numbers. In fact, once you've dealt with 3, you need only test numbers of the form 6k-1 and 6k+1. That should triple the speed of your loop. -- Richard Heathfield <http://www.cpax.org.uk> Email: -http://www. +rjh@ "Usenet is a strange place" - dmr 29 July 1999 Sig line vacant - apply within
From: websnarf on 4 Feb 2010 20:31 On Feb 4, 12:36 pm, mukesh tiwari <mukeshtiwari.ii...(a)gmail.com> wrote: > Hello Everyone. I am trying to solve a problem on SPOJ (https://www.spoj.pl/problems/FACT0/)buti am getting time limit exceed.I > implemented pollard rho algorithm given on wikipedia.(http:// > en.wikipedia.org/wiki/Pollard%27s_rho_algorithm) Is there any other > faster version of this algorithm.Here is the full code but only major > part is Pollard and factor function. > Thank you The problem is that your mulmod() is too slow. You are doing multiplication by bit test, shift and add or something like that. That's a horrible way of doing multiplication to a long word with a modulo. You are doing one % per on-bit of the input. The / % operators are the slowest of the basic operations of your CPU. One thing you can do is implement a hard coded two-digit Karatsuba's method for the multiplication, to reduce your 128-bit output to calculating 3 64-bit outputs (with 32-bit inputs.) Then implement a "two-digit division" using successive 64 bit division estimations. I am sure you would never need to do more than 3 divisions, and can probably get away with 2 (but I have not tried it myself to be sure.) The other more standard method that is widely used is something called "Barrett modular reduction". I have not implemented this myself, but apparently that's the algorithm of choice for this. If you can guarantee that the integers are 32 bits or less, then you can use my implementation: http://www.azillionmonkeys.com/qed/primeat.zip which includes a primality test. Its for C, but easily adapted to Java or C++. Also note, that if your machine is actually running in 32 bit mode instead of 64 bit mode, but you are using these 64-bit primitives, then your results will be about 4 times slower. -- Paul Hsieh http://www.pobox.com/~qed/ http://bstring.sf.net/ > #include <vector> > #include <list> > #include <map> > #include <set> > #include <deque> > #include <queue> > #include <stack> > #include <bitset> > #include <algorithm> > #include <functional> > #include <numeric> > #include <utility> > #include <sstream> > #include <iostream> > #include <iomanip> > #include <cstdio> > #include <cmath> > #include <cstdlib> > #include <cctype> > #include <string> > #include <cstring> > #include <cstdio> > #include <cmath> > #include <cstdlib> > #include <ctime> > using namespace std; > typedef long long ll; > bool B[1000000]; > vector<long long>v,p; > > void isprime() > { > B[0]=B[1]=1; > for(int i=2;i*i<1000000;i++) > if(!B[i]) for(int j=i*i;j<1000000;j+=i) B[j]=1; > > for(int i=0;i<1000000;i++) if(!B[i]) p.push_back((long long)i); > //for(int i=0;i<10;i++) cout<<p[i]<<" "; > //cout<<endl; > } > > long long mulmod(long long a,long long b,long long c){ > long long x = 0,y=a%c; > while(b > 0){ > if(b%2 == 1){ > x = (x+y)%c; > } > y = (y*2)%c; > b /= 2; > } > return x%c; > > } > > ll modulo(ll a,ll b,ll c){ > ll x=1,y=a; // ll is taken to avoid overflow of intermediate > results > while(b > 0){ > if(b%2 == 1){ > x=mulmod(x,y,c); > } > y = mulmod(y,y,c); // squaring the base > b /= 2; > } > return x%c;} > > bool Miller(long long p,int iteration){ > if(p<2){ > return false; > } > if(p!=2 && p%2==0){ > return false; > } > long long s=p-1; > while(s%2==0){ > s/=2; > } > for(int i=0;i<iteration;i++){ > long long a=rand()%(p-1)+1,temp=s; > long long mod=modulo(a,temp,p); > while(temp!=p-1 && mod!=1 && mod!=p-1){ > mod=mulmod(mod,mod,p); > temp *= 2; > } > if(mod!=p-1 && temp%2==0){ > return false; > } > } > return true; > > } > > ll PollardRho(ll n) > { > //cout<<n<<endl; > if(n%2==0) return 2LL; > ll x=(ll) random(); > ll y=x; > ll c=(ll)random(); > ll d; > do > { > x=((x*x)%n+c)%n; > //if(x<0) cout<<x<<endl; > y=((y*y)%n+c)%n; > y=((y*y)%n+c)%n; > d=__gcd(abs(x-y),n); > //cout<<"d= "<<d<<endl; > }while(d==1); > return d; > } > > void factor(ll n) > { > if(n==1) return ; > if(Miller(n,1)) { v.push_back(n); return ;} > ll div=PollardRho(n); > factor(div); > factor(n/div); > } > > int main() > { > > ll n; > isprime(); > while(scanf("%lld",&n) && n!=0) > { > v.clear(); > for(int i=0;i<p.size() && p[i]*p[i]<=n;) > { > if(n%p[i]==0) > { > v.push_back(p[i]); > n/=p[i]; > } > else i++; > } > //cout<<n<<endl; > factor(n); > sort(v.begin(),v.end()); > v.push_back(0); > int p=1; > for(int i=0;i+1<v.size();i++) > { > if(v[i]==v[i+1]) p++; > else > { > cout<<v[i]<<"^"<<p; > if(i+2<v.size()) printf(" "); > else printf("\n"); > p=1; > } > } > > } > }
From: aslan on 12 Feb 2010 05:42 "mukesh tiwari" <mukeshtiwari.iiitm(a)gmail.com>, iletisinde sunu yazdi, news:d618772c-845c-47e5-9d94-e836c92fce3a(a)m24g2000prn.googlegroups.com... > Hello Everyone. I am trying to solve a problem on SPOJ (https:// > www.spoj.pl/problems/FACT0/)but i am getting time limit exceed.I > implemented pollard rho algorithm given on wikipedia.(http:// > en.wikipedia.org/wiki/Pollard%27s_rho_algorithm) Is there any other > faster version of this algorithm.Here is the full code but only major > part is Pollard and factor function. > Thank you > > #include <vector> > #include <list> > #include <map> > #include <set> > #include <deque> > #include <queue> > #include <stack> > #include <bitset> > #include <algorithm> > #include <functional> > #include <numeric> > #include <utility> > #include <sstream> > #include <iostream> > #include <iomanip> > #include <cstdio> > #include <cmath> > #include <cstdlib> > #include <cctype> > #include <string> > #include <cstring> > #include <cstdio> > #include <cmath> > #include <cstdlib> > #include <ctime> > using namespace std; > typedef long long ll; > bool B[1000000]; > vector<long long>v,p; > > void isprime() > { > B[0]=B[1]=1; > for(int i=2;i*i<1000000;i++) > if(!B[i]) for(int j=i*i;j<1000000;j+=i) B[j]=1; > > for(int i=0;i<1000000;i++) if(!B[i]) p.push_back((long long)i); > //for(int i=0;i<10;i++) cout<<p[i]<<" "; > //cout<<endl; > } > > > long long mulmod(long long a,long long b,long long c){ > long long x = 0,y=a%c; > while(b > 0){ > if(b%2 == 1){ > x = (x+y)%c; > } > y = (y*2)%c; > b /= 2; > } > return x%c; > } > > ll modulo(ll a,ll b,ll c){ > ll x=1,y=a; // ll is taken to avoid overflow of intermediate > results > while(b > 0){ > if(b%2 == 1){ > x=mulmod(x,y,c); > } > y = mulmod(y,y,c); // squaring the base > b /= 2; > } > return x%c; > } > bool Miller(long long p,int iteration){ > if(p<2){ > return false; > } > if(p!=2 && p%2==0){ > return false; > } > long long s=p-1; > while(s%2==0){ > s/=2; > } > for(int i=0;i<iteration;i++){ > long long a=rand()%(p-1)+1,temp=s; > long long mod=modulo(a,temp,p); > while(temp!=p-1 && mod!=1 && mod!=p-1){ > mod=mulmod(mod,mod,p); > temp *= 2; > } > if(mod!=p-1 && temp%2==0){ > return false; > } > } > return true; > } > > > ll PollardRho(ll n) > { > //cout<<n<<endl; > if(n%2==0) return 2LL; > ll x=(ll) random(); > ll y=x; > ll c=(ll)random(); > ll d; > do > { > x=((x*x)%n+c)%n; > //if(x<0) cout<<x<<endl; > y=((y*y)%n+c)%n; > y=((y*y)%n+c)%n; > d=__gcd(abs(x-y),n); > //cout<<"d= "<<d<<endl; > }while(d==1); > return d; > } > > void factor(ll n) > { > if(n==1) return ; > if(Miller(n,1)) { v.push_back(n); return ;} > ll div=PollardRho(n); > factor(div); > factor(n/div); > } > > int main() > { > > > ll n; > isprime(); > while(scanf("%lld",&n) && n!=0) > { > v.clear(); > for(int i=0;i<p.size() && p[i]*p[i]<=n;) > { > if(n%p[i]==0) > { > v.push_back(p[i]); > n/=p[i]; > } > else i++; > } > //cout<<n<<endl; > factor(n); > sort(v.begin(),v.end()); > v.push_back(0); > int p=1; > for(int i=0;i+1<v.size();i++) > { > if(v[i]==v[i+1]) p++; > else > { > cout<<v[i]<<"^"<<p; > if(i+2<v.size()) printf(" "); > else printf("\n"); > p=1; > } > } > > } > } This one outputs the prime factorisation for the given example: #include <stdlib.h> #include <stdio.h> #include <math.h> #include <memory.h> #include <string.h> #include <time.h> typedef __int64 int64; struct timer { clock_t t; timer() { t = clock(); } ~timer() { printf("runtime %.3f secs\n", getTime()); } double getTime() { return ((double)clock()-(double)t)/(double)CLOCKS_PER_SEC; } }; void get_divisors(int64 n) { printf("%I64d = ",n); if (!(n&1)) { int np=1; while (n/=2, !(n&1)) ++np; printf("2^%d ", np); } if (!(n%3)) { int np=1; while (n/=3, !(n%3)) ++np; printf("3^%d ", np); } int64 p=1, a=2, hi=sqrt(n)+1; while (a=6-a, p+=a, p<=hi) { if (!(n%p)) { int np=1; while (n/=p, !(n%p)) ++np; printf("%I64d^%d ", p, np); } } if (n>1) printf("%I64d^1 ", n); printf("\n"); } int main(int argc, char**argv) { timer t; get_divisors(3111989); get_divisors(13091989); get_divisors(2432902008176640000); get_divisors(77145199750673); return 0; } /* 3111989 = 317^1 9817^1 13091989 = 17^2 89^1 509^1 2432902008176640000 = 2^18 3^8 5^4 7^2 11^1 13^1 17^1 19^1 77145199750673 = 328439^1 234884407^1 runtime 0.048 secs Press any key to continue */
|
Pages: 1 Prev: Draft paper submission deadline is extended: SETP-10 Next: Blender's UI toolkit? |