commit 4151bc7ad1433a6ed9dec777600264fa5009af0c
parent 0552439452889c80c6299cc11c34257395784d5a
Author: olikru <olikru@tkruger.se>
Date: Mon, 27 May 2024 16:59:46 +0200
working on batchgcd
Diffstat:
6 files changed, 263 insertions(+), 6 deletions(-)
diff --git a/Makefile b/Makefile
@@ -29,14 +29,19 @@ cyclefind.o\
pierre.o\
fmpzio.o\
wiener.o\
-prodtree.o
+prodtree.o\
+batchgcd.o
TOOLS=\
tools/hnpsolve\
tools/lcgfloyd\
tools/dbtool\
tools/pierre\
tools/varmkorv\
-tools/small
+tools/small\
+tools/batchgcd\
+tools/pelofskegcd\
+tools/gen\
+tools/fileify
PRECOMPUTERS=\
precomputers/header
TEST=test
@@ -67,18 +72,30 @@ test: testbin
tools/hnpsolve: $(OBJ) $(@:S/$/.c/)
$(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/hnpsolve.c $(OBJ) $(LDFLAGS)
-tools/lcgfloyd: $(OBJ)
+tools/lcgfloyd: $(OBJ) $(@:S/$/.c/)
$(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/lcgfloyd.c $(OBJ) $(LDFLAGS)
-tools/pierre: $(OBJ)
+tools/pierre: $(OBJ) $(@:S/$/.c/)
$(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/pierre.c $(OBJ) $(LDFLAGS)
-tools/varmkorv: $(OBJ)
+tools/varmkorv: $(OBJ) $(@:S/$/.c/)
$(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/varmkorv.c $(OBJ) $(LDFLAGS)
-tools/small: $(OBJ)
+tools/small: $(OBJ) $(@:S/$/.c/)
$(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/small.c $(OBJ) $(LDFLAGS)
+tools/batchgcd: $(OBJ) $(@:S/$/.c/)
+ $(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/batchgcd.c $(OBJ) $(LDFLAGS)
+
+tools/pelofskegcd: $(OBJ) $(@:S/$/.c/)
+ $(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/pelofskegcd.c $(OBJ) $(LDFLAGS)
+
+tools/gen: $(OBJ) $(@:S/$/.c/)
+ $(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/gen.c $(OBJ) $(LDFLAGS)
+
+tools/fileify: $(OBJ) $(@:S/$/.c/)
+ $(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/fileify.c $(OBJ) $(LDFLAGS)
+
tools/dbtool: $(OBJ) $(@:S/$/.c/)
$(CC) -I. $(CFLAGS) -o $@ $(TOOLS_DIR)/dbtool.c $(OBJ) $(LDFLAGS) -lsqlite3
diff --git a/batchgcd.c b/batchgcd.c
@@ -0,0 +1,120 @@
+#include <err.h>
+#include <fmpz.h>
+#include <fmpz_vec.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "prodtree.h"
+#include "batchgcd.h"
+
+static inline void
+remtree_compute(fmpz *tree, const slong tree_size)
+{
+ slong i = tree_size - 2;
+ slong parent_distance = 1;
+
+ while (1) {
+ fmpz_mul(&tree[i], &tree[i], &tree[i]);
+ fmpz_mod(&tree[i], &tree[i + parent_distance], &tree[i]);
+ parent_distance++;
+ i--;
+
+ fmpz_mul(&tree[i], &tree[i], &tree[i]);
+ fmpz_mod(&tree[i], &tree[i + parent_distance], &tree[i]);
+ if(i == 0) break;
+ i--;
+ }
+}
+
+void
+batch_gcd(fmpz *r, const fmpz *v, slong n)
+{
+ fmpz *tree;
+ slong tree_size, i;
+
+ tree_size = prodtree_compute(&tree, v, n);
+ remtree_compute(tree, tree_size);
+
+ for(i = 0; i < n; i++) {
+ fmpz_fdiv_q(&r[i], &tree[i], &v[i]);
+ fmpz_gcd(&r[i], &r[i], &v[i]);
+ }
+
+ free(tree);
+}
+
+static inline slong
+upper_bound_2pow(const slong x)
+{
+ slong r = 1;
+
+ while (r < x) {
+ r <<= 1;
+
+ if (r == 0) {
+ err(EXIT_FAILURE, "upper bound 2pow failure");
+ }
+ }
+
+ return r;
+}
+
+
+static inline slong
+incomplete_prodtree_compute(fmpz** res, const fmpz *v, slong nv)
+{
+ slong ubound = upper_bound_2pow(nv);
+ slong tree_size = 2 * ubound - 1;
+ slong i, j;
+
+ *res = _fmpz_vec_init(tree_size);
+
+ // row 0, v and padded with ones
+ for (i = 0; i < nv; i++) {
+ fmpz_set(&(*res)[i], &v[i]);
+ }
+ for (i = nv; i < ubound; i++) {
+ fmpz_set_ui(&(*res)[i], 1);
+ }
+
+ // rest, not root
+ j = 0;
+ for (i = ubound; i < tree_size-1; i++) {
+ fmpz_mul(&(*res)[i], &(*res)[j], &(*res)[j + 1]);
+ j += 2;
+ }
+
+ return tree_size;
+}
+
+// see: https://github.com/epelofske65537/binary_tree_Batch_GCD
+void
+pelofske_gcd(fmpz *r, const fmpz *v, slong n)
+{
+ fmpz *tree;
+ fmpz *bvec;
+ fmpz *bbtree;
+ slong tree_size, i, ninb, bbtree_size;
+
+ tree_size = incomplete_prodtree_compute(&tree, v, n);
+
+ bvec = _fmpz_vec_init(tree_size);
+ ninb = 0;
+
+ for(i = 0; i < tree_size-1; i+=2) {
+ fmpz_gcd(&bvec[ninb], &tree[i], &tree[i+1]);
+ if (fmpz_cmp_ui(&bvec[ninb], 1) != 0) {
+ ninb++;
+ }
+ }
+
+ bbtree_size = prodtree_compute(&bbtree, bvec, ninb);
+
+ for(i = 0; i < n; i++) {
+ fmpz_gcd(&r[i], &bbtree[bbtree_size-1], &v[i]);
+ }
+
+ _fmpz_vec_clear(tree, tree_size);
+ _fmpz_vec_clear(bvec, tree_size);
+ _fmpz_vec_clear(bbtree, bbtree_size);
+}
diff --git a/batchgcd.h b/batchgcd.h
@@ -0,0 +1,7 @@
+#ifndef _BATCHGCD_H_
+#define _BATCHGCD_H_
+
+void batch_gcd(fmpz *r, const fmpz *v, slong n);
+void pelofske_gcd(fmpz *r, const fmpz *v, slong n);
+
+#endif
diff --git a/tools/batchgcd.c b/tools/batchgcd.c
@@ -0,0 +1,29 @@
+#include <flint.h>
+#include <fmpz_vec.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "batchgcd.h"
+#include "fmpzio.h"
+
+int
+main(int argc, char *argv[])
+{
+ fmpz *v, *res;
+ slong nv, i;
+
+ nv = read_hex_lines(&v);
+ res = _fmpz_vec_init(nv);
+
+ batch_gcd(res, v, nv);
+
+ for (i = 0; i < nv; i++) {
+ fmpz_print(&res[i]);
+ printf("\n");
+ }
+
+ free(v);
+ free(res);
+
+ return 0;
+}
diff --git a/tools/gen.c b/tools/gen.c
@@ -0,0 +1,55 @@
+#include <flint.h>
+#include <fmpz_vec.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "fmpzio.h"
+
+static void
+print_hex(fmpz_t n)
+{
+ char *s = fmpz_get_str(NULL, 16, n);
+ printf("%s\n", s);
+ free(s);
+}
+
+static inline void
+print_usage(const char *prog) {
+ printf("Usage: %s <number of moduli> <bit size of modulus>\n", prog);
+}
+
+int
+main(int argc, char *argv[])
+{
+ size_t i;
+ uint64_t nmoduli, moduli_size;
+ flint_rand_t state;
+ fmpz_t p;
+
+ if(argc != 3) {
+ print_usage(argv[0]);
+ exit(EXIT_FAILURE);
+ }
+
+ nmoduli = strtoull(argv[1], NULL, 10);
+ moduli_size = strtoull(argv[2], NULL, 10);
+
+ flint_randinit(state);
+ fmpz_init(p);
+
+ // fmpz_randprime is not uniform, but good enough for the
+ // current purpose
+
+ for (i = 0; i < nmoduli; i++) {
+ fmpz_randprime(p, state, moduli_size >> 1, 0);
+ while (fmpz_sizeinbase(p, 2) < (moduli_size >> 1)) {
+ fmpz_randprime(p, state, moduli_size >> 1, 0);
+ }
+ print_hex(p);
+ }
+
+ flint_randclear(state);
+ fmpz_clear(p);
+
+ return 0;
+}
diff --git a/tools/pelofskegcd.c b/tools/pelofskegcd.c
@@ -0,0 +1,29 @@
+#include <flint.h>
+#include <fmpz_vec.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "batchgcd.h"
+#include "fmpzio.h"
+
+int
+main(int argc, char *argv[])
+{
+ fmpz *v, *res;
+ slong nv, i;
+
+ nv = read_hex_lines(&v);
+ res = _fmpz_vec_init(nv);
+
+ pelofske_gcd(res, v, nv);
+
+ for (i = 0; i < nv; i++) {
+ fmpz_print(&res[i]);
+ printf("\n");
+ }
+
+ free(v);
+ free(res);
+
+ return 0;
+}
