diff --git a/test.c b/test.c
index 82611eb..c772950 100644
--- a/test.c
+++ b/test.c
@@ -214,6 +214,12 @@ void test_nans() {
         "1%(1%0)",
         "(1%0)%1",
         "fac(-1)",
+        "ncr(2, 4)",
+        "ncr(-2, 4)",
+        "ncr(2, -4)",
+        "npr(2, 4)",
+        "npr(-2, 4)",
+        "npr(2, -4)",
     };
 
     int i;
@@ -242,6 +248,12 @@ void test_infs() {
             "log(0)",
             "pow(2,10000000)",
             "fac(300)",
+            "ncr(300,100)",
+            "ncr(300000,100)",
+            "ncr(300000,100)*8",
+            "npr(3,2)*ncr(300000,100)",
+            "npr(100,90)",
+            "npr(30,25)",
     };
 
     int i;
@@ -618,13 +630,28 @@ void test_pow() {
 
 void test_combinatorics() {
     test_case cases[] = {
-            {"fac(0) ", 1},
-            {"fac(0.2) ", 1},
-            {"fac(1) ", 1},
-            {"fac(2) ", 2},
-            {"fac(3) ", 6},
-            {"fac(4.8) ", 24},
-            {"fac(10) ", 3628800},
+            {"fac(0)", 1},
+            {"fac(0.2)", 1},
+            {"fac(1)", 1},
+            {"fac(2)", 2},
+            {"fac(3)", 6},
+            {"fac(4.8)", 24},
+            {"fac(10)", 3628800},
+
+            {"ncr(0,0)", 1},
+            {"ncr(10,1)", 10},
+            {"ncr(10,0)", 1},
+            {"ncr(10,10)", 1},
+            {"ncr(16,7)", 11440},
+            {"ncr(16,9)", 11440},
+            {"ncr(100,95)", 75287520},
+
+            {"npr(0,0)", 1},
+            {"npr(10,1)", 10},
+            {"npr(10,0)", 1},
+            {"npr(10,10)", 3628800},
+            {"npr(20,5)", 1860480},
+            {"npr(100,4)", 94109400},
     };
 
 
diff --git a/tinyexpr.c b/tinyexpr.c
index f2a6175..115ea11 100755
--- a/tinyexpr.c
+++ b/tinyexpr.c
@@ -125,14 +125,29 @@ static double fac(double a) {/* simplest version of fac */
     if (a > UINT_MAX)
         return INFINITY;
     unsigned int ua = (unsigned int)(a);
-    unsigned long int result = 1, i = 1;
+    unsigned long int result = 1, i;
     for (i = 1; i <= ua; i++) {
-        if (ua > ULONG_MAX / result)
+        if (i > ULONG_MAX / result)
             return INFINITY;
         result *= i;
     }
     return (double)result;
 }
+static double ncr(double n, double r) {
+    if (n < 0.0 || r < 0.0 || n < r) return NAN;
+    if (n > UINT_MAX || r > UINT_MAX) return INFINITY;
+    unsigned long int un = (unsigned int)(n), ur = (unsigned int)(r), i;
+    unsigned long int result = 1;
+    if (ur > un / 2) ur = un - ur;
+    for (i = 1; i <= ur; i++) {
+        if (result > ULONG_MAX / (un - ur + i))
+            return INFINITY;
+        result *= un - ur + i;
+        result /= i;
+    }
+    return result;
+}
+static double npr(double n, double r) {return ncr(n, r) * fac(r);}
 
 static const te_variable functions[] = {
     /* must be in alphabetical order */
@@ -155,6 +170,8 @@ static const te_variable functions[] = {
     {"log", log10,    TE_FUNCTION1 | TE_FLAG_PURE, 0},
 #endif
     {"log10", log10,  TE_FUNCTION1 | TE_FLAG_PURE, 0},
+    {"ncr", ncr,      TE_FUNCTION2 | TE_FLAG_PURE, 0},
+    {"npr", npr,      TE_FUNCTION2 | TE_FLAG_PURE, 0},
     {"pi", pi,        TE_FUNCTION0 | TE_FLAG_PURE, 0},
     {"pow", pow,      TE_FUNCTION2 | TE_FLAG_PURE, 0},
     {"sin", sin,      TE_FUNCTION1 | TE_FLAG_PURE, 0},