Return to squfof.gp CVS log

Up to [local] / OpenXM_contrib / pari / examples

Annotation of OpenXM_contrib/pari/examples/squfof.gp, Revision 1.1.1.1

1.1       maekawa     1: squfof(n)=
                      2: {
                      3:
                      4:   if (isprime(n), return(n));
                      5:   if (issquare(n), return(sqrtint(n)));
                      6:
                      7:   p=factor(n,0)[1,1];
                      8:   if (p!=n, return(p));
                      9:
                     10:   if (n%4==1,
                     11:     dd=n; d=sqrtint(dd); b=(((d-1)\2)<<1) + 1
                     12:   ,
                     13:     dd=n<<2; d=sqrtint(dd); b=(d\2)<<1
                     14:   );
                     15:   f=Qfb(1,b,(b^2-dd)>>2,0.);
                     16:   q=[]; lq=0; ii=0; l=sqrtint(d);
                     17:
                     18:   while (1,
                     19:     f=qfbred(f,3,dd,d);
                     20:     ii++; a=component(f,1);
                     21:     if (!(ii%2),
                     22:       if (issquare(a),
                     23:        as=sqrtint(a); j=1;
                     24:        while (j<=lq,
                     25:          if (as==q[j], break);
                     26:          j++
                     27:        );
                     28:        if (j>lq, break)
                     29:       );
                     30:     );
                     31:
                     32:     if (abs(a)<=l,
                     33:        q=vecconcat(q,abs(a));
                     34:        print(q); lq++;
                     35:     )
                     36:   );
                     37:
                     38:   gs=gcd(as,dd);
                     39:   print("i = ",ii); print(f);
                     40:   gs=gcd(gs,bb=component(f,2));
                     41:   if (gs>1, return (gs));
                     42:
                     43:   g=qfbred(Qfb(as,-bb,as*component(f,3),0.),3,dd,d);
                     44:   b=component(g,2);
                     45:   until (b1==b,
                     46:     b1=b; g=qfbred(g,3,dd,d);
                     47:     b=component(g,2);
                     48:   );
                     49:   a=abs(component(g,1));
                     50:   if (!(a%2),a>>=1);
                     51:   a
                     52: }