34 files changed, 1008 insertions, 1078 deletions

diff --git a/README.md b/README.md
index a94bb9b..d63b733 100644
--- a/README.md
+++ b/README.md
@@ -23,7 +23,7 @@ Tool for generating Elliptic curve domain parameters.
 #### Generation options
 
  - `-c / --count=COUNT`		Generate multiple curves.
- - `-k / --cofactor=BOUND`	Generate a curve with cofactor up to `BOUND`.
+ - `-k / --cofactor=VALUE`	Generate a curve with cofactor of `VALUE`.
  - `-K / --koblitz[=A]`		Generate a Koblitz curve (a in {0, 1}, b = 1).
  - `-p / --prime`			Generate a curve with prime order.
  - `--points=TYPE`			Generate points of given `TYPE` (random/prime/all/nonprime/none).
@@ -73,19 +73,19 @@ Generate a prime field, uniquely generated, prime order curve, don't ask for inp
         "a": "0x9c083973bdca36ea71078bbaabab4947",
         "b": "0x3d986a0206bfbe1ba62c858df54385e9",
         "order": "0xa5393890f26881d9394aece3bc2d9b47",
-        "generators": [
+        "subgroups": [
             {
                 "x": "0x5acc17d6a44e8f8d30e877f4fef8712f",
                 "y": "0x6864dd64e80609abd1797c8de1febb9f",
                 "order": "0xa5393890f26881d9394aece3bc2d9b47",
-                "cofactor": "0x1"
-            }
-        ],
-        "points": [
-            {
-                "x": "0x9c7878930ddf5bfb705102f652754e7",
-                "y": "0x4b15a7bb808cb3579fd4c2ce42f628de",
-                "order": "0xa5393890f26881d9394aece3bc2d9b47"
+                "cofactor": "0x1",
+                "points": [
+                    {
+                        "x": "0x9c7878930ddf5bfb705102f652754e7",
+                        "y": "0x4b15a7bb808cb3579fd4c2ce42f628de",
+                        "order": "0xa5393890f26881d9394aece3bc2d9b47"
+                    }
+                ]
             }
         ]
     }
diff --git a/docs/output.md b/docs/output.md
index ec18419..8ec18d8 100644
--- a/docs/output.md
+++ b/docs/output.md
@@ -1,6 +1,6 @@
 # Output
 
-There are only output format currently supported in ecgen is JSON.
+There only output format currently supported in ecgen is JSON.
 
 ## JSON
 
@@ -16,34 +16,34 @@ Self-explanatory format. The curve dictionaries are enclosed in an array as you
 	    "a": "0x223fcc1306c21406",
 	    "b": "0x4114b86128071651",
 	    "order": "0xfa4101a5b65111b6",
-	    "generators": [
+	    "subgroups": [
 	        {
 	            "x": "0x41b71e83794c614a",
 	            "y": "0xc43c15e114f16ba1",
 	            "order": "0xfa4101a5b65111b6",
-	            "cofactor": "0x1"
-	        }
-	    ],
-	    "points": [
-	        {
-	            "x": "0x7427f55c615d0c60",
-	            "y": "0x0000000000000000",
-	            "order": "0x2"
-	        },
-	        {
-	            "x": "0x357cbff05dacc66c",
-	            "y": "0x4f5d5e523a38a35a",
-	            "order": "0xd"
-	        },
-	        {
-	            "x": "0x4ef832dbb406dac6",
-	            "y": "0x684eb7a227fc23c3",
-	            "order": "0x8b60007"
-	        },
-	        {
-	            "x": "0xd82bdb55db6bb3ef",
-	            "y": "0xc103b77986a1c2e3",
-	            "order": "0x11ade0181"
+	            "cofactor": "0x1",
+                "points": [
+                    {
+                        "x": "0x7427f55c615d0c60",
+                        "y": "0x0000000000000000",
+                        "order": "0x2"
+                    },
+                    {
+                        "x": "0x357cbff05dacc66c",
+                        "y": "0x4f5d5e523a38a35a",
+                        "order": "0xd"
+                    },
+                    {
+                        "x": "0x4ef832dbb406dac6",
+                        "y": "0x684eb7a227fc23c3",
+                        "order": "0x8b60007"
+                    },
+                    {
+                        "x": "0xd82bdb55db6bb3ef",
+                        "y": "0xc103b77986a1c2e3",
+                        "order": "0x11ade0181"
+                    }
+                ]
 	        }
 	    ]
 	}]
@@ -61,44 +61,44 @@ Self-explanatory format. The curve dictionaries are enclosed in an array as you
 	    "a": "0x3869c3f4bbeb4087",
 	    "b": "0xba60557c283d94cd",
 	    "order": "0xffffffffea115f5a",
