Return to rho.gp CVS log

Up to [local] / OpenXM_contrib / pari / examples

Annotation of OpenXM_contrib/pari/examples/rho.gp, Revision 1.1.1.1

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