diff --git a/arrays.scad b/arrays.scad
index e4de121..7d85d8c 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -53,13 +53,12 @@ function _same_type(a,b, depth) =
 
 
 // Function: select()
+// Topics: List Handling
 // Description:
 //   Returns a portion of a list, wrapping around past the beginning, if end<start. 
 //   The first item is index 0. Negative indexes are counted back from the end.
 //   The last item is -1.  If only the `start` index is given, returns just the value
 //   at that position.
-// Topics: List Handling
-// See Also: slice(), subindex(), last()
 // Usage:
 //   item = select(list,start);
 //   list = select(list,start,end);
@@ -67,6 +66,7 @@ function _same_type(a,b, depth) =
 //   list = The list to get the portion of.
 //   start = The index of the first item.
 //   end = The index of the last item.
+// See Also: slice(), subindex(), last()
 // Example:
 //   l = [3,4,5,6,7,8,9];
 //   a = select(l, 5, 6);   // Returns [8,9]
@@ -97,33 +97,36 @@ function select(list, start, end) =
 
 
 // Function: slice()
-// Usage:
-//   list = slice(list,start,end);
 // Topics: List Handling
-// See Also: select(), subindex(), last()
+// Usage:
+//   list = slice(list,s,e);
 // Description:
-//   Returns a slice of a list.  The first item is index 0.
-//   Negative indexes are counted back from the end.  The last item is -1.
+//   Returns a slice of a list, from the first position `s` up to and including the last position `e`.
+//   The first item in the list is at index 0.  Negative indexes are counted back from the end.
+//   An index of -1 refers to the last list item.
 // Arguments:
-//   list = The array/list to get the slice of.
-//   start = The index of the first item to return.
-//   end = The index after the last item to return, unless negative, in which case the last item to return.
+//   list = The list to get the slice of.
+//   s = The index of the first item to return.
+//   e = The index of the last item to return.
+// See Also: select(), subindex(), last()
 // Example:
-//   a = slice([3,4,5,6,7,8,9], 3, 5);   // Returns [6,7]
+//   a = slice([3,4,5,6,7,8,9], 3, 5);   // Returns [6,7,8]
 //   b = slice([3,4,5,6,7,8,9], 2, -1);  // Returns [5,6,7,8,9]
-//   c = slice([3,4,5,6,7,8,9], 1, 1);   // Returns []
-//   d = slice([3,4,5,6,7,8,9], 6, -1);  // Returns [9]
+//   c = slice([3,4,5,6,7,8,9], 1, 1);   // Returns [4]
+//   d = slice([3,4,5,6,7,8,9], 5);      // Returns [8,9]
 //   e = slice([3,4,5,6,7,8,9], 2, -2);  // Returns [5,6,7,8]
-function slice(list,start,end) =
-    assert( is_list(list), "Invalid list" )
-    assert( is_finite(start) && is_finite(end), "Invalid number(s)" )
-    let( l = len(list) )
-    l==0
-      ? []
-      : let(
-            s = start<0? (l+start) : start,
-            e = end<0? (l+end+1) : end
-        ) [for (i=[s:1:e-1]) if (e>s) list[i]];
+//   f = slice([3,4,5,6,7,8,9], 4, 3;    // Returns []
+function slice(list,s=0,e=-1) =
+    assert(is_list(list))
+    assert(is_int(s))
+    assert(is_int(e))
+    !list? [] :
+    let(
+        l = len(list),
+        s = constrain(s + (s<0? l : 0), 0, l-1),
+        e = constrain(e + (e<0? l : 0), 0, l-1)
+    )
+    [if (e>=s) for (i=[s:1:e]) list[i]];
 
 
 // Function: last()
@@ -367,7 +370,7 @@ function list_decreasing(list) =
 // Usage:
 //   list = repeat(val, n);
 // Topics: List Handling
-// See Also: list_range()
+// See Also: range(), rangex()
 // Description:
 //   Generates a list or array of `n` copies of the given value `val`.
 //   If the count `n` is given as a list of counts, then this creates a
@@ -387,37 +390,36 @@ function repeat(val, n, i=0) =
     [for (j=[1:1:n[i]]) repeat(val, n, i+1)];
 
 
-// Function: list_range()
+// Function: range()
 // Usage:
-//   list = list_range(n=, <s=>, <e=>);
-//   list = list_range(n=, <s=>, <step=>);
-//   list = list_range(e=, <step=>);
-//   list = list_range(s=, e=, <step=>);
+//   list = range(n, <s=>, <e=>);
+//   list = range(n, <s=>, <step=>);
+//   list = range(e=, <step=>);
+//   list = range(s=, e=, <step=>);
 // Topics: List Handling
-// See Also: repeat()
 // Description:
 //   Returns a list, counting up from starting value `s`, by `step` increments,
 //   until either `n` values are in the list, or it reaches the end value `e`.
 //   If both `n` and `e` are given, returns `n` values evenly spread from `s`
 //   to `e`, and `step` is ignored.
 // Arguments:
-//   ---
 //   n = Desired number of values in returned list, if given.
+//   ---
 //   s = Starting value.  Default: 0
 //   e = Ending value to stop at, if given.
 //   step = Amount to increment each value.  Default: 1
+// See Also: repeat(), rangex()
 // Example:
-//   a = list_range(4);                  // Returns [0,1,2,3]
-//   b = list_range(n=4, step=2);        // Returns [0,2,4,6]
-//   c = list_range(n=4, s=3, step=3);   // Returns [3,6,9,12]
-//   d = list_range(n=5, s=0, e=10);     // Returns [0, 2.5, 5, 7.5, 10]
-//   e = list_range(e=3);                // Returns [0,1,2,3]
-//   f = list_range(e=7, step=2);        // Returns [0,2,4,6]
-//   g = list_range(s=3, e=5);           // Returns [3,4,5]
-//   h = list_range(s=3, e=8, step=2);   // Returns [3,5,7]
-//   i = list_range(s=4, e=8.3, step=2); // Returns [4,6,8]
-//   j = list_range(n=4, s=[3,4], step=[2,3]);  // Returns [[3,4], [5,7], [7,10], [9,13]]
-function list_range(n, s=0, e, step) =
+//   a = range(4);                  // Returns [0,1,2,3]
+//   b = range(n=4, step=2);        // Returns [0,2,4,6]
+//   c = range(n=4, s=3, step=3);   // Returns [3,6,9,12]
+//   d = range(n=5, s=0, e=10);     // Returns [0, 2.5, 5, 7.5, 10]
+//   e = range(e=3);                // Returns [0,1,2,3]
+//   f = range(e=7, step=2);        // Returns [0,2,4,6]
+//   g = range(s=3, e=5);           // Returns [3,4,5]
+//   h = range(s=3, e=8, step=2);   // Returns [3,5,7]
+//   i = range(s=4, e=8.3, step=2); // Returns [4,6,8]
+function range(n, s=0, e, step) =
     assert( is_undef(n) || is_finite(n), "Parameter `n` must be a number.")
     assert( is_undef(n) || is_undef(e) || is_undef(step), "At most 2 of n, e, and step can be given.")
     let( step = (n!=undef && e!=undef)? (e-s)/(n-1) : default(step,1) )
@@ -428,6 +430,50 @@ function list_range(n, s=0, e, step) =
         [for (v=[s:step:e]) v] ;
     
 