-	    "generators": [
+	    "subgroups": [
 	        {
 	            "x": "0x4dcc3b35abcc13e4",
 	            "y": "0x4380aa919232c5b1",
 	            "order": "0xffffffffea115f5a",
-	            "cofactor": "0x1"
-	        }
-	    ],
-	    "points": [
-	        {
-	            "x": "0x0000000000000000",
-	            "y": "0x0000000000000000",
-	            "order": "0x2"
-	        },
-	        {
-	            "x": "0xde23b127982c9db7",
-	            "y": "0x1c63aa4b52327ca0",
-	            "order": "0x13"
-	        },
-	        {
-	            "x": "0xcb00f2cf13a8bab9",
-	            "y": "0x192a159ad11df21b",
-	            "order": "0xf1"
-	        },
-	        {
-	            "x": "0x3b9cd618b48e9a73",
-	            "y": "0xdacc98d7c32d8c0b",
-	            "order": "0x649"
-	        },
-	        {
-	            "x": "0xdbdf22552a10abbe",
-	            "y": "0x0dd6e81435a11566",
-	            "order": "0x2f94b"
-	        },
-	        {
-	            "x": "0x963ae5e637a9917a",
-	            "y": "0x1974294368f7f6d1",
-	            "order": "0x6203c5"
+	            "cofactor": "0x1",
+	            "points": [
+                    {
+                        "x": "0x0000000000000000",
+                        "y": "0x0000000000000000",
+                        "order": "0x2"
+                    },
+                    {
+                        "x": "0xde23b127982c9db7",
+                        "y": "0x1c63aa4b52327ca0",
+                        "order": "0x13"
+                    },
+                    {
+                        "x": "0xcb00f2cf13a8bab9",
+                        "y": "0x192a159ad11df21b",
+                        "order": "0xf1"
+                    },
+                    {
+                        "x": "0x3b9cd618b48e9a73",
+                        "y": "0xdacc98d7c32d8c0b",
+                        "order": "0x649"
+                    },
+                    {
+                        "x": "0xdbdf22552a10abbe",
+                        "y": "0x0dd6e81435a11566",
+                        "order": "0x2f94b"
+                    },
+                    {
+                        "x": "0x963ae5e637a9917a",
+                        "y": "0x1974294368f7f6d1",
+                        "order": "0x6203c5"
+                    }
+                ]
 	        }
 	    ]
 	}]
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
diff --git a/test/src/gen/test_gens.c b/test/src/gen/test_gens.c
index 49f82e2..db11e3d 100644
--- a/test/src/gen/test_gens.c
+++ b/test/src/gen/test_gens.c
@@ -79,7 +79,8 @@ Test(gens, test_gens_check_anomalous) {
 	};
 	point_t gen = {.point = mkvec2(mkintmodu(4, 19), mkintmodu(14, 19)),
 	               .order = stoi(16)};
-	point_t *generators[1] = {&gen};
+	subgroup_t gen_subgroup = {.generator = &gen};
+	subgroup_t *generators[1] = {&gen_subgroup};
 	curve.generators = generators;
 	curve.ngens = 1;
 
@@ -116,7 +117,9 @@ Test(gens, test_gens_check_embedding) {
 	                .order = stoi(6)};
 	point_t gen2 = {.point = mkvec2(mkintmodu(10, 19), mkintmodu(4, 19)),
 	                .order = stoi(3)};
-	point_t *generators[2] = {&gen1, &gen2};
+	subgroup_t gen1_subgroup = {.generator = &gen1};
+	subgroup_t gen2_subgroup = {.generator = &gen2};
+	subgroup_t *generators[2] = {&gen1_subgroup, &gen2_subgroup};
 	curve.generators = generators;
 	curve.ngens = 2;
 
