SEGV in factorint (fwd)

[Gerhard Niklasch <nikl@mathematik.tu-muenchen.de>]

With Karim's blessing, here's the patch that cures the occasional
segfaults in factorint(_,1) in 2.0.9.alpha.  (It is still possible
after applying this to cause memory corruption and other funny things
by hitting Ctrl-C in the middle of an MPQS run;  this is known and
fixed for 2.0.10.  If you want to see the factorizers at work, type
first
gp > default(debug,4);default(primelimit,66000);
and then a series of factorint(_,2)s.  Higher debug levels are even
more verbose;  at level 7, MPQS will create nice visual effects by
printing the huge binary matrices before and after Gauss over F_2...)

Enjoy, Gerhard


Forwarded message:
From: Karim BELABAS <Karim.Belabas@math.u-psud.fr>
Date: 	Wed, 17 Jun 1998 13:33:14 +0200
Message-Id: <199806171133.NAA01347@geo.math.u-psud.fr>
To: nikl@mathematik.tu-muenchen.de
Subject: SEGV in factorint


should be cured by the following patch  [GN: remainder of commentary
censored -- Karim asks you NOT to read the code being replaced by
this ;^]

*** src/basemath/ifactor1.c.orig	Tue Jun 16 16:44:49 1998
--- src/basemath/ifactor1.c	Wed Jun 17 13:27:24 1998
***************
*** 1,6 ****
  /* $Id: ifactor1.c,v 2.0.0.9 1998/06/09 11:47:19 belabas Exp belabas $ */
  #include "pari.h"
- static GEN N;
  
  /*********************************************************************/
  /**                                                                 **/
--- 1,5 ----
***************
*** 10,47 ****
  static GEN sqrt1, sqrt2, t1, t;
  static long r1;
  
! static void
  init_miller(GEN n)
  {
-   N = (signe(n) < 0)? absi(n): n;
    t=subis(n,1); r1=vali(t); t1 = shifti(t,-r1);
    sqrt1=cgeti(lg(t)); sqrt1[1]=evalsigne(0)|evallgefint(2);
    sqrt2=cgeti(lg(t)); sqrt2[1]=evalsigne(0)|evallgefint(2);
  }
  
! /* is N strong pseudo-prime for base a ? `End matching' (check for square
   * roots of -1) added by GN */
  static int
! bad_for_base(GEN a)
  {
!   GEN c2, c = powmodulo(a,t1,N);
    long r;
  
    if (!is_pm1(c) && !egalii(t,c)) /* go fishing for -1, not for 1 */
    {
      for (r=r1-1; r; r--)	/* (this saves one squaring/reduction) */
      {
!       c2=c; c=modii(sqri(c),N);
        if (egalii(t,c)) break;
      }
      if (!r) return 1;
      /* sqrt(-1) seen, compare or remember */
      if (signe(sqrt1))		/* we saw one earlier: compare */
      {
!       /* check if too many sqrt(-1)s mod N */
        if (!egalii(c2,sqrt1) && !egalii(c2,sqrt2)) return 1; 
      }
!     else { affii(c2,sqrt1); affii(subii(N,c2),sqrt2); } /* remember */
    }
    return 0;
  }
--- 9,46 ----
  static GEN sqrt1, sqrt2, t1, t;
  static long r1;
  
! static GEN
  init_miller(GEN n)
  {
    t=subis(n,1); r1=vali(t); t1 = shifti(t,-r1);
    sqrt1=cgeti(lg(t)); sqrt1[1]=evalsigne(0)|evallgefint(2);
    sqrt2=cgeti(lg(t)); sqrt2[1]=evalsigne(0)|evallgefint(2);
+   return (signe(n) < 0)? absi(n): n;
  }
  
! /* is n strong pseudo-prime for base a ? `End matching' (check for square
   * roots of -1) added by GN */
  static int
! bad_for_base(GEN n, GEN a)
  {
!   GEN c2, c = powmodulo(a,t1,n);
    long r;
  
    if (!is_pm1(c) && !egalii(t,c)) /* go fishing for -1, not for 1 */
    {
      for (r=r1-1; r; r--)	/* (this saves one squaring/reduction) */
      {
!       c2=c; c=modii(sqri(c),n);
        if (egalii(t,c)) break;
      }
      if (!r) return 1;
      /* sqrt(-1) seen, compare or remember */
      if (signe(sqrt1))		/* we saw one earlier: compare */
      {
!       /* check if too many sqrt(-1)s mod n */
        if (!egalii(c2,sqrt1) && !egalii(c2,sqrt2)) return 1; 
      }
!     else { affii(c2,sqrt1); affii(subii(n,c2),sqrt2); } /* remember */
    }
    return 0;
  }
***************
*** 57,67 ****
    if (lgefint(n)==3 && (ulong)(n[2])<=3) return (n[2] != 1); 
  
    if (!mod2(n)) return 0;
!   init_miller(n); av2=avma;
    for (i=1; i<=k; i++)
    {
      do r = smodsi(mymyrand(),n); while (!r);
!     if (bad_for_base(stoi(r))) { avma=av; return 0; }
      avma=av2;
    }
    avma=av; return 1;
--- 56,66 ----
    if (lgefint(n)==3 && (ulong)(n[2])<=3) return (n[2] != 1); 
  
    if (!mod2(n)) return 0;
!   n = init_miller(n); av2=avma;
    for (i=1; i<=k; i++)
    {
      do r = smodsi(mymyrand(),n); while (!r);
!     if (bad_for_base(n, stoi(r))) { avma=av; return 0; }
      avma=av2;
    }
    avma=av; return 1;
***************
*** 78,88 ****
    static long pr[] = {0, 2,3,5,7,11,13,17,19,23,29};
  
    if (!mod2(n)) return 0;
!   init_miller(n); av2=avma;
    for (i=1; i<=k; i++)
    {
      r = smodsi(pr[i],n); if (!r) break;
!     if (bad_for_base(stoi(r))) { avma = av; return 0; }
      avma=av2;
    }
    avma=av; return 1;
--- 77,87 ----
    static long pr[] = {0, 2,3,5,7,11,13,17,19,23,29};
  
    if (!mod2(n)) return 0;
!   n = init_miller(n); av2=avma;
    for (i=1; i<=k; i++)
    {
      r = smodsi(pr[i],n); if (!r) break;
!     if (bad_for_base(n, stoi(r))) { avma = av; return 0; }
      avma=av2;
    }
    avma=av; return 1;
***************
*** 163,169 ****
  /**   is composite, and has no small prime divisor (P^-(N) > 65.536). **/
  /**                                                                   **/
  /***********************************************************************/
! static GEN gl, *W;
  #define nbc 10
  
  static int
--- 162,168 ----
  /**   is composite, and has no small prime divisor (P^-(N) > 65.536). **/
  /**                                                                   **/
  /***********************************************************************/
! static GEN N, gl, *W;
  #define nbc 10
  
  static int
--
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (33 1) 01 69 15 57 48
F-91405 Orsay (France)           Fax: (33 1) 01 69 15 60 19
--
PARI/GP Home Page: http://pari.home.ml.org