[BACK]Return to squfof.gp CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / examples

Diff for /OpenXM_contrib/pari-2.2/examples/Attic/squfof.gp between version 1.1 and 1.2

version 1.1, 2001/10/02 11:16:59 version 1.2, 2002/09/11 07:26:45
Line 1 
Line 1 
 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));    p = factor(n,0)[1,1];
   if (issquare(n), return(sqrtint(n)));    if (p != n, return(p));
   
   p=factor(n,0)[1,1];  
   if (p!=n, return(p));  
   
   if (n%4==1,    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.);    f = Qfb(1, b, (b^2-D)>>2);
   q=[]; lq=0; ii=0; l=sqrtint(d);    l = sqrtint(d);
   
     q = []; lq = 0; i = 0;
   while (1,    while (1,
     f=qfbred(f,3,dd,d);      i++;
     ii++; a=component(f,1);      f = qfbred(f, 3, D, d);
     if (!(ii%2),      a = component(f, 1);
       if (issquare(a),      if (!(i%2) && issquare(a, &as),
         as=sqrtint(a); j=1;        j = 1;
         while (j<=lq,        while (j<=lq,
           if (as==q[j], break);          if (as == q[j], break);
           j++          j++
         );  
         if (j>lq, break)  
       );        );
         if (j > lq, break)
     );      );
   
     if (abs(a)<=l,      if (abs(a) <= l,
         q=concat(q,abs(a));        q = concat(q, abs(a));
         print(q); lq++;        print(q); lq++;
     )      )
   );    );
   
   gs=gcd(as,dd);    print("i = ", i); print(f);
   print("i = ",ii); print(f);    bb = component(f, 2);
   gs=gcd(gs,bb=component(f,2));    gs = gcd([as, bb, D]);
   if (gs>1, return (gs));    if (gs > 1, return (gs));
   
   g=qfbred(Qfb(as,-bb,as*component(f,3),0.),3,dd,d);    f = Qfb(as, -bb, as*component(f,3));
   b=component(g,2);    g = qfbred(f, 3, D, d);
   until (b1==b,    b = component(g, 2);
     b1=b; g=qfbred(g,3,dd,d);    until (b1 == b,
     b=component(g,2);      b1 = b; g = qfbred(g, 3, D, d);
       b = component(g, 2);
   );    );
   a=abs(component(g,1));    a = abs(component(g, 1));
   if (!(a%2),a>>=1);    if (a % 2, a, a>>1);
   a  
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>