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