Annotation of OpenXM_contrib/pari-2.2/examples/squfof.gp, Revision 1.2
1.2 ! noro 1: \\ return one non-trivial divisor of n > 1 using Shanks's SQUFOF
! 2: squfof(n) =
1.1 noro 3: {
1.2 ! noro 4: if (isprime(n) || issquare(n, &n), return(n));
1.1 noro 5:
1.2 ! noro 6: p = factor(n,0)[1,1];
! 7: if (p != n, return(p));
1.1 noro 8:
9: if (n%4==1,
1.2 ! noro 10: D = n; d = sqrtint(D); b = (((d-1)\2) << 1) + 1
1.1 noro 11: ,
1.2 ! noro 12: D = n << 2; d = sqrtint(D); b = (d\2) << 1
1.1 noro 13: );
1.2 ! noro 14: f = Qfb(1, b, (b^2-D)>>2);
! 15: l = sqrtint(d);
1.1 noro 16:
1.2 ! noro 17: q = []; lq = 0; i = 0;
1.1 noro 18: while (1,
1.2 ! noro 19: i++;
! 20: f = qfbred(f, 3, D, d);
! 21: a = component(f, 1);
! 22: if (!(i%2) && issquare(a, &as),
! 23: j = 1;
! 24: while (j<=lq,
! 25: if (as == q[j], break);
! 26: j++
1.1 noro 27: );
1.2 ! noro 28: if (j > lq, break)
1.1 noro 29: );
30:
1.2 ! noro 31: if (abs(a) <= l,
! 32: q = concat(q, abs(a));
! 33: print(q); lq++;
1.1 noro 34: )
35: );
36:
1.2 ! noro 37: print("i = ", i); print(f);
! 38: bb = component(f, 2);
! 39: gs = gcd([as, bb, D]);
! 40: if (gs > 1, return (gs));
! 41:
! 42: f = Qfb(as, -bb, as*component(f,3));
! 43: g = qfbred(f, 3, D, d);
! 44: b = component(g, 2);
! 45: until (b1 == b,
! 46: b1 = b; g = qfbred(g, 3, D, d);
! 47: b = component(g, 2);
1.1 noro 48: );
1.2 ! noro 49: a = abs(component(g, 1));
! 50: if (a % 2, a, a>>1);
1.1 noro 51: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>