+// Function: rangex()
+// Usage:
+//   list = rangex(n, <s=>, <e=>);
+//   list = rangex(n, <s=>, <step=>);
+//   list = rangex(e=, <step=>);
+//   list = rangex(s=, e=, <step=>);
+// Topics: List Handling
+// Description:
+//   Returns a list, counting up from starting value `s`, by `step` increments, until
+//   either `n` values are in the list, or it reaches the value just before the end value `e`.
+//   If both `n` and `e` are given, returns `n` values evenly spread from `s` to the value
+//   just before `e`, and `step` is ignored.
+// Arguments:
+//   n = Desired number of values in returned list, if given.
+//   ---
+//   s = Starting value.  Default: 0
+//   e = Ending value to stop at, if given.
+//   step = Amount to increment each value.  Default: 1
+// See Also: repeat(), range()
+// Example:
+//   a = rangex(4);                  // Returns [0,1,2,3]
+//   b = rangex(5,e=1);              // Returns [0, 0.2, 0.4, 0.6, 0.8]
+//   c = rangex(n=4, step=2);        // Returns [0,2,4,6]
+//   d = rangex(n=4, step=0.25);     // Returns [0, 0.25, 0.5, 0.75]
+//   e = rangex(n=4, s=3, step=3);   // Returns [3,6,9,12]
+//   f = rangex(n=5, s=0, e=10);     // Returns [0, 2, 4, 6, 8]
+//   g = rangex(e=3);                // Returns [0,1,2]
+//   h = rangex(e=7, step=2);        // Returns [0,2,4,6]
+//   i = rangex(s=3, e=5);           // Returns [3,4]
+//   j = rangex(s=3, e=8, step=2);   // Returns [3,5,7]
+//   k = rangex(s=2, e=8, step=2);   // Returns [2,4,6]
+//   l = rangex(s=2, e=8.1, step=2); // Returns [2,4,6,8]
+function rangex(n, s=0, e, step) =
+    assert( is_undef(n) || is_finite(n), "Parameter `n` must be a number.")
+    assert( is_undef(n) || is_undef(e) || is_undef(step), "At most 2 of n, e, and step can be given.")
+    let( step = (n!=undef && e!=undef)? (e-s)/n : default(step,1) )
+    is_undef(e)
+      ? assert( is_consistent([s, step]), "Incompatible data.")
+        [for (i=[0:1:n-1]) s+step*i ]
+      : assert( is_vector([s,step,e]), "Start `s`, step `step` and end `e` must be numbers.")
+        let(steps=floor((e-s)/step+0.5))
+        [for (i=[0:1:steps-1]) s+step*i ];
+    
+
 
 // Section: List Manipulation
 
diff --git a/attachments.scad b/attachments.scad
index 04d8a17..0a283d1 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -442,7 +442,7 @@ function find_anchor(anchor, geom) =
     let(
         cp = select(geom,-3),
         offset = anchor==CENTER? CENTER : select(geom,-2),
-        anchors = select(geom,-1),
+        anchors = last(geom),
         type = geom[0]
     )
     is_string(anchor)? (
diff --git a/beziers.scad b/beziers.scad
index 79fe650..fb49adc 100644
--- a/beziers.scad
+++ b/beziers.scad
@@ -769,7 +769,7 @@ function path_to_bezier(path, closed=false, tangents, uniform=false, size, relsi
     )
     assert(min(sizevect)>0, "Size and relsize must be greater than zero")
     [
-        for(i=[0:lastpt-1])
+        for(i=[0:1:lastpt-1])
             let(
                 first = path[i],
                 second = select(path,i+1),
diff --git a/debug.scad b/debug.scad
index 6e4bdae..fccc6c9 100644
--- a/debug.scad
+++ b/debug.scad
@@ -377,7 +377,7 @@ module show_anchors(s=10, std=true, custom=true) {
         }
     }
     if (custom) {
-        for (anchor=select($parent_geom,-1)) {
+        for (anchor=last($parent_geom)) {
             attach(anchor[0]) {
                 if(two_d) {
                     anchor_arrow2d(s, color="cyan");
diff --git a/gears.scad b/gears.scad
index 5bff0fc..52ea9bd 100644
--- a/gears.scad
+++ b/gears.scad
@@ -960,11 +960,11 @@ function bevel_gear(
                 each apply(xflip() * zrot(360*tooth/teeth) * m, path3d(profile))
             ]
         ],
-        thickness = abs(verts1[0][0].z - select(verts1,-1)[0].z),
+        thickness = abs(verts1[0][0].z - last(verts1)[0].z),
         vertices = [for (x=verts1) down(thickness/2, p=reverse(x))],
         sides_vnf = vnf_vertex_array(vertices, caps=false, col_wrap=true, reverse=true),
-        top_verts = select(vertices,-1),
-        bot_verts = select(vertices,0),
+        top_verts = last(vertices),
+        bot_verts = vertices[0],
         gear_pts = len(top_verts),
         face_pts = gear_pts / teeth,
         top_faces =[
@@ -1431,8 +1431,8 @@ function worm_gear(
                 )
             ]
         ],
-        top_verts = select(profiles,-1),
-        bot_verts = select(profiles,0),
+        top_verts = last(profiles),
+        bot_verts = profiles[0],
         face_pts = len(tooth_profile),
         gear_pts = face_pts * teeth,
         top_faces =[
diff --git a/hull.scad b/hull.scad
index 5658bec..a7b73c7 100644
--- a/hull.scad
+++ b/hull.scad
@@ -158,7 +158,7 @@ function _hull_collinear(points) =
 //   %polyhedron(points=pts, faces=faces);
 function hull3d_faces(points) =
     assert(is_path(points,3),"Invalid input to hull3d_faces")
-    len(points) < 3 ? list_range(len(points))
+    len(points) < 3 ? range(len(points))
   : let ( // start with a single non-collinear triangle
           tri = noncollinear_triple(points, error=false)
         )
diff --git a/joiners.scad b/joiners.scad
index 8c0625d..2dd8175 100644
--- a/joiners.scad
+++ b/joiners.scad
@@ -968,7 +968,7 @@ module rabbit_clip(type, length, width,  snap, thickness, depth, compression=0.1
                       );
     assert(fullpath[4].y < fullpath[3].y, "Pin is too wide for its length");
     
-    snapmargin = -snap + select(sidepath,-1).x;// - compression;
+    snapmargin = -snap + last(sidepath).x;// - compression;
     if (is_pin){
       if (snapmargin<0) echo("WARNING: The snap is too large for the clip to squeeze to fit its socket")
       echo(snapmargin=snapmargin);
diff --git a/math.scad b/math.scad
index c63efb8..2297a29 100644
--- a/math.scad
+++ b/math.scad
@@ -495,7 +495,7 @@ function rand_int(minval, maxval, N, seed=undef) =
 function gaussian_rands(mean, stddev, N=1, seed=undef) =
     assert( is_finite(mean+stddev+N) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.")
     let(nums = is_undef(seed)? rands(0,1,N*2) : rands(0,1,N*2,seed))
-    [for (i = list_range(N)) mean + stddev*sqrt(-2*ln(nums[i*2]))*cos(360*nums[i*2+1])];
+    [for (i = range(N)) mean + stddev*sqrt(-2*ln(nums[i*2]))*cos(360*nums[i*2+1])];
 
 
 // Function: log_rands()
@@ -612,7 +612,7 @@ function _cumsum(v,_i=0,_acc=[]) =
         v, _i+1,
         concat(
             _acc,
-            [_i==0 ? v[_i] : select(_acc,-1)+v[_i]]
+            [_i==0 ? v[_i] : last(_acc) + v[_i]]
         )
     );
 
@@ -905,7 +905,7 @@ function _back_substitute(R, b, x=[]) =
     : let(
           newvalue = len(x)==0
             ? b[ind]/R[ind][ind]
-            : (b[ind]-select(R[ind],ind+1,-1) * x)/R[ind][ind]
+            : (b[ind]-list_tail(R[ind],ind+1) * x)/R[ind][ind]
       )
       _back_substitute(R, b, concat([newvalue],x));
 
@@ -1635,7 +1635,7 @@ function poly_mult(p,q) =
   is_undef(q) ?
     len(p)==2 
         ? poly_mult(p[0],p[1]) 
-    : poly_mult(p[0], poly_mult(select(p,1,-1)))
+    : poly_mult(p[0], poly_mult(list_tail(p)))
   :
   assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
     p*p==0 || q*q==0
@@ -1678,7 +1678,7 @@ function _poly_div(n,d,q) =
 ///    or give epsilon for approximate zeros.  
 function _poly_trim(p,eps=0) =
     let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i])
-    len(nz)==0 ? [0] : select(p,nz[0],-1);
+    len(nz)==0 ? [0] : list_tail(p,nz[0]);
 
 
 // Function: poly_add()
@@ -1719,7 +1719,7 @@ function poly_roots(p,tol=1e-14,error_bound=false) =
     let( p = _poly_trim(p,eps=0) )
     assert( p!=[0], "Input polynomial cannot be zero." )
     p[len(p)-1] == 0 ?                                       // Strip trailing zero coefficients
-        let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound))
+        let( solutions = poly_roots(list_head(p),tol=tol, error_bound=error_bound))
         (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
                     : [[0,0], each solutions]) :
     len(p)==1 ? (error_bound ? [[],[]] : []) :               // Nonzero constant case has no solutions
diff --git a/modular_hose.scad b/modular_hose.scad
index 365bab3..0a878a1 100644
--- a/modular_hose.scad
+++ b/modular_hose.scad
@@ -169,9 +169,9 @@ module modular_hose(size, type, clearance=0, waist_len, anchor=BOTTOM, spin=0,or
     type=="segment"? concat(back(midlength,p=smallend),yflip(p=bigend))
   : type=="small" || type=="ball" ?
           concat(back(midlength,p=smallend),
-                 [[select(smallend,-1).x,0],[ smallend[0].x,0]])
+                 [[last(smallend).x,0],[ smallend[0].x,0]])
   : concat( back(midlength,p=bigend), 
-                  [[select(bigend,-1).x,0],[ bigend[0].x,0]]);
+                  [[last(bigend).x,0],[ bigend[0].x,0]]);
   bounds = pointlist_bounds(shape);
   center = mean(bounds);
   attachable(anchor,spin,orient,l=bounds[1].y-bounds[0].y, r=bounds[1].x)
diff --git a/paths.scad b/paths.scad
index 28858fa..d5ebae7 100644
--- a/paths.scad
+++ b/paths.scad
@@ -183,7 +183,7 @@ function path_length(path,closed=false) =
 function path_segment_lengths(path, closed=false) =
     [
         for (i=[0:1:len(path)-2]) norm(path[i+1]-path[i]),
-        if (closed) norm(path[0]-select(path,-1))
+        if (closed) norm(path[0]-last(path))
     ]; 
 
 
@@ -482,13 +482,13 @@ function path_chamfer_and_rounding(path, closed, chamfer, rounding) =
 				p2 = select(path,i),
 				crn1 = select(corner_paths,i-1),
 				crn2 = corner_paths[i],
-				l1 = norm(select(crn1,-1)-p1),
+				l1 = norm(last(crn1)-p1),
 				l2 = norm(crn2[0]-p2),
 				needed = l1 + l2,
 				seglen = norm(p2-p1),
 				check = assert(seglen >= needed, str("Path segment ",i," is too short to fulfill rounding/chamfering for the adjacent corners."))
 			) each crn2,
-			if (!closed) select(path,-1)
+			if (!closed) last(path)
 		]
 	) deduplicate(out);
 
