diff --git a/Makefile b/Makefile
index a49b7ef..13edf45 100644
--- a/Makefile
+++ b/Makefile
@@ -39,6 +39,7 @@ BUILD_LIST+=exprtk_simple_example_12
 BUILD_LIST+=exprtk_simple_example_13
 BUILD_LIST+=exprtk_simple_example_14
 BUILD_LIST+=exprtk_simple_example_15
+BUILD_LIST+=exprtk_simple_example_16
 
 all: $(BUILD_LIST)
 
@@ -93,7 +94,10 @@ exprtk_simple_example_14: exprtk_simple_example_14.cpp exprtk.hpp
 exprtk_simple_example_15: exprtk_simple_example_15.cpp exprtk.hpp
 	$(COMPILER) $(OPTIONS) exprtk_simple_example_15 exprtk_simple_example_15.cpp $(LINKER_OPT)
 
-pgo: exprtk_test.cpp exprtk_benchmark.cpp exprtk.hpp
+exprtk_simple_example_16: exprtk_simple_example_16.cpp exprtk.hpp
+	$(COMPILER) $(OPTIONS) exprtk_simple_example_16 exprtk_simple_example_16.cpp $(LINKER_OPT)
+
+pgo: exprtk_benchmark.cpp exprtk.hpp
 	$(COMPILER) $(BASE_OPTIONS) -O3 -march=native -fprofile-generate -o exprtk_benchmark exprtk_benchmark.cpp $(LINKER_OPT)
 	./exprtk_benchmark
 	$(COMPILER) $(BASE_OPTIONS) -O3 -march=native -fprofile-use -o exprtk_benchmark exprtk_benchmark.cpp $(LINKER_OPT)
@@ -116,6 +120,7 @@ strip_bin:
 	strip -s exprtk_simple_example_13
 	strip -s exprtk_simple_example_14
 	strip -s exprtk_simple_example_15
+	strip -s exprtk_simple_example_16
 
 valgrind_check:
 	valgrind --leak-check=full --show-reachable=yes --track-origins=yes --log-file=exprtk_test_valgrind.log -v ./exprtk_test
@@ -135,6 +140,7 @@ valgrind_check:
 	valgrind --leak-check=full --show-reachable=yes --track-origins=yes --log-file=exprtk_simple_example_13_valgrind.log -v ./exprtk_simple_example_13
 	valgrind --leak-check=full --show-reachable=yes --track-origins=yes --log-file=exprtk_simple_example_14_valgrind.log -v ./exprtk_simple_example_14
 	valgrind --leak-check=full --show-reachable=yes --track-origins=yes --log-file=exprtk_simple_example_15_valgrind.log -v ./exprtk_simple_example_15	
+	valgrind --leak-check=full --show-reachable=yes --track-origins=yes --log-file=exprtk_simple_example_16_valgrind.log -v ./exprtk_simple_example_16	
 
 clean:
 	rm -f core.* *~ *.o *.bak *stackdump gmon.out *.gcda *.gcno *.gcnor *.gch
diff --git a/exprtk.hpp b/exprtk.hpp
index 360f017..dbf43c4 100644
--- a/exprtk.hpp
+++ b/exprtk.hpp
@@ -3280,7 +3280,7 @@ namespace exprtk
          e_sf4ext44 = 2044, e_sf4ext45 = 2045, e_sf4ext46 = 2046, e_sf4ext47 = 2047,
          e_sf4ext48 = 2048, e_sf4ext49 = 2049, e_sf4ext50 = 2050, e_sf4ext51 = 2051,
          e_sf4ext52 = 2052, e_sf4ext53 = 2053, e_sf4ext54 = 2054, e_sf4ext55 = 2055,
-         e_sf4ext56 = 2056, e_sf4ext57 = 2057, e_sf4ext58 = 2058
+         e_sf4ext56 = 2056, e_sf4ext57 = 2057, e_sf4ext58 = 2058, e_sf4ext59 = 2059
       };
 
       struct base_operation_t
