Mark Dickinson on Wed, 18 Jul 2001 18:26:20 -0400 (EDT)


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

polrootspadic bug resolved


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
 


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/07/18 22:09:51
***************
*** 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,2168 ----
    return gerepileuptoint(av, r); /* = resii(x, y) */
  }
  
! /* x % (2^n), assuming x, n >= 0 */
  GEN
  resmod2n(GEN x, long n)
  {
!   long l,k,lx,ly,top;
    GEN y, z, t;
  
!   if (!signe(x) || !n) return gzero;
! 
!   l = (n-1) % BITS_IN_LONG + 1;
!   k = (n-1) / BITS_IN_LONG;
    lx = lgefint(x);
!   if (lx < k+3) return icopy(x);
! 
!   top = x[lx-k-1] & ((1<<l)-1);
!   if (!top) {   /* strip leading zeros from the result */
!     while(k && !x[lx-k]) k--;
!     if (!k) return gzero;
!     ly = k+2;
!   } else {
!     ly = k+3;
!   }
! 
    z = cgeti(ly);
    z[1] = evalsigne(1)|evallgefint(ly);
    y = z + ly;
    t = x + lx;
    for ( ;k; k--) *--y = *--t;
!   if(top) *--y = top;
    return z;
  }