summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorH. Peter Anvin <hpa@zytor.com>2007-10-12 01:05:50 (GMT)
committerH. Peter Anvin <hpa@zytor.com>2007-10-12 01:05:50 (GMT)
commit3a01eee61c1ada7fcc83ae0c0cd5f8f25b8d5001 (patch)
treef86b9b1f51901ff89ca6287a30e1b7c116789306
downloadpbn-3a01eee61c1ada7fcc83ae0c0cd5f8f25b8d5001.zip
pbn-3a01eee61c1ada7fcc83ae0c0cd5f8f25b8d5001.tar.gz
pbn-3a01eee61c1ada7fcc83ae0c0cd5f8f25b8d5001.tar.bz2
pbn-3a01eee61c1ada7fcc83ae0c0cd5f8f25b8d5001.tar.xz
Simple multiprecision library
-rw-r--r--Makefile18
-rw-r--r--mpn.h78
-rw-r--r--mpn_add.c118
-rw-r--r--mpn_cmp.c45
-rw-r--r--mpn_dump.c35
-rw-r--r--mpn_init.c164
-rw-r--r--mpn_mul.c53
-rw-r--r--mpnint.h23
-rw-r--r--test.c57
9 files changed, 591 insertions, 0 deletions
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..af7d553
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,18 @@
+CC = gcc
+CFLAGS = -g -D_FORTIFY_SOURCE=2 -W -Wall
+LDFLAGS =
+
+LIBOBJ = mpn_add.o mpn_cmp.o mpn_dump.o mpn_init.o mpn_mul.o
+
+TESTS = test
+
+.c.o:
+ $(CC) $(CFLAGS) -c -o $@ $<
+
+all: $(TESTS)
+
+test: test.o $(LIBOBJ)
+ $(CC) $(LDFLAGS) -o $@ $^
+
+clean:
+ rm -f $(TESTS) *.o
diff --git a/mpn.h b/mpn.h
new file mode 100644
index 0000000..98f4883
--- /dev/null
+++ b/mpn.h
@@ -0,0 +1,78 @@
+/* ----------------------------------------------------------------------- *
+ *
+ * Copyright 2007 H. Peter Anvin - All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, Inc.,
+ * 59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
+ * incorporated herein by reference.
+ *
+ * ----------------------------------------------------------------------- */
+
+/*
+ * mpn.h
+ *
+ * Simple multi-precision signed integer library
+ */
+
+#ifndef MPN_H
+#define MPN_H
+
+#include <stddef.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <inttypes.h>
+
+/* Configurables */
+
+#define MPN_ATOM_BITS 32
+#define MPN_ATOM_xFMT "08"PRIx32
+typedef uint32_t mpn_atom_t;
+typedef uint64_t mpn_2atom_t;
+typedef int32_t mpn_satom_t;
+typedef int64_t mpn_2satom_t;
+
+struct mpn {
+ int ref;
+ int len;
+ int bits;
+ int minus;
+ mpn_atom_t num[1];
+};
+
+
+#define _mpn_malloc(x) malloc(x)
+#define _mpn_zalloc(x) calloc(x,1)
+#define _mpn_free(x) free(x)
+#define _mpn_realloc(x,y) realloc(x,y)
+
+#define mpn_ref(x) ((x)->ref++, (x))
+
+#if defined(__GNUC__) && !defined(MPN_INIT)
+extern inline void mpn_free(struct mpn *__mpn)
+{
+ if (! --__mpn->ref)
+ _mpn_free(__mpn);
+}
+#else
+void mpn_free(struct mpn *);
+#endif
+struct mpn *mpn_addsub(struct mpn *, struct mpn *, int issub);
+#define mpn_add(x,y) mpn_addsub(x,y,0)
+#define mpn_sub(x,y) mpn_addsub(x,y,1)
+int mpn_cmp(struct mpn *, struct mpn *);
+struct mpn *mpn_mul(struct mpn *, struct mpn *);
+struct mpn *mpn_shr(struct mpn *, int);
+struct mpn *mpn_shl(struct mpn *, int);
+void mpn_dump(FILE *, struct mpn *);
+struct mpn *mpn_int(mpn_satom_t);
+struct mpn *mpn_new(int);
+struct mpn *mpn_dup(const struct mpn *);
+struct mpn *mpn_dupx(const struct mpn *, int);
+struct mpn *mpn_cow(struct mpn *, int);
+
+/* Internal functions */
+void _mpn_adjust_bits(struct mpn *mpn);
+
+#endif /* MPN_H */
diff --git a/mpn_add.c b/mpn_add.c
new file mode 100644
index 0000000..d5c732a
--- /dev/null
+++ b/mpn_add.c
@@ -0,0 +1,118 @@
+/* ----------------------------------------------------------------------- *
+ *
+ * Copyright 2007 H. Peter Anvin - All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, Inc.,
+ * 59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
+ * incorporated herein by reference.
+ *
+ * ----------------------------------------------------------------------- */
+
+/*
+ * mpn_add.c
+ */
+
+#include "mpnint.h"
+
+struct mpn *mpn_addsub(struct mpn *s1, struct mpn *s2, int issub)
+{
+ int bits, len;
+ struct mpn *d;
+ int i;
+ int s1minus = s1->minus;
+ int s2minus = s2->minus ^ issub;
+ int dminus;
+
+ bits = s2->bits > s1->bits ? s2->bits+1 : s1->bits+1;
+ len = (bits+MPN_ATOM_BITS-1)/MPN_ATOM_BITS;
+
+ if (s1 == s2) {
+ s1 = mpn_dupx(s1, len);
+ mpn_free(s2); /* Drop one reference */
+ }
+
+ if (s1minus == s2minus) {
+ mpn_2atom_t c;
+ /* Sign is the same, addition */
+
+ if (s1->ref > s2->ref) {
+ struct mpn *t = s1;
+ s1 = s2;
+ s2 = t;
+ }
+
+ d = mpn_cow(s1, len);
+ d->bits = bits;
+
+ c = 0;
+ for (i = 0; i < len; i++) {
+ c += d->num[i];
+ if (i < s2->len)
+ c += s2->num[i];
+
+ d->num[i] = c;
+ c >>= MPN_ATOM_BITS;
+ }
+ mpn_free(s2);
+
+ /* d->minus already copied from s1->minus */
+ } else {
+ mpn_2satom_t c; /* Relies on 2's complement with
+ proper right arithmetic shift... */
+ /* We know the sign is different */
+ dminus = 0;
+
+ if (s2minus) {
+ struct mpn *t = s1;
+ s1 = s2;
+ s2 = t;
+ dminus = !dminus;
+ }
+
+ /* Now s1 is the addend and s2 is the subtrahend */
+ if (mpn_cmp(s1, s2) < 0) {
+ /* addend is smaller than subtrahend, reverse and
+ invert minus */
+ struct mpn *t = s1;
+ s1 = s2;
+ s2 = t;
+ dminus = !dminus;
+ }
+
+ if (s1->ref > s2->ref) {
+ d = mpn_cow(s2, len);
+
+ c = 0;
+ for (i = 0; i < len; i++) {
+ c -= d->num[i];
+ if (i < s1->len)
+ c += s1->num[i];
+
+ d->num[i] = c;
+ c >>= MPN_ATOM_BITS;
+ }
+ mpn_free(s1);
+ } else {
+ d = mpn_cow(s1, len);
+
+ c = 0;
+ for (i = 0; i < len; i++) {
+ c += d->num[i];
+ if (i < s2->len)
+ c -= s2->num[i];
+
+ d->num[i] = c;
+ c >>= MPN_ATOM_BITS;
+ }
+ mpn_free(s2);
+ }
+
+ d->minus = dminus;
+ }
+
+ _mpn_adjust_bits(d);
+
+ return d;
+}
diff --git a/mpn_cmp.c b/mpn_cmp.c
new file mode 100644
index 0000000..59ffcf4
--- /dev/null
+++ b/mpn_cmp.c
@@ -0,0 +1,45 @@
+/* ----------------------------------------------------------------------- *
+ *
+ * Copyright 2007 H. Peter Anvin - All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, Inc.,
+ * 59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
+ * incorporated herein by reference.
+ *
+ * ----------------------------------------------------------------------- */
+
+/*
+ * mpn_cmp.c
+ */
+
+#include "mpnint.h"
+
+int mpn_cmp(struct mpn *s1, struct mpn *s2)
+{
+ int lt = s1->minus ? 1 : -1;
+ int gt = -lt;
+ int i;
+
+ if (s1->minus != s2->minus)
+ return gt;
+
+ if (s1->bits < s2->bits)
+ return lt;
+ else if (s1->bits > s2->bits)
+ return gt;
+
+ i = (s1->bits+MPN_ATOM_BITS-1)/MPN_ATOM_BITS;
+ while (i >= 0) {
+ if (s1->num[i] > s2->num[i])
+ return gt;
+ else if (s1->num[i] < s2->num[i])
+ return lt;
+ i--;
+ }
+
+ return 0;
+}
+
+
diff --git a/mpn_dump.c b/mpn_dump.c
new file mode 100644
index 0000000..bffb165
--- /dev/null
+++ b/mpn_dump.c
@@ -0,0 +1,35 @@
+/* ----------------------------------------------------------------------- *
+ *
+ * Copyright 2007 H. Peter Anvin - All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, Inc.,
+ * 59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
+ * incorporated herein by reference.
+ *
+ * ----------------------------------------------------------------------- */
+
+/*
+ * mpn_dump.c
+ *
+ * Output an MPN in hexadecimal form, mainly intended for debugging
+ */
+
+#include <stdio.h>
+#include "mpnint.h"
+
+void mpn_dump(FILE *f, struct mpn *mpn)
+{
+ int i;
+
+ putc(mpn->minus ? '-' : '+', f);
+
+ for (i = mpn->len-1; i >= 0; i--) {
+ fprintf(f, "%" MPN_ATOM_xFMT, mpn->num[i]);
+ if (i)
+ putc('_', f);
+ }
+
+ mpn_free(mpn);
+}
diff --git a/mpn_init.c b/mpn_init.c
new file mode 100644
index 0000000..2770b9f
--- /dev/null
+++ b/mpn_init.c
@@ -0,0 +1,164 @@
+/* ----------------------------------------------------------------------- *
+ *
+ * Copyright 2007 H. Peter Anvin - All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, Inc.,
+ * 59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
+ * incorporated herein by reference.
+ *
+ * ----------------------------------------------------------------------- */
+
+/*
+ * mpn_int.c
+ */
+
+#define MPN_INIT 1
+#include <string.h>
+#include "mpnint.h"
+
+struct mpn *mpn_new(int len)
+{
+ struct mpn *mpn;
+
+ mpn = _mpn_zalloc(sizeof(struct mpn)+(len-1)*sizeof(mpn_atom_t));
+ if (!mpn)
+ return NULL;
+
+ mpn->ref = 1;
+ mpn->len = len;
+
+ return mpn;
+}
+
+void mpn_free(struct mpn *mpn)
+{
+ if (! --mpn->ref)
+ _mpn_free(mpn);
+}
+
+#if MPN_ATOM_BITS != 32 && MPN_ATOM_BITS != 64
+# error "Fix mpn_ilog2p1 for new MPN_ATOM_BITS size"
+#endif
+
+static int mpn_ilog2p1(mpn_atom_t v)
+{
+ int p = 1;
+
+#if MPN_ATOM_BITS >= 64
+ if (v & 0xffffffff00000000) {
+ p += 32;
+ v >>= 32;
+ }
+#endif
+
+ if (v & 0xffff0000) {
+ p += 16;
+ v >>= 16;
+ }
+ if (v & 0xff00) {
+ p += 8;
+ v >>= 8;
+ }
+ if (v & 0xf0) {
+ p += 4;
+ v >>= 4;
+ }
+ if (v & 0xc) {
+ p += 2;
+ v >>= 2;
+ }
+ if (v & 0x2) {
+ p += 1;
+ v >>= 1;
+ }
+
+ return p;
+}
+
+struct mpn *mpn_int(mpn_satom_t v)
+{
+ struct mpn *mpn = mpn_new(1);
+ if (!mpn)
+ return NULL;
+
+ if (v < 0) {
+ v = -v;
+ mpn->minus = 1;
+ } else {
+ mpn->minus = 0;
+ }
+
+ mpn->num[0] = v;
+ mpn->bits = mpn_ilog2p1(v);
+ return mpn;
+}
+
+/* Produce a duplicate */
+struct mpn *mpn_dup(const struct mpn *src)
+{
+ int bytes = sizeof(struct mpn)+(src->len-1)*sizeof(mpn_atom_t);
+ struct mpn *dst = _mpn_malloc(bytes);
+
+ memcpy(dst, src, bytes);
+ dst->ref = 1;
+
+ return dst;
+}
+
+/* Produce a duplicate, possibly extended */
+struct mpn *mpn_dupx(const struct mpn *src, int len)
+{
+ int slack = src->len >= len ? 0 : len-src->len;
+ int bytes = sizeof(struct mpn)+(src->len-1)*sizeof(mpn_atom_t);
+ struct mpn *dst = _mpn_malloc(bytes+slack*sizeof(mpn_atom_t));
+
+ if (!dst)
+ return NULL;
+
+ memcpy(dst, src, bytes);
+ dst->ref = 1;
+ memset(&dst->num[slack], 0, slack*sizeof(mpn_atom_t));
+
+ return dst;
+}
+
+/* Returns a duplicate or clears src for reuse, consumes a reference */
+struct mpn *mpn_cow(struct mpn *src, int len)
+{
+ int olen = src->len;
+
+ if (src->ref <= 1) {
+ src->ref = 1;
+ if (len > olen) {
+ src = _mpn_realloc(src, sizeof(struct mpn)+
+ (len-1)*sizeof(mpn_atom_t));
+ src->len = len;
+ }
+ memset(&src->num[olen], 0, (len-olen)*sizeof(mpn_atom_t));
+ return src;
+ }
+
+ src->ref--;
+ return mpn_dupx(src, len);
+}
+
+/* Adjust the bits field of a specific mpn, and clear the minus flag
+ if the value is actually zero. */
+void _mpn_adjust_bits(struct mpn *mpn)
+{
+ int i;
+ mpn_atom_t v;
+
+ for (i = mpn->len-1; i >= 0; i--) {
+ v = mpn->num[i];
+ if (v) {
+ mpn->bits = i*MPN_ATOM_BITS + mpn_ilog2p1(v);
+ return;
+ }
+ }
+
+ mpn->bits = 0;
+ mpn->minus = 0;
+}
diff --git a/mpn_mul.c b/mpn_mul.c
new file mode 100644
index 0000000..5d4477b
--- /dev/null
+++ b/mpn_mul.c
@@ -0,0 +1,53 @@
+/* ----------------------------------------------------------------------- *
+ *
+ * Copyright 2007 H. Peter Anvin - All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, Inc.,
+ * 59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
+ * incorporated herein by reference.
+ *
+ * ----------------------------------------------------------------------- */
+
+/*
+ * mpn_mul.c
+ */
+
+#include "mpnint.h"
+
+struct mpn *mpn_mul(struct mpn *s1, struct mpn *s2)
+{
+ size_t bits, len;
+ struct mpn *d;
+ int i, j;
+ mpn_atom_t x1;
+ mpn_2atom_t c;
+
+ bits = s1->bits + s2->bits;
+ len = (bits+MPN_ATOM_BITS-1)/MPN_ATOM_BITS;
+
+ d = mpn_new(len);
+ d->bits = bits;
+
+ for (i = 0; i < s1->len; i++) {
+ x1 = s1->num[i];
+ c = 0;
+ for (j = 0; j < s2->len; j++) {
+ c += (mpn_2atom_t)x1 * s2->num[j];
+ d->num[i+j] += (mpn_atom_t)c;
+ c >>= MPN_ATOM_BITS;
+ }
+ if (c)
+ d->num[i+j] += (mpn_atom_t)c;
+ }
+ d->minus = s1->minus ^ s2->minus;
+
+ mpn_free(s1);
+ mpn_free(s2);
+
+ _mpn_adjust_bits(d);
+
+ return d;
+}
+
diff --git a/mpnint.h b/mpnint.h
new file mode 100644
index 0000000..d97529d
--- /dev/null
+++ b/mpnint.h
@@ -0,0 +1,23 @@
+/* ----------------------------------------------------------------------- *
+ *
+ * Copyright 2007 H. Peter Anvin - All Rights Reserved
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, Inc.,
+ * 59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
+ * incorporated herein by reference.
+ *
+ * ----------------------------------------------------------------------- */
+
+/*
+ * mpn library internal header file
+ */
+
+#ifndef MPNINT_H
+#define MPNINT_H
+
+#include "mpn.h"
+
+
+#endif /* MPNINT_H */
diff --git a/test.c b/test.c
new file mode 100644
index 0000000..e2ecc40
--- /dev/null
+++ b/test.c
@@ -0,0 +1,57 @@
+#include <stdio.h>
+#include <inttypes.h>
+#include <stdlib.h>
+#include "mpn.h"
+
+int main(int argc, char *argv[])
+{
+ uint64_t i, j;
+ struct mpn *m, *mj;
+ int nx = (argc < 2) ? 64 : atoi(argv[1]);
+ int n;
+
+ m = mpn_int(5);
+ i = 5;
+
+ n = nx;
+ while (n--) {
+ printf("%016"PRIx64" : ", i);
+ mpn_dump(stdout, mpn_ref(m));
+ putchar('\n');
+
+ i *= 5;
+ m = mpn_mul(m, mpn_int(5));
+ }
+ mpn_free(m);
+
+ m = mpn_int(5);
+ i = 5;
+ n = nx;
+ while (n--) {
+ printf("%016"PRIx64" : ", i);
+ mpn_dump(stdout, mpn_ref(m));
+ putchar('\n');
+
+ i += i;
+ m = mpn_add(m, mpn_ref(m));
+ }
+
+ n = nx;
+ j = 20;
+ mj = mpn_int(20);
+ while (n--) {
+ printf("%016"PRIx64" : ", i);
+ mpn_dump(stdout, mpn_ref(m));
+ printf(" : ");
+ mpn_dump(stdout, mpn_ref(mj));
+ putchar('\n');
+
+ i -= j;
+ j += j;
+ m = mpn_sub(m, mpn_ref(mj));
+ mj = mpn_add(mj, mpn_ref(mj));
+ }
+
+ return 0;
+}
+