[BACK]Return to squfof.gp CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari / examples
File: [local] / OpenXM_contrib / pari / examples / Attic / squfof.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.

squfof(n)=
{

  if (isprime(n), return(n));
  if (issquare(n), return(sqrtint(n)));

  p=factor(n,0)[1,1];
  if (p!=n, return(p));

  if (n%4==1,
    dd=n; d=sqrtint(dd); b=(((d-1)\2)<<1) + 1
  ,
    dd=n<<2; d=sqrtint(dd); b=(d\2)<<1
  );
  f=Qfb(1,b,(b^2-dd)>>2,0.);
  q=[]; lq=0; ii=0; l=sqrtint(d);

  while (1,
    f=qfbred(f,3,dd,d);
    ii++; a=component(f,1);
    if (!(ii%2),
      if (issquare(a),
	as=sqrtint(a); j=1;
	while (j<=lq,
	  if (as==q[j], break);
	  j++
	);
	if (j>lq, break)
      );
    );

    if (abs(a)<=l,
	q=vecconcat(q,abs(a));
	print(q); lq++;
    )
  );

  gs=gcd(as,dd);
  print("i = ",ii); print(f);
  gs=gcd(gs,bb=component(f,2));
  if (gs>1, return (gs));

  g=qfbred(Qfb(as,-bb,as*component(f,3),0.),3,dd,d);
  b=component(g,2);
  until (b1==b,
    b1=b; g=qfbred(g,3,dd,d);
    b=component(g,2);
  );
  a=abs(component(g,1));
  if (!(a%2),a>>=1);
  a
}