| Ilya Zakharevich on Sun, 10 Oct 1999 16:42:11 -0400 (EDT) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| [PATCH 2.0.15] Bitops in PARI |
[This patch was created on 2.0.15, but bulk of testsuite was run on
2.0.17 as well.]
Finally, a lot of work and more than 80000 test cases behind, I
managed to squeeze the full suite of bitops into PARI. Note that
bitneg(), bitor(), bitxor() and bitand() are required to make PARI into a
complete replacement of Perl's numeric subsystem.
Here is why bitnegimply() may be wished too: there are 16 different
binary ops with 2 0-or-1 arguments and 0-or-1 output (there are 4
possible inputs, thus 2^4 possible maps). Out of them 6 do not depend
on one of the arguments. Out of remaining 10 consider those 5 which
send (0,0) to 0.
First of all, the other 5 can be described as 1-op(x,y) with op(x,y)
in the above list.
Second, if we extend op() to act on sequences of 0-or-1 componentwise,
then these 5 ops sends sequences with finite support to sequences with
finite support. Thus it makes sense to consider these 5 ops as
bitwise operations on integers. And these 5 ops become
bitor(x,y), bitxor(x,y), bitand(x,y), bitnegimply(x,y), bitnegimply(y,x)
As I explained already, the other 5 may be obtained as, say,
bitneg(bitxor(x,y))
Here enters bitneg: it sends a sequence of all 0s into a sequence of
all 1s. There may be different ways to encode such a sequence in
PARI's type system. I chose the simplest one: encode them as negative
numbers. This has many different names 2-complement, 2-adic integers, etc.
Enjoy,
Ilya
--- ./src/language/helpmsg.c-ini Fri May 21 13:10:30 1999
+++ ./src/language/helpmsg.c Sat Oct 9 20:47:04 1999
@@ -38,7 +38,12 @@ char *helpmessages_basic[]={
"bigomega(x): number of repeated prime divisors of x",
"binary(x): gives the vector formed by the binary digits of x (x integer)",
"binomial(x,y): binomial coefficient x*(x-1)...*(x-y+1)/y! defined for y in Z and any x",
+ "bitand(x,y): bitwise \"and\" of two integers x and y. Negative numbers behave as if modulo big power of 2",
+ "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",
+ "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",
"bnfclgp(P,{tech=[]}): compute the class group of the number field defined by the polynomial P. If P is a non-zero integer, it is interpreted as a quadratic discriminant. See manual for details about tech",
--- ./src/language/init.c-ini Fri Jun 4 03:42:28 1999
+++ ./src/language/init.c Sat Oct 9 20:46:20 1999
@@ -1291,7 +1291,12 @@ entree functions_basic[]={
{"bigomega",18,(void*)gbigomega,4,"G"},
{"binary",18,(void*)binaire,2,"G"},
{"binomial",21,(void*)binome,4,"GL"},
+{"bitand",2,(void*)gbitand,2,"GG"},
+{"bitneg",2,(void*)gbitneg,2,"GD-1,L,"},
+{"bitnegimply",2,(void*)gbitnegimply,2,"GG"},
+{"bitor",2,(void*)gbitor,2,"GG"},
{"bittest",2,(void*)gbittest,2,"GG"},
+{"bitxor",2,(void*)gbitxor,2,"GG"},
{"bnfcertify",10,(void*)certifybuchall,6,"lG"},
{"bnfclassunit",99,(void*)bnfclassunit0,6,"GD0,L,DGp"},
{"bnfclgp",99,(void*)classgrouponly,6,"GDGp"},
--- ./src/basemath/gen1.c-ini Thu May 20 12:04:38 1999
+++ ./src/basemath/gen1.c Sun Oct 10 16:10:42 1999
@@ -1977,3 +1977,422 @@ gdiv(GEN x, GEN y)
err(typeer,"division");
return NULL; /* not reached */
}
+
+/********************************************************************/
+/** **/
+/** BITMAP OPS **/
+/** **/
+/********************************************************************/
+
+/* Normalize a non-negative integer. */
+static void
+inormalize(GEN x, long known_zero_words)
+{
+ int xl = lgefint(x);
+ int i, j;
+
+ /* Normalize */
+ i = 2 + known_zero_words;
+ while (i < xl) {
+ if (x[i])
+ break;
+ i++;
+ }
+ j = 2;
+ while (i < xl)
+ x[j++] = x[i++];
+ xl -= i - j;
+ setlgefint(x, xl);
+ if (xl == 2)
+ setsigne(x,0);
+}
+
+/* Truncate a non-negative integer to a number of bits. */
+static void
+ibittrunc(GEN x, long bits, long normalized)
+{
+ int xl = lgefint(x);
+ int len_out = ((bits + BITS_IN_LONG - 1) >> TWOPOTBITS_IN_LONG) + 2;
+ int known_zero_words, i = 2 + xl - len_out;
+
+ if (xl < len_out && normalized)
+ return;
+ /* Check whether mask is trivial */
+ if (!(bits & (BITS_IN_LONG - 1))) {
+ if (xl == len_out && normalized)
+ return;
+ } else if (len_out <= xl) {
+ /* Non-trival mask is given by a formula, if x is not
+ normalized, this works even in the exceptional case */
+ x[i] = x[i] & ((1 << (bits & (BITS_IN_LONG - 1))) - 1);
+ if (x[i] && xl == len_out)
+ return;
+ }
+ /* Normalize */
+ if (xl <= len_out) /* Not normalized */
+ known_zero_words = 0;
+ else
+ known_zero_words = xl - len_out;
+ inormalize(x, known_zero_words);
+}
+
+/* Increment/decrement absolute value of non-zero integer in place.
+ It is assumed that the increment will not increase the length.
+ Decrement may result in a non-normalized value.
+ Returns 1 if the increment overflows (thus the result is 0). */
+static int
+incdec(GEN x, long incdec)
+{
+ long *xf = x + 2, *xl;
+ long len = lgefint(x);
+ const ulong uzero = 0;
+
+ xl = x + len;
+ if (incdec == 1) {
+ while (--xl >= xf) {
+ if (*xl != ~uzero) {
+ (*xl)++;
+ return 0;
+ }
+ *xl = 0;
+ }
+ return 1;
+ } else {
+ while (--xl >= xf) {
+ if (*xl != 0) {
+ (*xl)--;
+ return 0;
+ }
+ *xl = (long)~uzero;
+ }
+ return 0;
+ }
+}
+
+GEN
+gbitneg(GEN x, long bits)
+{
+ long xl, len_out, i, j;
+ const ulong uzero = 0;
+
+ if (typ(x) != t_INT)
+ err(typeer, "bitwise negation");
+ if (bits < -1)
+ err(talker, "negative exponent in bitwise negation");
+ if (bits == -1)
+ return gsub(gneg(gun),x);
+ if (bits == 0)
+ return gzero;
+ if (signe(x) == -1) { /* Consider as if mod big power of 2 */
+ x = gcopy(x);
+ setsigne(x, 1);
+ incdec(x, -1);
+ /* Now truncate this! */
+ ibittrunc(x, bits, x[2]);
+ return x;
+ }
+ xl = lgefint(x);
+ len_out = ((bits + BITS_IN_LONG - 1) >> TWOPOTBITS_IN_LONG) + 2;
+ if (len_out > xl) { /* Need to grow */
+ GEN out = cgeti(len_out);
+ int j = 2;
+
+ if (!(bits & (BITS_IN_LONG - 1)))
+ out[2] = ~uzero;
+ else
+ out[2] = (1 << (bits & (BITS_IN_LONG - 1))) - 1;
+ for (i = 3; i < len_out - xl + 2; i++)
+ out[i] = ~uzero;
+ while (i < len_out)
+ out[i++] = ~x[j++];
+ setlgefint(out, len_out);
+ setsigne(out,1);
+ return out;
+ }
+ x = gcopy(x);
+ for (i = 2; i < xl; i++)
+ x[i] = ~x[i];
+ ibittrunc(x, bits, x[2]);
+ return x;
+}
+
+/* bitwise 'and' of two positive integers (any integers, but we ignore sign).
+ * Inputs are not necessary normalized. */
+static GEN
+ibitand(GEN x, GEN y)
+{
+ long lx, ly, lout;
+ long *xp, *yp, *outp, *xlim;
+ GEN out;
+
+ lx=lgefint(x); ly=lgefint(y);
+ if (lx > ly)
+ lout = ly;
+ else
+ lout = lx;
+ xlim = x + lx;
+ xp = xlim + 2 - lout;
+ yp = y + 2 + ly - lout;
+ out = cgeti(lout);
+ outp = out + 2;
+ while (xp < xlim)
+ *outp++ = (*xp++) & (*yp++);
+ setsigne(out,1);
+ setlgefint(out,lout);
+ if (lout == 2)
+ setsigne(out,0);
+ else if ( !out[2] )
+ inormalize(out, 1);
+ return out;
+}
+
+#define swaplen(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
+
+/* bitwise 'or' of two positive integers (any integers, but we ignore sign).
+ * Inputs are not necessary normalized. */
+static GEN
+ibitor(GEN x, GEN y, long exclusive)
+{
+ long lx, ly, lout;
+ long *xp, *yp, *outp, *xlim, *xprep;
+ GEN out;
+
+ lx=lgefint(x); ly=lgefint(y);
+ if (lx < ly)
+ swaplen(x,y,lx,ly);
+ lout = lx;
+ xlim = x + lx;
+ xp = xlim + 2 - ly;
+ yp = y + 2;
+ out = cgeti(lout);
+ outp = out + 2;
+ if (lx > ly) {
+ xprep = x + 2;
+ while (xprep < xp)
+ *outp++ = *xprep++;
+ }
+ if (exclusive) {
+ while (xp < xlim)
+ *outp++ = (*xp++) ^ (*yp++);
+ } else {
+ while (xp < xlim)
+ *outp++ = (*xp++) | (*yp++);
+ }
+ setsigne(out,1);
+ setlgefint(out,lout);
+ if (lout == 2)
+ setsigne(out,0);
+ else if ( !out[2] )
+ inormalize(out, 1);
+ return out;
+}
+
+/* bitwise negated 'implies' of two positive integers (any integers, but we
+ * ignore sign). "Neg-Implies" is x & ~y unless "negated".
+ * Inputs are not necessary normalized. */
+static GEN
+ibitnegimply(GEN x, GEN y)
+{
+ long lx, ly, lout, inverted = 0;
+ long *xp, *yp, *outp, *xlim, *xprep;
+ GEN out;
+
+ lx=lgefint(x); ly=lgefint(y);
+ if (lx < ly) {
+ inverted = 1;
+ swaplen(x,y,lx,ly);
+ }
+ /* x is longer than y */
+ lout = lx;
+ xlim = x + lx;
+ xp = xlim + 2 - ly;
+ yp = y + 2;
+ out = cgeti(lout);
+ outp = out + 2;
+ if (lx > ly) {
+ xprep = x + 2;
+ if (!inverted) { /* x & ~y */
+ while (xprep < xp)
+ *outp++ = *xprep++;
+ } else { /* ~x & y */
+ while (xprep++ < xp)
+ *outp++ = 0;
+ }
+ }
+ if (inverted) { /* ~x & y */
+ while (xp < xlim)
+ *outp++ = ~(*xp++) & (*yp++);
+ } else {
+ while (xp < xlim)
+ *outp++ = (*xp++) & ~(*yp++);
+ }
+ setsigne(out,1);
+ setlgefint(out,lout);
+ if (lout == 2)
+ setsigne(out,0);
+ else if ( !out[2] )
+ inormalize(out, 1);
+ return out;
+}
+
+static GEN
+inegate_inplace(GEN z, long ltop)
+{
+ long lbot, o;
+
+ /* Negate z */
+ o = incdec(z, 1); /* Can overflow z... */
+ setsigne(z, -1);
+ if (!o)
+ return z;
+ else if (lgefint(z) == 2)
+ setsigne(z, 0);
+ incdec(z,-1); /* Restore a normalized value */
+ lbot = avma;
+ z = gsub(z,gun);
+ return gerepile(ltop,lbot,z);
+}
+
+GEN
+gbitor(GEN x, GEN y)
+{
+ long sx,sy;
+ long ltop, o;
+ GEN z;
+
+ if (typ(x) != t_INT || typ(y) != t_INT)
+ err(typeer, "bitwise or");
+ sx=signe(x); if (!sx) return icopy(y);
+ sy=signe(y); if (!sy) return icopy(x);
+ if (sx == 1) {
+ GEN t;
+
+ if (sy == 1)
+ return ibitor(x,y,0);
+ goto posneg;
+ } else if (sy == -1) {
+ ltop = avma;
+ incdec(x, -1); incdec(y, -1);
+ z = ibitand(x,y);
+ incdec(x, 1); incdec(y, 1); /* Restore x and y... */
+ } else {
+ z = x; x = y; y = z;
+ /* x is positive, y is negative. The result will be negative. */
+ posneg:
+ ltop = avma;
+ incdec(y, -1);
+ z = ibitnegimply(y,x); /* ~x & y */
+ incdec(y, 1);
+ }
+ return inegate_inplace(z, ltop);
+}
+
+GEN
+gbitand(GEN x, GEN y)
+{
+ long sx,sy;
+ long ltop, o;
+ GEN z;
+
+ if (typ(x) != t_INT || typ(y) != t_INT)
+ err(typeer, "bitwise and");
+ sx=signe(x); if (!sx) return gzero;
+ sy=signe(y); if (!sy) return gzero;
+ if (sx == 1) {
+ GEN t;
+
+ if (sy == 1)
+ return ibitand(x,y);
+ goto posneg;
+ } else if (sy == -1) {
+ ltop = avma;
+ incdec(x, -1); incdec(y, -1);
+ z = ibitor(x,y,0);
+ incdec(x, 1); incdec(y, 1); /* Restore x and y... */
+ return inegate_inplace(z, ltop);
+ } else {
+ z = x; x = y; y = z;
+ /* x is positive, y is negative. The result will be positive. */
+ posneg:
+ ltop = avma;
+ incdec(y, -1);
+ /* x & ~y */
+ z = ibitnegimply(x,y); /* x & ~y */
+ incdec(y, 1);
+ return z;
+ }
+}
+
+GEN
+gbitxor(GEN x, GEN y)
+{
+ long sx,sy;
+ long ltop, o;
+ GEN z;
+
+ if (typ(x) != t_INT || typ(y) != t_INT)
+ err(typeer, "bitwise xor");
+ sx=signe(x); if (!sx) return icopy(y);
+ sy=signe(y); if (!sy) return icopy(x);
+ if (sx == 1) {
+ GEN t;
+
+ if (sy == 1)
+ return ibitor(x,y,1);
+ goto posneg;
+ } else if (sy == -1) {
+ incdec(x, -1); incdec(y, -1);
+ z = ibitor(x,y,1);
+ incdec(x, 1); incdec(y, 1); /* Restore x and y... */
+ return z;
+ } else {
+ z = x; x = y; y = z;
+ /* x is positive, y is negative. The result will be negative. */
+ posneg:
+ ltop = avma;
+ incdec(y, -1);
+ /* ~(x ^ ~y) == x ^ y */
+ z = ibitor(x,y,1);
+ incdec(y, 1);
+ }
+ return inegate_inplace(z, ltop);
+}
+
+GEN
+gbitnegimply(GEN x, GEN y) /* x & ~y */
+{
+ long sx,sy;
+ long ltop, o;
+ GEN z;
+
+ if (typ(x) != t_INT || typ(y) != t_INT)
+ err(typeer, "bitwise negated imply");
+ sx=signe(x); if (!sx) return gzero;
+ sy=signe(y); if (!sy) return icopy(x);
+ if (sx == 1) {
+ GEN t;
+
+ if (sy == 1)
+ return ibitnegimply(x,y);
+ /* x is positive, y is negative. The result will be positive. */
+ incdec(y, -1);
+ z = ibitand(x,y);
+ incdec(y, 1);
+ return z;
+ } else if (sy == -1) {
+ /* both negative. The result will be positive. */
+ incdec(x, -1); incdec(y, -1);
+ /* ((~x) & ~(~y)) == ~x & y */
+ z = ibitnegimply(y,x);
+ incdec(x, 1); incdec(y, 1); /* Restore x and y... */
+ return z;
+ } else {
+ /* x is negative, y is positive. The result will be negative. */
+ ltop = avma;
+ incdec(x, -1);
+ /* ~((~x) & ~y) == x | y */
+ z = ibitor(x,y,0);
+ incdec(x, 1);
+ }
+ return inegate_inplace(z, ltop);
+}
--- ./src/headers/paridecl.h-ini Thu May 20 12:05:10 1999
+++ ./src/headers/paridecl.h Sat Oct 9 20:49:48 1999
@@ -131,7 +131,12 @@ GEN classno(GEN x);
GEN classno2(GEN x);
GEN fibo(long n);
GEN fundunit(GEN x);
+GEN gbitand(GEN x, GEN y);
+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 gbitxor(GEN x, GEN y);
GEN gboundcf(GEN x, long k);
GEN gcarrecomplet(GEN x, GEN *pt);
GEN gcarreparfait(GEN x);