| Ilya Zakharevich on Sun, 7 Oct 2001 13:48:21 -0400 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| [PATCH CVS+] multi-bit and 2-adic bittest() |
This patch enables
a) extraction of several bits from a bitmap so that the result is
not a vector of bits, but a bitmap too;
b) Optional support of negative bitmaps in bittest() meaning the
same as for other bitmap operations (2s-complement).
All this via the optional argument to bittest(bitmap, offset, {bits=1}).
Enjoy,
Ilya
P.S. To support all the boundary cases turned out to be extremely messy.
Here is the test script (I ran it only to i ~= 30, running to
completion now, each i takes around 1min, and there are ~700 of them...):
{ bt(x,n,c)=
if(c<0,c=-c;x=2^1301+x);
if(n<0,x=shift(x,-n); n=0);
shift( bitand(x,2^n*(2^c-1)), -n )
}
/* Make a "random" number: 72 zero bits, 72 true bits, ~72 bits "random",
then bits with regular increments: 5, 10, 15, 20, 40, 80, 160,
then 72 true bits. */
bits_0=72; bits_1=72; bits_1a=72;
N=5^31;
random_bits = floor(log(N)/log(2)) + 1;
N = (2^bits_1 - 1)*2^bits_0 + 2^(bits_1+bits_0)*N;
N = N + 2^(bits_1 + bits_0 + random_bits)*(2^5 + 2^15 + 2^30 + 2^50 + 2^90 + 2^170 + 2^330*(2^bits_1a - 1));
bits = floor(log(N)/log(2)) + 1;
lim = bits + 70;
{ for(i=0,bits,print1(";");
for(n=-70,lim,
rest = lim - n;
for(c=0,rest,
a = bt(N,n,c); b = bittest(N,n,c);
if( a != b,
print("Error: N="N", n="n", c="c"; got "b", expect "a));
a = bt(-N,n,-c); b = bittest(-N,n,-c);
if( a != b,
print("Error: N="-N", n="n", c="-c"; got "b", expect "a))
));
N=floor(N/2);
lim = lim - 1) \\ shift the pattern down
}
==================================================================
--- ./src/basemath/arith1.c.orig Sat Oct 6 15:36:23 2001
+++ ./src/basemath/arith1.c Sat Oct 6 16:08:51 2001
@@ -115,6 +115,30 @@ garith_proto2gs(GEN f(GEN,long), GEN x,
}
GEN
+garith_proto3ggs(GEN f(GEN,GEN,long), GEN x, GEN y, long z)
+{
+ long l,i,tx = typ(x);
+ GEN t;
+
+ if (is_matvec_t(tx))
+ {
+ l=lg(x); t=cgetg(l,tx);
+ for (i=1; i<l; i++) t[i]= (long) garith_proto3ggs(f,(GEN) x[i],y,z);
+ return t;
+ }
+ if (tx != t_INT) err(arither1);
+ tx = typ(y);
+ if (is_matvec_t(tx))
+ {
+ l=lg(y); t=cgetg(l,tx);
+ for (i=1; i<l; i++) t[i]= (long) garith_proto3ggs(f,x,(GEN) y[i],z);
+ return t;
+ }
+ if (tx != t_INT) err(arither1);
+ return f(x,y,z);
+}
+
+GEN
gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y)
{
int tx=typ(x);
--- ./src/basemath/arith2.c.orig Sat Oct 6 16:14:08 2001
+++ ./src/basemath/arith2.c Sun Oct 7 13:20:30 2001
@@ -25,6 +25,7 @@ extern GEN arith_proto(long f(GEN), GEN
extern GEN arith_proto2(long f(GEN,GEN), GEN x, GEN n);
extern GEN garith_proto(GEN f(GEN), GEN x, int do_error);
extern GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y);
+extern GEN garith_proto3ggs(GEN f(GEN,GEN,long), GEN x, GEN y, long z);
#define sqru(i) (muluu((i),(i)))
@@ -1944,6 +1945,100 @@ bittest(GEN x, long n)
return n? 1: 0;
}
+GEN
+bittest_many(GEN x, GEN gn, long c)
+{
+ long lx = lgefint(x), l1, l2, b1, b2, two_adic = 0, l_res, skip;
+ ulong *p, *pnew;
+ long n = itos(gn), extra_words = 0, partial_bits;
+ GEN res;
+
+ if (c == 1) /* Shortcut */
+ return bittest(x,n) ? gun : gzero;
+ /* Negative n with n+c>0 are too hairy to implement... */
+ if (!signe(x) || c == 0)
+ return gzero;
+ if (c < 0) { /* Negative x means 2s complemenent */
+ c = -c;
+ if (signe(x) < 0)
+ two_adic = 1; /* treat x 2-adically... */
+ }
+ if (n < 0) {
+ long oa, na;
+
+ if (n + c <= 0)
+ return gzero;
+ oa = avma;
+ res = bittest_many(x, gzero, two_adic? -(n+c) : n+c);
+ na = avma;
+ res = shifti3(res,-n,0);
+ return gerepile(oa,na,res);
+ }
+ partial_bits = (c & (BITS_IN_LONG-1));
+ l2 = lx-1 - (n>>TWOPOTBITS_IN_LONG); /* Last=least-significant word */
+ l1 = lx-1 - ((n + c - 1)>>TWOPOTBITS_IN_LONG); /* First word */
+ b2 = (n & (BITS_IN_LONG-1)); /* Last bit, skip less-signifant */
+ b1 = ((n + c - 1) & (BITS_IN_LONG-1)); /* First bit, skip more-significant */
+ if (l2 < 2) { /* always: l1 <= l2 */
+ if (!two_adic)
+ return gzero;
+ /* x is non-zero, so it behaves as -1. */
+ return gbitneg(gzero,c);
+ }
+ if (l1 < 2) {
+ if (two_adic) /* If b1 < b2, bits are set by prepend in shift_r */
+ extra_words = 2 - l1 - (b1 < b2);
+ else
+ partial_bits = b2 ? BITS_IN_LONG - b2 : 0;
+ l1 = 2;
+ b1 = (BITS_IN_LONG-1); /* Include all bits in this word */
+ }
+ skip = (b1 < b2); /* Skip the first word of the shift */
+ l_res = l2 - l1 + 1 + 2 - skip + extra_words; /* A coarse estimate */
+ p = new_chunk(l_res) + 2 + extra_words;
+ shift_r(p - skip, x + l1, x + l2 + 1, 0, b2);
+ if (two_adic) { /* Check the low bits of x */
+ int i = lx, nonzero = 0;
+
+ /* Flip the bits */
+ pnew = p + l_res - 2 - extra_words;
+ while (--pnew >= p)
+ *pnew = MAXULONG - *pnew;
+ /* Fill the extra words */
+ while (extra_words--)
+ *--p = MAXULONG;
+ /* The result is one less than 2s-complement if the lower-bits
+ of x were all 0. */
+ while (--i > l2) {
+ if (x[i]) {
+ nonzero = 1;
+ break;
+ }
+ }
+ if (!nonzero && b2) /* Check the partial word too */
+ nonzero = x[l2] & ((1<<b2) - 1);
+ if (!nonzero) { /* Increment res. Do not underflow to before p */
+ pnew = p + l_res - 2;
+ while (--pnew >= p)
+ if (++*pnew)
+ break;
+ }
+ }
+ if (partial_bits)
+ *p &= (1<<partial_bits) - 1;
+ pnew = p;
+ while ((pnew < p + l_res - 2) && !*pnew) /* Skip 0 words */
+ pnew++;
+ avma = (long)(pnew - 2);
+ if (pnew >= p - 2 + l_res)
+ return gzero;
+ l_res -= (pnew - p);
+ res = pnew - 2;
+ res[1] = evalsigne(1)|evallgefint(l_res);
+ res[0] = evaltyp(t_INT)|evallg(l_res);
+ return res;
+}
+
static long
bittestg(GEN x, GEN n)
{
@@ -1954,6 +2049,12 @@ GEN
gbittest(GEN x, GEN n)
{
return arith_proto2(bittestg,x,n);
+}
+
+GEN
+gbittest3(GEN x, GEN n, long c)
+{
+ return garith_proto3ggs(bittest_many,x,n,c);
}
/***********************************************************************/
--- ./src/headers/paridecl.h.orig Sat Oct 6 14:41:21 2001
+++ ./src/headers/paridecl.h Sat Oct 6 16:03:47 2001
@@ -235,6 +235,7 @@ GEN gbitneg(GEN x, long n);
GEN gbitnegimply(GEN x, GEN y);
GEN gbitor(GEN x, GEN y);
GEN gbittest(GEN x, GEN n);
+GEN gbittest3(GEN x, GEN n, long c);
GEN gbitxor(GEN x, GEN y);
GEN gboundfact(GEN n, long lim);
GEN gissquarefree(GEN x);
@@ -1030,6 +1031,7 @@ int ratlift(GEN x, GEN m, GEN *a, GE
GEN resss(long x, long y);
double rtodbl(GEN x);
GEN shifti(GEN x, long n);
+void shift_r(ulong *target, ulong *source, ulong *source_end, ulong prepend, ulong sh);
GEN shifti3(GEN x, long n, long flag);
long smodsi(long x, GEN y);
GEN sqri(GEN x);
--- ./src/kernel/none/mp.c.orig Sat Oct 6 16:14:27 2001
+++ ./src/kernel/none/mp.c Sat Oct 6 18:32:20 2001
@@ -30,6 +30,7 @@ Foundation, Inc., 59 Temple Place - Suit
*/
/* z2 := z1[imin..imax].f shifted left sh bits (feeding f from the right) */
+/* These macros work only for sh != 0 !!! */
#define shift_left2(z2,z1,imin,imax,f, sh,m) {\
register ulong _l,_k = ((ulong)f)>>m;\
GEN t1 = z1 + imax, t2 = z2 + imax, T = z1 + imin;\
@@ -45,18 +46,21 @@ Foundation, Inc., 59 Temple Place - Suit
shift_left2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
}
-/* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */
-#define shift_right2(z2,z1,imin,imax,f, sh,m) {\
- register GEN t1 = z1 + imin, t2 = z2 + imin, T = z1 + imax;\
- register ulong _k,_l = *t1++;\
- *t2++ = (_l>>(ulong)sh) | ((ulong)f<<(ulong)m);\
- while (t1 < T) {\
- _k = _l<<(ulong)m; _l = *t1++;\
- *t2++ = (_l>>(ulong)sh) | _k;\
+#define shift_words_r(target,source,source_end,prepend, sh, sh_complement) {\
+ register ulong _k,_l = *source++;\
+ *target++ = (_l>>(ulong)sh) | ((ulong)prepend<<(ulong)sh_complement);\
+ while (source < source_end) {\
+ _k = _l<<(ulong)sh_complement; _l = *source++;\
+ *target++ = (_l>>(ulong)sh) | _k;\
}\
}
+#define shift_right2(z2,z1,imin,imax,f, sh,m) {\
+ register GEN s = (z1) + (imin), ta = (z2) + (imin), se = (z1) + (imax);\
+ shift_words_r(ta,s,se,(f),(sh),(m)); \
+}
+/* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */
#define shift_right(z2,z1,imin,imax,f, sh) {\
- register const ulong _m = BITS_IN_LONG - sh;\
+ register const ulong _m = BITS_IN_LONG - (sh);\
shift_right2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
}
@@ -311,31 +315,31 @@ shifti3(GEN x, long n, long flag)
}
/* With FLAG: round up instead of rounding down */
if (flag) { /* Check low bits of x */
- i = lx;
- flag = 0;
- while (--i >= lyorig) {
- if (x[i]) {
- flag = 1; /* Need to increment y by 1 */
- break;
+ i = lx;
+ flag = 0;
+ while (--i >= lyorig) {
+ if (x[i]) {
+ flag = 1; /* Need to increment y by 1 */
+ break;
+ }
}
- }
- if (!flag && m)
- flag = x[lyorig - 1] & ((1<<m) - 1);
- }
- if (flag) { /* Increment i */
- i = ly;
- while (1) {
- if (--i < 2) { /* Need to extend y on the left */
- if (avma <= bot)
- err(errpile);
- avma = (long)(--y); ly++;
- y[2] = 1;
- break;
+ if (!flag && m) /* Check the partial word too */
+ flag = x[lyorig - 1] & ((1<<m) - 1);
+ if (flag) { /* Increment i */
+ i = ly;
+ while (1) {
+ if (--i < 2) { /* Extend y on the left */
+ if (avma <= bot)
+ err(errpile);
+ avma = (long)(--y); ly++;
+ y[2] = 1;
+ break;
+ }
+ if (++y[i])
+ break;
+ /* On the next cycle: increment the next longword... */
+ }
}
- if (++y[i])
- break;
- /* Now we need to propagate the bit into the next longword... */
- }
}
}
y[1] = evalsigne(s)|evallgefint(ly);
@@ -4279,4 +4283,29 @@ mpsqrtl(GEN a)
}
while (x < y);
return y;
+}
+
+/* target should point to a buffer of source_end - source + 1 ulongs.
+
+ fills this buffer by bits of ulongs in source..source_end-1 shifted
+ right sh units; the "most significant" sh bits of the result are
+ set to be the least significant sh bits of prepend.
+
+ The ordering of bits in this bitmap is the same as for t_INT.
+
+ sh should not exceed BITS_IN_LONG.
+ */
+void
+shift_r(ulong *target, ulong *source, ulong *source_end, ulong prepend, ulong sh)
+{
+ if (sh) { /* shift_words_r() works */
+ register ulong sh_complement = BITS_IN_LONG - sh;
+
+ shift_words_r(target, source, source_end, prepend, sh, sh_complement);
+ } else {
+ int i;
+
+ for (i=0; i < source_end - source; i++)
+ target[i] = source[i];
+ }
}
--- ./src/language/init.c.orig Sat Oct 6 16:12:14 2001
+++ ./src/language/init.c Sat Oct 6 16:51:34 2001
@@ -1756,7 +1756,7 @@ entree functions_basic[]={
{"bitneg",99,(void*)gbitneg,2,"GD-1,L,"},
{"bitnegimply",2,(void*)gbitnegimply,2,"GG"},
{"bitor",2,(void*)gbitor,2,"GG"},
-{"bittest",2,(void*)gbittest,2,"GG"},
+{"bittest",99,(void*)gbittest3,2,"GGD1,L,"},
{"bitxor",2,(void*)gbitxor,2,"GG"},
{"bnfcertify",10,(void*)certifybuchall,6,"lG"},
{"bnfclassunit",99,(void*)bnfclassunit0,6,"GD0,L,DGp"},
--- ./src/language/helpmsg.c.orig Sat Oct 6 16:46:44 2001
+++ ./src/language/helpmsg.c Sun Oct 7 13:16:07 2001
@@ -55,7 +55,7 @@ char *helpmessages_basic[]={
"bitneg(x,{n=-1}): bitwise negation of an integers x truncated to n bits. n=-1 means represent infinite sequences of bit 1 as negative numbers. Negative numbers behave as if modulo big power of 2",
"bitnegimpy(x,y): bitwise \"negated imply\" of two integers x and y, in other words, x BITAND BITNEG(y). Negative numbers behave as if modulo big power of 2",
"bitor(x,y): bitwise \"or\" of two integers x and y. Negative numbers behave as if modulo big power of 2",
- "bittest(x,n): gives bit number n (coefficient of 2^n) of the integer x",
+ "bittest(x,n,{c=1}): extracts |c| bits starting from number n (coefficient of 2^n) of the integer x, returning the bits as an integer bitmap; negative c means: treat negative values of x as if modulo big power of 2; bits at negative offsets are zeros",
"bitxor(x,y): bitwise \"exclusive or\" of two integers x and y. Negative numbers behave as if modulo big power of 2",
"bnfcertify(bnf): certify the correctness (i.e. remove the GRH) of the bnf data output by bnfclassunit or bnfinit",
"bnfclassunit(P,{flag=0},{tech=[]}): compute the class group, regulator of the number field defined by the polynomial P, and also the fundamental units if they are not too large. flag and tech are both optional. flag can be any of 0: default, 1: insist on having fundamental units, 2: do not compute units. See manual for details about tech. P may also be a non-zero integer, and is then considered as the discriminant of a quadratic order",