diff --git a/M2/Macaulay2/e/CMakeLists.txt b/M2/Macaulay2/e/CMakeLists.txt
index fc65c439a56..ce702082fef 100644
--- a/M2/Macaulay2/e/CMakeLists.txt
+++ b/M2/Macaulay2/e/CMakeLists.txt
@@ -89,7 +89,6 @@ set(HPPONLYFILES
   newdelete.hpp
   ExponentVector.hpp
   res-a0-pair.hpp
-  ringelem.hpp
   schur-poly-heap.hpp
   smat.hpp
   util.hpp
@@ -266,6 +265,7 @@ set(SRCLIST
   res-a1
   res-a2
   ring
+  ringelem
   ringmap
   sagbi
   schorder
diff --git a/M2/Macaulay2/e/GF.cpp b/M2/Macaulay2/e/GF.cpp
index 1081e94263c..2357e28388a 100644
--- a/M2/Macaulay2/e/GF.cpp
+++ b/M2/Macaulay2/e/GF.cpp
@@ -256,15 +256,16 @@ bool GF::promote(const Ring *Rf, const ring_elem f, ring_elem &result) const
 
   if (Rf != _originalR) return false;
 
+  const Ring *K = _originalR->getCoefficients();
+
   result = from_long(0);
   int exp[1];
-  for (Nterm *t = f; t != NULL; t = t->next)
+  for (Nterm& t : f)
     {
-      std::pair<bool, long> b =
-          _originalR->getCoefficients()->coerceToLongInteger(t->coeff);
+      std::pair<bool, long> b = K->coerceToLongInteger(t.coeff);
       assert(b.first);
       ring_elem coef = from_long(b.second);
-      _originalR->getMonoid()->to_expvector(t->monom, exp);
+      _originalR->getMonoid()->to_expvector(t.monom, exp);
       // exp[0] is the variable we want.  Notice that since the ring is a
       // quotient,
       // this degree is < n (where Q_ = P^n).
diff --git a/M2/Macaulay2/e/Makefile.files b/M2/Macaulay2/e/Makefile.files
index f5ffec84798..5f9566315bf 100644
--- a/M2/Macaulay2/e/Makefile.files
+++ b/M2/Macaulay2/e/Makefile.files
@@ -159,6 +159,7 @@ INTERFACE = \
 	pfaff \
 	relem \
 	ring \
+	ringelem \
 	ringmap \
 	schur \
 	skew \
@@ -205,7 +206,6 @@ NAMES_H = \
 	VectorArithmetic \
 	MemoryBlock \
 	NCAlgebras/Range \
-	ringelem \
 	style \
 	util \
 	hash \
diff --git a/M2/Macaulay2/e/bibasis/involutive.hpp b/M2/Macaulay2/e/bibasis/involutive.hpp
index eabb40bc5b9..2fa1141905d 100644
--- a/M2/Macaulay2/e/bibasis/involutive.hpp
+++ b/M2/Macaulay2/e/bibasis/involutive.hpp
@@ -168,10 +168,10 @@ namespace BIBasis
 
             Polynom<MonomType>* currentPolynom = new Polynom<MonomType>();
 
-            for (Nterm* currentTerm = polynomVector->coeff; currentTerm; currentTerm = currentTerm->next)
+            for (Nterm& currentTerm : polynomVector->coeff)
             {
                 exponents_t monomVector = newarray_atomic(int, independ);
-                monoid->to_expvector(currentTerm->monom, monomVector);
+                monoid->to_expvector(currentTerm.monom, monomVector);
                 
                 //construct Monom for every term
                 MonomType* currentMonom = new MonomType();
diff --git a/M2/Macaulay2/e/franzi-interface.cpp b/M2/Macaulay2/e/franzi-interface.cpp
index 4dd45d1bc57..6b9a7978ea8 100644
--- a/M2/Macaulay2/e/franzi-interface.cpp
+++ b/M2/Macaulay2/e/franzi-interface.cpp
@@ -51,9 +51,9 @@ IntermediateBasis BRPSfromMatrix(const Matrix *m)
       vec v = m->elem(i);
       if (v == 0) continue;
       BRP temp;
-      for (Nterm *f = v->coeff; f != 0; f = f->next)
+      for (Nterm& f : v->coeff)
         {
-          M->to_expvector(f->monom, exp);
+          M->to_expvector(f.monom, exp);
           brMonomial mono = exponentsToLong(n, exp);
           temp = temp + BRP(mono);
         }
diff --git a/M2/Macaulay2/e/hilb.cpp b/M2/Macaulay2/e/hilb.cpp
index 375f918f63b..b4bc0c3f540 100644
--- a/M2/Macaulay2/e/hilb.cpp
+++ b/M2/Macaulay2/e/hilb.cpp
@@ -711,9 +711,9 @@ int hilb_comp::coeff_of(const RingElement *h, int deg)
 
   exponents_t exp = newarray_atomic(int, P->n_vars());
   int result = 0;
-  for (Nterm *f = h->get_value(); f != NULL; f = f->next)
+  for (Nterm& f : h->get_value())
     {
-      P->getMonoid()->to_expvector(f->monom, exp);
+      P->getMonoid()->to_expvector(f.monom, exp);
       if (exp[0] < deg)
         {
           ERROR("incorrect Hilbert function given");
@@ -725,7 +725,7 @@ int hilb_comp::coeff_of(const RingElement *h, int deg)
       else if (exp[0] == deg)
         {
           std::pair<bool, long> res =
-              P->getCoefficientRing()->coerceToLongInteger(f->coeff);
+              P->getCoefficientRing()->coerceToLongInteger(f.coeff);
           assert(res.first &&
                  std::abs(res.second) < std::numeric_limits<int>::max());
           int n = static_cast<int>(res.second);
diff --git a/M2/Macaulay2/e/matrix-stream.hpp b/M2/Macaulay2/e/matrix-stream.hpp
index cd22f7c78d9..3fb33e4127b 100644
--- a/M2/Macaulay2/e/matrix-stream.hpp
+++ b/M2/Macaulay2/e/matrix-stream.hpp
@@ -101,7 +101,7 @@ void matrixToStream(const Matrix* M, T& stream)
       for (; i.valid(); i.next())
         {
           Nterm* t = i.entry();
-          for (Nterm* s = t; s != 0; s = s->next) nterms++;
+          for ([[maybe_unused]] Nterm& s : t) nterms++;
         }
       stream.appendPolynomialBegin(nterms);
 
@@ -110,13 +110,13 @@ void matrixToStream(const Matrix* M, T& stream)
       for (; i.valid(); i.next())
         {
           Nterm* t = i.entry();
-          for (Nterm* s = t; s != 0; s = s->next)
+          for (Nterm& s : t)
             {
-              P->getMonoid()->to_expvector(s->monom, exp);
+              P->getMonoid()->to_expvector(s.monom, exp);
               stream.appendTermBegin(i.row());
               for (size_t j = 0; j < nvars; j++)
                 if (exp[j] != 0) stream.appendExponent(j, exp[j]);
-              std::pair<bool, long> b = KK->coerceToLongInteger(s->coeff);
+              std::pair<bool, long> b = KK->coerceToLongInteger(s.coeff);
               assert(b.first);
               int a = static_cast<int>(
                   b.second);  // This will fit, as the charac fits into an int
diff --git a/M2/Macaulay2/e/matrix.cpp b/M2/Macaulay2/e/matrix.cpp
index be2119e1dbc..b1aeb5419e6 100644
--- a/M2/Macaulay2/e/matrix.cpp
+++ b/M2/Macaulay2/e/matrix.cpp
@@ -904,9 +904,9 @@ M2_arrayintOrNull Matrix::support() const
     for (int i = 0; i < R->n_vars(); i++) exp[i] = exp2[i] = 0;
     for (int j = 0; j < n_cols(); j++)
       for (vec v = elem(j); v != nullptr; v = v->next)
-        for (const Nterm *f = v->coeff; f != nullptr; f = f->next)
+        for (Nterm& f : v->coeff)
         {
-          R->getMonoid()->to_expvector(f->monom, exp2);
+          R->getMonoid()->to_expvector(f.monom, exp2);
           for (int k = 0; k < n; k++)
             if (exp2[k] != 0 && exp[k] == 0)
             {
@@ -1718,10 +1718,10 @@ Matrix /* or null */ *Matrix::monomials(M2_arrayint vars) const
       vec v = elem(c);
       for (; v != nullptr; v = v->next)
         {
-          for (Nterm *t = v->coeff; t != nullptr; t = t->next)
+          for (Nterm& t : v->coeff)
             {
               exponents_t exp1 = newarray_atomic(int, vars->len + 1);
-              M->to_expvector(t->monom, exp);
+              M->to_expvector(t.monom, exp);
               for (unsigned int i = 0; i < vars->len; i++)
                 exp1[i] = exp[vars->array[i]];
               exp1[vars->len] = v->comp;
@@ -1798,15 +1798,15 @@ static vec coeffs_of_vec(exponent_table *E,
   vec result = nullptr;
   for (vec g = f; g != nullptr; g = g->next)
     {
-      for (Nterm *h = g->coeff; h != nullptr; h = h->next)
+      for (Nterm& h : g->coeff)
         {
-          M->to_expvector(h->monom, exp);
+          M->to_expvector(h.monom, exp);
           get_part_of_expvector(vars, exp, g->comp, scratch_exp);
           int val = static_cast<int>(exponent_table_get(E, scratch_exp));
           if (val > 0)
             {
               M->from_expvector(exp, mon);
-              ring_elem t = P->make_flat_term(h->coeff, mon);
+              ring_elem t = P->make_flat_term(h.coeff, mon);
               vec v = P->make_vec(val - 1, t);
               v->next = result;
               result = v;
diff --git a/M2/Macaulay2/e/poly.cpp b/M2/Macaulay2/e/poly.cpp
index 419ab09aeec..1140c01bc61 100644
--- a/M2/Macaulay2/e/poly.cpp
+++ b/M2/Macaulay2/e/poly.cpp
@@ -323,16 +323,16 @@ int PolyRing::index_of_var(const ring_elem a) const
   return result;
 }
 
-M2_arrayint PolyRing::support(const ring_elem a) const
+// TODO: remove only M2_arrayint in this file
+M2_arrayint PolyRing::support(const ring_elem f) const
 {
   exponents_t EXP1 = ALLOCATE_EXPONENTS(exp_size);
   exponents_t EXP2 = ALLOCATE_EXPONENTS(exp_size);
   for (int i = 0; i < n_vars(); i++) EXP1[i] = 0;
-  for (const Nterm *f = a; f != 0; f = f->next)
+  for (Nterm& t : f)
     {
-      M_->to_expvector(f->monom, EXP2);
-      for (int j = 0; j < n_vars(); j++)
-        if (EXP2[j] != 0) EXP1[j] = 1;
+      M_->to_expvector(t.monom, EXP2);
+      exponents::lcm(n_vars(), EXP1, EXP2, EXP1);
     }
   int nelems = 0;
   for (int i = 0; i < n_vars(); i++)
@@ -444,8 +444,8 @@ ring_elem PolyRing::preferred_associate(ring_elem ff) const
   return t;
 }
 
-ring_elem PolyRing::preferred_associate_divisor(ring_elem ff) const
-// ff is an element of 'this'.
+ring_elem PolyRing::preferred_associate_divisor(ring_elem f) const
+// f is an element of 'this'.
 // result is in the coefficient ring
 // If the coefficient ring of this is
 //   ZZ -- gcd of all, same sign as lead coeff
@@ -456,20 +456,18 @@ ring_elem PolyRing::preferred_associate_divisor(ring_elem ff) const
 {
   ring_elem result = getCoefficients()->zero();
   assert(getCoefficients()->has_associate_divisors());
-  for (Nterm *f = ff; f != NULL; f = f->next)
-    {
-      if (!getCoefficients()->lower_associate_divisor(result, f->coeff))
-        // ie it cannot change, no matter what next coeff is
-        return result;
-    }
+  for (Nterm& t : f)
+    if (!getCoefficients()->lower_associate_divisor(result, t.coeff))
+      // ie it cannot change, no matter what next coeff is
+      return result;
   return result;
 }
 
-bool PolyRing::is_unit(const ring_elem ff) const
+bool PolyRing::is_unit(const ring_elem f) const
 {
-  Nterm *f = ff;
-  if (f == NULL) return false;
-  if (f->next == NULL && M_->is_one(f->monom) && K_->is_unit(f->coeff))
+  Nterm *t = f;
+  if (begin(f) == end(f)) return false;
+  if (t->next == NULL && M_->is_one(t->monom) && K_->is_unit(t->coeff))
     return true;
   return false;
 }
@@ -523,22 +521,19 @@ int PolyRing::compare_elems(const ring_elem f, const ring_elem g) const
 bool PolyRing::is_homogeneous(const ring_elem f) const
 {
   if (!is_graded()) return false;
-  Nterm *t = f;
-  if (t == 0) return true;
-  bool result = true;
-  monomial e = degree_monoid()->make_one();
-  monomial degf = degree_monoid()->make_one();
-  M_->multi_degree(t->monom, degf);
-  for (t = t->next; t != NULL; t = t->next)
-    {
-      M_->multi_degree(t->monom, e);
-      if (0 != degree_monoid()->compare(degf, e))
-        {
-          result = false;
-          break;
-        }
+  if (begin(f) == end(f)) return true;
+  auto DM = degree_monoid();
+  auto degf = DM->make_one();
+  Nterm& t = *begin(f);
+  M_->multi_degree(t.monom, degf);
+  for (Nterm& t : f)
+    {
+      auto e = DM->make_one();
+      M_->multi_degree(t.monom, e);
+      if (EQ != DM->compare(degf, e))
+        return false;
     }
-  return result;
+  return true;
 }
 
 void PolyRing::degree(const ring_elem f, monomial degf) const
@@ -548,25 +543,28 @@ void PolyRing::degree(const ring_elem f, monomial degf) const
 
 bool PolyRing::multi_degree(const ring_elem f, monomial degf) const
 {
-  Nterm *t = f;
-  monomial e = degree_monoid()->make_one();
-  if (t == 0 || M_->n_vars() == 0)
+  auto DM = degree_monoid();
+  if (begin(f) == end(f) || M_->n_vars() == 0)
     {
-      degree_monoid()->one(degf);
+      DM->one(degf);
       return true;
     }
-  M_->multi_degree(t->monom, degf);
   bool result = true;
-  for (t = t->next; t != NULL; t = t->next)
+  Nterm& t = *begin(f);
+  M_->multi_degree(t.monom, degf);
+  for (Nterm& t : f)
     {
-      M_->multi_degree(t->monom, e);
-      if (0 != degree_monoid()->compare(degf, e))
+      auto e = DM->make_one();
+      M_->multi_degree(t.monom, e);
+      if (EQ != DM->compare(degf, e))
         {
           result = false;
-          ;
-          degree_monoid()->lcm(degf, e, degf);
+          DM->lcm(degf, e, degf);
         }
     }
+  // for (int i = 0; i < DM->monomial_size(); i++)
+  //   std::cout << degf[i] << " , ";
+  // std::cout << std::endl;
   return result;
 }
 
@@ -601,15 +599,13 @@ ring_elem PolyRing::homogenize(const ring_elem f,
   // assert(wts[v] != 0);
   // If an error occurs, then return 0, and set gError.
 
-  monomial exp = newarray_atomic(int, nvars_);
-  int maxlen = (wts->len < nvars_ ? wts->len : nvars_);
+  exponents_t exp = newarray_atomic(int, nvars_);
   Nterm head;
   Nterm *result = &head;
-  for (Nterm *a = f; a != NULL; a = a->next)
+  for (Nterm& a : f)
     {
-      M_->to_expvector(a->monom, exp);
-      int e = 0;
-      for (int i = 0; i < maxlen; i++) e += wts->array[i] * exp[i];
+      M_->to_expvector(a.monom, exp);
+      auto e = exponents::weight(n_vars(), exp, wts);
       if (((d - e) % wts->array[v]) != 0)
         {
           // We cannot homogenize, so clean up and exit.
@@ -622,7 +618,7 @@ ring_elem PolyRing::homogenize(const ring_elem f,
       if (is_skew_ && skew_.is_skew_var(v) && exp[v] > 1) continue;
       result->next = new_term();
       result = result->next;
-      result->coeff = K_->copy(a->coeff);
+      result->coeff = K_->copy(a.coeff);
       M_->from_expvector(exp, result->monom);
     }
   result->next = NULL;
@@ -744,12 +740,12 @@ ring_elem PolyRing::negate(const ring_elem f) const
 {
   Nterm head;
   Nterm *result = &head;
-  for (Nterm *a = f; a != NULL; a = a->next)
+  for (Nterm& a : f)
     {
       result->next = new_term();
       result = result->next;
-      result->coeff = K_->negate(a->coeff);
-      M_->copy(a->monom, result->monom);
+      result->coeff = K_->negate(a.coeff);
+      M_->copy(a.monom, result->monom);
     }
   result->next = NULL;
   return head.next;
@@ -778,12 +774,12 @@ ring_elem PolyRing::mult_by_term(const ring_elem f,
 {
   Nterm head;
   Nterm *result = &head;
-  for (Nterm *a = f; a != NULL; a = a->next)
+  for (Nterm& a : f)
     {
       result->next = new_term();
       result = result->next;
-      result->coeff = K_->mult(a->coeff, c);
-      M_->mult(m, a->monom, result->monom);
+      result->coeff = K_->mult(a.coeff, c);
+      M_->mult(m, a.monom, result->monom);
     }
   result->next = NULL;
   return head.next;
@@ -814,11 +810,8 @@ void PolyRing::divide_coeff_to(ring_elem &f, ring_elem a) const
 ring_elem PolyRing::mult(const ring_elem f, const ring_elem g) const
 {
   polyheap H(this);
-  for (Nterm *a = f; a != NULL; a = a->next)
-    {
-      ring_elem h = mult_by_term(g, a->coeff, a->monom);
-      H.add(h);
-    }
+  for (Nterm& a : f)
+    H.add(mult_by_term(g, a.coeff, a.monom));
   return H.value();
 }
 
@@ -1346,17 +1339,14 @@ ring_elem PolyRing::eval(const RingMap *map,
   // The way we collect the result depends on whether the target ring
   // is a polynomial ring: if so, use a heap structure.  If not, just add to the
   // result.
-
-  gc_vector<int> vp;
   const Ring *target = map->get_ring();
   SumCollector *H = target->make_SumCollector();
 
-  for (Nterm *t = f; t != NULL; t = t->next)
+  for (Nterm& t : f)
     {
-      vp.resize(0);
-      M_->to_varpower(t->monom, vp);
-      ring_elem g = map->eval_term(K_, t->coeff, vp.data(), first_var, n_vars());
-      H->add(g);
+      gc_vector<int> vp;
+      M_->to_varpower(t.monom, vp);
+      H->add(map->eval_term(K_, t.coeff, vp.data(), first_var, n_vars()));
     }
   ring_elem result = H->getValue();
   delete H;
@@ -1367,15 +1357,15 @@ ring_elem PolyRing::zeroize_tiny(gmp_RR epsilon, const ring_elem f) const
 {
   Nterm head;
   Nterm *result = &head;
-  for (Nterm *a = f; a != NULL; a = a->next)
+  for (Nterm& a : f)
     {
-      ring_elem c = K_->zeroize_tiny(epsilon, a->coeff);
+      ring_elem c = K_->zeroize_tiny(epsilon, a.coeff);
       if (!K_->is_zero(c))
         {
           result->next = new_term();
           result = result->next;
           result->coeff = c;
-          M_->copy(a->monom, result->monom);
+          M_->copy(a.monom, result->monom);
         }
     }
   result->next = NULL;
@@ -1383,8 +1373,7 @@ ring_elem PolyRing::zeroize_tiny(gmp_RR epsilon, const ring_elem f) const
 }
 void PolyRing::increase_maxnorm(gmp_RRmutable norm, const ring_elem f) const
 {
-  for (Nterm *a = f; a != NULL; a = a->next)
-    K_->increase_maxnorm(norm, a->coeff);
+  for (Nterm& a : f) K_->increase_maxnorm(norm, a.coeff);
 }
 
 /////////////////////////////////////////
@@ -1861,13 +1850,13 @@ ring_elem PolyRing::get_part(const M2_arrayint wts,
 
   exponents_t exp = newarray_atomic(int, M_->n_vars());
 
-  for (Nterm *t = f; t != 0; t = t->next)
+  for (Nterm& t : f)
     {
-      M_->to_expvector(t->monom, exp);
+      M_->to_expvector(t.monom, exp);
       long wt = exponents::weight(M_->n_vars(), exp, wts);
       if (lobound_given && wt < lobound) continue;
       if (hibound_given && wt > hibound) continue;
-      inresult->next = copy_term(t);
+      inresult->next = copy_term(&t);
       inresult = inresult->next;
     }
 
@@ -1900,13 +1889,13 @@ ring_elem PolyRing::make_logical_term(const Ring *coeffR,
   Nterm *inresult = &head;
   exponents_t exp = newarray_atomic(int, M_->n_vars());
   exponents::copy(nvars0, exp0, exp);  // Sets the first part of exp
-  for (Nterm *f = a; f != 0; f = f->next)
+  for (Nterm& f : a)
     {
       Nterm *t = new_term();
       inresult->next = t;
       inresult = t;
-      t->coeff = f->coeff;
-      logicalK->getMonoid()->to_expvector(f->monom, exp + nvars0);
+      t->coeff = f.coeff;
+      logicalK->getMonoid()->to_expvector(f.monom, exp + nvars0);
       M_->from_expvector(exp, t->monom);
     }
   inresult->next = 0;
@@ -1955,7 +1944,7 @@ ring_elem PolyRing::get_terms(int nvars0,
 int PolyRing::n_flat_terms(const ring_elem f) const
 {
   int result = 0;
-  for (Nterm *a = f; a != NULL; a = a->next) result++;
+  for ([[maybe_unused]] Nterm& a : f) result++;
   return result;
 }
 
@@ -2012,15 +2001,15 @@ ring_elem PolyRing::diff(ring_elem a, ring_elem b, int use_coeff) const
 {
   polyheap H(this);
   Nterm *d = new_term();
-  for (Nterm *s = a; s != 0; s = s->next)
+  for (Nterm& s : a)
     {
-      for (Nterm *t = b; t != 0; t = t->next)
+      for (Nterm& t : b)
         {
-          d->coeff = diff_term(s->monom, t->monom, d->monom, use_coeff);
+          d->coeff = diff_term(s.monom, t.monom, d->monom, use_coeff);
           if (!K_->is_zero(d->coeff))
             {
-              K_->mult_to(d->coeff, s->coeff);
-              K_->mult_to(d->coeff, t->coeff);
+              K_->mult_to(d->coeff, s.coeff);
+              K_->mult_to(d->coeff, t.coeff);
               d->next = 0;
               H.add(d);
               d = new_term();
@@ -2104,8 +2093,8 @@ void PolyRing::sort(Nterm *&f) const
 
 bool PolyRing::in_subring(int nslots, const ring_elem a) const
 {
-  for (Nterm *t = a; t != 0; t = t->next)
-    if (!M_->in_subring(nslots, t->monom)) return false;
+  for (Nterm& t : a)
+    if (!M_->in_subring(nslots, t.monom)) return false;
   return true;
 }
 
@@ -2136,9 +2125,9 @@ void PolyRing::monomial_divisor(const ring_elem a, exponents_t exp) const
 // monomials of 'a'.
 {
   exponents_t exp1 = newarray_atomic(int, n_vars());
-  for (const Nterm *t = a; t != 0; t = t->next)
+  for (Nterm& t : a)
     {
-      M_->to_expvector(t->monom, exp1);
+      M_->to_expvector(t.monom, exp1);
       exponents::gcd(n_vars(), exp1, exp, exp);
     }
 }
@@ -2152,16 +2141,16 @@ ring_elem PolyRing::divide_by_var(int n, int d, const ring_elem a) const
   Nterm head;
   Nterm *result = &head;
   exponents_t exp = newarray_atomic(int, n_vars());
-  for (Nterm *t = a; t != 0; t = t->next)
+  for (Nterm& t : a)
     {
-      M_->to_expvector(t->monom, exp);
+      M_->to_expvector(t.monom, exp);
       if (exp[n] >= d)
         exp[n] -= d;
       else
         continue;
       result->next = new_term();
       result = result->next;
-      result->coeff = t->coeff;
+      result->coeff = t.coeff;
       M_->from_expvector(exp, result->monom);
     }
   freemem(exp);
@@ -2173,12 +2162,12 @@ ring_elem PolyRing::divide_by_expvector(const_exponents exp, const ring_elem a)
 {
   Nterm *result = 0;
   exponents_t exp0 = newarray_atomic(int, n_vars());
-  for (Nterm *t = a; t != 0; t = t->next)
+  for (Nterm& t : a)
     {
-      M_->to_expvector(t->monom, exp0);
+      M_->to_expvector(t.monom, exp0);
       exponents::quotient(n_vars(), exp0, exp, exp0);
       Nterm *u = new_term();
-      u->coeff = t->coeff;
+      u->coeff = t.coeff;
       M_->from_expvector(exp0, u->monom);
       u->next = result;
       result = u;
@@ -2191,7 +2180,7 @@ ring_elem PolyRing::divide_by_expvector(const_exponents exp, const ring_elem a)
 void PolyRing::lower_content(ring_elem &c, ring_elem g) const
 // c is a content elem, g is in ring
 {
-  for (Nterm *t = g; t != 0; t = t->next) K_->lower_content(c, t->coeff);
+  for (Nterm& t : g) K_->lower_content(c, t.coeff);
 }
 
 ring_elem PolyRing::content(ring_elem f) const
@@ -2277,13 +2266,13 @@ ring_elem PolyRing::lead_term(int nparts, const ring_elem f) const
   Nterm head;
   Nterm *result = &head;
   int nslots = M_->n_slots(nparts);
-  for (Nterm *a = f; a != NULL; a = a->next)
+  for (Nterm& a : f)
     {
-      if (M_->compare(nslots, lead->monom, a->monom) != EQ) break;
+      if (M_->compare(nslots, lead->monom, a.monom) != EQ) break;
       result->next = new_term();
       result = result->next;
-      result->coeff = a->coeff;
-      M_->copy(a->monom, result->monom);
+      result->coeff = a.coeff;
+      M_->copy(a.monom, result->monom);
     }
   result->next = NULL;
   return head.next;
@@ -2349,14 +2338,14 @@ vec PolyRing::vec_coefficient_of_var(vec v, int x, int e) const
     {
       Nterm head;
       Nterm *result = &head;
-      for (Nterm *f = t->coeff; f != NULL; f = f->next)
+      for (Nterm& f : t->coeff)
         {
-          M_->to_expvector(f->monom, exp);
+          M_->to_expvector(f.monom, exp);
           if (exp[x] != e) continue;
           exp[x] = 0;
           result->next = new_term();
           result = result->next;
-          result->coeff = f->coeff;
+          result->coeff = f.coeff;
           M_->from_expvector(exp, result->monom);
         }
       result->next = NULL;
@@ -2386,9 +2375,9 @@ vec PolyRing::vec_top_coefficient(const vec v, int &x, int &e) const
 
   exponents_t exp = newarray_atomic(int, n_vars());
   for (vec t = v; t != 0; t = t->next)
-    for (Nterm *f = t->coeff; f != 0; f = f->next)
+    for (Nterm& f : t->coeff)
       {
-        M_->to_expvector(f->monom, exp);
+        M_->to_expvector(f.monom, exp);
         for (int i = 0; i < x; i++)
           {
             if (exp[i] > 0)
@@ -2418,9 +2407,9 @@ void PolyRing::determine_common_denominator_QQ(ring_elem f,
   // We assume that 'this' is QQ[M].
   if (getCoefficients() != globalQQ) return;
 
-  for (Nterm *t = f; t != 0; t = t->next)
+  for (Nterm& t : f)
     {
-      mpq_srcptr a = MPQ_VAL(t->coeff);
+      mpq_srcptr a = MPQ_VAL(t.coeff);
       mpz_lcm(denom_so_far, denom_so_far, mpq_denref(a));
     }
 }
@@ -2470,11 +2459,11 @@ gbvector *PolyRing::translate_gbvector_from_ringelem_QQ(ring_elem coeff) const
   GBRing *GR = get_gb_ring();
   gbvector head;
   gbvector *inresult = &head;
-  for (Nterm *t = coeff; t != 0; t = t->next)
+  for (Nterm& t : coeff)
     {
       // make a gbvector node.
-      ring_elem a = globalZZ->RingZZ::from_int(mpq_numref(MPQ_VAL(t->coeff)));
-      gbvector *g = GR->gbvector_term(0, a, t->monom, 0);
+      ring_elem a = globalZZ->RingZZ::from_int(mpq_numref(MPQ_VAL(t.coeff)));
+      gbvector *g = GR->gbvector_term(0, a, t.monom, 0);
       inresult->next = g;
       inresult = inresult->next;
     }
@@ -2488,10 +2477,10 @@ gbvector *PolyRing::translate_gbvector_from_ringelem(ring_elem coeff) const
   GBRing *GR = get_gb_ring();
   gbvector head;
   gbvector *inresult = &head;
-  for (Nterm *t = coeff; t != 0; t = t->next)
+  for (Nterm& t : coeff)
     {
       // make a gbvector node.
-      gbvector *g = GR->gbvector_term(0, t->coeff, t->monom, 0);
+      gbvector *g = GR->gbvector_term(0, t.coeff, t.monom, 0);
       inresult->next = g;
       inresult = inresult->next;
     }
@@ -2519,14 +2508,14 @@ gbvector *PolyRing::translate_gbvector_from_vec_QQ(
     {
       inresult = &head;
       int comp = w->comp + 1;
-      for (Nterm *t = w->coeff; t != 0; t = t->next)
+      for (Nterm& t : w->coeff)
         {
           // make a gbvector node.
-          mpq_srcptr b = MPQ_VAL(t->coeff);
+          mpq_srcptr b = MPQ_VAL(t.coeff);
           mpz_mul(a, result_denominator.get_mpz(), mpq_numref(b));
           mpz_divexact(a, a, mpq_denref(b));
           gbvector *g = GR->gbvector_term(
-              F, globalZZ->RingZZ::from_int(a), t->monom, comp);
+              F, globalZZ->RingZZ::from_int(a), t.monom, comp);
           inresult->next = g;
           inresult = inresult->next;
         }
@@ -2553,10 +2542,10 @@ gbvector *PolyRing::translate_gbvector_from_vec(
     {
       inresult = &head;
       int comp = w->comp + 1;
-      for (Nterm *t = w->coeff; t != 0; t = t->next)
+      for (Nterm& t : w->coeff)
         {
           // make a gbvector node.
-          gbvector *g = GR->gbvector_term(F, t->coeff, t->monom, comp);
+          gbvector *g = GR->gbvector_term(F, t.coeff, t.monom, comp);
           inresult->next = g;
           inresult = inresult->next;
         }
diff --git a/M2/Macaulay2/e/polyring.cpp b/M2/Macaulay2/e/polyring.cpp
index f0bf8c82a53..2f0f2a99bc8 100644
--- a/M2/Macaulay2/e/polyring.cpp
+++ b/M2/Macaulay2/e/polyring.cpp
@@ -203,10 +203,10 @@ unsigned int PolynomialRing::computeHashValue(const ring_elem a) const
   unsigned int hash = 0;
   unsigned int seed1 = 103;
   unsigned int seed2 = 347654;
-  for (const Nterm *t = a.poly_val; t != 0; t = t->next)
+  for (Nterm& t : a.poly_val)
     {
-      unsigned int hash1 = getCoefficientRing()->computeHashValue(t->coeff);
-      unsigned int hash2 = getMonoid()->computeHashValue(t->monom);
+      unsigned int hash1 = getCoefficientRing()->computeHashValue(t.coeff);
+      unsigned int hash2 = getMonoid()->computeHashValue(t.monom);
       hash += seed1 * hash1 + seed2 * hash2;
       seed1 += 463633;
       seed2 += 7858565;
diff --git a/M2/Macaulay2/e/relem.cpp b/M2/Macaulay2/e/relem.cpp
index 3de4727bbeb..993aeb61434 100644
--- a/M2/Macaulay2/e/relem.cpp
+++ b/M2/Macaulay2/e/relem.cpp
@@ -455,6 +455,7 @@ bool RingElement::getSmallIntegerCoefficients(
     std::vector<long> &result_coeffs) const
 {
   const PolynomialRing *R = get_ring()->cast_to_PolynomialRing();
+  const Ring *K = R->getCoefficientRing();
   if (R == 0 || R->n_vars() != 1)
     {
       throw exc::engine_error(
@@ -472,10 +473,9 @@ bool RingElement::getSmallIntegerCoefficients(
   result_coeffs.resize(deg + 1);
   for (int i = 0; i <= deg; i++) result_coeffs[i] = 0;
   int exp[1];
-  for (Nterm *t = get_value(); t != NULL; t = t->next)
+  for (Nterm& t : get_value())
     {
-      std::pair<bool, long> res =
-          R->getCoefficientRing()->coerceToLongInteger(t->coeff);
+      std::pair<bool, long> res = K->coerceToLongInteger(t.coeff);
       if (not res.first)
         {
           // At this point, the answer is meaningless
@@ -484,7 +484,7 @@ bool RingElement::getSmallIntegerCoefficients(
         }
       long coeff = res.second;
 
-      R->getMonoid()->to_expvector(t->monom, exp);
+      R->getMonoid()->to_expvector(t.monom, exp);
       assert(exp[0] >= 0);
       assert(exp[0] <= deg);
       result_coeffs[exp[0]] = coeff;
diff --git a/M2/Macaulay2/e/res-a0-poly.cpp b/M2/Macaulay2/e/res-a0-poly.cpp
index 8727408dbfa..cf588f3868e 100644
--- a/M2/Macaulay2/e/res-a0-poly.cpp
+++ b/M2/Macaulay2/e/res-a0-poly.cpp
@@ -124,13 +124,13 @@ res2term *res2_poly::ring_mult_by_term(const ring_elem f,
 {
   res2term head;
   res2term *result = &head;
-  for (const Nterm *tm = f; tm != NULL; tm = tm->next)
+  for (Nterm& tm : f)
     {
       result->next = new_term();
       result = result->next;
       result->comp = x;
-      result->coeff = K->mult(c, tm->coeff);
-      M->mult(tm->monom, m, result->monom);
+      result->coeff = K->mult(c, tm.coeff);
+      M->mult(tm.monom, m, result->monom);
     }
   result->next = NULL;
   return head.next;
@@ -299,13 +299,13 @@ res2term *res2_poly::from_vector(const VECTOR(res2_pair *)& base,
   res2term *result = &head;
 
   for (vecterm *w = v; w != NULL; w = w->next)
-    for (Nterm *t = w->coeff; t != 0; t = t->next)
+    for (Nterm& t : w->coeff)
       {
         result->next = new_term();
         result = result->next;
         result->comp = base[w->comp];
-        result->coeff = t->coeff;
-        M->copy(t->monom, result->monom);
+        result->coeff = t.coeff;
+        M->copy(t.monom, result->monom);
         M->mult(result->monom, result->comp->syz->monom, result->monom);
       }
   result->next = NULL;
diff --git a/M2/Macaulay2/e/res-a1-poly.cpp b/M2/Macaulay2/e/res-a1-poly.cpp
index 87bf4308891..3d5cf9b1ac3 100644
--- a/M2/Macaulay2/e/res-a1-poly.cpp
+++ b/M2/Macaulay2/e/res-a1-poly.cpp
@@ -107,13 +107,13 @@ resterm *res_poly::ring_mult_by_term(const ring_elem f,
 {
   resterm head;
   resterm *result = &head;
-  for (const Nterm *tm = f; tm != NULL; tm = tm->next)
+  for (Nterm& tm : f)
     {
       result->next = new_term();
       result = result->next;
       result->comp = x;
-      result->coeff = K->mult(c, tm->coeff);
-      M->mult(tm->monom, m, result->monom);
+      result->coeff = K->mult(c, tm.coeff);
+      M->mult(tm.monom, m, result->monom);
     }
   result->next = NULL;
   return head.next;
@@ -281,13 +281,13 @@ resterm *res_poly::from_vector(const VECTOR(res_pair *)& base, const vec v) cons
   resterm head;
   resterm *result = &head;
   for (vecterm *w = v; w != NULL; w = w->next)
-    for (Nterm *t = w->coeff; t != 0; t = t->next)
+    for (Nterm& t : w->coeff)
       {
         result->next = new_term();
         result = result->next;
         result->comp = base[w->comp];
-        result->coeff = t->coeff;
-        M->copy(t->monom, result->monom);
+        result->coeff = t.coeff;
+        M->copy(t.monom, result->monom);
         M->mult(result->monom, result->comp->base_monom, result->monom);
       }
   result->next = NULL;
diff --git a/M2/Macaulay2/e/ringelem.cpp b/M2/Macaulay2/e/ringelem.cpp
new file mode 100644
index 00000000000..e12b608494a
--- /dev/null
+++ b/M2/Macaulay2/e/ringelem.cpp
@@ -0,0 +1,8 @@
+#include <ringelem.hpp>
+
+// TODO: don't use NULL or 0 to check end of Nterm
+TermIterator<Nterm> begin(Nterm* ptr) { return TermIterator<Nterm>(ptr); }
+TermIterator<Nterm> end(Nterm*)       { return TermIterator<Nterm>{}; }
+
+TermIterator<vecterm> begin(vecterm* ptr) { return TermIterator<vecterm>(ptr); }
+TermIterator<vecterm> end(vecterm*)       { return TermIterator<vecterm>{}; }
diff --git a/M2/Macaulay2/e/ringelem.hpp b/M2/Macaulay2/e/ringelem.hpp
index 38f45bfc482..1153ec1e319 100644
--- a/M2/Macaulay2/e/ringelem.hpp
+++ b/M2/Macaulay2/e/ringelem.hpp
@@ -3,12 +3,9 @@
 #ifndef _ringelem_hh_
 #define _ringelem_hh_
 
-#include <stddef.h>
-#if !defined(SAFEC_EXPORTS)
-#include <engine-exports.h>
-#endif
-#include "interface/gmp-util.h"
-#include "newdelete.hpp"
+#include "M2/math-include.h"  // for mpfi_srcptr, mpfr_srcptr, mpq_srcptr
+#include "monoid.hpp"         // for monomial
+#include "newdelete.hpp"      // for our_new_delete
 
 using ZZ = mpz_srcptr;
 using ZZmutable = mpz_ptr;
@@ -94,14 +91,18 @@ union ring_elem
   const schur_poly* get_schur_poly() const { return schur_poly_val; }
 };
 
+/* Implements a linked list of ring monomials along with coefficients */
 struct Nterm
 {
   Nterm *next;
   ring_elem coeff;
+  // TODO: should this have type monomial?
   int monom[1];
 };
 
 typedef struct vecterm *vec;
+/* Implements a linked list of module monomials along with coefficients */
+// TODO: why is this garbage collected?
 struct vecterm : public our_new_delete
 {
   vec next;
@@ -109,6 +110,37 @@ struct vecterm : public our_new_delete
   ring_elem coeff;
 };
 
+/* Implements an iterator for linked list-based multi-termed structs
+ *
+ * For example, the functions begin(Nterm*) and end(Nterm*) return
+ * a TermIterator<Nterm> object which makes for(Nterm& t : f) work. */
+template<typename T>
+struct TermIterator
+{
+  T* p;
+
+  TermIterator():              p(nullptr) {}
+  TermIterator(T* ptr):        p(ptr)     {}
+  TermIterator(ring_elem ptr): p(ptr)     {}
+
+  TermIterator& operator++() { p = p->next; return *this; }
+
+  T const& operator*() const  { return *p; }
+  T&       operator*()        { return *p; }
+  T const* operator->() const { return p; }
+  T*       operator->()       { return p; }
+
+  bool operator==(TermIterator const& rhs) const { return p == rhs.p; }
+  bool operator!=(TermIterator const& rhs) const { return p != rhs.p; }
+};
+
+TermIterator<Nterm> begin(Nterm* ptr);
+TermIterator<Nterm> end(Nterm*);
+
+TermIterator<vecterm> begin(vecterm* ptr);
+TermIterator<vecterm> end(vecterm*);
+
+
 #define MPQ_VAL(f) ((f).get_mpq())
 
 // these should only be used as temporary const.  Do not store results!
diff --git a/M2/Macaulay2/e/schur.cpp b/M2/Macaulay2/e/schur.cpp
index 854c8d5dd72..4259777da58 100644
--- a/M2/Macaulay2/e/schur.cpp
+++ b/M2/Macaulay2/e/schur.cpp
@@ -266,14 +266,14 @@ ring_elem SchurRing::mult_by_term(const ring_elem f,
 {
   // return c*m*f
   ring_elem result = ZERO_RINGELEM;
-  for (Nterm *t = f; t != NULL; t = t->next)
+  for (Nterm& t : f)
     {
-      ring_elem a = K_->mult(c, t->coeff);
-      ring_elem g = const_cast<SchurRing *>(this)->mult_monomials(t->monom, m);
-      for (Nterm *s = g; s != NULL; s = s->next)
+      ring_elem a = K_->mult(c, t.coeff);
+      ring_elem g = const_cast<SchurRing *>(this)->mult_monomials(t.monom, m);
+      for (Nterm& s : g)
         {
-          ring_elem b = K_->mult(a, s->coeff);
-          s->coeff = b;
+          ring_elem b = K_->mult(a, s.coeff);
+          s.coeff = b;
         }
       Nterm *gt = g;
       sort(gt);
@@ -332,12 +332,12 @@ ring_elem SchurRing::dimension(const ring_elem f) const
   ring_elem result = K_->from_long(0);
   mpz_t dim;
   mpz_init(dim);
-  for (Nterm *t = f; t != NULL; t = t->next)
+  for (Nterm& t : f)
     {
-      to_partition(t->monom, EXP);
+      to_partition(t.monom, EXP);
       dimension(EXP, dim);
       ring_elem h = K_->from_int(dim);
-      ring_elem h2 = K_->mult(t->coeff, h);
+      ring_elem h2 = K_->mult(t.coeff, h);
       K_->add_to(result, h2);
       K_->remove(h);
     }
@@ -363,12 +363,12 @@ void SchurRing::elem_text_out(buffer &o,
     }
 
   p_one = false;
-  for (Nterm *t = f; t != NULL; t = t->next)
+  for (Nterm& t : f)
     {
-      int isone = M_->is_one(t->monom);
+      int isone = M_->is_one(t.monom);
       p_parens = !isone;
-      K_->elem_text_out(o, t->coeff, p_one, p_plus, p_parens);
-      to_partition(t->monom, EXP);
+      K_->elem_text_out(o, t.coeff, p_one, p_plus, p_parens);
+      to_partition(t.monom, EXP);
       o << "{" << EXP[1];
       for (int i = 2; i <= nvars_ && EXP[i] != 0; i++) o << "," << EXP[i];
       o << "}";
diff --git a/M2/Macaulay2/e/skewpoly.cpp b/M2/Macaulay2/e/skewpoly.cpp
index 0ec4caa58d9..d882c53c791 100644
--- a/M2/Macaulay2/e/skewpoly.cpp
+++ b/M2/Macaulay2/e/skewpoly.cpp
@@ -41,9 +41,9 @@ ring_elem SkewPolynomialRing::antipode(const ring_elem f) const
 
   exponents_t EXP = ALLOCATE_EXPONENTS(exp_size);
 
-  for (Nterm *s = f; s != NULL; s = s->next)
+  for (Nterm& s : f)
     {
-      M_->to_expvector(s->monom, EXP);
+      M_->to_expvector(s.monom, EXP);
       int deg = skew_.skew_degree(EXP);
       // sign is (-1)^ (binomial(deg,2))
       // deg = 0: sign is 1
@@ -55,8 +55,8 @@ ring_elem SkewPolynomialRing::antipode(const ring_elem f) const
       int sign = sign4[mod4];
       Nterm *t = new_term();
       t->next = 0;
-      t->coeff = (sign == 1 ? s->coeff : K_->negate(s->coeff));
-      M_->copy(s->monom, t->monom);
+      t->coeff = (sign == 1 ? s.coeff : K_->negate(s.coeff));
+      M_->copy(s.monom, t->monom);
       inresult->next = t;
       inresult = inresult->next;
     }
@@ -76,18 +76,18 @@ ring_elem SkewPolynomialRing::mult_by_term(const ring_elem f,
   exponents_t EXP2 = ALLOCATE_EXPONENTS(exp_size);
   M_->to_expvector(m, EXP1);
 
-  for (Nterm *s = f; s != NULL; s = s->next)
+  for (Nterm& s : f)
     {
-      M_->to_expvector(s->monom, EXP2);
+      M_->to_expvector(s.monom, EXP2);
       int sign = skew_.mult_sign(EXP1, EXP2);
       if (sign == 0) continue;
 
       Nterm *t = new_term();
       t->next = 0;
-      t->coeff = K_->mult(c, s->coeff);
+      t->coeff = K_->mult(c, s.coeff);
       if (sign < 0) K_->negate_to(t->coeff);
 
-      M_->mult(m, s->monom, t->monom);
+      M_->mult(m, s.monom, t->monom);
       inresult->next = t;
       inresult = inresult->next;
     }
diff --git a/M2/Macaulay2/e/tower.cpp b/M2/Macaulay2/e/tower.cpp
index 46d02469d13..34689737cdc 100644
--- a/M2/Macaulay2/e/tower.cpp
+++ b/M2/Macaulay2/e/tower.cpp
@@ -450,11 +450,10 @@ ring_elem Tower::translate(const PolynomialRing *R, ring_elem fR) const
   int nvars = R->n_vars();
   TowerPolynomial result = 0;
   exponents_t exp = new int[nvars];
-  for (Nterm *t = fR; t != 0; t = t->next)
+  for (Nterm& t : fR)
     {
-      M->to_expvector(t->monom, exp);
-
-      std::pair<bool, long> res = K->coerceToLongInteger(t->coeff);
+      M->to_expvector(t.monom, exp);
+      std::pair<bool, long> res = K->coerceToLongInteger(t.coeff);
       assert(res.first);
       int c1 = static_cast<int>(res.second);
 
diff --git a/M2/Macaulay2/e/weylalg.cpp b/M2/Macaulay2/e/weylalg.cpp
index 0f6c4e1b046..7f3046b55bb 100644
--- a/M2/Macaulay2/e/weylalg.cpp
+++ b/M2/Macaulay2/e/weylalg.cpp
@@ -377,12 +377,12 @@ Nterm *WeylAlgebra::weyl_diff(const ring_elem c,
         deriv_exp[_commutative[i]] = derivatives[i];
       }
 
-  for (const Nterm *t = g; t != 0; t = t->next)
+  for (Nterm& t : ring_elem(g))
     {
       // This first part checks whether the x-part of t->monom is divisible by
       // 'derivatives'.  If so, true is returned, and the resulting monomial is
       // set.
-      M_->to_expvector(t->monom, result_exp);
+      M_->to_expvector(t.monom, result_exp);
       extractCommutativePart(result_exp, exp);
       if (divides(derivatives, exp))
         {
@@ -392,7 +392,7 @@ Nterm *WeylAlgebra::weyl_diff(const ring_elem c,
               K_->remove(a);
               continue;
             }
-          ring_elem b = K_->mult(a, t->coeff);
+          ring_elem b = K_->mult(a, t.coeff);
           K_->remove(a);
           if (K_->is_zero(b))
             {