diff --git a/test/src/gen/test_gp.c b/test/src/gen/test_gp.c
deleted file mode 100644
index 9e3e1a9..0000000
--- a/test/src/gen/test_gp.c
+++ /dev/null
@@ -1,249 +0,0 @@
-/*
- * ecgen, tool for generating Elliptic curve domain parameters
- * Copyright (C) 2017 J08nY
- */
-#include <criterion/criterion.h>
-#include "gen/gp.h"
-#include "test/default.h"
-#include "util/bits.h"
-
-TestSuite(gp, .init = default_setup, .fini = default_teardown);
-
-Test(gp, test_gp_gen_seed) {
-	curve_t curve = {0};
-	arg_t arg = {.args = "() -> { return(Vecsmall([1,0])); }", .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_SEED);
-	cr_assert_eq(ret, 1, );
-	cr_assert_not_null(curve.seed, );
-	cr_assert_not_null(curve.seed->seed, );
-	cr_assert(bits_eq(curve.seed->seed, bits_from_bin("10")), );
-}
-
-Test(gp, test_gp_gen_field) {
-	curve_t curve = {0};
-	arg_t arg = {.args = "(seed) -> { return(19); }", .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_FIELD);
-	cr_assert_eq(ret, 1, );
-	cr_assert(gequal(curve.field, stoi(19)), );
-}
-
-Test(gp, test_gp_gen_a) {
-	curve_t curve = {.field = stoi(19)};
-	arg_t arg = {.args = "(seed, field) -> { return(Mod(3,field)); }",
-	             .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_A);
-	cr_assert_eq(ret, 1, );
-	cr_assert(gequal(curve.a, mkintmodu(3, 19)), );
-}
-
-Test(gp, test_gp_gen_b) {
-	curve_t curve = {.field = stoi(19), .a = mkintmodu(3, 19)};
-	arg_t arg = {.args = "(seed, field, a) -> { return(a * 2); }", .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_B);
-	cr_assert_eq(ret, 1, );
-	cr_assert(gequal(curve.b, mkintmodu(6, 19)), );
-}
-
-Test(gp, test_gp_gen_curve) {
-	curve_t curve = {
-	    .field = stoi(19), .a = mkintmodu(3, 19), .b = mkintmodu(6, 19)};
-	arg_t arg = {
-	    .args = "(seed, field, a, b) -> { return(ellinit([a,b], field)); }",
-	    .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_CURVE);
-	cr_assert_eq(ret, 1, );
-	cr_assert(gequal(curve.curve,
-	                 ellinit(mkvec2(curve.a, curve.b), curve.field, 0)), );
-}
-
-Test(gp, test_gp_gen_order) {
-	curve_t curve = {.field = stoi(19),
-	                 .a = mkintmodu(3, 19),
-	                 .b = mkintmodu(6, 19),
-	                 .curve = ellinit(mkvec2(stoi(3), stoi(6)), stoi(19), 0)};
-	arg_t arg = {
-	    .args = "(seed, field, a, b, curve) -> { return(ellsea(curve)); }",
-	    .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_ORDER);
-	cr_assert_eq(ret, 1, );
-	cr_assert(gequal(ellsea(curve.curve, 0), curve.order), );
-}
-
-Test(gp, test_gp_gen_generators) {
-	curve_t curve = {.field = stoi(19),
-	                 .a = mkintmodu(3, 19),
-	                 .b = mkintmodu(6, 19),
-	                 .curve = ellinit(mkvec2(stoi(3), stoi(6)), stoi(19), 0),
-	                 .order = stoi(16)};
-	arg_t arg = {.args =
-	                 "(seed, field, a, b, curve, order) -> { "
-	                 "return(ellgenerators(curve)); }",
-	             .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_GENERATORS);
-	cr_assert_eq(ret, 1, );
-
-	GEN ellgens = ellgenerators(curve.curve);
-	long len = glength(ellgens);
-	cr_assert_eq(len, curve.ngens, );
-	for (long i = 1; i <= len; ++i) {
-		cr_assert(gequal(gel(ellgens, i), curve.generators[i - 1]->point), );
-	}
-}
-
-Test(gp, test_gp_gen_points) {
-	curve_t curve = {
-	    .field = stoi(19),
-	    .a = mkintmodu(3, 19),
-	    .b = mkintmodu(6, 19),
-	    .curve = ellinit(mkvec2(stoi(3), stoi(6)), stoi(19), 0),
-	    .order = stoi(16),
-	};
-	point_t gen = {.point = mkvec2(mkintmodu(4, 19), mkintmodu(14, 19))};
-	point_t *generators[1] = {&gen};
-	curve.generators = generators;
-	curve.ngens = 1;
-
-	arg_t arg = {.args =
-	                 "(seed, field, a, b, curve, order, gens) -> { "
-	                 "return([ellmul(curve,gens[1],2)]); }",
-	             .nargs = 1};
-
-	int ret = gp_gen(&curve, &arg, OFFSET_POINTS);
-	cr_assert_eq(ret, 1, );
-	cr_assert_eq(curve.npoints, 1, );
-	cr_assert(gequal(curve.points[0]->point,
-	                 ellmul(curve.curve, gen.point, stoi(2))), );
-}
-
-Test(gp, test_gp_check_seed) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {.seed = &seed};
-
-	arg_t arg = {.args = "(seed) -> { return(1);}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_SEED);
-	cr_assert_eq(ret, 1, );
-}
-
-Test(gp, test_gp_check_field) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {.seed = &seed, .field = stoi(19)};
-
-	arg_t arg = {.args = "(seed, field) -> { if(field == 19, return(1));}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_FIELD);
-	cr_assert_eq(ret, 1, );
-}
-
-Test(gp, test_gp_check_a) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {.seed = &seed, .field = stoi(19), .a = mkintmodu(3, 19)};
-
-	arg_t arg = {.args =
-	                 "(seed, field, a) -> { if(a == Mod(3,19), return(1));}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_A);
-	cr_assert_eq(ret, 1, );
-}
-
-Test(gp, test_gp_check_b) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {.seed = &seed,
-	                 .field = stoi(19),
-	                 .a = mkintmodu(3, 19),
-	                 .b = mkintmodu(5, 19)};
-
-	arg_t arg = {
-	    .args = "(seed, field, a, b) -> { if(b == Mod(5,19), return(1));}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_B);
-	cr_assert_eq(ret, 1, );
-}
-
-Test(gp, test_gp_check_curve) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {.seed = &seed,
-	                 .field = stoi(19),
-	                 .a = mkintmodu(3, 19),
-	                 .b = mkintmodu(5, 19),
-	                 .curve = ellinit(mkvec2(stoi(3), stoi(5)), stoi(19), 0)};
-
-	arg_t arg = {.args =
-	                 "(seed, field, a, b, curve) -> { if(curve == ellinit([3, "
-	                 "5], 19), return(1));}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_CURVE);
-	cr_assert_eq(ret, 1, );
-}
-
-Test(gp, test_gp_check_order) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {.seed = &seed,
-	                 .field = stoi(19),
-	                 .a = mkintmodu(3, 19),
-	                 .b = mkintmodu(5, 19),
-	                 .curve = ellinit(mkvec2(stoi(3), stoi(5)), stoi(19), 0),
-	                 .order = stoi(16)};
-
-	arg_t arg = {.args =
-	                 "(seed, field, a, b, curve, order) -> { if(order == 16, "
-	                 "return(1));}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_ORDER);
-	cr_assert_eq(ret, 1, );
-}
-
-Test(gp, test_gp_check_generators) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {
-	    .seed = &seed,
-	    .field = stoi(19),
-	    .a = mkintmodu(3, 19),
-	    .b = mkintmodu(6, 19),
-	    .curve = ellinit(mkvec2(stoi(3), stoi(6)), stoi(19), 0),
-	    .order = stoi(16),
-	};
-	point_t gen = {.point = mkvec2(mkintmodu(4, 19), mkintmodu(14, 19))};
-	point_t *generators[1] = {&gen};
-	curve.generators = generators;
-	curve.ngens = 1;
-
-	arg_t arg = {.args =
-	                 "(seed, field, a, b, curve, order, gens) -> { if(gens == "
-	                 "ellgenerators(curve), return(1));}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_GENERATORS);
-	cr_assert_eq(ret, 1, );
-}
-
-Test(gp, test_gp_check_points) {
-	seed_t seed = {.seed = bits_from_hex("ff")};
-	curve_t curve = {
-	    .seed = &seed,
-	    .field = stoi(19),
-	    .a = mkintmodu(3, 19),
-	    .b = mkintmodu(6, 19),
-	    .curve = ellinit(mkvec2(stoi(3), stoi(6)), stoi(19), 0),
-	    .order = stoi(16),
-	};
-	point_t gen = {.point = mkvec2(mkintmodu(4, 19), mkintmodu(14, 19))};
-	point_t *generators[1] = {&gen};
-	curve.generators = generators;
-	curve.ngens = 1;
-	curve.points = generators;
-	curve.npoints = 1;
-
-	arg_t arg = {.args =
-	                 "(seed, field, a, b, curve, order, gens, points) -> { "
-	                 "if(points == ellgenerators(curve), return(1));}"};
-
-	int ret = gp_check(&curve, &arg, OFFSET_POINTS);
-	cr_assert_eq(ret, 1, );
-}
\ No newline at end of file
diff --git a/test/src/gen/test_hex.c b/test/src/gen/test_hex.c
index aa558ec..d1f1931 100644
--- a/test/src/gen/test_hex.c
+++ b/test/src/gen/test_hex.c
@@ -91,7 +91,8 @@ Test(hex, test_hex_generators) {
 	char *what = "ABCDE";
 	point_t gen = {.point = mkvec2(mkintmod(strtoi(where), strtoi(field)),
 	                               mkintmod(stoi(52), strtoi(field)))};
-	point_t *generators[1] = {&gen};
+	subgroup_t gen_subgroup = {.generator = &gen};
+	subgroup_t *subgroups[1] = {&gen_subgroup};
 
 	seed_t seed = {.seed = bits_from_hex("12345")};
 	curve_t curve = {.seed = &seed,
@@ -99,7 +100,7 @@ Test(hex, test_hex_generators) {
 	                 .a = mkintmod(stoi(15), strtoi(field)),
 	                 .b = mkintmod(stoi(20), strtoi(field)),
 	                 .order = stoi(22),
-	                 .generators = generators,
+	                 .generators = subgroups,
 	                 .ngens = 1};
 
 	arg_t arg = {.args = what, .nargs = 1};
@@ -115,10 +116,13 @@ Test(hex, test_hex_points) {
 	char *what = "ABCDE";
 	point_t gen = {.point = mkvec2(mkintmod(stoi(23), strtoi(field)),
 	                               mkintmod(stoi(52), strtoi(field)))};
-	point_t *generators[1] = {&gen};
+	subgroup_t gen_subgroup = {.generator = &gen};
+	subgroup_t *subgroups[1] = {&gen_subgroup};
 	point_t point = {.point = mkvec2(mkintmod(strtoi(where), strtoi(field)),
 	                                 mkintmod(stoi(31), strtoi(field)))};
 	point_t *points[1] = {&point};
+	gen_subgroup.npoints = 1;
+	gen_subgroup.points = points;
 
 	seed_t seed = {.seed = bits_from_hex("12345")};
 	curve_t curve = {.seed = &seed,
@@ -126,10 +130,8 @@ Test(hex, test_hex_points) {
 	                 .a = mkintmod(stoi(15), strtoi(field)),
 	                 .b = mkintmod(stoi(20), strtoi(field)),
 	                 .order = stoi(22),
-	                 .generators = generators,
-	                 .ngens = 1,
-	                 .points = points,
-	                 .npoints = 1};
+	                 .generators = subgroups,
+	                 .ngens = 1};
 
 	arg_t arg = {.args = what, .nargs = 1};
 
diff --git a/test/src/gen/test_point.c b/test/src/gen/test_point.c
index ca23ba6..cb27495 100644
--- a/test/src/gen/test_point.c
+++ b/test/src/gen/test_point.c
@@ -5,6 +5,7 @@
 
 #include <criterion/criterion.h>
 #include "gen/point.h"
+#include "math/subgroup.h"
 #include "test/io.h"
 
 TestSuite(point, .init = io_setup, .fini = io_teardown);
@@ -31,118 +32,160 @@ Test(point, test_points_clone) {
 	points_free_deep(&others, 1);
 }
 
+// TODO: add utility functions to setup the two example curves used in these tests.
+
 Test(point, test_point_random) {
 	// curve = ellinit([1, 3], 23), order = 27
 	GEN e = ellinit(mkvec2s(1, 3), stoi(23), -1);
-	curve_t curve = {.order = stoi(27), .curve = e};
+	GEN gen = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
+	point_t gen_point = {.point = gen, .order = stoi(27), .cofactor = stoi(1)};
+	subgroup_t gen_subgroup = {.generator = &gen_point};
+	subgroup_t *generators[1] = {&gen_subgroup};
+	curve_t curve = {.order = stoi(27), .curve = e, .generators = generators};
 	int ret = point_gen_random(&curve, NULL, OFFSET_POINTS);
 
 	cr_assert_eq(ret, 1, "Point wasn't generated.");
-	cr_assert_eq(curve.npoints, 1, "Incorrect number of points.");
-	cr_assert_not_null(curve.points, "Points are null.");
-	cr_assert(ellisoncurve(e, curve.points[0]->point), "Point not on curve.");
-	cr_assert(gequal(ellorder(e, curve.points[0]->point, NULL),
-	                 curve.points[0]->order),
+	cr_assert_eq(gen_subgroup.npoints, 1, "Incorrect number of points.");
+	cr_assert_not_null(gen_subgroup.points, "Points are null.");
+	cr_assert(ellisoncurve(e, gen_subgroup.points[0]->point), "Point not on curve.");
+	cr_assert(gequal(ellorder(e, gen_subgroup.points[0]->point, NULL),
+					 gen_subgroup.points[0]->order),
 	          "Point has wrong order set.");
 
-	points_free_deep(&curve.points, 1);
+	points_free_deep(&gen_subgroup.points, gen_subgroup.npoints);
+}
+
+Test(point, test_point_random_more_gens) {
+	GEN e = ellinit(mkvec2s(2,3), stoi(23), -1);
+	GEN one_gen = mkvec2(mkintmodu(6, 23), mkintmodu(1,23));
+	point_t *one_gen_point = point_new();
+	one_gen_point->point = one_gen;
+	one_gen_point->order = stoi(12);
+	one_gen_point->cofactor = stoi(2);
+	subgroup_t *one_subgroup = subgroup_new();
+	one_subgroup->generator = one_gen_point;
+	GEN two_gen = mkvec2(mkintmodu(20, 23), mkintmodu(19, 23));
+	point_t *two_gen_point = point_new();
+	two_gen_point->point = two_gen;
+	two_gen_point->order =  stoi(6);
+	two_gen_point->cofactor = stoi(4);
+	subgroup_t *two_subgroup = subgroup_new();
+	two_subgroup->generator = two_gen_point;
+	subgroup_t **subgroups = subgroups_new(2);
+	subgroups[0] = one_subgroup;
+	subgroups[1] = two_subgroup;
+	curve_t curve = {.order = stoi(24), .curve = e, .generators = subgroups, .ngens = 2};
+	int ret = point_gen_random(&curve, NULL, OFFSET_POINTS);
+
+	cr_assert_eq(ret, 1, "Point wasn't generated.");
+	size_t generated = 0;
+	for (size_t i = 0; i < 2; ++i) {
+		subgroup_t *subgroup = curve.generators[i];
+		if (subgroup->npoints > 0) {
+			generated += subgroup->npoints;
+			cr_assert_not_null(subgroup->points, "Points are null.");
+			cr_assert(ellisoncurve(e, subgroup->points[0]->point), "Point not on curve.");
+			cr_assert(gequal(ellorder(e, subgroup->points[0]->point, NULL),
+							 subgroup->points[0]->order),
+					  "Point has wrong order set.");
+		}
+	}
+	cr_assert_eq(generated, 1, "Point wasn't saved.");
+
+	subgroups_free_deep(&subgroups, 2);
 }
 
 Test(point, test_points_random) {
 	// curve = ellinit([1, 3], 23), order = 27
 	GEN e = ellinit(mkvec2s(1, 3), stoi(23), -1);
-	curve_t curve = {.order = stoi(27), .curve = e};
+	GEN gen = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
+	point_t gen_point = {.point = gen, .order = stoi(27), .cofactor = stoi(1)};
+	subgroup_t gen_subgroup = {.generator = &gen_point};
+	subgroup_t *generators[1] = {&gen_subgroup};
+	curve_t curve = {.order = stoi(27), .curve = e, .generators = generators, .ngens = 1};
 	size_t npoints = 3;
 	arg_t arg = {.args = &npoints, .nargs = 1};
 	int ret = points_gen_random(&curve, &arg, OFFSET_POINTS);
 
 	cr_assert_eq(ret, 1, "Points weren't generated.");
-	cr_assert_eq(curve.npoints, npoints, "Incorrect number of points.");
-	cr_assert_not_null(curve.points, "Points are null.");
+	cr_assert_eq(gen_subgroup.npoints, npoints, "Incorrect number of points.");
+	cr_assert_not_null(gen_subgroup.points, "Points are null.");
 	for (size_t i = 0; i < npoints; i++) {
-		point_t *point = curve.points[i];
+		point_t *point = gen_subgroup.points[i];
 		cr_assert(ellisoncurve(e, point->point), "Point not on curve.");
 		cr_assert(gequal(ellorder(e, point->point, NULL), point->order),
 		          "Point has wrong order set.");
 	}
 
-	points_free_deep(&curve.points, npoints);
+	points_free_deep(&gen_subgroup.points, npoints);
 }
 
 Test(point, test_points_trial) {
 	// curve = ellinit([1, 3], 23), order = 27
 	GEN e = ellinit(mkvec2s(1, 3), stoi(23), -1);
-	point_t **gens = points_new(1);
-	gens[0] = point_new();
-	gens[0]->order = stoi(27);
-	gens[0]->cofactor = stoi(1);
-	gens[0]->point = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
-	curve_t curve = {
-	    .order = stoi(27), .curve = e, .ngens = 1, .generators = gens};
+	GEN gen = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
+	point_t gen_point = {.point = gen, .order = stoi(27), .cofactor = stoi(1)};
+	subgroup_t gen_subgroup = {.generator = &gen_point};
+	subgroup_t *generators[1] = {&gen_subgroup};
+	curve_t curve = {.order = stoi(27), .curve = e, .generators = generators, .ngens = 1};
 	pari_ulong prime = 3;
 	arg_t arg = {.args = &prime, .nargs = 1};
 	int ret = points_gen_trial(&curve, &arg, OFFSET_POINTS);
 
 	cr_assert_eq(ret, 1, "Points weren't generated.");
-	cr_assert_eq(curve.npoints, 1, "Incorrect number of points.");
-	cr_assert_not_null(curve.points, "Points are null.");
-	cr_assert(ellisoncurve(e, curve.points[0]->point), "Point not on curve.");
-	cr_assert(gequal(ellorder(e, curve.points[0]->point, NULL),
-	                 curve.points[0]->order),
+	cr_assert_eq(gen_subgroup.npoints, 1, "Incorrect number of points.");
+	cr_assert_not_null(gen_subgroup.points, "Points are null.");
+	cr_assert(ellisoncurve(e, gen_subgroup.points[0]->point), "Point not on curve.");
+	cr_assert(gequal(ellorder(e, gen_subgroup.points[0]->point, NULL),
+					 gen_subgroup.points[0]->order),
 	          "Point has wrong order set.");
-	cr_assert(gequal(curve.points[0]->order, utoi(prime)),
+	cr_assert(gequal(gen_subgroup.points[0]->order, utoi(prime)),
 	          "Point has wrong order.");
 
-	points_free_deep(&curve.points, 1);
-	points_free_deep(&curve.generators, 1);
+	points_free_deep(&gen_subgroup.points, 1);
 }
 
 Test(point, test_points_prime) {
 	// curve = ellinit([1, 3], 23), order = 27
 	GEN e = ellinit(mkvec2s(1, 3), stoi(23), -1);
-	point_t **gens = points_new(1);
-	gens[0] = point_new();
-	gens[0]->order = stoi(27);
-	gens[0]->cofactor = stoi(1);
-	gens[0]->point = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
-	curve_t curve = {
-	    .order = stoi(27), .curve = e, .ngens = 1, .generators = gens};
+	GEN gen = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
+	point_t gen_point = {.point = gen, .order = stoi(27), .cofactor = stoi(1)};
+	subgroup_t gen_subgroup = {.generator = &gen_point};
+	subgroup_t *generators[1] = {&gen_subgroup};
+	curve_t curve = {.order = stoi(27), .curve = e, .generators = generators, .ngens = 1};
 	pari_ulong prime = 3;
 	int ret = points_gen_prime(&curve, NULL, OFFSET_POINTS);
 
 	cr_assert_eq(ret, 1, "Points weren't generated.");
-	cr_assert_eq(curve.npoints, 1, "Incorrect number of points.");
-	cr_assert_not_null(curve.points, "Points are null.");
-	cr_assert(ellisoncurve(e, curve.points[0]->point), "Point not on curve.");
-	cr_assert(gequal(ellorder(e, curve.points[0]->point, NULL),
-	                 curve.points[0]->order),
+	cr_assert_eq(gen_subgroup.npoints, 1, "Incorrect number of points.");
+	cr_assert_not_null(gen_subgroup.points, "Points are null.");
+	cr_assert(ellisoncurve(e, gen_subgroup.points[0]->point), "Point not on curve.");
+	cr_assert(gequal(ellorder(e, gen_subgroup.points[0]->point, NULL),
+					 gen_subgroup.points[0]->order),
 	          "Point has wrong order set.");
-	cr_assert(gequal(curve.points[0]->order, utoi(prime)),
+	cr_assert(gequal(gen_subgroup.points[0]->order, utoi(prime)),
 	          "Point has wrong order.");
 
-	points_free_deep(&curve.points, 1);
-	points_free_deep(&curve.generators, 1);
+	points_free_deep(&gen_subgroup.points, 1);
 }
 
 Test(point, test_points_all) {
 	// curve = ellinit([1, 3], 23), order = 27
 	GEN e = ellinit(mkvec2s(1, 3), stoi(23), -1);
-	point_t **gens = points_new(1);
-	gens[0] = point_new();
-	gens[0]->order = stoi(27);
-	gens[0]->cofactor = stoi(1);
-	gens[0]->point = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
-	curve_t curve = {
-	    .order = stoi(27), .curve = e, .ngens = 1, .generators = gens};
+	GEN gen = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
+	point_t gen_point = {.point = gen, .order = stoi(27), .cofactor = stoi(1)};
+	subgroup_t gen_subgroup = {.generator = &gen_point};
+	subgroup_t *generators[1] = {&gen_subgroup};
+	curve_t curve = {.order = stoi(27), .curve = e, .generators = generators, .ngens = 1};
 	GEN orders = mkvec3s(3, 9, 27);
 	size_t npoints = 3;
 	int ret = points_gen_allgroups(&curve, NULL, OFFSET_POINTS);
 
 	cr_assert_eq(ret, 1, "Points weren't generated.");
-	cr_assert_eq(curve.npoints, npoints, "Incorrect number of points.");
-	cr_assert_not_null(curve.points, "Points are null.");
+	cr_assert_eq(gen_subgroup.npoints, npoints, "Incorrect number of points.");
+	cr_assert_not_null(gen_subgroup.points, "Points are null.");
 	for (size_t i = 0; i < npoints; i++) {
-		point_t *point = curve.points[i];
+		point_t *point = gen_subgroup.points[i];
 		cr_assert(ellisoncurve(e, point->point), "Point not on curve.");
 		cr_assert(gequal(ellorder(e, point->point, NULL), point->order),
 		          "Point has wrong order set.");
@@ -150,29 +193,26 @@ Test(point, test_points_all) {
 		          "Point has wrong order.");
 	}
 
-	points_free_deep(&curve.points, 1);
-	points_free_deep(&curve.generators, 1);
+	points_free_deep(&gen_subgroup.points, 1);
 }
 
 Test(point, test_points_nonprime) {
 	// curve = ellinit([1, 3], 23), order = 27
 	GEN e = ellinit(mkvec2s(1, 3), stoi(23), -1);
-	point_t **gens = points_new(1);
-	gens[0] = point_new();
-	gens[0]->order = stoi(27);
-	gens[0]->cofactor = stoi(1);
-	gens[0]->point = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
-	curve_t curve = {
-	    .order = stoi(27), .curve = e, .ngens = 1, .generators = gens};
+	GEN gen = mkvec2(mkintmodu(15, 23), mkintmodu(14, 23));
+	point_t gen_point = {.point = gen, .order = stoi(27), .cofactor = stoi(1)};
+	subgroup_t gen_subgroup = {.generator = &gen_point};
+	subgroup_t *generators[1] = {&gen_subgroup};
+	curve_t curve = {.order = stoi(27), .curve = e, .generators = generators, .ngens = 1};
 	GEN orders = mkvec2s(9, 27);
 	size_t npoints = 2;
 	int ret = points_gen_nonprime(&curve, NULL, OFFSET_POINTS);
 
 	cr_assert_eq(ret, 1, "Points weren't generated.");
-	cr_assert_eq(curve.npoints, npoints, "Incorrect number of points.");
-	cr_assert_not_null(curve.points, "Points are null.");
+	cr_assert_eq(gen_subgroup.npoints, npoints, "Incorrect number of points.");
+	cr_assert_not_null(gen_subgroup.points, "Points are null.");
 	for (size_t i = 0; i < npoints; i++) {
-		point_t *point = curve.points[i];
+		point_t *point = gen_subgroup.points[i];
 		cr_assert(ellisoncurve(e, point->point), "Point not on curve.");
 		cr_assert(gequal(ellorder(e, point->point, NULL), point->order),
 		          "Point has wrong order set.");
@@ -180,6 +220,5 @@ Test(point, test_points_nonprime) {
 		          "Point has wrong order.");
 	}
 
-	points_free_deep(&curve.points, 1);
-	points_free_deep(&curve.generators, 1);
+	points_free_deep(&gen_subgroup.points, 1);
 }
diff --git a/test/src/math/test_subgroup.c b/test/src/math/test_subgroup.c
new file mode 100644
index 0000000..1acb468
--- /dev/null
+++ b/test/src/math/test_subgroup.c
@@ -0,0 +1,40 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017 J08nY
+ */
+#include <criterion/criterion.h>
+#include "gen/point.h"
+#include "math/subgroup.h"
+#include "test/default.h"
+
+TestSuite(subgroup, .init = default_setup, .fini = default_teardown);
+
+Test(subgroup, test_prime_factors) {
+	GEN divs = subgroups_prime(stoi(12));
+	GEN vec = mkvec2s(2, 3);
+	cr_assert(gequal(divs, vec), "Factors not equal!");
+}
+
+Test(subgroup, test_prime_factors_other) {
+	GEN divs = subgroups_prime(stoi(27));
+	GEN vec = gtovec(stoi(3));
+	cr_assert(gequal(divs, vec), "Factors not equal!");
+}
+
+Test(subgroup, test_prime_prime) {
+	GEN divs = subgroups_prime(stoi(5));
+	GEN vec = gtovec(stoi(5));
+	cr_assert(gequal(divs, vec), "Factors not equal!");
+}
+
+Test(subgroup, test_nonprime_factors) {
+	GEN divs = subgroups_nonprime(stoi(27));
+	GEN vec = mkvec2s(9, 27);
+	cr_assert(gequal(divs, vec), "Factors not equal!");
+}
+
+Test(subgroup, test_all_factors) {
+	GEN divs = subgroups_all(stoi(27));
+	GEN vec = mkvec3s(3, 9, 27);
+	cr_assert(gequal(divs, vec), "Factors not equal!");
+}
\ No newline at end of file
diff --git a/test/src/math/test_subgroups.c b/test/src/math/test_subgroups.c
deleted file mode 100644
index f944a4f..0000000
--- a/test/src/math/test_subgroups.c
+++ /dev/null
@@ -1,73 +0,0 @@
-/*
- * ecgen, tool for generating Elliptic curve domain parameters
- * Copyright (C) 2017 J08nY
- */
-#include <criterion/criterion.h>
-#include "gen/point.h"
-#include "math/subgroups.h"
-#include "test/default.h"
-
-TestSuite(subgroups, .init = default_setup, .fini = default_teardown);
-
-Test(subgroups, test_prime_factors) {
-	curve_t curve = {.order = stoi(12)};
-	cfg->prime = false;
-	GEN divs = subgroups_prime(&curve);
-	GEN vec = mkvec2s(2, 3);
-	cr_assert(gequal(divs, vec), "Factors not equal!");
-}
-
-Test(subgroups, test_prime_factors_other) {
-	curve_t curve = {.order = stoi(27)};
-	cfg->prime = false;
-	GEN divs = subgroups_prime(&curve);
-	GEN vec = gtovec(stoi(3));
-	cr_assert(gequal(divs, vec), "Factors not equal!");
-}
-
-Test(subgroups, test_prime_prime) {
-	curve_t curve = {.order = stoi(5)};
-	cfg->prime = true;
-	GEN divs = subgroups_prime(&curve);
-	GEN vec = gtovec(stoi(5));
-	cr_assert(gequal(divs, vec), "Factors not equal!");
-}
-
-Test(subgroups, test_nonprime_factors) {
-	// curve = ellinit([1, 3], 23), order = 27
-	curve_t curve = {.order = stoi(27), .ngens = 1};
-	cfg->prime = false;
-	GEN divs = subgroups_nonprime(&curve);
-	GEN vec = mkvec2s(9, 27);
-	cr_assert(gequal(divs, vec), "Factors not equal!");
-}
-
-Test(subgroups, test_all_factors) {
-	// curve = ellinit([1, 3], 23), order = 27
-	curve_t curve = {.order = stoi(27), .ngens = 1};
-	cfg->prime = false;
-	GEN divs = subgroups_all(&curve);
-	GEN vec = mkvec3s(3, 9, 27);
-	cr_assert(gequal(divs, vec), "Factors not equal!");
-}
-
-Test(subgroups, test_all_factors_two_gens) {
-	// curve = ellinit([2, 3], 23), order = 24
-	point_t **gens = points_new(2);
-	gens[0] = point_new();
-	gens[0]->order = stoi(12);
-	gens[0]->cofactor = stoi(2);
-	gens[0]->point = mkvec2(mkintmodu(6, 23), mkintmodu(1, 23));
-	gens[1] = point_new();
-	gens[1]->order = stoi(6);
-	gens[1]->cofactor = stoi(4);
-	gens[1]->point = mkvec2(mkintmodu(20, 23), mkintmodu(19, 23));
-
-	curve_t curve = {.order = stoi(24), .ngens = 2, .generators = gens};
-	cfg->prime = false;
-	GEN divs = subgroups_all(&curve);
-	GEN vec = mkvecn(8, stoi(2), stoi(3), stoi(4), stoi(6), stoi(12), stoi(2),
-	                 stoi(3), stoi(6));
-	cr_assert(gequal(divs, vec), "Factors not equal!");
-	points_free_deep(&gens, 2);
-}
diff --git a/test/src/util/test_random.c b/test/src/util/test_random.c
index be45f95..8f10453 100644
--- a/test/src/util/test_random.c
+++ b/test/src/util/test_random.c
@@ -32,6 +32,14 @@ Test(random, test_random_int) {
 	}
 }
 
+Test(random, test_random_range) {
+	for (size_t i = 0; i < 100; ++i) {
+		GEN j = random_range(stoi(100), stoi(5000));
+		cr_assert_lt(cmpii(j, stoi(5000)), 0, );
+		cr_assert_geq(cmpii(j, stoi(100)), 0, );
+	}
+}
+
 Test(random, test_random_field_element_fp) {
 	GEN fp = random_int(25);
 	for (size_t i = 0; i < 100; ++i) {