diff --git a/exprtk.hpp b/exprtk.hpp
index 6aa8536..1e33de8 100644
--- a/exprtk.hpp
+++ b/exprtk.hpp
@@ -21334,6 +21334,7 @@
          if (result)
          {
             brkcnt_list_.push_front(false);
+
             if (0 == (loop_body = parse_multi_sequence("for-loop")))
             {
                set_error(
@@ -36257,6 +36258,45 @@
    };
 
    template <typename T>
+   class nthelement : public exprtk::igeneric_function<T>
+   {
+   public:
+
+      typedef typename exprtk::igeneric_function<T> igfun_t;
+      typedef typename igfun_t::parameter_list_t    parameter_list_t;
+      typedef typename igfun_t::generic_type        generic_type;
+      typedef typename generic_type::scalar_view    scalar_t;
+      typedef typename generic_type::vector_view    vector_t;
+
+      using exprtk::igeneric_function<T>::operator();
+
+      nthelement()
+      : exprtk::igeneric_function<T>("VT|VTTT")
+      {}
+
+      inline T operator()(const std::size_t& ps_index, parameter_list_t parameters)
+      {
+         vector_t vec(parameters[0]);
+
+         std::size_t n  = 0;
+         std::size_t r0 = 0;
+         std::size_t r1 = vec.size() - 1;
+
+         if (!scalar_t(parameters[1]).to_uint(n))
+            return T(0);
+
+         if ((1 == ps_index) && !details::load_vector_range<T>::process(parameters,r0,r1,2,3))
+            return std::numeric_limits<T>::quiet_NaN();
+         else if (details::invalid_range(vec,r0,r1))
+            return std::numeric_limits<T>::quiet_NaN();
+
+         std::nth_element(vec.begin() + r0, vec.begin() + r0 + n , vec.begin() + r1 + 1);
+
+         return T(1);
+      }
+   };
+
+   template <typename T>
    class sumk : public exprtk::igeneric_function<T>
    {
    public:
@@ -36666,6 +36706,7 @@
       shift_left <T> sl;
       shift_right<T> sr;
       sort       <T> st;
+      nthelement <T> ne;
       sumk       <T> sk;
       axpy       <T> b1_axpy;
       axpby      <T> b1_axpby;
@@ -36703,6 +36744,8 @@
             return false;
          else if (!symtab.add_function("sort"         ,st))
             return false;
+         else if (!symtab.add_function("nth_element"  ,ne))
+            return false;
          else if (!symtab.add_function("sumk"         ,sk))
             return false;
          else if (!symtab.add_function("axpy"    ,b1_axpy))
diff --git a/exprtk_simple_example_14.cpp b/exprtk_simple_example_14.cpp
index 3be1065..b7e71f1 100644
--- a/exprtk_simple_example_14.cpp
+++ b/exprtk_simple_example_14.cpp
@@ -41,7 +41,6 @@
    expression_t expression;
 
    parser_t parser;
-
    parser.compile(stddev_program,expression);
 
    T stddev = expression.value();
diff --git a/exprtk_test.cpp b/exprtk_test.cpp
index 733a0bf..5a9f423 100644
--- a/exprtk_test.cpp
+++ b/exprtk_test.cpp
@@ -6532,6 +6532,7 @@
                      "var v[5] := {1,3,5,2,4}; var r[5] := {1,2,3,4,5}; sort(v); sum(v == r) == v[]",
                      "var v[5] := {5,4,3,2,1}; var r[5] := {1,2,3,4,5}; sort(v); sum(v == r) == v[]",
                      "var v[5] := {1,4,2,3,5}; var r[5] := {1,2,3,4,5}; sort(v,1,3); sum(v == r) == v[]",
+                     "var v[5] := {5,4,2,3,1}; var r[5] := {5,2,3,4,1}; sort(v,1,3); sum(v == r) == v[]",
                      "var v[5] := {3,1,2,4,5}; var r[5] := {1,2,3,4,5}; sort(v,0,2); sum(v == r) == v[]",
                      "var v[5] := {1,2,5,3,4}; var r[5] := {1,2,3,4,5}; sort(v,2,4); sum(v == r) == v[]",
 
@@ -6549,12 +6550,18 @@
                      "var v[5] := {3,1,2,4,5}; var r[5] := {3,2,1,4,5}; sort(v,'descending',0,2); sum(v == r) == v[]",
                      "var v[5] := {1,2,5,3,4}; var r[5] := {1,2,5,4,3}; sort(v,'descending',2,4); sum(v == r) == v[]",
 
-                     " var a := 2; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var r[3] := [0]; r := a*x+y; axpy(a,x,y); sum(y == r) == y[]",
-                     " var a := 2; var b := 3; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var r[3] := [0]; r := a*x+b*y; axpby(a,x,b,y); sum(y == r) == y[]",
+                     "var v[9] := {7,8,9,1,2,3,4,5,6}; nth_element(v,trunc(v[] / 2)); v[v[] / 2] == 5",
+                     "var v[9] := {7,8,9,1,2,3,4,5,6}; nth_element(v,trunc(v[] / 3)); v[v[] / 3] == 4",
 
-                     " var a := 2; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var z[3] := [0]; var r[3] := [0]; r := a*x+y; axpyz(a,x,y,z); sum(z == r) == z[]",
-                     " var a := 2; var b := 3; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var z[3] := [0]; var r[3] := [0]; r := a*x+b*y; axpbyz(a,x,b,y,z); sum(z == r) == z[]",
-                     " var a := 2; var b := 3; var x[3] := {1,2,3}; var z[3] := [0]; var r[3] := [0]; r := a*x+b; axpbz(a,x,b,z); sum(z == r) == z[]",
+                     "var v[9] := {7,8,9,1,2,3,4,5,6}; nth_element(v,trunc(v[] / 2)); sort(v,0,trunc(v[] / 2)); (v[v[] / 2] == 5) and (v[0] == 1)",
+                     "var v[9] := {7,8,9,1,2,3,4,5,6}; nth_element(v,trunc(v[] / 3)); sort(v,0,trunc(v[] / 3)); (v[v[] / 3] == 4) and (v[0] == 1)",
+
+                     " var a := 2; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var r[3] := [0]; r := a * x + y; axpy(a,x,y); sum(y == r) == y[]",
+                     " var a := 2; var b := 3; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var r[3] := [0]; r := a * x + b * y; axpby(a,x,b,y); sum(y == r) == y[]",
+
+                     " var a := 2; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var z[3] := [0]; var r[3] := [0]; r := a * x + y; axpyz(a,x,y,z); sum(z == r) == z[]",
+                     " var a := 2; var b := 3; var x[3] := {1,2,3}; var y[3] := {1,2,3}; var z[3] := [0]; var r[3] := [0]; r := a * x + b * y; axpbyz(a,x,b,y,z); sum(z == r) == z[]",
+                     " var a := 2; var b := 3; var x[3] := {1,2,3}; var z[3] := [0]; var r[3] := [0]; r := a * x + b; axpbz(a,x,b,z); sum(z == r) == z[]",
                   };
 
       const std::size_t expr_str_list_size = sizeof(expr_str_list) / sizeof(std::string);
diff --git a/readme.txt b/readme.txt
index 4831116..56d1aef 100644
--- a/readme.txt
+++ b/readme.txt
@@ -1301,7 +1301,7 @@
        avg, max, min, mul,  dot, dotk, sum, sumk,  count, all_true,
        all_false, any_true, any_false
    (h) Transformation operations:
-       copy, rotate-left/right, shift-left/right, sort
+       copy, rotate-left/right, shift-left/right, sort, nth_element
    (i) BLAS-L1:
        axpy, axpby, axpyz, axpbyz, axpbz
 
@@ -2734,11 +2734,11 @@
    (e) count       (f) copy
    (g) rotate-left (h) rotate-right
    (i) shift-left  (j) shift-right
-   (k) sort        (l) sumk
-   (m) axpy        (n) axpby
-   (o) axpyz       (p) axpbyz
-   (q) axpbz       (r) dot
-   (s) dotk
+   (k) sort        (l) nth_element
+   (m) sumk        (n) axpy
+   (o) axpby       (p) axpyz
+   (q) axpbyz      (r) axpbz
+   (s) dot         (t) dotk
 
      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~