25 files changed, 772 insertions, 612 deletions

diff --git a/src/exhaustive/brainpool.c b/src/exhaustive/brainpool.c
index 17c2c65..62a3352 100644
--- a/src/exhaustive/brainpool.c
+++ b/src/exhaustive/brainpool.c
@@ -3,11 +3,13 @@
  * Copyright (C) 2017-2018 J08nY
  */
 
+#include <misc/types.h>
 #include "brainpool.h"
 #include "gen/gens.h"
 #include "gen/point.h"
 #include "gen/seed.h"
 #include "io/output.h"
+#include "math/subgroup.h"
 #include "util/bits.h"
 #include "util/str.h"
 
@@ -237,10 +239,12 @@ GENERATOR(brainpool_gen_gens) {
 		return INT_MIN;
 	}
 
-	curve->generators = points_new(1);
+	curve->generators = subgroups_new(1);
 	curve->ngens = 1;
+	subgroup_t *sub = subgroup_new();
+	curve->generators[0] = sub;
 	point_t *G = point_new();
-	curve->generators[0] = G;
+	sub->generator = G;
 	G->point = ellmul(curve->curve, P, k);
 	G->order = ellorder(curve->curve, G->point, NULL);
 	G->cofactor = divii(curve->order, G->order);
@@ -252,7 +256,7 @@ GENERATOR(brainpool_gen_gens) {
 
 CHECK(brainpool_check_gens) {
 	pari_sp ltop = avma;
-	point_t *G = curve->generators[0];
+	point_t *G = curve->generators[0]->generator;
 	GEN min_degree = divis(subii(G->order, gen_1), 100);
 	if (mpcmp(min_degree, gens_get_embedding(curve->field, G->order)) >= 0) {
 		avma = ltop;
diff --git a/src/exhaustive/exhaustive.c b/src/exhaustive/exhaustive.c
index e7d6350..7ad4705 100644
--- a/src/exhaustive/exhaustive.c
+++ b/src/exhaustive/exhaustive.c
@@ -43,13 +43,15 @@ static void exhaustive_ginit(gen_f *generators) {
 		if (cfg->prime) {
 			generators[OFFSET_ORDER] = &order_gen_prime;
 		} else if (cfg->cofactor) {
-			generators[OFFSET_ORDER] = &order_gen_smallfact;
+			generators[OFFSET_ORDER] = &order_gen_cofactor;
 		} else {
 			generators[OFFSET_ORDER] = &order_gen_any;
 		}
 
 		if (cfg->unique) {
 			generators[OFFSET_GENERATORS] = &gens_gen_one;
+		} else if (cfg->cofactor) {
+			generators[OFFSET_GENERATORS] = &gens_gen_cofactor;
 		} else {
 			generators[OFFSET_GENERATORS] = &gens_gen_any;
 		}
@@ -144,7 +146,7 @@ static void exhaustive_ginit(gen_f *generators) {
 		if (cfg->prime) {
 			generators[OFFSET_ORDER] = &order_gen_prime;
 		} else if (cfg->cofactor) {
-			generators[OFFSET_ORDER] = &order_gen_smallfact;
+			generators[OFFSET_ORDER] = &order_gen_cofactor;
 		} else if (cfg->method == METHOD_ANOMALOUS) {
 			generators[OFFSET_ORDER] = &anomalous_gen_order;
 		} else {
@@ -161,6 +163,8 @@ static void exhaustive_ginit(gen_f *generators) {
 
 		if (cfg->unique) {
 			generators[OFFSET_GENERATORS] = &gens_gen_one;
+		} else if (cfg->cofactor) {
+			generators[OFFSET_GENERATORS] = &gens_gen_cofactor;
 		} else {
 			generators[OFFSET_GENERATORS] = &gens_gen_any;
 		}
@@ -249,9 +253,9 @@ static void exhaustive_ainit(arg_t **gen_argss, arg_t **check_argss) {
 	if (cfg->cofactor) {
 		arg_t *order_arg = arg_new();
 		arg_t *gens_arg = arg_new();
-		order_arg->args = &cfg->cofactor_bound;
+		order_arg->args = &cfg->cofactor_value;
 		order_arg->nargs = 1;
-		gens_arg->args = &cfg->cofactor_bound;
+		gens_arg->args = &cfg->cofactor_value;
 		gens_arg->nargs = 1;
 		gen_argss[OFFSET_ORDER] = order_arg;
 		gen_argss[OFFSET_GENERATORS] = gens_arg;
diff --git a/src/gen/curve.c b/src/gen/curve.c
index f89044a..8c7d4c9 100644
--- a/src/gen/curve.c
+++ b/src/gen/curve.c
@@ -4,7 +4,7 @@
  */
 #include "curve.h"
 #include "math/twists.h"
-#include "point.h"
+#include "math/subgroup.h"
 #include "seed.h"
 #include "util/memory.h"
 
@@ -18,13 +18,9 @@ curve_t *curve_copy(const curve_t *src, curve_t *dest) {
 	if (src->curve) dest->curve = gcopy(src->curve);
 	if (src->order) dest->order = gcopy(src->order);
 	if (src->generators) {
-		dest->generators = points_new_copy(src->generators, src->ngens);
+		dest->generators = subgroups_new_copy(src->generators, src->ngens);
 		dest->ngens = src->ngens;
 	}
-	if (src->points) {
-		dest->points = points_new_copy(src->points, src->npoints);
-		dest->npoints = src->npoints;
-	}
 	return dest;
 }
 
@@ -41,13 +37,9 @@ curve_t *curve_clone(const curve_t *src, curve_t *dest) {
 	if (src->curve) dest->curve = gclone(src->curve);
 	if (src->order) dest->order = gclone(src->order);
 	if (src->generators) {
-		dest->generators = points_new_clone(src->generators, src->ngens);
+		dest->generators = subgroups_new_clone(src->generators, src->ngens);
 		dest->ngens = src->ngens;
 	}
-	if (src->points) {
-		dest->points = points_new_clone(src->points, src->npoints);
-		dest->npoints = src->npoints;
-	}
 	return dest;
 }
 
@@ -59,8 +51,7 @@ curve_t *curve_new_clone(const curve_t *src) {
 void curve_free(curve_t **curve) {
 	if (*curve) {
 		seed_free(&(*curve)->seed);
-		points_free_deep(&(*curve)->generators, (*curve)->ngens);
-		points_free_deep(&(*curve)->points, (*curve)->npoints);
+		subgroups_free_deep(&(*curve)->generators, (*curve)->ngens);
 
 		if ((*curve)->curve) {
 			// TODO, this is possibly dangerous...
diff --git a/src/gen/gens.c b/src/gen/gens.c
index ba94c64..9d870e7 100644
--- a/src/gen/gens.c
+++ b/src/gen/gens.c
@@ -3,21 +3,28 @@
  * Copyright (C) 2017-2018 J08nY
  */
 #include "gens.h"
+#include "math/subgroup.h"
 #include "exhaustive/arg.h"
 #include "point.h"
 
+
+static subgroup_t *gens_point(GEN point, const curve_t *curve) {
+	subgroup_t *sub = subgroup_new();
+	point_t *p = point_new();
+	sub->generator = p;
+	p->point = gcopy(point);
+	p->order = ellorder(curve->curve, p->point, NULL);
+	p->cofactor = divii(curve->order, p->order);
+	return sub;
+}
+
 static int gens_put(curve_t *curve, GEN generators, long len) {
-	curve->generators = points_new((size_t)len);
-	curve->ngens = (size_t)len;
+	curve->generators = subgroups_new((size_t) len);
+	curve->ngens = (size_t) len;
 
 	for (long i = 1; i <= len; ++i) {
-		point_t *p = point_new();
-		p->point = gcopy(gel(generators, i));
-		p->order = ellorder(curve->curve, p->point, NULL);
-		p->cofactor = divii(curve->order, p->order);
-		curve->generators[i - 1] = p;
+		curve->generators[i - 1] = gens_point(gel(generators, i), curve);
 	}
-
 	return 1;
 }
 
@@ -38,10 +45,56 @@ GENERATOR(gens_gen_one) {
 	return gens_put(curve, generators, len);
 }
 
+GENERATOR(gens_gen_cofactor) {
+	HAS_ARG(args);
+	pari_ulong cofactor = *(pari_ulong *) args->args;
+	pari_sp ltop = avma;
+	GEN order = diviuexact(curve->order, cofactor);
+
+	GEN generators = ellff_get_gens(curve->curve);
+	long len = glength(generators);
+	point_t *p = NULL;
+	for (long i = 1; i <= len; ++i) {
+		GEN gen = gel(generators, i);
+		GEN gen_order = ellorder(curve->curve, gen, NULL);
+
+		if (equalii(order, gen_order)) {
+			p = point_new();
+			p->point = gcopy(gen);
+			p->order = gen_order;
+			p->cofactor = utoi(cofactor);
+			break;
+		}
+		GEN res = cgeti(DEFAULTPREC);
+		if (dvdiuz(gen_order, cofactor, res)) {
+			p = point_new();
+			p->point = gcopy(ellmul(curve->curve, gen, utoi(cofactor)));
+			p->order = res;
+			p->cofactor = utoi(cofactor);
+			break;
+		}
+	}
+
+	if (p) {
+		curve->ngens = (size_t) (len + 1);
+		curve->generators = subgroups_new(curve->ngens);
+		for (long i = 1; i <= len; ++i) {
+			curve->generators[i] = gens_point(gel(generators, i), curve);
+		}
+		subgroup_t *sub = subgroup_new();
+		sub->generator = p;
+		curve->generators[0] = sub;
+		return 1;
+	} else {
+		avma = ltop;
+		return -5;
+	}
+}
+
 CHECK(gens_check_anomalous) {
 	if (cfg->field == FIELD_BINARY) return 1;
 	for (size_t i = 0; i < curve->ngens; ++i) {
-		if (mpcmp(curve->field, curve->generators[i]->order) == 0) {
+		if (mpcmp(curve->field, curve->generators[i]->generator->order) == 0) {
 			return -5;
 		}
 	}
@@ -62,7 +115,7 @@ CHECK(gens_check_embedding) {
 
 	for (size_t i = 0; i < curve->ngens; ++i) {
 		GEN power =
-		    gens_get_embedding(curve->field, curve->generators[i]->order);
+				gens_get_embedding(curve->field, curve->generators[i]->generator->order);
 
 		if (mpcmp(power, mind) <= 0) {
 			avma = ltop;
@@ -75,7 +128,7 @@ CHECK(gens_check_embedding) {
 
 UNROLL(gens_unroll) {
 	if (curve->generators) {
-		points_free_deep(&curve->generators, curve->ngens);
+		subgroups_free_deep(&curve->generators, curve->ngens);
 	}
 	return -1;
 }
diff --git a/src/gen/gens.h b/src/gen/gens.h
index 0929bbf..4d1ea88 100644
--- a/src/gen/gens.h
+++ b/src/gen/gens.h
@@ -16,6 +16,7 @@
  *
  * @param curve A curve_t being generated
  * @param args unused
+ * @param state
  * @return state diff
  */
 GENERATOR(gens_gen_any);
@@ -25,10 +26,22 @@ GENERATOR(gens_gen_any);
  *
  * @param curve A curve_t being generated
  * @param args unused
+ * @param state
  * @return state diff
  */
 GENERATOR(gens_gen_one);
 
+
+/**
+ * GENERATOR(gen_f)
+ *
+ * @param curve
+ * @param args
+ * @param state
+ * @return
+ */
+GENERATOR(gens_gen_cofactor);
+
 /**
  * CHECK(check_f)
  *
diff --git a/src/gen/gp.c b/src/gen/gp.c
deleted file mode 100644
index 1800f26..0000000
--- a/src/gen/gp.c
+++ /dev/null
@@ -1,176 +0,0 @@
-/*
- * ecgen, tool for generating Elliptic curve domain parameters
- * Copyright (C) 2017-2018 J08nY
- */
-#include "gp.h"
-#include "exhaustive/arg.h"
-#include "io/output.h"
-#include "point.h"
-#include "seed.h"
-#include "util/bits.h"
-
-static point_t **gp_points(const curve_t *curve, GEN point_vec) {
-	long len = glength(point_vec);
-	point_t **result = points_new((size_t)len);
-
-	for (long i = 1; i <= len; ++i) {
-		point_t *point = point_new();
-		point->point = gel(point_vec, i);
-		point->order = ellorder(curve->curve, point->point, NULL);
-		result[i - 1] = point;
-	}
-	return result;
-}
-
-static point_t **gp_gens(const curve_t *curve, GEN gens_vec) {
-	point_t **result = gp_points(curve, gens_vec);
-
-	long len = glength(gens_vec);
-	for (long i = 1; i <= len; ++i) {
-		point_t *gen = result[i - 1];
-		gen->cofactor = divii(curve->order, gen->order);
-	}
-	return result;
-}
-
-GENERATOR(gp_gen) {
-	HAS_ARG(args);
-	pari_sp ltop = avma;
-	GEN closure = compile_str(args->args);
-	GEN params = zerovec(state - OFFSET_SEED);
-
-	if (state > OFFSET_SEED) {
-		if (curve->seed && curve->seed->seed) {
-			gel(params, 1) = bits_to_bitvec(curve->seed->seed);
-		}
-	}
-
-	if (state > OFFSET_FIELD) {
-		gel(params, 2) = curve->field;
-	}
-
-	if (state > OFFSET_A) {
-		gel(params, 3) = curve->a;
-	}
-
-	if (state > OFFSET_B) {
-		gel(params, 4) = curve->b;
-	}
-
-	if (state > OFFSET_CURVE) {
-		gel(params, 5) = curve->curve;
-	}
-
-	if (state > OFFSET_ORDER) {
-		gel(params, 6) = curve->order;
-	}
-
-	if (state > OFFSET_GENERATORS) {
-		GEN gens = zerovec(curve->ngens);
-		for (size_t i = 0; i < curve->ngens; ++i) {
-			gel(gens, i + 1) = curve->generators[i]->point;
-		}
-		gel(params, 7) = gens;
-	}
-
-	if (state > OFFSET_POINTS) {
-		GEN points = zerovec(curve->npoints);
-		for (size_t i = 0; i < curve->npoints; ++i) {
-			gel(points, i + 1) = curve->points[i]->point;
-		}
-		gel(params, 8) = points;
-	}
-
-	GEN res = call0(closure, zerovec(0));
-	res = call0(res, params);
-
-	res = gerepileupto(ltop, res);
-	switch (state) {
-		case OFFSET_SEED:
-			curve->seed = seed_new();
-			curve->seed->seed = bits_from_bitvec(res);
-			break;
-		case OFFSET_FIELD:
-			curve->field = res;
-			break;
-		case OFFSET_A:
-			curve->a = res;
-			break;
-		case OFFSET_B:
-			curve->b = res;
-			break;
-		case OFFSET_CURVE:
-			curve->curve = res;
-			break;
-		case OFFSET_ORDER:
-			curve->order = res;
-			break;
-		case OFFSET_GENERATORS:
-			curve->ngens = (size_t)glength(res);
-			curve->generators = gp_gens(curve, res);
-			break;
-		case OFFSET_POINTS:
-			curve->npoints = (size_t)glength(res);
-			curve->points = gp_points(curve, res);
-			break;
-		case OFFSET_END:
-			break;
-	}
-	return 1;
-}
-
-CHECK(gp_check) {
-	HAS_ARG(args);
-	pari_sp ltop = avma;
-	GEN closure = compile_str(args->args);
-	GEN params = zerovec(state - OFFSET_SEED + 1);
-
-	if (state >= OFFSET_SEED) {
-		if (curve->seed && curve->seed->seed) {
-			gel(params, 1) = bits_to_bitvec(curve->seed->seed);
-		}
-	}
-
-	if (state >= OFFSET_FIELD) {
-		gel(params, 2) = curve->field;
-	}
-
-	if (state >= OFFSET_A) {
-		gel(params, 3) = curve->a;
-	}
-
-	if (state >= OFFSET_B) {
-		gel(params, 4) = curve->b;
-	}
-
-	if (state >= OFFSET_CURVE) {
-		gel(params, 5) = curve->curve;
-	}
-
-	if (state >= OFFSET_ORDER) {
-		gel(params, 6) = curve->order;
-	}
-
-	if (state >= OFFSET_GENERATORS) {
-		GEN gens = zerovec(curve->ngens);
-		for (size_t i = 0; i < curve->ngens; ++i) {
-			gel(gens, i + 1) = curve->generators[i]->point;
-		}
-		gel(params, 7) = gens;
-	}
-
-	if (state >= OFFSET_POINTS) {
-		GEN points = zerovec(curve->npoints);
-		for (size_t i = 0; i < curve->npoints; ++i) {
-			gel(points, i + 1) = curve->points[i]->point;
-		}
-		gel(params, 8) = points;
-	}
-
-	GEN res = call0(closure, zerovec(0));
-	res = call0(res, params);
-
-	int result = (int)itos(res);
-	avma = ltop;
-	return result;
-}
\ No newline at end of file
diff --git a/src/gen/gp.h b/src/gen/gp.h
deleted file mode 100644
index 4506816..0000000
--- a/src/gen/gp.h
+++ /dev/null
@@ -1,30 +0,0 @@
-/*
- * ecgen, tool for generating Elliptic curve domain parameters
- * Copyright (C) 2017-2018 J08nY
- */
-/**
- * @file gp.h
- */
-#ifndef ECGEN_GP_H
-#define ECGEN_GP_H
-
-#include "misc/types.h"
-
-/**
- * @brief
- * @param curve
- * @param args
- * @return
- */
-GENERATOR(gp_gen);
-
-/**
- * @brief
- * @param curve
- * @param args
- * @param state
- * @return
- */
-CHECK(gp_check);
-
-#endif  // ECGEN_GP_H
diff --git a/src/gen/hex.c b/src/gen/hex.c
index 39200a8..08c280f 100644
--- a/src/gen/hex.c
+++ b/src/gen/hex.c
@@ -2,20 +2,30 @@
  * ecgen, tool for generating Elliptic curve domain parameters
  * Copyright (C) 2017-2018 J08nY
  */
+#include <misc/types.h>
 #include "hex.h"
 #include "exhaustive/arg.h"
 #include "field.h"
 #include "util/bits.h"
 #include "util/memory.h"
+#include "util/str.h"
+
+static char *hex_point(point_t *point) {
+	GEN fx = field_elementi(gel(point->point, 1));
+	GEN fy = field_elementi(gel(point->point, 2));
+	char *fxs = pari_sprintf("%P0#*x", cfg->hex_digits, fx);
+	char *fxy = pari_sprintf("%P0#*x", cfg->hex_digits, fy);
+	char *result = str_joinv(",", fxs, fxy, NULL);
+	pari_free(fxs);
+	pari_free(fxy);
+	return result;
+}
 
 static char *hex_points(point_t *points[], size_t len) {
 	char *p[len];
 	for (size_t i = 0; i < len; ++i) {
 		point_t *pt = points[i];
-		GEN fx = field_elementi(gel(pt->point, 1));
-		GEN fy = field_elementi(gel(pt->point, 2));
-		p[i] = pari_sprintf("%P0#*x,%P0#*x,", cfg->hex_digits, fx,
-		                    cfg->hex_digits, fy);
+		p[i] = hex_point(pt);
 	}
 
 	size_t total = 1;
@@ -26,7 +36,7 @@ static char *hex_points(point_t *points[], size_t len) {
 	char *result = try_calloc(total);
 	for (size_t i = 0; i < len; ++i) {
 		strcat(result, p[i]);
-		pari_free(p[i]);
+		try_free(p[i]);
 	}
 	return result;
 }
@@ -36,7 +46,7 @@ CHECK(hex_check_param) {
 
 	char *search_hex = try_strdup(args->args);
 	char *p = search_hex;
-	for (; *p; ++p) *p = (char)tolower(*p);
+	for (; *p; ++p) *p = (char) tolower(*p);
 
 	char *params[OFFSET_END] = {NULL};
 	bool pari[OFFSET_END] = {false};
@@ -50,7 +60,7 @@ CHECK(hex_check_param) {
 	if (state >= OFFSET_FIELD) {
 		if (cfg->field == FIELD_PRIME) {
 			params[OFFSET_FIELD] =
-			    pari_sprintf("%P0#*x", cfg->hex_digits, curve->field);
+					pari_sprintf("%P0#*x", cfg->hex_digits, curve->field);
 			pari[OFFSET_FIELD] = true;
 		} else if (cfg->field == FIELD_BINARY) {
 		}
@@ -58,28 +68,42 @@ CHECK(hex_check_param) {
 
 	if (state >= OFFSET_A) {
 		params[OFFSET_A] =
-		    pari_sprintf("%P0#*x", cfg->hex_digits, field_elementi(curve->a));
+				pari_sprintf("%P0#*x", cfg->hex_digits, field_elementi(curve->a));
 		pari[OFFSET_A] = true;
 	}
 
 	if (state >= OFFSET_B) {
 		params[OFFSET_B] =
-		    pari_sprintf("%P0#*x", cfg->hex_digits, field_elementi(curve->b));
+				pari_sprintf("%P0#*x", cfg->hex_digits, field_elementi(curve->b));
 		pari[OFFSET_B] = true;
 	}
 
 	if (state >= OFFSET_ORDER) {
 		params[OFFSET_ORDER] =
-		    pari_sprintf("%P0#*x", cfg->hex_digits, curve->order);
+				pari_sprintf("%P0#*x", cfg->hex_digits, curve->order);
 		pari[OFFSET_ORDER] = true;
 	}
 
 	if (state >= OFFSET_GENERATORS) {
-		params[OFFSET_GENERATORS] = hex_points(curve->generators, curve->ngens);
+		char *subgroups[curve->ngens];
+		for (size_t i = 0; i < curve->ngens; ++i) {
+			subgroups[i] = hex_point(curve->generators[i]->generator);
+		}
+		params[OFFSET_GENERATORS] = str_join(",", subgroups, curve->ngens);
+		for (size_t i = 0; i < curve->ngens; ++i) {
+			try_free(subgroups[i]);
+		}
 	}
 
 	if (state >= OFFSET_POINTS) {
-		params[OFFSET_POINTS] = hex_points(curve->points, curve->npoints);
+		char *subgroups[curve->ngens];
+		for (size_t i = 0; i < curve->ngens; ++i) {
+			subgroups[i] = hex_points(curve->generators[i]->points, curve->generators[i]->npoints);
+		}
+		params[OFFSET_POINTS] = str_join(",", subgroups, curve->ngens);
+		for (size_t i = 0; i < curve->ngens; ++i) {
+			try_free(subgroups[i]);
+		}
 	}
 
 	int result = OFFSET_FIELD - state;
diff --git a/src/gen/order.c b/src/gen/order.c
index 16d597a..0c0779c 100644
--- a/src/gen/order.c
+++ b/src/gen/order.c
@@ -44,39 +44,23 @@ GENERATOR(order_gen_sea) {
 	}
 }
 
-GENERATOR(order_gen_smallfact) {
+GENERATOR(order_gen_cofactor) {
 	HAS_ARG(args);
-	pari_ulong smallfact = *(pari_ulong *)args->args;
+	pari_ulong cofactor = *(pari_ulong *)args->args;
 	pari_sp ltop = avma;
-	GEN fact = mpfact(smallfact);
-	pari_ulong lfact = 0;
-	if (lgefint(fact) > 3) {
-		lfact = 0;
-	} else {
-		lfact = itou(fact);
-	}
-
-	GEN order = ellsea(curve->curve, lfact);
-	if (gequal0(order)) {
+	GEN order = ellff_get_card(curve->curve);
+	GEN res = cgeti(DEFAULTPREC);
+	if (!dvdiiz(order, utoi(cofactor), res)) {
 		avma = ltop;
-		return -4;
+		return  -4;
 	}
-
-	GEN factors = factor(order);
-	GEN primes = gel(factors, 1);
-	GEN powers = gel(factors, 2);
-	long len = glength(primes);
-	GEN total = gen_1;
-	for (long i = 1; i < len; ++i) {
-		GEN pow = powii(gel(primes, i), gel(powers, i));
-		total = mulii(total, pow);
-		if (abscmpiu(total, smallfact) > 0) {
-			avma = ltop;
-			return -4;
-		}
+	if (!isprime(res)) {
+		avma = ltop;
+		return  -4;
 	}
+	verbose_log("cofactor");
 
-	curve->order = gerepileupto(ltop, order);
+	curve->order = gerepilecopy(ltop, order);
 	obj_insert_shallow(curve->curve, 1, curve->order);
 	return 1;
 }
diff --git a/src/gen/order.h b/src/gen/order.h
index 97d31ab..e77937f 100644
--- a/src/gen/order.h
+++ b/src/gen/order.h
@@ -48,11 +48,11 @@ GENERATOR(order_gen_sea);
  * GENERATOR(gen_f)
  *
  * @param curve A curve_t being generated
- * @param args pari_ulong passed to ellsea(curve, smallfact)
+ * @param args pari_ulong the desired cofactor
  * @param state
  * @return state diff
  */
-GENERATOR(order_gen_smallfact);
+GENERATOR(order_gen_cofactor);
 
 /**
  * GENERATOR(gen_f)
diff --git a/src/gen/point.c b/src/gen/point.c
index 4d716e8..1f17d0f 100644
--- a/src/gen/point.c
+++ b/src/gen/point.c
@@ -4,9 +4,9 @@
  */
 #include "point.h"
 #include "exhaustive/arg.h"
-#include "io/output.h"
-#include "math/subgroups.h"
+#include "math/subgroup.h"
 #include "util/memory.h"
+#include "util/random.h"
 
 point_t *point_new(void) { return try_calloc(sizeof(point_t)); }
 
@@ -93,74 +93,85 @@ void points_free_deep(point_t ***points, size_t npoints) {
 }
 
 GENERATOR(point_gen_random) {
-	point_t *p = point_new();
-	p->point = genrand(curve->curve);
-	p->order = ellorder(curve->curve, p->point, NULL);
+	long which_gen = itos(random_range(gen_0, stoi(curve->ngens)));
 
-	curve->points = points_new(1);
-	curve->points[0] = p;
-	curve->npoints = 1;
+	subgroup_t *subgroup = curve->generators[which_gen];
+	GEN mul = random_range(gen_0, subgroup->generator->order);
+	GEN p = ellmul(curve->curve, subgroup->generator->point, mul);
+	point_t *point = point_new();
+	point->point = p;
+	point->order = ellorder(curve->curve, p, NULL);
+	subgroup->npoints = 1;
+	subgroup->points = points_new(1);
+	subgroup->points[0] = point;
 	return 1;
 }
 
 GENERATOR(points_gen_random) {
 	HAS_ARG(args);
 
-	size_t npoints = *(size_t *)args->args;
+	size_t npoints = *(size_t *) args->args;
+	size_t npoints_per_gen[curve->ngens];
+
+	for (size_t i = 0; i < curve->ngens; ++i) {
+		npoints_per_gen[i] = 0;
+	}
 
-	curve->points = points_new(npoints);
-	curve->npoints = npoints;
 	for (size_t i = 0; i < npoints; ++i) {
-		point_t *p = point_new();
-		p->point = genrand(curve->curve);
-		p->order = ellorder(curve->curve, p->point, NULL);
-		curve->points[i] = p;
+		long which_gen = itos(random_range(gen_0, stoi(curve->ngens)));
+		npoints_per_gen[which_gen]++;
 	}
-	return 1;
-}
 
-static int points_from_orders(curve_t *curve, GEN orders) {
-	// TODO better stack code
-	size_t norders = (size_t)glength(orders);
+	for (size_t i = 0; i < curve->ngens; ++i) {
+		subgroup_t *subgroup = curve->generators[i];
+		subgroup->npoints = npoints_per_gen[i];
+		subgroup->points = points_new(npoints_per_gen[i]);
 
-	curve->points = points_new(norders);
-	curve->npoints = norders;
+		for (size_t j = 0; j < npoints_per_gen[i]; ++j) {
+			GEN mul = random_range(gen_0, subgroup->generator->order);
+			GEN p = ellmul(curve->curve, subgroup->generator->point, mul);
+			point_t *point = point_new();
+			point->point = p;
+			point->order = ellorder(curve->curve, p, NULL);
+			subgroup->points[j] = point;
+		}
+	}
+	return 1;
+}
 
-	for (size_t ngen = 0; ngen < curve->ngens; ++ngen) {
-		point_t *gen = curve->generators[ngen];
+static point_t **points_from_orders(GEN curve, point_t *generator, GEN orders) {
+	size_t norders = (size_t) glength(orders);
+	point_t **result = points_new(norders);
 
-		for (long i = 0; i < norders; ++i) {
-			GEN num = gel(orders, i + 1);
-			if (curve->points[i] == NULL) {
-				pari_sp ftop = avma;
+	for (long i = 0; i < norders; ++i) {
+		pari_sp ftop = avma;
+		GEN num = gel(orders, i + 1);
 
-				GEN point = NULL;
-				if (equalii(gen->order, num)) {
-					point = gcopy(gen->point);
-				} else if (dvdii(gen->order, num)) {
-					GEN mul = divii(gen->order, num);
-					point = ellmul(curve->curve, gen->point, mul);
-				}
+		GEN point = NULL;
+		if (equalii(generator->order, num)) {
+			point = gcopy(generator->point);
+		} else if (dvdii(generator->order, num)) {
+			GEN mul = divii(generator->order, num);
+			point = ellmul(curve, generator->point, mul);
+		}
 
-				if (point) {
-					debug_log("VERIFY %Ps %Ps", num,
-					          ellorder(curve->curve, point, NULL));
-					curve->points[i] = point_new();
-					gerepileall(ftop, 1, &point);
-					curve->points[i]->point = point;
-					curve->points[i]->order = gcopy(num);
-				}
-			}
+		if (point) {
+			debug_log("VERIFY %Ps %Ps", num,
+					  ellorder(curve, point, NULL));
+			result[i] = point_new();
+			gerepileall(ftop, 1, &point);
+			result[i]->point = point;
+			result[i]->order = gcopy(num);
 		}
 	}
 
-	return 1;
+	return result;
 }
 
 GENERATOR(points_gen_trial) {
 	HAS_ARG(args);
 
-	pari_ulong *primes = (pari_ulong *)args->args;
+	pari_ulong *primes = (pari_ulong *) args->args;
 	size_t nprimes = args->nargs;
 
 	GEN orders = gtovec0(gen_0, nprimes);
@@ -168,32 +179,78 @@ GENERATOR(points_gen_trial) {
 		gel(orders, i) = utoi(primes[i - 1]);
 	}
 
-	return points_from_orders(curve, orders);
+	GEN orders_per_gen[curve->ngens];
+
+	for (size_t i = 0; i < curve->ngens; ++i) {
+		orders_per_gen[i] = gen_0;
+	}
+
+	for (size_t j = 0; j < nprimes; ++j) {
+		GEN num = gel(orders, j + 1);
+
+		for (size_t i = 0; i < curve->ngens; ++i) {
+			point_t *gen = curve->generators[i]->generator;
+			if (equalii(gen->order, num) || dvdii(gen->order, num)) {
+				if (orders_per_gen[i] == gen_0) {
+					orders_per_gen[i] = gtovec(num);
+				} else {
+					vec_append(orders_per_gen[i], num);
+				}
+				break;
+			}
+		}
+		debug_log("Should not happen.");
+	}
+
+	for (size_t i = 0; i < curve->ngens; ++i) {
+		subgroup_t *subgroup = curve->generators[i];
+		if (orders_per_gen[i] != gen_0) {
+			subgroup->npoints = (size_t) glength(orders_per_gen[i]);
+			subgroup->points = points_from_orders(curve->curve, subgroup->generator, orders_per_gen[i]);
+		}
+	}
+
+	return 1;
 }
 
 GENERATOR(points_gen_prime) {
-	GEN primes = subgroups_prime(curve);
-	return points_from_orders(curve, primes);
+	for (size_t i = 0; i < curve->ngens; ++i) {
+		GEN primes = subgroups_prime(curve->generators[i]->generator->order);
+		curve->generators[i]->npoints = (size_t) glength(primes);
+		curve->generators[i]->points = points_from_orders(curve->curve, curve->generators[i]->generator, primes);
+	}
+
+	return 1;
 }
 
 GENERATOR(points_gen_allgroups) {
-	GEN groups = subgroups_all(curve);
-	return points_from_orders(curve, groups);
+	for (size_t i = 0; i < curve->ngens; ++i) {
+		GEN primes = subgroups_all(curve->generators[i]->generator->order);
+		curve->generators[i]->npoints = (size_t) glength(primes);
+		curve->generators[i]->points = points_from_orders(curve->curve, curve->generators[i]->generator, primes);
+	}
+
+	return 1;
 }
 
 GENERATOR(points_gen_nonprime) {
-	GEN groups = subgroups_nonprime(curve);
-	if (!groups) {
-		// No non-prime order points on curve.
-		return 1;
-	} else {
-		return points_from_orders(curve, groups);
+	for (size_t i = 0; i < curve->ngens; ++i) {
+		GEN primes = subgroups_nonprime(curve->generators[i]->generator->order);
+		if (primes) {
+			curve->generators[i]->npoints = (size_t) glength(primes);
+			curve->generators[i]->points = points_from_orders(curve->curve, curve->generators[i]->generator, primes);
+		}
 	}
+
+	return 1;
 }
 
 UNROLL(points_unroll) {
-	if (curve->points) {
-		points_free_deep(&curve->points, curve->npoints);
+	if (curve->generators) {
+		for (size_t i = 0; i < curve->ngens; ++i) {
+			points_free_deep(&curve->generators[i]->points, curve->generators[i]->npoints);
+		}
+
 	}
 	return -1;
 }
diff --git a/src/gen/point.h b/src/gen/point.h
index c411dee..c97e738 100644
--- a/src/gen/point.h
+++ b/src/gen/point.h
@@ -66,7 +66,7 @@ point_t **points_new(size_t num);
  * @param num
  * @return
  */
-point_t **points_copy(point_t **const src, point_t **dest, size_t num);
+point_t **points_copy(point_t **src, point_t **dest, size_t num);
 
 /**
  *
@@ -74,7 +74,7 @@ point_t **points_copy(point_t **const src, point_t **dest, size_t num);
  * @param num
  * @return
  */
-point_t **points_new_copy(point_t **const src, size_t num);
+point_t **points_new_copy(point_t **src, size_t num);
 
 /**
  *
@@ -83,7 +83,7 @@ point_t **points_new_copy(point_t **const src, size_t num);
  * @param num
  * @return
  */
-point_t **points_clone(point_t **const src, point_t **dest, size_t num);
+point_t **points_clone(point_t **src, point_t **dest, size_t num);
 
 /**
  *
@@ -91,7 +91,7 @@ point_t **points_clone(point_t **const src, point_t **dest, size_t num);
  * @param num
  * @return
  */
-point_t **points_new_clone(point_t **const src, size_t num);
+point_t **points_new_clone(point_t **src, size_t num);
 
 /**
  *
diff --git a/src/io/cli.c b/src/io/cli.c
index c192021..ae11836 100644
--- a/src/io/cli.c
+++ b/src/io/cli.c
@@ -24,7 +24,6 @@ enum opt_keys {
 	OPT_ORDER = 'n',
 	OPT_KOBLITZ = 'K',
 	OPT_UNIQUE = 'u',
-	OPT_FORMAT = 't',
 	OPT_OUTPUT = 'o',
 	OPT_INPUT = 'f',
 	OPT_APPEND = 'a',
@@ -37,8 +36,6 @@ enum opt_keys {
 	OPT_TSTACK,
 	OPT_TIMEOUT,
 	OPT_ANOMALOUS,
-	OPT_GPGEN,
-	OPT_GPCHECK,
 	OPT_HEXCHECK,
 	OPT_BRAINPOOL_RFC,
 	OPT_TWIST,
@@ -62,11 +59,9 @@ struct argp_option cli_options[] = {
 		{0,               0,                 0,        0,                 "Generation options:",                                                                   3},
 		{"random",        OPT_RANDOM,        0,        0,                 "Generate a random curve (using Random approach).",                                      3},
 		{"prime",         OPT_PRIME,         0,        0,                 "Generate a curve with prime order.",                                                    3},
-		{"cofactor",      OPT_COFACTOR,      "BOUND",  0,                 "Generate a curve with cofactor up to BOUND.",                                           3},
+		{"cofactor",      OPT_COFACTOR,      "VALUE",  0,                 "Generate a curve with cofactor of VALUE.",                                              3},
 		{"koblitz",       OPT_KOBLITZ,       "A",    OPTION_ARG_OPTIONAL, "Generate a Koblitz curve (a in {0, 1}, b = 1).",                                        3},
 		{"unique",        OPT_UNIQUE,        0,        0,                 "Generate a curve with only one generator.",                                             3},
-		{"gp-gen",        OPT_GPGEN,         "FUNC",   0,                 "Generate a curve param using a GP function. **NOT IMPLEMENTED**",                       3},
-		{"gp-check",      OPT_GPCHECK,       "FUNC",   0,                 "Check a generated curve param using a GP function. **NOT IMPLEMENTED**",                3},
 		{"hex-check",     OPT_HEXCHECK,      "HEX",    0,                 "Check a generated curve param hex expansion for the HEX string.",                       3},
 		{"points",        OPT_POINTS,        "TYPE",   0,                 "Generate points of given type (random/prime/all/nonprime/none).",                       3},
 		{"count",         OPT_COUNT,         "COUNT",  0,                 "Generate multiple curves.",                                                             3},
@@ -275,7 +270,7 @@ error_t cli_parse(int key, char *arg, struct argp_state *state) {
 			break;
 		case OPT_COFACTOR:
 			cfg->cofactor = true;
-			cfg->cofactor_bound = strtol(arg, NULL, 10);
+			cfg->cofactor_value = strtol(arg, NULL, 10);
 			break;
 		case OPT_KOBLITZ:
 			cfg->koblitz = true;
@@ -290,12 +285,6 @@ error_t cli_parse(int key, char *arg, struct argp_state *state) {
 		case OPT_UNIQUE:
 			cfg->unique = true;
 			break;
-		case OPT_GPGEN:
-			cfg->gp_gens[cfg->gp_gens_size++] = arg;
-			break;
-		case OPT_GPCHECK:
-			cfg->gp_checks[cfg->gp_checks_size++] = arg;
-			break;
 		case OPT_HEXCHECK: {
 			char *str_start = arg;
 			if (strlen(arg) > 2) {
diff --git a/src/io/output.c b/src/io/output.c
index 4d8fca0..6fb82d3 100644
--- a/src/io/output.c
+++ b/src/io/output.c
@@ -19,6 +19,34 @@ char *output_malloc(const char *what) {
 	return s;
 }
 
+
+static JSON_Value *output_json_point(point_t *point) {
+	JSON_Value *point_value = json_value_init_object();
+	JSON_Object *point_object = json_value_get_object(point_value);
+
+	char *x = pari_sprintf(
+			"%P0#*x", cfg->hex_digits,
+			field_elementi(gel(point->point, 1)));
+	json_object_set_string(point_object, "x", x);
+	pari_free(x);
+	char *y = pari_sprintf(
+			"%P0#*x", cfg->hex_digits,
+			field_elementi(gel(point->point, 2)));
+	json_object_set_string(point_object, "y", y);
+	pari_free(y);
+	char *p_order = pari_sprintf("%P#x", point->order);
+	json_object_set_string(point_object, "order", p_order);
+	pari_free(p_order);
+	if (point->cofactor) {
+		char *cofactor =
+				pari_sprintf("%P#x", point->cofactor);
+		json_object_set_string(point_object, "cofactor", cofactor);
+		pari_free(cofactor);
+	}
+
+	return point_value;
+}
+
 static JSON_Value *output_jjson(curve_t *curve) {
 	pari_sp ltop = avma;
 	// root object/value is curve
@@ -77,68 +105,22 @@ static JSON_Value *output_jjson(curve_t *curve) {
 		JSON_Array *gens_array = json_value_get_array(gens_value);
 
 		for (size_t i = 0; i < curve->ngens; ++i) {
-			JSON_Value *point_value = json_value_init_object();
-			JSON_Object *point_object = json_value_get_object(point_value);
+			JSON_Value *gen_value = output_json_point(curve->generators[i]->generator);
+			JSON_Object *gen_object = json_value_get_object(gen_value);
 
-			char *x = pari_sprintf(
-			    "%P0#*x", cfg->hex_digits,
-			    field_elementi(gel(curve->generators[i]->point, 1)));
-			json_object_set_string(point_object, "x", x);
-			pari_free(x);
-			char *y = pari_sprintf(
-			    "%P0#*x", cfg->hex_digits,
-			    field_elementi(gel(curve->generators[i]->point, 2)));
-			json_object_set_string(point_object, "y", y);
-			pari_free(y);
-			char *p_order = pari_sprintf("%P#x", curve->generators[i]->order);
-			json_object_set_string(point_object, "order", p_order);
-			pari_free(p_order);
-			if (curve->generators[i]->cofactor) {
-				char *cofactor =
-				    pari_sprintf("%P#x", curve->generators[i]->cofactor);
-				json_object_set_string(point_object, "cofactor", cofactor);
-				pari_free(cofactor);
+			if (curve->generators[i]->npoints) {
+				JSON_Value *gens_points_value = json_value_init_array();
+				JSON_Array *gens_points_array = json_value_get_array(gens_points_value);
+				for (size_t j = 0; j < curve->generators[i]->npoints; ++j) {
+					json_array_append_value(gens_points_array, output_json_point(curve->generators[i]->points[j]));
+				}
+				json_object_set_value(gen_object, "points", gens_points_value);
 			}
-
-			json_array_append_value(gens_array, point_value);
+			json_array_append_value(gens_array, gen_value);
 		}
-
-		json_object_set_value(root_object, "generators", gens_value);
+		json_object_set_value(root_object, "subgroups", gens_value);
 	}
 
-	if (curve->npoints) {
-		JSON_Value *points_value = json_value_init_array();
-		JSON_Array *points_array = json_value_get_array(points_value);
-
-		for (size_t i = 0; i < curve->npoints; ++i) {
-			JSON_Value *point_value = json_value_init_object();
-			JSON_Object *point_object = json_value_get_object(point_value);
-
-			char *x =
-			    pari_sprintf("%P0#*x", cfg->hex_digits,
-			                 field_elementi(gel(curve->points[i]->point, 1)));
-			json_object_set_string(point_object, "x", x);
-			pari_free(x);
-			char *y =
-			    pari_sprintf("%P0#*x", cfg->hex_digits,
-			                 field_elementi(gel(curve->points[i]->point, 2)));
-			json_object_set_string(point_object, "y", y);
-			pari_free(y);
-			char *p_order = pari_sprintf("%P#x", curve->points[i]->order);
-			json_object_set_string(point_object, "order", p_order);
-			pari_free(p_order);
-			if (curve->points[i]->cofactor) {
-				char *cofactor =
-				    pari_sprintf("%P#x", curve->points[i]->cofactor);
-				json_object_set_string(point_object, "cofactor", cofactor);
-				pari_free(cofactor);
-			}
-
-			json_array_append_value(points_array, point_value);
-		}
-
-		json_object_set_value(root_object, "points", points_value);
-	}
 	avma = ltop;
 	return root_value;
 }
diff --git a/src/math/subgroup.c b/src/math/subgroup.c
new file mode 100644
index 0000000..3d78db5
--- /dev/null
+++ b/src/math/subgroup.c
@@ -0,0 +1,237 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017-2018 J08nY
+ */
+#include "gen/point.h"
+#include "subgroup.h"
+#include "util/memory.h"
+
+subgroup_t *subgroup_new(void) {
+	return try_calloc(sizeof(subgroup_t));
+}
+
+subgroup_t *subgroup_copy(const subgroup_t *src, subgroup_t *dst) {
+	if (src->generator)
+		dst->generator = point_new_copy(src->generator);
+	if (src->points) {
+		dst->points = points_new_copy(src->points, src->npoints);
+		dst->npoints = src->npoints;
+	}
+	return dst;
+}
+
+subgroup_t *subgroup_new_copy(const subgroup_t *src) {
+	subgroup_t *result = subgroup_new();
+	return subgroup_copy(src, result);
+}
+
+subgroup_t *subgroup_clone(const subgroup_t *src, subgroup_t *dst) {
+	if (src->generator)
+		dst->generator = point_new_clone(src->generator);
+	if (src->points) {
+		dst->points = points_new_clone(src->points, src->npoints);
+		dst->npoints = src->npoints;
+	}
+	return dst;
+}
+
+subgroup_t *subgroup_new_clone(const subgroup_t *src) {
+	subgroup_t *result = subgroup_new();
+	return subgroup_clone(src, result);
+}
+
+void subgroup_free(subgroup_t **subgroup) {
+	if (*subgroup) {
+		if ((*subgroup)->generator) {
+			point_free(&(*subgroup)->generator);
+		}
+		try_free(*subgroup);
+		*subgroup = NULL;
+	}
+}
+
+void subgroup_free_deep(subgroup_t **subgroup) {
+	if (*subgroup) {
+		points_free_deep(&(*subgroup)->points, (*subgroup)->npoints);
+		subgroup_free(subgroup);
+	}
+}
+
+subgroup_t **subgroups_new(size_t num) {
+	return try_calloc(num * sizeof(subgroup_t *));
+}
+
+subgroup_t **subgroups_copy(subgroup_t **const src, subgroup_t **dest, size_t num) {
+	for (size_t i = 0; i < num; ++i) {
+		dest[i] = subgroup_new_copy(src[i]);
+	}
+	return dest;
+}
+
+subgroup_t **subgroups_new_copy(subgroup_t **const src, size_t num) {
+	subgroup_t **result = subgroups_new(num);
+	return subgroups_copy(src, result, num);
+}
+
+subgroup_t **subgroups_clone(subgroup_t **const src, subgroup_t **dest, size_t num) {
+	for (size_t i = 0; i < num; ++i) {
+		dest[i] = subgroup_new_clone(src[i]);
+	}
+	return dest;
+}
+
+subgroup_t **subgroups_new_clone(subgroup_t **const src, size_t num) {
+	subgroup_t **result = subgroups_new(num);
+	return subgroups_clone(src, result, num);
+}
+
+void subgroups_free(subgroup_t ***subgroups) {
+	if (*subgroups) {
+		try_free(*subgroups);
+		*subgroups = NULL;
+	}
+}
+
+void subgroups_free_deep(subgroup_t ***subgroups, size_t num) {
+	if (*subgroups) {
+		for (size_t i = 0; i < num; ++i) {
+			subgroup_free(&(*subgroups)[i]);
+		}
+		subgroups_free(subgroups);
+	}
+}
+
+/**
+ * @brief All prime divisors of a given integer with multiplicity.
+ *
+ * subgroups_divisors(27) = [3, 3, 3]
+ * @param order
+ * @return a t_VEC of prime divisors.
+ */
+static GEN subgroups_divisors(GEN order) {
+	GEN factors = Z_factor(order);
+	GEN primes = gel(factors, 1);
+	GEN multiples = gel(factors, 2);
+	long uniqs = glength(primes);
+
+	long size = 0;
+	for (long i = 1; i <= uniqs; ++i) {
+		size += itos(gel(multiples, i));
+	}
+	GEN result = gtovec0(gen_0, size);
+
+	long count = 0;
+	for (long i = 1; i <= uniqs; ++i) {
+		long multiple = itos(gel(multiples, i));
+		for (long j = 1; j <= multiple; ++j) {
+			gel(result, ++count) = gel(primes, i);
+		}
+	}
+	return result;
+}
+
+/**
+ * @brief All factors consisting of at least <code>min_bits</code> prime
+ * <code>factors</code>.
+ *
+ * @param factors
+ * @param min_bits
+ * @return a t_VEC of factors
+ */
+static GEN subgroups_2n_factors(GEN factors, size_t min_bits) {
+	pari_sp ltop = avma;
+	long nprimes = glength(factors);
+	if (nprimes == min_bits) return NULL;
+	GEN amount = int2n(nprimes);
+	GEN groups = gtovec0(gen_0, itos(amount) - (min_bits * nprimes) - 1);
+
+	size_t i = 0;
+	for (size_t count = 1; count < (size_t) (1) << nprimes; ++count) {
+		pari_sp btop = avma;
+		GEN result = gen_1;
+		size_t bits = 0;
+		for (long bit = 0; bit < nprimes; ++bit) {
+			size_t mask = (size_t) (1) << bit;
+			if (count & mask) {
+				result = mulii(result, gel(factors, bit + 1));
+				bits++;
+			}
+		}
+		if (bits > min_bits) {
+			gel(groups, ++i) = result;
+		} else {
+			avma = btop;
+		}
+	}
+	GEN ret = gtoset(groups);
+	return gerepilecopy(ltop, ret);
+}
+
+GEN subgroups_prime(GEN order) {
+	GEN factors = Z_factor(order);
+	return gtovec(gel(factors, 1));
+}
+
+GEN subgroups_nonprime(GEN order) {
+	if (isprime(order)) {
+		return NULL;
+	}
+	GEN factors = subgroups_divisors(order);
+	return subgroups_2n_factors(factors, 1);
+}
+
+GEN subgroups_all(GEN order) {
+	if (isprime(order)) {
+		return gtovec(order);
+	}
+	GEN factors = subgroups_divisors(order);
+	return subgroups_2n_factors(factors, 0);
+}
+
+/**
+ * @brief
+ * @param curve
+ * @param min_bits
+ * @return
+ */
+/*
+static GEN subgroups_2n_gens(const curve_t *curve, size_t min_bits) {
+	GEN one_factors = subgroups_divisors(curve->generators[0]->order);
+	GEN one = subgroups_2n_factors(one_factors, min_bits);
+	GEN other_factors = subgroups_divisors(curve->generators[1]->order);
+	GEN other = subgroups_2n_factors(other_factors, min_bits);
+	if (!one) {
+		return other;
+	}
+	if (!other) {
+		return one;
+	}
+	GEN result = gtovec0(gen_0, glength(one) + glength(other));
+	for (long i = 1; i <= glength(result); ++i) {
+		if (i <= glength(one)) {
+			gel(result, i) = gel(one, i);
+		} else {
+			gel(result, i) = gel(other, i - glength(one));
+		}
+	}
+	return result;
+}
+ */
+
+/**
+ * @brief
+ * @param curve
+ * @param min_bits
+ * @return
+ */
+/*
+static GEN subgroups_2n(const curve_t *curve, size_t min_bits) {
+   if (curve->ngens == 1) {
+	   GEN factors = subgroups_divisors(curve->order);
+	   return subgroups_2n_factors(factors, min_bits);
+   }
+
+   return subgroups_2n_gens(curve, min_bits);
+}
+
+*/
\ No newline at end of file
diff --git a/src/math/subgroup.h b/src/math/subgroup.h
new file mode 100644
index 0000000..ad36826
--- /dev/null
+++ b/src/math/subgroup.h
@@ -0,0 +1,139 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017-2018 J08nY
+ */
+/**
+ * @file subgroup.h
+ */
+#ifndef ECGEN_SUBGROUPS_H
+#define ECGEN_SUBGROUPS_H
+
+#include <pari/pari.h>
+#include "misc/types.h"
+
+/**
+ * @brief
+ * @return
+ */
+subgroup_t *subgroup_new(void);
+
+/**
+ * @brief
+ * @param src
+ * @param dst
+ * @return
+ */
+subgroup_t *subgroup_copy(const subgroup_t *src, subgroup_t *dst);
+
+/**
+ * @brief
+ * @param src
+ * @return
+ */
+subgroup_t *subgroup_new_copy(const subgroup_t *src);
+
+/**
+ * @brief
+ * @param src
+ * @param dst
+ * @return
+ */
+subgroup_t *subgroup_clone(const subgroup_t *src, subgroup_t *dst);
+
+/**
+ * @brief
+ * @param src
+ * @return
+ */
+subgroup_t *subgroup_new_clone(const subgroup_t *src);
+
+/**
+ * @brief
+ * @param subgroup
+ */
+void subgroup_free(subgroup_t **subgroup);
+
+/**
+ * @brief
+ * @param subgroup
+ */
+void subgroup_free_deep(subgroup_t **subgroup);
+
+/**
+ * @brief
+ * @param num
+ * @return
+ */
+subgroup_t **subgroups_new(size_t num);
+
+/**
+ *
+ * @param src
+ * @param dest
+ * @param num
+ * @return
+ */
+subgroup_t **subgroups_copy(subgroup_t **src, subgroup_t **dest, size_t num);
+
+/**
+ *
+ * @param src
+ * @param num
+ * @return
+ */
+subgroup_t **subgroups_new_copy(subgroup_t **src, size_t num);
+
+/**
+ *
+ * @param src
+ * @param dest
+ * @param num
+ * @return
+ */
+subgroup_t **subgroups_clone(subgroup_t **src, subgroup_t **dest, size_t num);
+
+/**
+ *
+ * @param src
+ * @param num
+ * @return
+ */
+subgroup_t **subgroups_new_clone(subgroup_t **src, size_t num);
+
+/**
+ * @brief
+ * @param subgroups
+ */
+void subgroups_free(subgroup_t ***subgroups);
+
+/**
+ * @brief
+ * @param subgroups
+ * @param num
+ */
+void subgroups_free_deep(subgroup_t ***subgroups, size_t num);
+
+/**
+ * @brief All prime factors of a given integer, without multipliticity.
+ *
+ * subgroups_factors(27) = [3]
+ * @param order
+ * @return a t_VEC of prime factors.
+ */
+GEN subgroups_prime(GEN order);
+
+/**
+ * @brief All nonprime subgroup orders of a given integer.
+ * @param order
+ * @return a t_VEC of nonprime factors.
+ */
+GEN subgroups_nonprime(GEN order);
+
+/**
+ * @brief All all subgroup orders of a given integer.
+ * @param order
+ * @return a t_VEC of all factors.
+ */
+GEN subgroups_all(GEN order);
+
+#endif  // ECGEN_SUBGROUPS_H
diff --git a/src/math/subgroups.c b/src/math/subgroups.c
deleted file mode 100644
index 8fe3918..0000000
--- a/src/math/subgroups.c
+++ /dev/null
@@ -1,150 +0,0 @@
-/*
- * ecgen, tool for generating Elliptic curve domain parameters
- * Copyright (C) 2017-2018 J08nY
- */
-#include "subgroups.h"
-
-/**
- * @brief All prime factors of a given integer.
- *
- * subgroups_factors(27) = [3]
- * @param order
- * @return a t_VEC of prime factors.
- */
-static GEN subgroups_factors(GEN order) {
-	GEN factors = Z_factor(order);
-	return gtovec(gel(factors, 1));
-}
-
-/**
- * @brief All prime divisors of a given integer with multiplicity.
- *
- * subgroups_divisors(27) = [3, 3, 3]
- * @param order
- * @return a t_VEC of prime divisors.
- */
-static GEN subgroups_divisors(GEN order) {
-	GEN factors = Z_factor(order);
-	GEN primes = gel(factors, 1);
-	GEN multiples = gel(factors, 2);
-	long uniqs = glength(primes);
-
-	long size = 0;
-	for (long i = 1; i <= uniqs; ++i) {
-		size += itos(gel(multiples, i));
-	}
-	GEN result = gtovec0(gen_0, size);
-
-	long count = 0;
-	for (long i = 1; i <= uniqs; ++i) {
-		long multiple = itos(gel(multiples, i));
-		for (long j = 1; j <= multiple; ++j) {
-			gel(result, ++count) = gel(primes, i);
-		}
-	}
-	return result;
-}
-
-/**
- * @brief All factors consisting of at least <code>min_bits</code> prime
- * <code>factors</code>.
- *
- * @param factors
- * @param min_bits
- * @return a t_VEC of factors
- */
-static GEN subgroups_2n_factors(GEN factors, size_t min_bits) {
-	pari_sp ltop = avma;
-	long nprimes = glength(factors);
-	if (nprimes == min_bits) return NULL;
-	GEN amount = int2n(nprimes);
-	GEN groups = gtovec0(gen_0, itos(amount) - (min_bits * nprimes) - 1);
-
-	size_t i = 0;
-	for (size_t count = 1; count < (size_t)(1) << nprimes; ++count) {
-		pari_sp btop = avma;
-		GEN result = gen_1;
-		size_t bits = 0;
-		for (long bit = 0; bit < nprimes; ++bit) {
-			size_t mask = (size_t)(1) << bit;
-			if (count & mask) {
-				result = mulii(result, gel(factors, bit + 1));
-				bits++;
-			}
-		}
-		if (bits > min_bits) {
-			gel(groups, ++i) = result;
-		} else {
-			avma = btop;
-		}
-	}
-	GEN ret = gtoset(groups);
-	return gerepilecopy(ltop, ret);
-}
-
-/**
- * @brief
- * @param curve
- * @param min_bits
- * @return
- */
-static GEN subgroups_2n_gens(const curve_t *curve, size_t min_bits) {
-	GEN one_factors = subgroups_divisors(curve->generators[0]->order);
-	GEN one = subgroups_2n_factors(one_factors, min_bits);
-	GEN other_factors = subgroups_divisors(curve->generators[1]->order);
-	GEN other = subgroups_2n_factors(other_factors, min_bits);
-	if (!one) {
-		return other;
-	}
-	if (!other) {
-		return one;
-	}
-	GEN result = gtovec0(gen_0, glength(one) + glength(other));
-	for (long i = 1; i <= glength(result); ++i) {
-		if (i <= glength(one)) {
-			gel(result, i) = gel(one, i);
-		} else {
-			gel(result, i) = gel(other, i - glength(one));
-		}
-	}
-	return result;
-}
-
-/**
- * @brief
- * @param curve
- * @param min_bits
- * @return
- */
-static GEN subgroups_2n(const curve_t *curve, size_t min_bits) {
-	if (curve->ngens == 1) {
-		GEN factors = subgroups_divisors(curve->order);
-		return subgroups_2n_factors(factors, min_bits);
-	}
-
-	return subgroups_2n_gens(curve, min_bits);
-}
-
-GEN subgroups_prime(const curve_t *curve) {
-	if (cfg->prime || isprime(curve->order)) {
-		return gtovec(curve->order);
-	}
-
-	return subgroups_factors(curve->order);
-}
-
-GEN subgroups_nonprime(const curve_t *curve) {
-	if (cfg->prime || isprime(curve->order)) {
-		return NULL;
-	}
-
-	return subgroups_2n(curve, 1);
-}
-
-GEN subgroups_all(const curve_t *curve) {
-	if (cfg->prime || isprime(curve->order)) {
-		return gtovec(curve->order);
-	}
-
-	return subgroups_2n(curve, 0);
-}
diff --git a/src/math/subgroups.h b/src/math/subgroups.h
deleted file mode 100644
index e08bfa5..0000000
--- a/src/math/subgroups.h
+++ /dev/null
@@ -1,35 +0,0 @@
-/*
- * ecgen, tool for generating Elliptic curve domain parameters
- * Copyright (C) 2017-2018 J08nY
- */
-/**
- * @file subgroups.h
- */
-#ifndef ECGEN_SUBGROUPS_H
-#define ECGEN_SUBGROUPS_H
-
-#include <pari/pari.h>
-#include "misc/types.h"
-
-/**
- * @brief Enumerates prime subgroup orders of a given curve.
- * @param curve
- * @return
- */
-GEN subgroups_prime(const curve_t *curve);
-
-/**
- * @brief Enumerates nonprime subgroup orders of a given curve.
- * @param curve
- * @return
- */
-GEN subgroups_nonprime(const curve_t *curve);
-
-/**
- * @brief Enumerates all subgroup orders of a given curve.
- * @param curve
- * @return
- */
-GEN subgroups_all(const curve_t *curve);
-
-#endif  // ECGEN_SUBGROUPS_H
diff --git a/src/math/twists.c b/src/math/twists.c
index 043594f..a6bb06d 100644
--- a/src/math/twists.c
+++ b/src/math/twists.c
@@ -3,7 +3,7 @@
  * Copyright (C) 2017-2018 J08nY
  */
 #include "twists.h"
-#include "gen/point.h"
+#include "math/subgroup.h"
 #include "gen/seed.h"
 
 void twist_rand_to(curve_t *to, const curve_t *of) {
@@ -32,6 +32,5 @@ void twist_rand_to(curve_t *to, const curve_t *of) {
 void twist_rand(curve_t *what) {
 	twist_rand_to(what, what);
 	seed_free(&what->seed);
-	points_free_deep(&what->points, what->npoints);
-	points_free_deep(&what->generators, what->ngens);
+	subgroups_free_deep(&what->generators, what->ngens);
 }
\ No newline at end of file
diff --git a/src/misc/config.h b/src/misc/config.h
index 9a8f8af..ab5b0ba 100644
--- a/src/misc/config.h
+++ b/src/misc/config.h
@@ -71,7 +71,7 @@ typedef struct {
 	long koblitz_value;
 	/** @brief Whether the curves should have a bound on the cofactor value. */
 	bool cofactor;
-	long cofactor_bound;
+	long cofactor_value;
 	/** @brief What seed algorithm, if any, to use to generate the curves. */
 	seed_e seed_algo;
 	/** @brief What seed to use, if any, to generate the curves. */
diff --git a/src/misc/types.h b/src/misc/types.h
index 1e74e04..f2e5ae2 100644
--- a/src/misc/types.h
+++ b/src/misc/types.h
@@ -67,6 +67,18 @@ typedef struct {
 } point_t;
 
 /**
+ * @brief A subgroup type.
+ * @param generator a point generating the subgroup
+ * @param npoints number of points stored in the subgroup
+ * @param points the stored points
+ */
+typedef struct {
+	point_t *generator;
+	size_t npoints;
+	point_t **points;
+} subgroup_t;
+
+/**
  * @brief A curve type.
  * @param seed a seed_t
  * @param field a t_INT or t_FFELT
@@ -75,9 +87,7 @@ typedef struct {
  * @param curve a t_ELL, curve object
  * @param order a t_INT, curve order
  * @param generators generators saved
- * @param ngens   numver of generators saved in the curve type
- * @param points  points saved
- * @param npoints number of points saved in the curve type
+ * @param ngens number of generators saved in the curve type
  */
 typedef struct {
 	seed_t *seed;
@@ -86,10 +96,8 @@ typedef struct {
 	GEN b;
 	GEN curve;
 	GEN order;
-	point_t **generators;
+	subgroup_t **generators;
 	size_t ngens;
-	point_t **points;
-	size_t npoints;
 } curve_t;
 
 /**
@@ -126,7 +134,7 @@ typedef struct {
  * @return state diff
  */
 #define GENERATOR(gen_name) \
-	int gen_name(curve_t *curve, arg_t *args, offset_e state)
+    int gen_name(curve_t *curve, arg_t *args, offset_e state)
 
 typedef GENERATOR((*gen_f));
 
@@ -138,7 +146,7 @@ typedef GENERATOR((*gen_f));
  * @return
  */
 #define UNROLL(unroll_name) \
-	int unroll_name(curve_t *curve, pari_sp from, pari_sp to)
+    int unroll_name(curve_t *curve, pari_sp from, pari_sp to)
 
 typedef UNROLL((*unroll_f));
 
diff --git a/src/util/random.c b/src/util/random.c
index 860cff0..eaeceaf 100644
--- a/src/util/random.c
+++ b/src/util/random.c
@@ -62,6 +62,16 @@ GEN random_int(unsigned long bits) {
 	return gerepilecopy(ltop, genrand(range));
 }
 
+GEN random_range(GEN lower, GEN upper) {
+	pari_sp ltop = avma;
+	if (gequal(lower, upper)) {
+		return gcopy(lower);
+	}
+
+	GEN range = mkvec2(lower, subis(upper, 1));
+	return gerepilecopy(ltop, genrand(range));
+}
+
 GEN random_field_element(GEN field) {
 	switch (typ(field)) {
 		case t_INT:
diff --git a/src/util/random.h b/src/util/random.h
index 56e2efd..f5844da 100644
--- a/src/util/random.h
+++ b/src/util/random.h
@@ -30,6 +30,14 @@ bool random_init(void);
 GEN random_prime(unsigned long bits);
 
 /**
+ * @brief
+ * @param lower
+ * @param upper
+ * @return
+ */
+GEN random_range(GEN lower, GEN upper);
+
+/**
  * @brief Generate random <code>bits</code> sized integer.
  * @param bits the size of the integer to generate
  * @return a random integer in range [2^(bits - 1), 2^bits]
diff --git a/src/util/str.c b/src/util/str.c
index f213c99..8b1bf92 100644
--- a/src/util/str.c
+++ b/src/util/str.c
@@ -5,6 +5,7 @@
 #include "str.h"
 #include <ctype.h>
 #include <string.h>
+#include <stdarg.h>
 #include "util/memory.h"
 
 const char *str_is_hex(const char *hex_str) {
@@ -22,16 +23,47 @@ const char *str_is_hex(const char *hex_str) {
 	return str_start;
 }
 
-char *str_join(char *strings[], size_t len) {
-	size_t total = 0;
+char *str_join(char *separator, char **strings, size_t len) {
+	size_t total = 1;
 	for (size_t i = 0; i < len; ++i) {
 		if (strings[i]) total += strlen(strings[i]);
 	}
+	if (separator)
+		total += (len - 1) * strlen(separator);
 	char *result = try_calloc(total);
 	for (size_t i = 0; i < len; ++i) {
 		if (strings[i]) {
 			strcat(result, strings[i]);
 		}
+		if (separator && i != len - 1) {
+			strcat(result, separator);
+		}
+	}
+	return result;
+}
+
+char *str_joinv(char *separator, ...) {
+	va_list valist;
+	va_start(valist, separator);
+
+	size_t len = 5;
+	char **strings = try_malloc(len * sizeof(char *));
+	char *next;
+	size_t i = 0;
+	while ((next = va_arg(valist, char *)) != NULL) {
+		if (i == len) {
+			strings = try_realloc(strings, (len * 2) * sizeof(char *));
+		}
+		strings[i++] = next;
 	}
+
+	va_end(valist);
+	char *result = str_join(separator, strings, i);
+	try_free(strings);
 	return result;
-}
\ No newline at end of file
+}
+
+char *str_concat(char **strings, size_t len) {
+	return str_join(NULL, strings, len);
+}
+
diff --git a/src/util/str.h b/src/util/str.h
index 8fd2540..fb828da 100644
--- a/src/util/str.h
+++ b/src/util/str.h
@@ -19,10 +19,27 @@ const char *str_is_hex(const char *hex_str);
 
 /**
  * @brief
+ * @param separator
  * @param strings
  * @param len
  * @return
  */
-char *str_join(char *strings[], size_t len);
+char *str_join(char *separator, char **strings, size_t len);
+
+/**
+ * @brief
+ * @param separator
+ * @param ...
+ * @return
+ */
+char *str_joinv(char *separator, ...);
+
+/**
+ * @brief
+ * @param strings
+ * @param len
+ * @return
+ */
+char *str_concat(char **strings, size_t len);
 
 #endif  // ECGEN_STR_H