@@ -783,16 +783,16 @@ function _path_fast_defragment(fragments, eps=EPSILON, _done=[]) =
     len(fragments)==0? _done :
     let(
         path = fragments[0],
-        endpt = select(path,-1),
+        endpt = last(path),
         extenders = [
             for (i = [1:1:len(fragments)-1]) let(
                 test1 = approx(endpt,fragments[i][0],eps=eps),
-                test2 = approx(endpt,select(fragments[i],-1),eps=eps)
+                test2 = approx(endpt,last(fragments[i]),eps=eps)
             ) if (test1 || test2) (test1? i : -1)
         ]
     ) len(extenders) == 1 && extenders[0] >= 0? _path_fast_defragment(
         fragments=[
-            concat(select(path,0,-2),fragments[extenders[0]]),
+            concat(list_head(path),fragments[extenders[0]]),
             for (i = [1:1:len(fragments)-1])
                 if (i != extenders[0]) fragments[i]
         ],
@@ -814,7 +814,7 @@ function _extreme_angle_fragment(seg, fragments, rightmost=true, eps=EPSILON) =
             for (i = idx(fragments)) let(
                 fragment = fragments[i],
                 fwdmatch = approx(seg[1], fragment[0], eps=eps),
-                bakmatch =  approx(seg[1], select(fragment,-1), eps=eps)
+                bakmatch =  approx(seg[1], last(fragment), eps=eps)
             ) [
                 fwdmatch,
                 bakmatch,
@@ -872,12 +872,12 @@ function assemble_a_path_from_fragments(fragments, rightmost=true, startfrag=0,
         // Found fragment is already closed
         [foundfrag, concat([path], remainder)]
     ) : let(
-        fragend = select(foundfrag,-1),
+        fragend = last(foundfrag),
         hits = [for (i = idx(path,e=-2)) if(approx(path[i],fragend,eps=eps)) i]
     ) hits? (
         let(
             // Found fragment intersects with initial path
-            hitidx = select(hits,-1),
+            hitidx = last(hits),
             newpath = list_head(path,hitidx),
             newfrags = concat(len(newpath)>1? [newpath] : [], remainder),
             outpath = concat(slice(path,hitidx,-2), foundfrag)
@@ -1239,17 +1239,17 @@ module path_spread(path, n, spacing, sp=undef, rotate_children=true, closed=fals
 {
     length = path_length(path,closed);
     distances = is_def(sp)? (
-        is_def(n) && is_def(spacing)? list_range(s=sp, step=spacing, n=n) :
-        is_def(n)? list_range(s=sp, e=length, n=n) :
-        list_range(s=sp, step=spacing, e=length)
+        is_def(n) && is_def(spacing)? range(s=sp, step=spacing, n=n) :
+        is_def(n)? range(s=sp, e=length, n=n) :
+        range(s=sp, step=spacing, e=length)
     ) : is_def(n) && is_undef(spacing)? (
         closed?
-            let(range=list_range(s=0,e=length, n=n+1)) list_head(range) :
-            list_range(s=0, e=length, n=n)
+            let(range=range(s=0,e=length, n=n+1)) list_head(range) :
+            range(s=0, e=length, n=n)
     ) : (
         let(
             n = is_def(n)? n : floor(length/spacing)+(closed?0:1),
-            ptlist = list_range(s=0,step=spacing,n=n),
+            ptlist = range(s=0,step=spacing,n=n),
             listcenter = mean(ptlist)
         ) closed?
             sort([for(entry=ptlist) posmod(entry-listcenter,length)]) :
@@ -1333,7 +1333,7 @@ function path_cut_points(path, dists, closed=false, direction=false) =
 function _path_cut_points(path, dists, closed=false, pind=0, dtotal=0, dind=0, result=[]) =
     dind == len(dists) ? result :
     let(
-        lastpt = len(result)==0? [] : select(result,-1)[0],       // location of last cut point
+        lastpt = len(result)==0? [] : last(result)[0],       // location of last cut point
         dpartial = len(result)==0? 0 : norm(lastpt-select(path,pind)),  // remaining length in segment
         nextpoint = dists[dind] < dpartial+dtotal  // Do we have enough length left on the current segment?
            ? [lerp(lastpt,select(path,pind),(dists[dind]-dtotal)/dpartial),pind] 
@@ -1419,7 +1419,7 @@ function _path_cuts_dir(path, cuts, closed=false, eps=1e-2) =
 function path_cut(path,cutdist,closed) =
   is_num(cutdist) ? path_cut(path,[cutdist],closed) :
   assert(is_vector(cutdist))
-  assert(select(cutdist,-1)<path_length(path,closed=closed),"Cut distances must be smaller than the path length")
+  assert(last(cutdist)<path_length(path,closed=closed),"Cut distances must be smaller than the path length")
   assert(cutdist[0]>0, "Cut distances must be strictly positive")
   let(
       cutlist = path_cut_points(path,cutdist,closed=closed),
@@ -1433,7 +1433,7 @@ function path_cut(path,cutdist,closed) =
           cutlist[i][0]==cutlist[i+1][0] ? []
           :
           [ if (!approx(cutlist[i][0], select(path,cutlist[i][1]))) cutlist[i][0],
-            each slice(path, cutlist[i][1], cutlist[i+1][1]),
+            each slice(path, cutlist[i][1], cutlist[i+1][1]-1),
             if (!approx(cutlist[i+1][0], select(path,cutlist[i+1][1]-1))) cutlist[i+1][0],
           ],
       [
@@ -1553,7 +1553,7 @@ function subdivide_path(path, N, refine, closed=true, exact=true, method="length
                 lerp(path[i],select(path,i+1), j/(add[i]+1))
             ]
         ],
-        closed? [] : [select(path,-1)]
+        closed? [] : [last(path)]
     );
 
 
@@ -1575,7 +1575,7 @@ function path_length_fractions(path, closed=false) =
                 norm(select(path,i+1)-path[i])
         ],
         partial_len = cumsum(lengths),
-        total_len = select(partial_len,-1)
+        total_len = last(partial_len)
     ) partial_len / total_len;
 
 
@@ -1603,10 +1603,10 @@ function resample_path(path, N, spacing, closed=false) =
     length = path_length(path,closed),
     N = is_def(N) ? N : round(length/spacing) + (closed?0:1), 
         spacing = length/(closed?N:N-1),     // Note: worried about round-off error, so don't include 
-        distlist = list_range(closed?N:N-1,step=spacing),  // last point when closed=false
+        distlist = range(closed?N:N-1,step=spacing),  // last point when closed=false
     cuts = path_cut_points(path, distlist, closed=closed)
    )
-   concat(subindex(cuts,0),closed?[]:[select(path,-1)]);  // Then add last point here
+   concat(subindex(cuts,0),closed?[]:[last(path)]);  // Then add last point here
 
 
 
diff --git a/quaternions.scad b/quaternions.scad
index eaf9e16..1dc7192 100644
--- a/quaternions.scad
+++ b/quaternions.scad
@@ -223,7 +223,7 @@ function Q_Cumulative(v, _i=0, _acc=[]) =
         v, _i+1,
         concat(
             _acc,
-            [_i==0 ? v[_i] : Q_Mul(v[_i], select(_acc,-1))]
+            [_i==0 ? v[_i] : Q_Mul(v[_i], last(_acc))]
         )
     );
 
diff --git a/regions.scad b/regions.scad
index 8b5ca6b..b76785c 100644
--- a/regions.scad
+++ b/regions.scad
@@ -758,7 +758,7 @@ function offset(
         // Note if !closed the last corner doesn't matter, so exclude it
         parallelcheck =
             (len(sharpcorners)==2 && !closed) ||
-            all_defined(select(sharpcorners,closed?0:1,-1))
+            all_defined(closed? sharpcorners : list_tail(sharpcorners))
     )
     assert(parallelcheck, "Path contains sequential parallel segments (either 180 deg turn or 0 deg turn")
     let(
diff --git a/rounding.scad b/rounding.scad
index cf95e7b..315986a 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -194,7 +194,7 @@ include <structs.scad>
 //   square = [[0,0],[1,0],[1,1],[0,1]];
 //   spiral = flatten(repeat(concat(square,reverse(square)),5));  // Squares repeat 10 times, forward and backward
 //   squareind = [for(i=[0:9]) each [i,i,i,i]];                   // Index of the square for each point
-//   z = list_range(40)*.2+squareind;
+//   z = range(40)*.2+squareind;
 //   path3d = hstack(spiral,z);                                   // 3D spiral
 //   rounding = squareind/20;
 //       // Setting k=1 means curvature won't be continuous, but curves are as round as possible
@@ -423,13 +423,12 @@ function _rounding_offsets(edgespec,z_dir=1) =
                         ) :
                         edgetype == "circle"? radius==0? [] : [for(i=[1:N]) [radius*(cos(i*90/N)-1), z_dir*abs(radius)*sin(i*90/N)]] :
                         /* smooth */ joint==0 ? [] :
-                        select(
-                                _bezcorner([[0,0],[0,z_dir*abs(joint)],[-joint,z_dir*abs(joint)]], k, $fn=N+2),
-                                1, -1
+                        list_tail(
+                                _bezcorner([[0,0],[0,z_dir*abs(joint)],[-joint,z_dir*abs(joint)]], k, $fn=N+2)
                         )
         )
   
-        quant(extra > 0? concat(offsets, [select(offsets,-1)+[0,z_dir*extra]]) : offsets, 1/1024);
+        quant(extra > 0? concat(offsets, [last(offsets)+[0,z_dir*extra]]) : offsets, 1/1024);
 
 
 
@@ -670,7 +669,7 @@ function _path_join(paths,joint,k=0.5,i=0,result=[],relocate=true,closed=false)
       new_result = [each select(result,loop?nextcut[1]:0,len(revresult)-1-firstcut[1]),
                     each bezpath,
                     nextcut[0],
-                    if (!loop) each select(nextpath,nextcut[1],-1)
+                    if (!loop) each list_tail(nextpath,nextcut[1])
                    ]
   )
   i==len(paths)-(closed?1:2)
@@ -900,7 +899,7 @@ function _make_offset_polyhedron(path,offsets, offset_type, flip_faces, quality,
                                  vertexcount=0, vertices=[], faces=[] )=
         offsetind==len(offsets)? (
                 let(
-                        bottom = list_range(n=len(path),s=vertexcount),
+                        bottom = range(n=len(path),s=vertexcount),
                         oriented_bottom = !flip_faces? bottom : reverse(bottom)
                 ) [vertices, concat(faces,[oriented_bottom])]
         ) : (
@@ -980,8 +979,8 @@ function offset_sweep(
                 : 0,
 
         // "Extra" height enlarges the result beyond the requested height, so subtract it
-        bottom_height = len(offsets_bot)==0 ? 0 : abs(select(offsets_bot,-1)[1]) - struct_val(bottom,"extra"),
-        top_height = len(offsets_top)==0 ? 0 : abs(select(offsets_top,-1)[1]) - struct_val(top,"extra"),
+        bottom_height = len(offsets_bot)==0 ? 0 : abs(last(offsets_bot)[1]) - struct_val(bottom,"extra"),
+        top_height = len(offsets_top)==0 ? 0 : abs(last(offsets_top)[1]) - struct_val(top,"extra"),
 
         height = one_defined([l,h,height], "l,h,height", dflt=u_add(bottom_height,top_height)),
         middle = height-bottom_height-top_height
@@ -1129,7 +1128,7 @@ function os_mask(mask, out=false, extra,check_valid, quality, offset_maxstep, of
   let(
       points = ([for(pt=polygon_shift(mask,origin_index[0])) [xfactor*max(pt.x,0),-max(pt.y,0)]])
   )
-  os_profile(deduplicate(move(-points[1],p=select(points,1,-1))), extra,check_valid,quality,offset_maxstep,offset);
+  os_profile(deduplicate(move(-points[1],p=list_tail(points))), extra,check_valid,quality,offset_maxstep,offset);
 
 
 // Module: convex_offset_extrude()
@@ -1248,8 +1247,8 @@ module convex_offset_extrude(
         offsets_top = _rounding_offsets(top, 1);
 
         // "Extra" height enlarges the result beyond the requested height, so subtract it
-        bottom_height = len(offsets_bot)==0 ? 0 : abs(select(offsets_bot,-1)[1]) - struct_val(bottom,"extra");
-        top_height = len(offsets_top)==0 ? 0 : abs(select(offsets_top,-1)[1]) - struct_val(top,"extra");
+        bottom_height = len(offsets_bot)==0 ? 0 : abs(last(offsets_bot)[1]) - struct_val(bottom,"extra");
+        top_height = len(offsets_top)==0 ? 0 : abs(last(offsets_top)[1]) - struct_val(top,"extra");
 
         height = one_defined([l,h,height], "l,h,height", dflt=u_add(bottom_height,top_height));
         assert(height>=0, "Height must be nonnegative");
@@ -1583,11 +1582,11 @@ function _stroke_end(width,left, right, spec) =
                         cutright = cut[1],
                         // Create updated paths taking into account clipping for end rotation
                         newright = intright?
-                                concat([pathclip[0]],select(right,pathclip[1],-1)) :
-                                concat([pathextend],select(right,1,-1)),
+                                concat([pathclip[0]],list_tail(right,pathclip[1])) :
+                                concat([pathextend],list_tail(right)),
                         newleft = !intright?
-                                concat([pathclip[0]],select(left,pathclip[1],-1)) :
-                                concat([pathextend],select(left,1,-1)),
+                                concat([pathclip[0]],list_tail(left,pathclip[1])) :
+                                concat([pathextend],list_tail(left)),
                         // calculate corner angles, which are different when the cut is negative (outside corner)
                         leftangle = cutleft>=0?
                                 vector_angle([newleft[1],newleft[0],newright[0]])/2 :
@@ -1953,11 +1952,11 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
 function bezier_patch_degenerate(patch, splinesteps=16, reverse=false) =
     assert(is_num(splinesteps), "splinesteps must be a number")
     let(
-        top_degen = patch[0][0] == select(patch[0],-1),
-        bot_degen = select(patch,-1)[0] == select(select(patch,-1),-1),
-        left_degen = patch[0][0] == select(patch,-1)[0],
-        right_degen = select(patch[0],-1) == select(select(patch,-1),-1),
-        samplepts = list_range(splinesteps+1)/splinesteps
+        top_degen = patch[0][0] == last(patch[0]),
+        bot_degen = last(patch)[0] == last(last(patch)),
+        left_degen = patch[0][0] == last(patch)[0],
+        right_degen = last(patch[0]) == last(last(patch)),
+        samplepts = range(splinesteps+1)/splinesteps
     )
     top_degen && bot_degen && left_degen && right_degen ?   // fully degenerate case
         [repeat([patch[0][0]],4), EMPTY_VNF] :
@@ -1965,19 +1964,19 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false) =
         let(  
             pts = bezier_points(subindex(patch,0), samplepts)
         )
-        [[pts,pts,[pts[0]],[select(pts,-1)]], EMPTY_VNF] :
+        [[pts,pts,[pts[0]],[last(pts)]], EMPTY_VNF] :
     left_degen && right_degen ?                             // double degenerate case (sides)
        let(
           pts = bezier_points(patch[0], samplepts)
        )
-       [[[pts[0]], [select(pts,-1)], pts, pts], EMPTY_VNF] :
+       [[[pts[0]], [last(pts)], pts, pts], EMPTY_VNF] :
     !top_degen && !bot_degen ?                             // non-degenerate case
        let(
            k=echo("non-degenerate case"),
            pts = bezier_patch_points(patch, samplepts, samplepts)
        )
        [
-        [subindex(pts,0), subindex(pts,len(pts)-1), pts[0], select(pts,-1)],
+        [subindex(pts,0), subindex(pts,len(pts)-1), pts[0], last(pts)],
         vnf_vertex_array(pts, reverse=reverse)
        ] :
     bot_degen ?                                           // only bottom is degenerate
@@ -1990,22 +1989,22 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false) =
        ] :
     // at this point top_degen is true                   // only top is degenerate
        let(
-           full_degen = patch[1][0] == select(patch[1],-1),
-           rowmax = full_degen ? list_range(splinesteps+1) :
+           full_degen = patch[1][0] == last(patch[1]),
+           rowmax = full_degen ? range(splinesteps+1) :
                                  [for(j=[0:splinesteps]) j<=splinesteps/2 ? 2*j : splinesteps],
            vbb=echo("single degenerate case"),
            bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(subindex(patch,i), samplepts)],
            pts = [
                   [bpatch[0][0]],
-                  for(j=[1:splinesteps]) bezier_points(subindex(bpatch,j), list_range(rowmax[j]+1)/rowmax[j])
+                  for(j=[1:splinesteps]) bezier_points(subindex(bpatch,j), range(rowmax[j]+1)/rowmax[j])
                  ],
            vnf = vnf_tri_array(pts, reverse=reverse)
         ) [
             [
              subindex(pts,0),
-             [for(row=pts) select(row,-1)],
+             [for(row=pts) last(row)],
              pts[0],
-             select(pts,-1),
+             last(pts),
             ],
             vnf
           ];
diff --git a/screws.scad b/screws.scad
index 27c6c5c..1ce6ec0 100644
--- a/screws.scad
+++ b/screws.scad
@@ -1085,7 +1085,7 @@ function _ISO_thread_tolerance(diameter, pitch, internal=false, tolerance=undef)
            ],
 
     rangepts = [0.99, 1.4, 2.8, 5.6, 11.2, 22.4, 45, 90, 180, 300],
-    d_ind = floor(lookup(diameter,hstack(rangepts,list_range(len(rangepts))))),
+    d_ind = floor(lookup(diameter,hstack(rangepts,range(len(rangepts))))),
     avgd = sqrt(rangepts[d_ind]* rangepts[d_ind+1]),
 
     T_d2_6 = 90*pow(P, 0.4)*pow(avgd,0.1),
diff --git a/shapes.scad b/shapes.scad
index 1a2ed21..c0eba7c 100644
--- a/shapes.scad
+++ b/shapes.scad
@@ -1339,7 +1339,7 @@ module spheroid(r, style="aligned", d, circum=false, anchor=CENTER, spin=0, orie
             path = [
                 let(a = merids[0]) [0, sin(a)],
                 for (a=merids) [cos(a), sin(a)],
-                let(a = select(merids,-1)) [0, sin(a)]
+                let(a = last(merids)) [0, sin(a)]
             ];
             scale(r) rotate(180) rotate_extrude(convexity=2,$fn=sides) polygon(path);
         } else {
@@ -1441,7 +1441,7 @@ function spheroid(r, style="aligned", d, circum=false, anchor=CENTER, spin=0, or
                 1,
             ],
             offs = cumsum(meridians),
-            pc = select(offs,-1)-1,
+            pc = last(offs)-1,
             os = octa_steps * 2
         ) [
             for (i=[0:1:3]) [0, 1+(i+1)%4, 1+i],
diff --git a/shapes2d.scad b/shapes2d.scad
index 34f92ae..5ba5c2a 100644
--- a/shapes2d.scad
+++ b/shapes2d.scad
@@ -209,10 +209,10 @@ module stroke(
     assert(is_undef(endcap_angle2)||is_num(endcap_angle2));
     assert(is_undef(joint_angle)||is_num(joint_angle));
 
-    endcap_shape1 = _shape_path(endcap1, select(width,0), endcap_width1, endcap_length1, endcap_extent1);
-    endcap_shape2 = _shape_path(endcap2, select(width,-1), endcap_width2, endcap_length2, endcap_extent2);
+    endcap_shape1 = _shape_path(endcap1, width[0], endcap_width1, endcap_length1, endcap_extent1);
+    endcap_shape2 = _shape_path(endcap2, last(width), endcap_width2, endcap_length2, endcap_extent2);
 
-    trim1 = select(width,0) * first_defined([
+    trim1 = width[0] * first_defined([
         trim1, trim,
         (endcap1=="arrow")? endcap_length1-0.01 :
         (endcap1=="arrow2")? endcap_length1*3/4 :
@@ -220,7 +220,7 @@ module stroke(
     ]);
     assert(is_num(trim1));
 
-    trim2 = select(width,-1) * first_defined([
+    trim2 = last(width) * first_defined([
         trim2, trim,
         (endcap2=="arrow")? endcap_length2-0.01 :
         (endcap2=="arrow2")? endcap_length2*3/4 :
@@ -244,8 +244,8 @@ module stroke(
             [lerp(width[epos[0]], width[(epos[0]+1)%len(width)], epos[1])]
         );
 
-        start_vec = select(path,0) - select(path,1);
-        end_vec = select(path,-1) - select(path,-2);
+        start_vec = path[0] - path[1];
+        end_vec = last(path) - select(path,-2);
 
         if (len(path[0]) == 2) {
             // Straight segments
@@ -301,7 +301,7 @@ module stroke(
             }
 
             // Endcap2
-            translate(select(path,-1)) {
+            translate(last(path)) {
                 mat = is_undef(endcap_angle2)? rot(from=BACK,to=end_vec) :
                     zrot(endcap_angle2);
                 multmatrix(mat) polygon(endcap_shape2);
@@ -402,9 +402,9 @@ module stroke(
             }
 
             // Endcap2
-            translate(select(path,-1)) {
-                multmatrix(select(rotmats,-1)) {
-                    $fn = select(sides,-1);
+            translate(last(path)) {
+                multmatrix(last(rotmats)) {
+                    $fn = last(sides);
                     if (is_undef(endcap_angle2)) {
                         rotate_extrude(convexity=convexity) {
                             right_half(planar=true) {
@@ -413,7 +413,7 @@ module stroke(
                         }
                     } else {
                         rotate([90,0,endcap_angle2]) {
-                            linear_extrude(height=max(select(widths,-1),0.001), center=true, convexity=convexity) {
+                            linear_extrude(height=max(last(widths),0.001), center=true, convexity=convexity) {
                                 polygon(endcap_shape2);
                             }
                         }
@@ -791,7 +791,7 @@ function _turtle_command(command, parm, parm2, state, index) =
         chvec = !in_list(command,needvec) || is_vector(parm,2),
         chnum = !in_list(command,neednum) || is_num(parm),
         vec_or_num = !in_list(command,needeither) || (is_num(parm) || is_vector(parm,2)),
-        lastpt = select(state[path],-1)
+        lastpt = last(state[path])
     )
     assert(chvec,str("\"",command,"\" requires a vector parameter at index ",index))
     assert(chnum,str("\"",command,"\" requires a numeric parameter at index ",index))
@@ -2339,8 +2339,8 @@ function mask2d_ogee(pattern, excess, anchor=CENTER, spin=0) =
                 0
             )
         ])),
-        tot_x = select(x,-1),
-        tot_y = select(y,-1),
+        tot_x = last(x),
+        tot_y = last(y),
         data = [
             for (i=idx(pattern,step=2)) let(
                 type = pattern[i],
diff --git a/skin.scad b/skin.scad
index e53def9..8275c51 100644
--- a/skin.scad
+++ b/skin.scad
@@ -444,8 +444,8 @@ function skin(profiles, slices, refine=1, method="direct", sampling, caps, close
     // Define this to be 1 if a profile is used on either side by a resampling method, zero otherwise.
     profile_resampled = [for(i=idx(profiles)) 
       1-(
-           i==0 ?  method_type[0] * (closed? select(method_type,-1) : 1) :
-           i==len(profiles)-1 ? select(method_type,-1) * (closed ? select(method_type,-2) : 1) :
+           i==0 ?  method_type[0] * (closed? last(method_type) : 1) :
+           i==len(profiles)-1 ? last(method_type) * (closed ? select(method_type,-2) : 1) :
          method_type[i] * method_type[i-1])],
     parts = search(1,[1,for(i=[0:1:len(profile_resampled)-2]) profile_resampled[i]!=profile_resampled[i+1] ? 1 : 0],0),
     plen = [for(i=idx(parts)) (i== len(parts)-1? len(refined_len) : parts[i+1]) - parts[i]],
@@ -534,8 +534,8 @@ function _skin_core(profiles, caps) =
                         prof1 = profiles[0]
                 ) [[for (i=idx(prof1)) plens[0]-1-i]],
                 secondcap = !caps[1] ? [] : let(
-            prof2 = select(profiles,-1),
-            eoff = sum(select(plens,0,-2))
+            prof2 = last(profiles),
+            eoff = sum(list_head(plens))
         ) [[for (i=idx(prof2)) eoff+i]]
     ) [vertices, concat(sidefaces,firstcap,secondcap)];
 
@@ -812,7 +812,7 @@ function _skin_tangent_match(poly1, poly2) =
         small = swap ? poly2 : poly1,
         curve_offset = centroid(small)-centroid(big),
         cutpts = [for(i=[0:len(small)-1]) _find_one_tangent(big, select(small,i,i+1),curve_offset=curve_offset)],
-        shift = select(cutpts,-1)+1,
+        shift = last(cutpts)+1,
         newbig = polygon_shift(big, shift),
         repeat_counts = [for(i=[0:len(small)-1]) posmod(cutpts[i]-select(cutpts,i-1),len(big))],
         newsmall = repeat_entries(small,repeat_counts)
@@ -894,7 +894,7 @@ function associate_vertices(polygons, split, curpoly=0) =
            str("Split ",cursplit," at polygon ",curpoly," has invalid vertices.  Must be in [0:",polylen-1,"]"))
     len(cursplit)==0 ? associate_vertices(polygons,split,curpoly+1) :
     let(
-      splitindex = sort(concat(list_range(polylen), cursplit)),
+      splitindex = sort(concat(range(polylen), cursplit)),
       newpoly = [for(i=[0:len(polygons)-1]) i<=curpoly ? select(polygons[i],splitindex) : polygons[i]]
     )
    associate_vertices(newpoly, split, curpoly+1);
@@ -970,10 +970,10 @@ function sweep(shape, transforms, closed=false, caps) =
         rtrans = reverse(transforms),
         vnfs = [
             for (rgn=regions) each [
-                for (path=select(rgn,0,-1))
+                for (path=rgn)
                     sweep(path, transforms, closed=closed, caps=false),
                 if (fullcaps[0]) region_faces(rgn, transform=transforms[0], reverse=true),
-                if (fullcaps[1]) region_faces(rgn, transform=select(transforms,-1)),
+                if (fullcaps[1]) region_faces(rgn, transform=last(transforms)),
             ],
         ],
         vnf = vnf_merge(vnfs)
@@ -1311,7 +1311,7 @@ function path_sweep(shape, path, method="incremental", normal, closed=false, twi
   assert(is_integer(symmetry) && symmetry>0, "symmetry must be a positive integer")
 //  let(shape = check_and_fix_path(shape,valid_dim=2,closed=true,name="shape"))
   assert(is_path(path), "input path is not a path")
-  assert(!closed || !approx(path[0],select(path,-1)), "Closed path includes start point at the end")
+  assert(!closed || !approx(path[0],last(path)), "Closed path includes start point at the end")
   let(
     path = path3d(path),
     caps = is_def(caps) ? caps :
@@ -1350,13 +1350,13 @@ function path_sweep(shape, path, method="incremental", normal, closed=false, twi
             // X axis.  Similarly, in the closed==false case the desired and actual transformations can only differ in the twist,
             // so we can need to calculate the twist angle so we can apply a correction, which we distribute uniformly over the whole path.
             reference_rot = closed ? rotations[0] :   
-                         is_undef(last_normal) ? select(rotations,-1) : 
+                         is_undef(last_normal) ? last(rotations) : 
                            let(
-                               last_tangent = select(tangents,-1),
+                               last_tangent = last(tangents),
                                lastynormal = last_normal - (last_normal * last_tangent) * last_tangent
                            )
                          affine3d_frame_map(y=lastynormal, z=last_tangent),
-            mismatch = transpose(select(rotations,-1)) * reference_rot, 
+            mismatch = transpose(last(rotations)) * reference_rot, 
             correction_twist = atan2(mismatch[1][0], mismatch[0][0]),
             // Spread out this extra twist over the whole sweep so that it doesn't occur
             // abruptly as an artifact at the last step.  
@@ -1492,7 +1492,7 @@ function _ofs_vmap(ofs,closed=false) =
     )
     [
      for(entry=ofs[1]) _ofs_face_edge(entry,firstlen),
-     if (!closed) _ofs_face_edge(select(ofs[1],-1),firstlen,second=true)
+     if (!closed) _ofs_face_edge(last(ofs[1]),firstlen,second=true)
     ];
 
 
diff --git a/strings.scad b/strings.scad
index be287de..147d558 100644
--- a/strings.scad
+++ b/strings.scad
@@ -335,7 +335,7 @@ function _str_find_last(str,pattern,sindex) =
         (sindex >=0 ? sindex : undef);
 
 function _str_find_all(str,pattern) =
-    pattern == "" ? list_range(len(str)) :
+    pattern == "" ? range(len(str)) :
     [for(i=[0:1:len(str)-len(pattern)]) if (_str_cmp(str,i,pattern)) i];
 
 
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 5b460e2..679f34c 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -29,12 +29,14 @@ test_select();
 
 
 module test_slice() {
-    assert(slice([3,4,5,6,7,8,9], 3, 5) == [6,7]);
-    assert(slice([3,4,5,6,7,8,9], 2, -1) == [5,6,7,8,9]);
-    assert(slice([3,4,5,6,7,8,9], 1, 1) == []);
-    assert(slice([3,4,5,6,7,8,9], 6, -1) == [9]);
-    assert(slice([3,4,5,6,7,8,9], 2, -2) == [5,6,7,8]);
-    assert(slice([], 2, -2) == []);
+    l = [3,4,5,6,7,8,9];
+    assert(slice(l, 5, 6) == [8,9]);
+    assert(slice(l, 5, 8) == [8,9]);
+    assert(slice(l, 5, 2) == []);
+    assert(slice(l, -3, -1) == [7,8,9]);
+    assert(slice(l, 3, 3) == [6]);
+    assert(slice(l, 4) == [7,8,9]);
+    assert(slice(l, -2) == [8,9]);
 }
 test_slice();
 
@@ -131,19 +133,36 @@ module test_repeat() {
 test_repeat();
 
 
-module test_list_range() {
-    assert(list_range(4) == [0,1,2,3]);
-    assert(list_range(n=4, step=2) == [0,2,4,6]);
-    assert(list_range(n=4, s=3, step=3) == [3,6,9,12]);
-    assert(list_range(e=3) == [0,1,2,3]);
-    assert(list_range(e=6, step=2) == [0,2,4,6]);
-    assert(list_range(s=3, e=5) == [3,4,5]);
-    assert(list_range(s=3, e=8, step=2) == [3,5,7]);
-    assert(list_range(s=4, e=8, step=2) == [4,6,8]);
-    assert(list_range(e=4, n=3) == [0,2,4]);
-    assert(list_range(n=4, s=[3,4], step=[2,3]) == [[3,4], [5,7], [7,10], [9,13]]);
+module test_range() {
+    assert(range(4) == [0,1,2,3]);
+    assert(range(n=4, step=2) == [0,2,4,6]);
+    assert(range(n=4, s=3, step=3) == [3,6,9,12]);
+    assert(range(e=3) == [0,1,2,3]);
+    assert(range(e=6, step=2) == [0,2,4,6]);
+    assert(range(s=3, e=5) == [3,4,5]);
+    assert(range(s=3, e=8, step=2) == [3,5,7]);
+    assert(range(s=4, e=8, step=2) == [4,6,8]);
+    assert(range(e=4, n=3) == [0,2,4]);
+    assert(range(n=4, s=[3,4], step=[2,3]) == [[3,4], [5,7], [7,10], [9,13]]);
 }
-test_list_range();
+test_range();
+
+
+module test_rangex() {
+    assert_equal(rangex(4), [0,1,2,3]);
+    assert_approx(rangex(5,e=1), [0, 0.2, 0.4, 0.6, 0.8]);
+    assert_equal(rangex(n=4, step=2), [0,2,4,6]);
+    assert_approx(rangex(n=4, step=0.25), [0,0.25,0.5,0.75]);
+    assert_equal(rangex(n=4, s=3, step=3), [3,6,9,12]);
+    assert_equal(rangex(e=3), [0,1,2]);
+    assert_equal(rangex(e=6, step=2), [0,2,4]);
+    assert_equal(rangex(s=3, e=5), [3,4]);
+    assert_equal(rangex(s=3, e=8, step=2), [3,5,7]);
+    assert_equal(rangex(s=4, e=8, step=2), [4,6]);
+    assert_equal(rangex(e=6, n=3), [0,2,4]);
+    assert_equal(rangex(n=4, s=[3,4], step=[2,3]), [[3,4], [5,7], [7,10], [9,13]]);
+}
+test_rangex();
 
 
 module test_reverse() {
@@ -303,7 +322,7 @@ test_enumerate();
 
 
 module test_shuffle() {
-    nums1 = [for (i=list_range(100)) i];
+    nums1 = [for (i=range(100)) i];
     nums2 = shuffle(nums1,33);
     nums3 = shuffle(nums2,99);
     assert(sort(nums2)==nums1);
diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad
index 817323d..5eb4b68 100644
--- a/tests/test_geometry.scad
+++ b/tests/test_geometry.scad
@@ -481,7 +481,7 @@ module test_circle_3points() {
     radii = rands(10,100,count,seed_value=390);
     angles = rands(0,360,count,seed_value=699);
     // 2D tests.
-    for(i = list_range(count)) {
+    for(i = range(count)) {
         cp = select(coords,i,i+1);
         r = radii[i];
         angs = sort(select(angles,i,i+2));
@@ -506,7 +506,7 @@ module test_circle_3points() {
             assert(approx(res[2], UP));
         }
     }
-    for(i = list_range(count)) {
+    for(i = range(count)) {
         cp = select(coords,i,i+1);
         r = radii[i];
         angs = sort(select(angles,i,i+2));
@@ -532,7 +532,7 @@ module test_circle_3points() {
         }
     }
     // 3D tests.
-    for(i = list_range(count)) {
+    for(i = range(count)) {
         cp = select(coords,i,i+2);
         r = radii[i];
         nrm = unit(select(coords,i+10,i+12));
@@ -559,7 +559,7 @@ module test_circle_3points() {
             assert(approx(res[2], n));
         }
     }
-    for(i = list_range(count)) {
+    for(i = range(count)) {
         cp = select(coords,i,i+2);
         r = radii[i];
         nrm = unit(select(coords,i+10,i+12));
@@ -984,8 +984,8 @@ module test_pointlist_bounds() {
 
 
 module test_closest_point() {
-    ptlist = [for (i=list_range(100)) rands(-100,100,2,seed_value=8463)];
-    testpts = [for (i=list_range(100)) rands(-100,100,2,seed_value=6834)];
+    ptlist = [for (i=range(100)) rands(-100,100,2,seed_value=8463)];
+    testpts = [for (i=range(100)) rands(-100,100,2,seed_value=6834)];
     for (pt = testpts) {
         pidx = closest_point(pt,ptlist);
         dists = [for (p=ptlist) norm(pt-p)];
@@ -997,8 +997,8 @@ module test_closest_point() {
 
 
 module test_furthest_point() {
-    ptlist = [for (i=list_range(100)) rands(-100,100,2,seed_value=8463)];
-    testpts = [for (i=list_range(100)) rands(-100,100,2,seed_value=6834)];
+    ptlist = [for (i=range(100)) rands(-100,100,2,seed_value=8463)];
+    testpts = [for (i=range(100)) rands(-100,100,2,seed_value=6834)];
     for (pt = testpts) {
         pidx = furthest_point(pt,ptlist);
         dists = [for (p=ptlist) norm(pt-p)];
diff --git a/tests/test_math.scad b/tests/test_math.scad
index 103be26..e943756 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -1122,7 +1122,7 @@ module test_real_roots(){
    // Wilkinson polynomial is a nasty test:
    assert_approx(
        sort(real_roots(poly_mult([[1,-1],[1,-2],[1,-3],[1,-4],[1,-5],[1,-6],[1,-7],[1,-8],[1,-9],[1,-10]]))),
-       list_range(n=10,s=1));
+       range(n=10,s=1));
    assert_equal(real_roots([3]), []);
    assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50]])),[]);
    assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50],[1,0,0]])),[0,0]);
diff --git a/threading.scad b/threading.scad
index ff5ee96..50e1336 100644
--- a/threading.scad
+++ b/threading.scad
@@ -228,13 +228,13 @@ module trapezoidal_threaded_rod(
             zrot(i*360/starts, p=vnf_vertex_array(thread_verts, reverse=left_handed)),
         for (i=[0:1:starts-1]) let(
             rmat = zrot(i*360/starts),
-            pts = deduplicate(select(thread_verts[0], 0, len(prof3d)+1)),
+            pts = deduplicate(list_head(thread_verts[0], len(prof3d)+1)),
             faces = [for (i=idx(pts,e=-2)) [0, i+1, i]],
             rfaces = left_handed? [for (x=faces) reverse(x)] : faces
         ) [apply(rmat,pts), rfaces],
         for (i=[0:1:starts-1]) let(
             rmat = zrot(i*360/starts),
-            pts = deduplicate(select(last(thread_verts), -len(prof3d)-2, -1)),
+            pts = deduplicate(list_tail(last(thread_verts), -len(prof3d)-2)),
             faces = [for (i=idx(pts,e=-2)) [len(pts)-1, i, i+1]],
             rfaces = left_handed? [for (x=faces) reverse(x)] : faces
         ) [apply(rmat,pts), rfaces]
diff --git a/turtle3d.scad b/turtle3d.scad
index c0631d0..7f1e4c5 100644
--- a/turtle3d.scad
+++ b/turtle3d.scad
@@ -479,7 +479,7 @@ function _tupdate(state, tran, pretran) =
     [
      concat(state[0],tran),
      concat(state[1],pretran),
-     each select(state,2,-1)
+     each list_tail(state,2)
     ];
 
 function _turtle3d_command(command, parm, parm2, state, index) =
@@ -504,8 +504,8 @@ function _turtle3d_command(command, parm, parm2, state, index) =
         chnum = (!in_list(command,neednum) || is_num(parm))
                 && (!in_list(command,numornothing) || (is_undef(parm) || is_num(parm))),
         chtran = !in_list(command,needtran) || is_matrix(parm,4,4),
-        lastT = select(state[trlist],-1),
-        lastPre = select(state[prelist],-1),
+        lastT = last(state[trlist]),
+        lastPre = last(state[prelist]),
         lastpt = apply(lastT,[0,0,0])
     )
     assert(chvec,str("\"",command,"\" requires a 3d vector parameter at index ",index))
diff --git a/vnf.scad b/vnf.scad
index 707f1c4..7635cf5 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -316,7 +316,7 @@ function vnf_vertex_array(
             )
         ]
     )
-    rows<=1 || cols<=1 ? vnf : 
+    rows<=1 || cols<=1 ? vnf :
     vnf_merge(cleanup=true, [
         vnf, [
             verts,
@@ -396,7 +396,7 @@ module vnf_polyhedron(vnf, convexity=2, extent=true, cp=[0,0,0], anchor="origin"
 // Usage:
 //   vnf_wireframe(vnf, <r|d>);
 // Description:
-//   Given a VNF, creates a wire frame ball-and-stick model of the polyhedron with a cylinder for each edge and a sphere at each vertex. 
+//   Given a VNF, creates a wire frame ball-and-stick model of the polyhedron with a cylinder for each edge and a sphere at each vertex.
 // Arguments:
 //   vnf = A vnf structure
 //   r|d = radius or diameter of the cylinders forming the wire frame.  Default: r=1
@@ -404,12 +404,12 @@ module vnf_polyhedron(vnf, convexity=2, extent=true, cp=[0,0,0], anchor="origin"
 //   $fn=32;
 //   ball = sphere(r=20, $fn=6);
 //   vnf_wireframe(ball,d=1);
-// Example: 
+// Example:
 //  include<BOSL2/polyhedra.scad>
 //  $fn=32;
 //  cube_oct = regular_polyhedron_info("vnf", name="cuboctahedron", or=20);
 //  vnf_wireframe(cube_oct);
-// Example: The spheres at the vertex are imperfect at aligning with the cylinders, so especially at low $fn things look prety ugly.  This is normal.  
+// Example: The spheres at the vertex are imperfect at aligning with the cylinders, so especially at low $fn things look prety ugly.  This is normal.
 //  include<BOSL2/polyhedra.scad>
 //  $fn=8;
 //  octahedron = regular_polyhedron_info("vnf", name="octahedron", or=20);
@@ -423,7 +423,7 @@ module vnf_wireframe(vnf, r, d)
                  ]);
   for (e=edges) extrude_from_to(vertex[e[0]],vertex[e[1]]) circle(r=r);
   move_copies(vertex) sphere(r=r);
-}  
+}
 
 
 // Function: vnf_volume()
@@ -451,7 +451,7 @@ function vnf_volume(vnf) =
 //   no holes; otherwise the results are undefined.
 
 // Divide the solid up into tetrahedra with the origin as one vertex.  The centroid of a tetrahedron is the average of its vertices.
-// The centroid of the total is the volume weighted average.  
+// The centroid of the total is the volume weighted average.
 function vnf_centroid(vnf) =
     let(
         verts = vnf[0],
@@ -481,7 +481,7 @@ function _triangulate_planar_convex_polygons(polys) =
         tris = [for (poly=polys) if (len(poly)==3) poly],
         bigs = [for (poly=polys) if (len(poly)>3) poly],
         newtris = [for (poly=bigs) select(poly,-2,0)],
-        newbigs = [for (poly=bigs) select(poly,0,-2)],
+        newbigs = [for (poly=bigs) list_head(poly)],
         newtris2 = _triangulate_planar_convex_polygons(newbigs),
         outtris = concat(tris, newtris, newtris2)
     ) outtris;
@@ -626,7 +626,7 @@ function vnf_bend(vnf,r,d,axis="Z") =
 //   bad edges and vertices, overlaid on a transparent gray polyhedron of the VNF.
 //   .
 //   Currently checks for these problems:
-//   Type    | Color    | Code         | Message 
+//   Type    | Color    | Code         | Message
 //   ------- | -------- | ------------ | ---------------------------------
 //   WARNING | Yellow   | BIG_FACE     | Face has more than 3 vertices, and may confuse CGAL.
 //   WARNING | Brown    | NULL_FACE    | Face has zero area.
@@ -695,7 +695,7 @@ function vnf_bend(vnf,r,d,axis="Z") =
 //       move([75,35,30],p=vnf_triangulate(linear_sweep(square(100,center=true), height=100)))
 //   ]);
 //   vnf_validate(vnf,size=2,check_isects=true);
-// Example: HOLE_EDGE Errors; Edges Adjacent to Holes.  
+// Example: HOLE_EDGE Errors; Edges Adjacent to Holes.
 //   vnf = skin([
 //       path3d(regular_ngon(n=4, d=100),0),
 //       path3d(regular_ngon(n=5, d=100),100)