diff --git a/sources/flintinterface.cc b/sources/flintinterface.cc
index 0c4fa78c..4c435877 100644
--- a/sources/flintinterface.cc
+++ b/sources/flintinterface.cc
@@ -37,6 +37,13 @@ static_assert(sizeof(slong) == sizeof(int64_t), "flint interface expects slong i
 	#] Types
 */
 
+/*
+ * FLINT's univariate poly has a dense representation. For sufficiently sparse polynomials it is
+ * faster to use mpoly instead, which is sparse. For a density <= this threshold we switch, where
+ * the density defined as is "number of terms" / "maximum degree".
+ */
+#define UNIVARIATE_DENSITY_THR 0.02f
+
 /*
 	#[ flint::cleanup :
 */
@@ -1061,7 +1068,9 @@ flint::var_map_t flint::get_variables(const vector <WORD *> &es, const bool with
 	const bool sort_vars) {
 
 	int32_t num_vars = 0;
-	// To be used if we sort by highest degree, as the polu code does.
+	// We count the total number of terms to determine "density".
+	uint32_t num_terms = 0;
+	// To be used if we sort by highest degree, as the poly code does.
 	vector<int32_t> degrees;
 	var_map_t var_map;
 
@@ -1071,8 +1080,10 @@ flint::var_map_t flint::get_variables(const vector <WORD *> &es, const bool with
 
 		// fast notation
 		if ( *e == -SNUMBER ) {
+			num_terms++;
 		}
 		else if ( *e == -SYMBOL ) {
+			num_terms++;
 			if ( !var_map.count(e[1]) ) {
 				var_map[e[1]] = num_vars++;
 				degrees.push_back(1);
@@ -1087,6 +1098,7 @@ flint::var_map_t flint::get_variables(const vector <WORD *> &es, const bool with
 		}
 		else {
 			for ( WORD i = with_arghead ? ARGHEAD:0; with_arghead ? i < e[0]:e[i] != 0; i += e[i] ) {
+				num_terms++;
 				if ( i+1 < i+e[i]-ABS(e[i+e[i]-1]) && e[i+1] != SYMBOL ) {
 					MLOCK(ErrorMessageLock);
 					MesPrint("ERROR: polynomials and polyratfuns must contain symbols only");
@@ -1152,6 +1164,16 @@ flint::var_map_t flint::get_variables(const vector <WORD *> &es, const bool with
 		}
 	}
 
+	if ( var_map.size() == 1 ) {
+		// In the univariate case, if the polynomials are sufficiently sparse force the use of the
+		// multivariate routines, which use a sparse representation, by adding a dummy map entry.
+		if ( (float)num_terms <= UNIVARIATE_DENSITY_THR * (float)degrees[0] ) {
+			// -1 will never be a symbol code. Built-in symbols from 0 to 19, and 20 is the first
+			// user symbol.
+			var_map[-1] = num_vars;
+		}
+	}
+
 	return var_map;
 }
 /*
diff --git a/sources/flintwrap.cc b/sources/flintwrap.cc
index 7a9b34a8..f211a196 100644
--- a/sources/flintwrap.cc
+++ b/sources/flintwrap.cc
@@ -36,10 +36,12 @@ WORD* flint_div(PHEAD WORD *a, WORD *b, const WORD must_fit_term) {
 	e.push_back(b);
 	const bool with_arghead = false;
 	const bool sort_vars = false;
-	const flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
+	flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
 
 	const bool return_rem = false;
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		return flint::divmod_mpoly(BHEAD a, b, return_rem, must_fit_term, var_map);
 	}
 	else {
@@ -54,11 +56,13 @@ int flint_factorize_argument(PHEAD WORD *argin, WORD *argout) {
 
 	const bool with_arghead = true;
 	const bool sort_vars = true;
-	const flint::var_map_t var_map = flint::get_variables(vector<WORD*>(1,argin), with_arghead,
+	flint::var_map_t var_map = flint::get_variables(vector<WORD*>(1,argin), with_arghead,
 		sort_vars);
 
 	const bool is_fun_arg = true;
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		flint::factorize_mpoly(BHEAD argin, argout, with_arghead, is_fun_arg, var_map);
 	}
 	else {
@@ -75,11 +79,13 @@ WORD* flint_factorize_dollar(PHEAD WORD *argin) {
 
 	const bool with_arghead = false;
 	const bool sort_vars = true;
-	const flint::var_map_t var_map = flint::get_variables(vector<WORD*>(1,argin), with_arghead,
+	flint::var_map_t var_map = flint::get_variables(vector<WORD*>(1,argin), with_arghead,
 		sort_vars);
 
 	const bool is_fun_arg = false;
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		return flint::factorize_mpoly(BHEAD argin, NULL, with_arghead, is_fun_arg, var_map);
 	}
 	else {
@@ -98,9 +104,11 @@ WORD* flint_gcd(PHEAD WORD *a, WORD *b, const WORD must_fit_term) {
 	e.push_back(b);
 	const bool with_arghead = false;
 	const bool sort_vars = true;
-	const flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
+	flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
 
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		return flint::gcd_mpoly(BHEAD a, b, must_fit_term, var_map);
 	}
 	else {
@@ -117,7 +125,9 @@ WORD* flint_inverse(PHEAD WORD *a, WORD *b) {
 	e.reserve(2);
 	e.push_back(a);
 	e.push_back(b);
-	const flint::var_map_t var_map = flint::get_variables(e, false, false);
+	flint::var_map_t var_map = flint::get_variables(e, false, false);
+	// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+	var_map.erase(-1);
 
 	if ( var_map.size() > 1 ) {
 		MLOCK(ErrorMessageLock);
@@ -138,9 +148,11 @@ WORD* flint_mul(PHEAD WORD *a, WORD *b) {
 	e.reserve(2);
 	e.push_back(a);
 	e.push_back(b);
-	const flint::var_map_t var_map = flint::get_variables(e, false, false);
+	flint::var_map_t var_map = flint::get_variables(e, false, false);
 
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		return flint::mul_mpoly(BHEAD a, b, var_map);
 	}
 	else {
@@ -175,9 +187,11 @@ WORD* flint_ratfun_add(PHEAD WORD *t1, WORD *t2) {
 	}
 	const bool with_arghead = true;
 	const bool sort_vars = true;
-	const flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
+	flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
 
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		flint::ratfun_add_mpoly(BHEAD t1, t2, oldworkpointer, var_map);
 	}
 	else {
@@ -241,9 +255,11 @@ int flint_ratfun_normalize(PHEAD WORD *term) {
 	}
 	const bool with_arghead = true;
 	const bool sort_vars = true;
-	const flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
+	flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
 
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		flint::ratfun_normalize_mpoly(BHEAD term, var_map);
 	}
 	else {
@@ -271,10 +287,12 @@ WORD* flint_rem(PHEAD WORD *a, WORD *b, const WORD must_fit_term) {
 	e.push_back(b);
 	const bool with_arghead = false;
 	const bool sort_vars = false;
-	const flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
+	flint::var_map_t var_map = flint::get_variables(e, with_arghead, sort_vars);
 
 	const bool return_rem = true;
 	if ( var_map.size() > 1 ) {
+		// a "-1" just variable forces mpoly routines for a univariate problem, remove it:
+		var_map.erase(-1);
 		return flint::divmod_mpoly(BHEAD a, b, return_rem, must_fit_term, var_map);
 	}
 	else {