Re: polrootspadic bug resolved

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

On Wed, 18 Jul 2001, Mark Dickinson wrote:
> The problem with polrootspadic that I reported earlier this week can be
> traced all the way back to a bug in resmod2n() in src/kernel/none/mp.c;
> this function computes x % 2^n for an integer x, but it turns out that it
> can sometimes fail to normalise the result properly (that is, the most
> significant word of the integer returned can sometimes be zero). The
> following patch seems to fix this for me; I hope it works for others as
> well  :)

Mark, I reworked your patch a little bit (in fact I had lost your patch and my net
access, so I had to start over from scratch :-) : I let l run in [0,BIL[ instead of
[1,BIL], used pointer arithmetic all the way through, and tried to help the
optimizer making the good choices.

I didn't run into discrepancies in any of the test suites or your original
examples; can you check that I didn't blunder ?

  Karim.

Index: src/kernel/none/mp.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/kernel/none/mp.c,v
retrieving revision 1.44
diff -c -r1.44 mp.c
*** src/kernel/none/mp.c	2001/07/02 20:54:10	1.44
--- src/kernel/none/mp.c	2001/08/26 17:21:28
***************
*** 2134,2163 ****
    return gerepileuptoint(av, r); /* = resii(x, y) */
  }

! /* x % (2^n) */
  GEN
  resmod2n(GEN x, long n)
  {
!   long l,k,lx,ly;
!   GEN y, z, t;

!   if (!signe(x)) return gzero;

!   l = n % BITS_IN_LONG;
!   k = n / BITS_IN_LONG;
!   ly = l? k+3: k+2;
    lx = lgefint(x);
!   if (lx < ly) return icopy(x);
!   z = cgeti(ly);
!   z[1] = evalsigne(1)|evallgefint(ly);
!   y = z + ly;
!   t = x + lx;
!   for ( ;k; k--) *--y = *--t;
!   if (l) *--y = *--t & ((1<<l) - 1);
! #if 0
!   if (!egalii(z, resii(x, shifti(gun,n))))
!     err(talker,"bug in resmod2n: n = %ld, x = %Z\n",n,x);
! #endif
    return z;
  }

--- 2134,2170 ----
    return gerepileuptoint(av, r); /* = resii(x, y) */
  }

! /* x % (2^n), assuming x, n >= 0 */
  GEN
  resmod2n(GEN x, long n)
  {
!   long hi,l,k,lx,ly;
!   GEN z, xd, zd;

!   if (!signe(x) || !n) return gzero;

!   l = n & (BITS_IN_LONG-1);    /* n % BITS_IN_LONG */
!   k = n >> TWOPOTBITS_IN_LONG; /* n / BITS_IN_LONG */
    lx = lgefint(x);
!   if (lx < k+3) return icopy(x);
!
!   xd = x + (lx-k-1);
!   /* x = |_|...|#|1|...|k| : copy the last l bits of # and the last k words
!    *            ^--- initial xd  */
!   hi = *xd & ((1<<l) - 1); /* last l bits of # = top bits of result */
!   if (!hi)
!   { /* strip leading zeroes from result */
!     xd++; while (k && !*xd) { k--; xd++; }
!     if (!k) return gzero;
!     ly = k+2;
!   }
!   else
!     ly = k+3;
!
!   zd = z = cgeti(ly);
!   *++zd = evalsigne(1) | evallgefint(ly);
!   if (hi) *++zd = hi;
!   for ( ;k; k--) *++zd = *++xd;
    return z;
  }


-- 
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://www.parigp-home.de/