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; }