[BACK]Return to rho.gp CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari / examples
File: [local] / OpenXM_contrib / pari / examples / Attic / rho.gp (download)

Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:30 2000 UTC (24 years, 5 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2
Changes since 1.1: +0 -0 lines 

Import PARI/GP 2.0.17 beta.

rho1(n)=
{
  local(x,y);

  x=2; y=5;
  while(gcd(y-x,n)==1,
    x=(x^2+1)%n;
    y=(y^2+1)%n; y=(y^2+1)%n
  );
  gcd(n,y-x)
}

rho2(n)=
{
  local(m);

  m=rho1(n);
  if (isprime(m),
    print(m)
  ,
    rho2(m)
  );
  if (isprime(n/m),
    print(n/m)
  ,
    rho2(n/m)
  );
}

rho(n)=
{
  local(m);

  m=factor(n,0); print(m);
  n=m[length(m[,1]),1];
  if (!isprime(n),rho2(n));
}

rhobrent(n)=
{
  x=y=x1=2; k=l=p=1; c=0;
  while (1,
    x=(x^2+1)%n; p=(p*(x1-x))%n; c++;
    if (c==20,
      if (gcd(p,n)>1, break);
      y=x;c=0
    );
    k--;
    if (!k,
      if (gcd(p,n)>1, break);

      x1=x; k=l; l <<= 1; \\ l = 2*l
      for (j=1,k, x=(x^2+1)%n);
      y=x; c=0
    )
  );
  g=1;
  while (g==1,
    y=(y^2+1)%n;
    g=gcd(x1-y,n)
  );
  if (g==n,error("algorithm fails"));
  g
}