Return to squfof.gp CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / examples

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