===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/examples/Attic/squfof.gp,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/examples/Attic/squfof.gp	2001/10/02 11:16:59	1.1
+++ OpenXM_contrib/pari-2.2/examples/Attic/squfof.gp	2002/09/11 07:26:45	1.2
@@ -1,52 +1,51 @@
-squfof(n)=
+\\ return one non-trivial divisor of n > 1 using Shanks's SQUFOF
+squfof(n) =
 {
+  if (isprime(n) || issquare(n, &n), return(n));
 
-  if (isprime(n), return(n));
-  if (issquare(n), return(sqrtint(n)));
+  p = factor(n,0)[1,1];
+  if (p != n, return(p));
 
-  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
+    D = n;      d = sqrtint(D); b = (((d-1)\2) << 1) + 1
   ,
-    dd=n<<2; d=sqrtint(dd); b=(d\2)<<1
+    D = n << 2; d = sqrtint(D); b = (d\2) << 1
   );
-  f=Qfb(1,b,(b^2-dd)>>2,0.);
-  q=[]; lq=0; ii=0; l=sqrtint(d);
+  f = Qfb(1, b, (b^2-D)>>2);
+  l = sqrtint(d);
 
+  q = []; lq = 0; i = 0;
   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)
+    i++;
+    f = qfbred(f, 3, D, d);
+    a = component(f, 1);
+    if (!(i%2) && issquare(a, &as),
+      j = 1;
+      while (j<=lq,
+        if (as == q[j], break);
+        j++
       );
+      if (j > lq, break)
     );
 
-    if (abs(a)<=l,
-	q=concat(q,abs(a));
-	print(q); lq++;
+    if (abs(a) <= l,
+      q = concat(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));
+  print("i = ", i); print(f);
+  bb = component(f, 2);
+  gs = gcd([as, bb, D]);
+  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);
+  f = Qfb(as, -bb, as*component(f,3));
+  g = qfbred(f, 3, D, d);
+  b = component(g, 2);
+  until (b1 == b,
+    b1 = b; g = qfbred(g, 3, D, d);
+    b = component(g, 2);
   );
-  a=abs(component(g,1));
-  if (!(a%2),a>>=1);
-  a
+  a = abs(component(g, 1));
+  if (a % 2, a, a>>1);
 }