@@ -5974,6 +5974,7 @@ namespace exprtk
       template <typename T> struct sfext56_op : public sf_base<T> { typedef typename sf_base<T>::Type Type; static inline T process(Type x, Type y, Type z, Type w) { return (x - y) * (z - w); } static inline std::string id() { return "(t-t)*(t-t)";} };
       template <typename T> struct sfext57_op : public sf_base<T> { typedef typename sf_base<T>::Type Type; static inline T process(Type x, Type y, Type z, Type w) { return (x - y) + (z - w); } static inline std::string id() { return "(t-t)+(t-t)";} };
       template <typename T> struct sfext58_op : public sf_base<T> { typedef typename sf_base<T>::Type Type; static inline T process(Type x, Type y, Type z, Type w) { return (x - y) - (z - w); } static inline std::string id() { return "(t-t)-(t-t)";} };
+      template <typename T> struct sfext59_op : public sf_base<T> { typedef typename sf_base<T>::Type Type; static inline T process(Type x, Type y, Type z, Type w) { return (x / y) + (z * w); } static inline std::string id() { return "(t/t)+(t*t)";} };
 
       template <typename T, typename SpecialFunction>
       class sf3_node : public trinary_node<T>
@@ -19723,7 +19724,7 @@ namespace exprtk
                   case_stmt(details::e_sf4ext52,details::sfext52_op) case_stmt(details::e_sf4ext53,details::sfext53_op)
                   case_stmt(details::e_sf4ext54,details::sfext54_op) case_stmt(details::e_sf4ext55,details::sfext55_op)
                   case_stmt(details::e_sf4ext56,details::sfext56_op) case_stmt(details::e_sf4ext57,details::sfext57_op)
-                  case_stmt(details::e_sf4ext58,details::sfext58_op)
+                  case_stmt(details::e_sf4ext58,details::sfext58_op) case_stmt(details::e_sf4ext59,details::sfext59_op)
 
                   #undef case_stmt
                   default : return error_node();
@@ -23920,7 +23921,7 @@ namespace exprtk
          register_sf4ext(44) register_sf4ext(45) register_sf4ext(46) register_sf4ext(47)
          register_sf4ext(48) register_sf4ext(49) register_sf4ext(50) register_sf4ext(51)
          register_sf4ext(52) register_sf4ext(53) register_sf4ext(54) register_sf4ext(55)
-         register_sf4ext(56) register_sf4ext(57) register_sf4ext(58)
+         register_sf4ext(56) register_sf4ext(57) register_sf4ext(58) register_sf4ext(59)
          #undef register_sf4ext
       }
 
diff --git a/exprtk_benchmark.cpp b/exprtk_benchmark.cpp
index cb3df48..130fec7 100644
--- a/exprtk_benchmark.cpp
+++ b/exprtk_benchmark.cpp
@@ -19,6 +19,7 @@
 #include <cstdio>
 #include <cmath>
 #include <iostream>
+#include <fstream>
 #include <string>
 #include <deque>
 
@@ -269,9 +270,20 @@ struct native
 };
 
 void pgo_primer();
+void perform_file_based_benchmark(const std::string& file_name, const std::size_t& rounds = 100000);
 
-int main()
+int main(int argc, char* argv[])
 {
+   if (argc >= 2)
+   {
+      const std::string file_name = argv[1];
+      if (argc == 2)
+         perform_file_based_benchmark(file_name);
+      else
+         perform_file_based_benchmark(file_name,atoi(argv[2]));
+      return 0;
+   }
+
    pgo_primer();
 
    double x = 0;
@@ -362,3 +374,158 @@ void pgo_primer()
       }
    }
 }
