Re: removeprimes()

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

[Igor:]
> I noticed this interesting behavior.
> 
> ? addprimes(2*primes(64));length(addprimes)
> 64
> ? removeprimes(addprimes);length(addprimes)
> 32
[...]
> So removeprimes() cuts the size of the list of
> declared primes in half.  I don't this this is how
> removeprimes() is described to behave.

Stupid bug (aren't they all :-). The argument to removeprimes() happens to be
the same (same pointer, not a copy) as the table we modify. This patch
special cases this situation (and cleans up the code while I was at it).

  Karim.

Index: src/basemath/arith2.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/basemath/arith2.c,v
retrieving revision 1.2
diff -c -r1.2 arith2.c
*** src/basemath/arith2.c	1999/10/15 11:30:40	1.2
--- src/basemath/arith2.c	1999/10/21 16:04:05
***************
*** 265,294 ****
    avma=av; return primetab;
  }
  
  GEN
  removeprimes(GEN prime)
  {
!   long i,tx, av = avma;
  
    if (!prime) return primetab;
    tx = typ(prime);
    if (is_vec_t(tx))
    {
!     for (i=1; i < lg(prime); i++) removeprimes((GEN) prime[i]);
!     return primetab;
!   }
!   if (tx != t_INT) err(typeer,"removeprimes");
!   if (is_pm1(prime)) return primetab;
!   if (signe(prime) < 0) prime=absi(prime);
! 
!   for (i=1; i < lg(primetab); i++)
!     if (egalii((GEN) primetab[i], prime))
      {
!       gunclone((GEN)primetab[i]); primetab[i]=0;
!       cleanprimetab(); avma=av; return primetab;
      }
!   err(talker,"prime is not in primetable in removeprimes");
!   return NULL; /* not reached */
  }
  
  /***********************************************************************/
--- 265,307 ----
    avma=av; return primetab;
  }
  
+ static GEN
+ removeprime(GEN prime)
+ {
+   long i;
+ 
+   if (typ(prime) != t_INT) err(typeer,"removeprime");
+   for (i=lg(primetab) - 1; i; i--)
+     if (absi_equal((GEN) primetab[i], prime))
+     {
+       gunclone((GEN)primetab[i]); primetab[i]=0;
+       cleanprimetab(); break;
+     }
+   if (!i) err(talker,"prime %Z is not in primetable", prime);
+   return primetab;
+ }
+ 
  GEN
  removeprimes(GEN prime)
  {
!   long i,tx;
  
    if (!prime) return primetab;
    tx = typ(prime);
    if (is_vec_t(tx))
    {
!     if (prime == primetab)
      {
!       for (i=1; i < lg(prime); i++) gunclone((GEN)prime[i]);
!       setlg(prime, 1);
      }
!     else
!     {
!       for (i=1; i < lg(prime); i++) removeprime((GEN) prime[i]);
!     }
!     return primetab;
!   }
!   return removeprime(prime);
  }
  
  /***********************************************************************/
__
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://hasse.mathematik.tu-muenchen.de/ntsw/pari/