summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorH. Peter Anvin <hpa@zytor.com>2007-10-12 23:30:17 (GMT)
committerH. Peter Anvin <hpa@zytor.com>2007-10-12 23:30:17 (GMT)
commitaf492b3c6ece9c9c1fa6b9b3bd2d61904d25cd26 (patch)
tree66e4775e8e4f22ecbc550b361e1fbb68b46ba09f
parent1981484ab028768afa3b2b0b746a00e6ed739891 (diff)
downloadpbn-af492b3c6ece9c9c1fa6b9b3bd2d61904d25cd26.zip
pbn-af492b3c6ece9c9c1fa6b9b3bd2d61904d25cd26.tar.gz
pbn-af492b3c6ece9c9c1fa6b9b3bd2d61904d25cd26.tar.bz2
pbn-af492b3c6ece9c9c1fa6b9b3bd2d61904d25cd26.tar.xz
Add working (but slow) division algorithm; fix shifts
- Add working (but slow) pbn_div(); - Make pbn_adjust_bits() a part of the official API - Fix bug in the simple-case shift code
-rw-r--r--Makefile2
-rw-r--r--pbn.h7
-rw-r--r--pbn_add.c2
-rw-r--r--pbn_and.c6
-rw-r--r--pbn_bit.c4
-rw-r--r--pbn_div.c85
-rw-r--r--pbn_dump.c7
-rw-r--r--pbn_init.c4
-rw-r--r--pbn_mul.c9
-rw-r--r--pbn_or.c2
-rw-r--r--pbn_shift.c8
-rw-r--r--pbn_xor.c6
-rw-r--r--pbnint.h2
-rw-r--r--test.c84
14 files changed, 176 insertions, 52 deletions
diff --git a/Makefile b/Makefile
index 5aa0d90..b97cd90 100644
--- a/Makefile
+++ b/Makefile
@@ -4,7 +4,7 @@ LDFLAGS =
LIBOBJ = pbn_add.o pbn_cmp.o pbn_dump.o pbn_init.o pbn_mul.o pbn_shift.o \
pbn_and.o pbn_or.o pbn_xor.o pbn_bit.o \
- pbn_abs.o
+ pbn_abs.o pbn_div.o
TESTS = test
diff --git a/pbn.h b/pbn.h
index aab69b6..42b0959 100644
--- a/pbn.h
+++ b/pbn.h
@@ -63,6 +63,7 @@ struct pbn *pbn_addsub(struct pbn *, struct pbn *, int issub);
int pbn_cmp(const struct pbn *, const struct pbn *);
int pbn_abscmp(const struct pbn *, const struct pbn *);
struct pbn *pbn_mul(struct pbn *, struct pbn *);
+int pbn_div(struct pbn **, struct pbn **, struct pbn *, struct pbn *);
struct pbn *pbn_abs(struct pbn *);
struct pbn *pbn_shr(struct pbn *, int);
@@ -78,13 +79,15 @@ int pbn_bit(const struct pbn *, int);
struct pbn *pbn_int(pbn_satom_t);
-void pbn_dump(FILE *, struct pbn *);
+/* Primarily for debugging; does NOT consume a reference */
+void pbn_dump(FILE *, const struct pbn *);
+
struct pbn *pbn_new(int);
struct pbn *pbn_dup(const struct pbn *);
struct pbn *pbn_dupx(const struct pbn *, int);
struct pbn *pbn_cow(struct pbn *, int);
/* Internal functions */
-struct pbn *_pbn_adjust_bits(struct pbn *pbn);
+struct pbn *pbn_adjust_bits(struct pbn *pbn);
#endif /* PBN_H */
diff --git a/pbn_add.c b/pbn_add.c
index 931a27b..1f8c1f8 100644
--- a/pbn_add.c
+++ b/pbn_add.c
@@ -120,5 +120,5 @@ struct pbn *pbn_addsub(struct pbn *s1, struct pbn *s2, int issub)
d->minus = dminus;
}
- return _pbn_adjust_bits(d);
+ return pbn_adjust_bits(d);
}
diff --git a/pbn_and.c b/pbn_and.c
index e0478e2..9e282e0 100644
--- a/pbn_and.c
+++ b/pbn_and.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -47,7 +47,5 @@ struct pbn *pbn_and(struct pbn *s1, struct pbn *s2)
d->num[i] &= s2->num[len];
pbn_free(s2);
- return _pbn_adjust_bits(d);
+ return pbn_adjust_bits(d);
}
-
-
diff --git a/pbn_bit.c b/pbn_bit.c
index 9049117..51b0677 100644
--- a/pbn_bit.c
+++ b/pbn_bit.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -44,7 +44,7 @@ struct pbn *pbn_clr_bit(struct pbn *pbn, int bit)
bit++;
if (bit == pbn->bits)
- return _pbn_adjust_bits(pbn); /* We just cleared the MSB */
+ return pbn_adjust_bits(pbn); /* We just cleared the MSB */
else
return pbn;
}
diff --git a/pbn_div.c b/pbn_div.c
index b9ee448..398c70e 100644
--- a/pbn_div.c
+++ b/pbn_div.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -20,23 +20,84 @@
Returns -1 on error (divide by zero) and frees the inputs.
- This is a "basecase" division and so is O(n^2) -- it assumes that
- the numbers we have to deal with aren't all that large. For really
- large numbers, we need a subquadratic algorithm like
- divide-and-conquer or Newton-Rhapson division. */
+ This is doing it the slow way, which is bitwise long division.
+ If this turns out to be a bottleneck, then it should be replaced
+ by smarter algorithms, but this algorithm should at least be
+ correct, so it should be usable for now. */
#include "pbn.h"
-int pbn_div(struct pbn **q, struct pbn **r, struct pbn *n, struct pbn *d)
+int pbn_div(struct pbn **qp, struct pbn **rp, struct pbn *n, struct pbn *d)
{
- if (d->bits == 0) {
- pbn_free(n);
+ struct pbn *q, *r;
+ int rv = d->bits == 0 ? -1 : 0;
+ int cm;
+ int shift, bits, len, i;
+
+ if (rv || (cm = pbn_abscmp(n, d)) < 0) {
+ /* We either had division by zero, or |n| < |d| */
+ int rv = d->bits == 0;
+
+ if (rp) {
+ if (d->minus && n->bits) {
+ n = pbn_cow(n, n->len);
+ n->minus ^= 1;
+ }
+ *rp = n;
+ } else {
+ pbn_free(n);
+ }
pbn_free(d);
- return -1; /* Division by zero */
+
+ if (qp)
+ *qp = pbn_int(0);
+
+ return rv; /* Division by zero? */
+ }
+
+ if (cm == 0) {
+ /* Not really a common case, but we got it anyway */
+ int quot = (n->minus ^ d->minus) ? -1 : 1;
+ pbn_free(d);
+ pbn_free(n);
+
+ if (rp)
+ *rp = pbn_int(0);
+ if (qp)
+ *qp = pbn_int(quot);
+
+ return 0;
}
-
+ shift = n->bits - d->bits;
+ bits = shift+1; /* Maximum possible number of quotient bits */
+ len = (bits+PBN_ATOM_BITS-1)/PBN_ATOM_BITS;
-
-}
+ r = pbn_cow(n, n->len);
+ q = pbn_new(len);
+ d = pbn_shl(d, shift);
+ q->minus = r->minus ^ d->minus;
+ r->minus = d->minus; /* Makes our lives easier for now... */
+
+ for (i = shift; i >= 0; i--) {
+ if (pbn_abscmp(r, d) >= 0) {
+ r = pbn_sub(r, d);
+ q = pbn_set_bit(q, i);
+ }
+ d = pbn_shr(d, 1);
+ }
+ r->minus = q->minus;
+
+ if (rp)
+ *rp = pbn_adjust_bits(r);
+ else
+ pbn_free(r);
+
+ if (qp)
+ *qp = pbn_adjust_bits(q);
+ else
+ pbn_free(q);
+
+ return 0;
+}
diff --git a/pbn_dump.c b/pbn_dump.c
index 5dbcde1..da2eb7d 100644
--- a/pbn_dump.c
+++ b/pbn_dump.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -19,7 +19,7 @@
#include <stdio.h>
#include "pbnint.h"
-void pbn_dump(FILE *f, struct pbn *pbn)
+void pbn_dump(FILE *f, const struct pbn *pbn)
{
int i;
@@ -30,6 +30,5 @@ void pbn_dump(FILE *f, struct pbn *pbn)
if (i)
putc('_', f);
}
-
- pbn_free(pbn);
+ fprintf(f, " (%d)", pbn->bits);
}
diff --git a/pbn_init.c b/pbn_init.c
index c7c053c..5bdebaf 100644
--- a/pbn_init.c
+++ b/pbn_init.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -144,7 +144,7 @@ struct pbn *pbn_cow(struct pbn *src, int len)
/* Adjust the bits field of a specific pbn, and clear the minus flag
if the value is actually zero. Returns its argument, to make
tailcalling possible. */
-struct pbn *_pbn_adjust_bits(struct pbn *pbn)
+struct pbn *pbn_adjust_bits(struct pbn *pbn)
{
int i;
pbn_atom_t v;
diff --git a/pbn_mul.c b/pbn_mul.c
index f7f35f3..2ec50d3 100644
--- a/pbn_mul.c
+++ b/pbn_mul.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -26,7 +26,7 @@ struct pbn *pbn_mul(struct pbn *s1, struct pbn *s2)
int i, j;
pbn_atom_t x1;
pbn_2atom_t c;
-
+
if (s1->bits == 0) {
pbn_free(s2);
return s1; /* s1 == 0 */
@@ -37,7 +37,7 @@ struct pbn *pbn_mul(struct pbn *s1, struct pbn *s2)
bits = s1->bits + s2->bits;
len = (bits+PBN_ATOM_BITS-1)/PBN_ATOM_BITS;
-
+
d = pbn_new(len);
d->bits = bits;
@@ -57,6 +57,5 @@ struct pbn *pbn_mul(struct pbn *s1, struct pbn *s2)
pbn_free(s1);
pbn_free(s2);
- return _pbn_adjust_bits(d);
+ return pbn_adjust_bits(d);
}
-
diff --git a/pbn_or.c b/pbn_or.c
index 294f401..387dbc8 100644
--- a/pbn_or.c
+++ b/pbn_or.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
diff --git a/pbn_shift.c b/pbn_shift.c
index 21d65cf..3a631bd 100644
--- a/pbn_shift.c
+++ b/pbn_shift.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -42,7 +42,7 @@ struct pbn *pbn_shl(struct pbn *src, int shift)
if (sb == 0) {
/* Particularly simple case... */
- memmove(&dst->num[sw], &dst->num[0], len-sw);
+ memmove(&dst->num[sw], &dst->num[0], (len-sw)*sizeof(pbn_atom_t));
} else {
b = dst->num[len-sw-1];
for (i = len-sw-1; i > 0; i--) {
@@ -85,10 +85,10 @@ struct pbn *pbn_shr(struct pbn *src, int shift)
sw = shift / PBN_ATOM_BITS;
sb = shift % PBN_ATOM_BITS;
-
+
if (sb == 0) {
/* Particularly simple case... */
- memmove(&dst->num[0], &dst->num[sw], len-sw);
+ memmove(&dst->num[0], &dst->num[sw], (len-sw)*sizeof(pbn_atom_t));
} else {
c = dst->num[sw];
for (i = sw; i < len-1; i++) {
diff --git a/pbn_xor.c b/pbn_xor.c
index a104e63..f5f0ff3 100644
--- a/pbn_xor.c
+++ b/pbn_xor.c
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
@@ -44,7 +44,5 @@ struct pbn *pbn_xor(struct pbn *s1, struct pbn *s2)
d->num[i] ^= s2->num[len];
pbn_free(s2);
- return _pbn_adjust_bits(d);
+ return pbn_adjust_bits(d);
}
-
-
diff --git a/pbnint.h b/pbnint.h
index 8915a1d..55d63f4 100644
--- a/pbnint.h
+++ b/pbnint.h
@@ -1,5 +1,5 @@
/* ----------------------------------------------------------------------- *
- *
+ *
* Copyright 2007 H. Peter Anvin - All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
diff --git a/test.c b/test.c
index b899741..17bd788 100644
--- a/test.c
+++ b/test.c
@@ -3,6 +3,62 @@
#include <stdlib.h>
#include "pbn.h"
+static void divide_test(void)
+{
+ struct pbn *n, *d, *r, *q;
+ int rv;
+
+ printf("\n*** Divide test ***\n\n");
+
+ /* 2^128/2^64 */
+
+ n = pbn_new(5);
+ n->num[4] = 1; /* 2^128 */
+ pbn_adjust_bits(n);
+
+ d = pbn_new(3);
+ d->num[2] = 1; /* 2^64 */
+ pbn_adjust_bits(d);
+
+ rv = pbn_div(&q, &r, n, d);
+
+ pbn_dump(stdout, q);
+ printf(" : ");
+ pbn_dump(stdout, r);
+ putchar('\n');
+
+ /* 2^128/3 */
+
+ n = pbn_new(5);
+ n->num[4] = 1; /* 2^128 */
+ pbn_adjust_bits(n);
+
+ d = pbn_int(3);
+ rv = pbn_div(&q, &r, n, d);
+
+ pbn_dump(stdout, q);
+ printf(" : ");
+ pbn_dump(stdout, r);
+ putchar('\n');
+
+ /* 2^128/(3*2^64) */
+
+ n = pbn_new(5);
+ n->num[4] = 1; /* 2^128 */
+ pbn_adjust_bits(n);
+
+ d = pbn_new(3);
+ d->num[2] = 3; /* 3*2^64 */
+ pbn_adjust_bits(d);
+
+ rv = pbn_div(&q, &r, n, d);
+
+ pbn_dump(stdout, q);
+ printf(" : ");
+ pbn_dump(stdout, r);
+ putchar('\n');
+}
+
int main(int argc, char *argv[])
{
uint64_t i, j;
@@ -13,10 +69,12 @@ int main(int argc, char *argv[])
m = pbn_int(5);
i = 5;
+ printf("\n*** Multiply test ***\n\n");
+
n = nx;
while (n--) {
printf("%016"PRIx64" : ", i);
- pbn_dump(stdout, pbn_ref(m));
+ pbn_dump(stdout, m);
putchar('\n');
i *= 5;
@@ -24,26 +82,30 @@ int main(int argc, char *argv[])
}
pbn_free(m);
+ printf("\n*** Add test ***\n\n");
+
m = pbn_int(5);
i = 5;
n = nx;
while (n--) {
printf("%016"PRIx64" : ", i);
- pbn_dump(stdout, pbn_ref(m));
+ pbn_dump(stdout, m);
putchar('\n');
i += i;
m = pbn_add(m, pbn_ref(m));
}
+ printf("\n*** Subtract test ***\n\n");
+
n = nx;
j = 20;
mj = pbn_int(20);
while (n--) {
printf("%016"PRIx64" : ", i);
- pbn_dump(stdout, pbn_ref(m));
+ pbn_dump(stdout, m);
printf(" : ");
- pbn_dump(stdout, pbn_ref(mj));
+ pbn_dump(stdout, mj);
putchar('\n');
i -= j;
@@ -51,16 +113,18 @@ int main(int argc, char *argv[])
m = pbn_sub(m, pbn_ref(mj));
mj = pbn_add(mj, pbn_ref(mj));
}
-
+
+ printf("\n*** Shift test ***\n\n");
+
printf("[ 0] ");
- pbn_dump(stdout, pbn_ref(m));
+ pbn_dump(stdout, m);
printf(" (%d)", m->bits);
putchar('\n');
for (n = 1; n <= 8; n++) {
m = pbn_shl(m, n);
printf("[%2d] ", n);
- pbn_dump(stdout, pbn_ref(m));
+ pbn_dump(stdout, m);
printf(" (%d)", m->bits);
putchar('\n');
}
@@ -68,11 +132,13 @@ int main(int argc, char *argv[])
for (n = 1; n <= 16; n++) {
m = pbn_shr(m, n);
printf("[%2d] ", n);
- pbn_dump(stdout, pbn_ref(m));
+ pbn_dump(stdout, m);
printf(" (%d)", m->bits);
putchar('\n');
}
+ pbn_free(m);
+
+ divide_test();
return 0;
}
-