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     ! 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: }