Return to rho.gp CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / examples

Annotation of OpenXM_contrib/pari-2.2/examples/rho.gp, Revision 1.2

1.1       noro        1: rho1(n)=
1.2     ! noro        2: { local(x,y);
1.1       noro        3:
                      4:   x=2; y=5;
                      5:   while(gcd(y-x,n)==1,
                      6:     x=(x^2+1)%n;
                      7:     y=(y^2+1)%n; y=(y^2+1)%n
                      8:   );
                      9:   gcd(n,y-x)
                     10: }
                     11:
                     12: rho2(n)=
1.2     ! noro       13: { local(m);
1.1       noro       14:
                     15:   m=rho1(n);
                     16:   if (isprime(m),
                     17:     print(m)
                     18:   ,
                     19:     rho2(m)
                     20:   );
                     21:   if (isprime(n/m),
                     22:     print(n/m)
                     23:   ,
                     24:     rho2(n/m)
                     25:   );
                     26: }
                     27:
                     28: rho(n)=
1.2     ! noro       29: { local(m = factor(n,0));
1.1       noro       30:
1.2     ! noro       31:   print(m); m = m[,1]; n = m[#m];
1.1       noro       32:   if (!isprime(n),rho2(n));
                     33: }
                     34:
                     35: rhobrent(n)=
1.2     ! noro       36: { local(x,y,x1,k,l,p,c,g);
        !            37:
        !            38:   x1=x=y=2; k=l=p=1; c=0;
1.1       noro       39:   while (1,
1.2     ! noro       40:     x=(x^2+1)%n; p=(p*(x1-x))%n;
        !            41:     c++;
1.1       noro       42:     if (c==20,
                     43:       if (gcd(p,n)>1, break);
1.2     ! noro       44:       y=x; c=0
1.1       noro       45:     );
                     46:     k--;
                     47:     if (!k,
                     48:       if (gcd(p,n)>1, break);
                     49:
1.2     ! noro       50:       x1=x; k=l; l <<= 1;
1.1       noro       51:       for (j=1,k, x=(x^2+1)%n);
                     52:       y=x; c=0
                     53:     )
                     54:   );
1.2     ! noro       55:   until (g!=1,
1.1       noro       56:     y=(y^2+1)%n;
                     57:     g=gcd(x1-y,n)
                     58:   );
1.2     ! noro       59:   if (g==n, error("algorithm fails"));
1.1       noro       60:   g
                     61: }