+
+std::size_t load_expression_file(const std::string& file_name, std::deque<std::string>& expression_list)
+{
+   std::ifstream stream(file_name.c_str());
+   if (!stream) return 0;
+   std::string buffer;
+   buffer.reserve(1024);
+   std::size_t line_count = 0;
+   while (std::getline(stream,buffer))
+   {
+      if (buffer.empty())
+         continue;
+      else if ('#' == buffer[0])
+         continue;
+      ++line_count;
+      expression_list.push_back(buffer);
+   }
+
+   return line_count;
+}
+
+void perform_file_based_benchmark(const std::string& file_name, const std::size_t& rounds)
+{
+   std::deque<std::string> expr_str_list;
+   if (0 == load_expression_file(file_name,expr_str_list))
+   {
+      std::cout << "Failed to load any expressions from: " << file_name << "\n";
+      return;
+   }
+
+   typedef exprtk::symbol_table<double> symbol_table_t;
+   typedef exprtk::expression<double>     expression_t;
+   typedef exprtk::parser<double>             parser_t;
+
+   std::deque<expression_t> expression_list;
+
+   symbol_table_t symbol_table;
+
+   double a = 1.1;
+   double b = 2.2;
+   double c = 3.3;
+   double x = 2.123456;
+   double y = 3.123456;
+   double z = 4.123456;
+   double w = 5.123456;
+
+   symbol_table.add_variable("a", a);
+   symbol_table.add_variable("b", b);
+   symbol_table.add_variable("c", c);
+
+   symbol_table.add_variable("x", x);
+   symbol_table.add_variable("y", y);
+   symbol_table.add_variable("z", z);
+   symbol_table.add_variable("w", w);
+
+   exprtk::polynomial<double, 1> poly01;
+   exprtk::polynomial<double, 2> poly02;
+   exprtk::polynomial<double, 3> poly03;
+   exprtk::polynomial<double, 4> poly04;
+   exprtk::polynomial<double, 5> poly05;
+   exprtk::polynomial<double, 6> poly06;
+   exprtk::polynomial<double, 7> poly07;
+   exprtk::polynomial<double, 8> poly08;
+   exprtk::polynomial<double, 9> poly09;
+   exprtk::polynomial<double,10> poly10;
+   exprtk::polynomial<double,11> poly11;
+   exprtk::polynomial<double,12> poly12;
+
+   symbol_table.add_function("poly01", poly01);
+   symbol_table.add_function("poly02", poly02);
+   symbol_table.add_function("poly03", poly03);
+   symbol_table.add_function("poly04", poly04);
+   symbol_table.add_function("poly05", poly05);
+   symbol_table.add_function("poly06", poly06);
+   symbol_table.add_function("poly07", poly07);
+   symbol_table.add_function("poly08", poly08);
+   symbol_table.add_function("poly09", poly09);
+   symbol_table.add_function("poly10", poly10);
+   symbol_table.add_function("poly11", poly11);
+   symbol_table.add_function("poly12", poly12);
+
+   static double e = exprtk::details::numeric::constant::e;
+   symbol_table.add_variable("e", e, true);
+
+   symbol_table.add_constants();
+
+   {
+      parser_t parser;
+      for (std::size_t i = 0; i < expr_str_list.size(); ++i)
+      {
+         expression_t expression;
+         expression.register_symbol_table(symbol_table);
+         if (!parser.compile(expr_str_list[i],expression))
+         {
+            printf("[perform_file_based_benchmark] - Parser Error: %s\tExpression: %s\n",
+                   parser.error().c_str(),
+                   expr_str_list[i].c_str());
+            return;
+         }
+         expression_list.push_back(expression);
+      }
+   }
+
+   exprtk::timer total_timer;
+
+   double single_eval_total_time = 0.0;
+
+   total_timer.start();
+   for (std::size_t i = 0; i < expression_list.size(); ++i)
+   {
+      expression_t& e = expression_list[i];
+
+      exprtk::timer timer;
+
+      a = 1.1;
+      b = 2.2;
+      c = 3.3;
+      x = 2.123456;
+      y = 3.123456;
+      z = 4.123456;
+      w = 5.123456;
+
+      timer.start();
+      double sum = 0.0;
+      for (std::size_t r = 0; r < rounds; ++r)
+      {
+         sum += e.value();
+         std::swap(a,b);
+         std::swap(x,y);
+      }
+      timer.stop();
+
+      printf("Expression %3d of %3d %9.3f ns\t%10d ns\t(%30.10f)  '%s'\n",
+             static_cast<int>(i + 1),
+             static_cast<int>(expression_list.size()),
+             (timer.time() * 1000000000.0) / (1.0 * rounds),
+             static_cast<unsigned int>(timer.time() * 1000000000.0),
+             sum,
+             expr_str_list[i].c_str());
+
+      fflush(stdout);
+
+      single_eval_total_time += (timer.time() * 1000000000.0) / (1.0 * rounds);
+   }
+   total_timer.stop();
+
+   printf("[*] Number Of Evals:        %15.0f\n",
+          rounds * (expression_list.size() * 1.0));
+
+   printf("[*] Total Time:             %9.3fsec\n",
+          total_timer.time());
+
+   printf("[*] Total Single Eval Time: %9.3fms\n",
+          single_eval_total_time / 1000000.0);
+}
diff --git a/exprtk_simple_example_01.cpp b/exprtk_simple_example_01.cpp
index 259eebc..3e04617 100644
--- a/exprtk_simple_example_01.cpp
+++ b/exprtk_simple_example_01.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void trig_function()
 {
    std::string expression_string = "clamp(-1.0,sin(2 * pi * x) + cos(x / 2 * pi),+1.0)";
diff --git a/exprtk_simple_example_02.cpp b/exprtk_simple_example_02.cpp
index bccc05b..fbc422c 100644
--- a/exprtk_simple_example_02.cpp
+++ b/exprtk_simple_example_02.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void square_wave()
 {
    std::string expr_string = "a*(4/pi)*"
diff --git a/exprtk_simple_example_03.cpp b/exprtk_simple_example_03.cpp
index 3c7cbba..54dc531 100644
--- a/exprtk_simple_example_03.cpp
+++ b/exprtk_simple_example_03.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void polynomial()
 {
    std::string expression_string = "25x^5 - 35x^4 - 15x^3 + 40x^2 - 15x + 1";
diff --git a/exprtk_simple_example_05.cpp b/exprtk_simple_example_05.cpp
index 35d706c..c9b3129 100644
--- a/exprtk_simple_example_05.cpp
+++ b/exprtk_simple_example_05.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 struct myfunc : public exprtk::ifunction<T>
 {
    myfunc()
diff --git a/exprtk_simple_example_06.cpp b/exprtk_simple_example_06.cpp
index 2a126e2..89faa9f 100644
--- a/exprtk_simple_example_06.cpp
+++ b/exprtk_simple_example_06.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void vector_function()
 {
    typedef exprtk::symbol_table<T> symbol_table_t;
diff --git a/exprtk_simple_example_11.cpp b/exprtk_simple_example_11.cpp
index e6d74a8..3f35bae 100644
--- a/exprtk_simple_example_11.cpp
+++ b/exprtk_simple_example_11.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void square_wave2()
 {
    typedef exprtk::symbol_table<T> symbol_table_t;
@@ -30,7 +30,7 @@ void square_wave2()
 
    std::string wave_program =
                   " var r := 0;                                         "
-                  " for(i := 0; i < 1000; i += 1)                       "
+                  " for(var i := 0; i < 1000; i += 1)                   "
                   " {                                                   "
                   "   r += (1 / (2i + 1)) * sin((4i + 2) * pi * f * t); "
                   " };                                                  "
diff --git a/exprtk_simple_example_12.cpp b/exprtk_simple_example_12.cpp
index d7a90f8..295a403 100644
--- a/exprtk_simple_example_12.cpp
+++ b/exprtk_simple_example_12.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void bubble_sort()
 {
    typedef exprtk::symbol_table<T> symbol_table_t;
diff --git a/exprtk_simple_example_13.cpp b/exprtk_simple_example_13.cpp
index 69b4523..3b12554 100644
--- a/exprtk_simple_example_13.cpp
+++ b/exprtk_simple_example_13.cpp
@@ -23,7 +23,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void savitzky_golay_filter()
 {
    typedef exprtk::symbol_table<T> symbol_table_t;
diff --git a/exprtk_simple_example_14.cpp b/exprtk_simple_example_14.cpp
index d5e427c..5ed5547 100644
--- a/exprtk_simple_example_14.cpp
+++ b/exprtk_simple_example_14.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void stddev_example()
 {
    typedef exprtk::expression<T> expression_t;
diff --git a/exprtk_simple_example_15.cpp b/exprtk_simple_example_15.cpp
index f735f72..9d24b5b 100644
--- a/exprtk_simple_example_15.cpp
+++ b/exprtk_simple_example_15.cpp
@@ -21,7 +21,7 @@
 #include "exprtk.hpp"
 
 
-template<typename T>
+template <typename T>
 void black_scholes_merton_model()
 {
    typedef exprtk::symbol_table<T> symbol_table_t;
diff --git a/exprtk_simple_example_16.cpp b/exprtk_simple_example_16.cpp
new file mode 100644
index 0000000..ac3d19b
--- /dev/null
+++ b/exprtk_simple_example_16.cpp
@@ -0,0 +1,84 @@
+/*
+ **************************************************************
+ *         C++ Mathematical Expression Toolkit Library        *
+ *                                                            *
+ * Simple Example 16                                          *
+ * Author: Arash Partow (1999-2014)                           *
+ * URL: http://www.partow.net/programming/exprtk/index.html   *
+ *                                                            *
+ * Copyright notice:                                          *
+ * Free use of the Mathematical Expression Toolkit Library is *
+ * permitted under the guidelines and in accordance with the  *
+ * most current version of the Common Public License.         *
+ * http://www.opensource.org/licenses/cpl1.0.php              *
+ *                                                            *
+ **************************************************************
+*/
+
+
+#include <cstdio>
+#include <cstdlib>
+#include <string>
+#include "exprtk.hpp"
+
+
+template <typename T>
+void linear_least_squares()
+{
+   typedef exprtk::symbol_table<T> symbol_table_t;
+   typedef exprtk::expression<T>     expression_t;
+   typedef exprtk::parser<T>             parser_t;
+
+   std::string linear_least_squares_program =
+                  " if (x[] == y[])                                       "
+                  " {                                                     "
+                  "   beta  := (sum(x * y) - sum(x) * sum(y) / x[]) /     "
+                  "            (sum(x^2) -  sum(x)^2 / x[]);              "
+                  "                                                       "
+                  "   alpha := avg(y) - beta * avg(x);                    "
+                  "                                                       "
+                  "   rmse  := sqrt(sum((beta * x + alpha - y)^2) / y[]); "
+                  " }                                                     "
+                  " else                                                  "
+                  " {                                                     "
+                  "   alpha := null;                                      "
+                  "   beta  := null;                                      "
+                  "   rmse  := null;                                      "
+                  " }                                                     ";
+
+   T x[] = {T(  1), T(  2), T(3), T(  4), T(  5), T(6), T(  7), T(  8), T(  9), T(10)};
+   T y[] = {T(8.7), T(6.8), T(6), T(5.6), T(3.8), T(3), T(2.4), T(1.7), T(0.4), T(-1)};
+
+   T alpha = T(0);
+   T beta  = T(0);
+   T rmse  = T(0);
+
+   symbol_table_t symbol_table;
+
+   symbol_table.add_variable("alpha",alpha);
+   symbol_table.add_variable("beta" ,beta );
+   symbol_table.add_variable("rmse" ,rmse );
+
+   symbol_table.add_vector("x",x);
+   symbol_table.add_vector("y",y);
+
+   expression_t expression;
+   expression.register_symbol_table(symbol_table);
+
+   parser_t parser;
+
+   parser.compile(linear_least_squares_program,expression);
+
+   expression.value();
+
+   printf("alpha: %15.12f\n",alpha);
+   printf("beta:  %15.12f\n",beta );
+   printf("rmse:  %15.12f\n",rmse );
+   printf("y = %15.12fx + %15.12f\n",beta,alpha);
+}
+
+int main()
+{
+   linear_least_squares<double>();
+   return 0;
+}
diff --git a/readme.txt b/readme.txt
index 44500ee..f8c3898 100644
--- a/readme.txt
+++ b/readme.txt
@@ -1346,6 +1346,7 @@ in a compilation failure.
 (17) exprtk_simple_example_13.cpp
 (18) exprtk_simple_example_14.cpp
 (19) exprtk_simple_example_15.cpp
+(20) exprtk_simple_example_16.cpp