diff --git a/README.md b/README.md
index e05f8fc..ad9f45a 100644
--- a/README.md
+++ b/README.md
@@ -3,7 +3,7 @@
 
  -#TinyExpr
+# TinyExpr
 
 TinyExpr is a very small recursive descent parser and evaluation engine for
 math expressions. It's handy when you want to add the ability to evaluation
@@ -12,7 +12,7 @@ math expressions at runtime without adding a bunch of cruft to you project.
 In addition to the standard math operators and precedence, TinyExpr also supports
 the standard C math functions and runtime binding of variables.
 
-##Features
+## Features
 
 - **ANSI C with no dependencies**.
 - Single source file and header file.
@@ -25,12 +25,12 @@ the standard C math functions and runtime binding of variables.
 - Easy to use and integrate with your code
 - Thread-safe, provided that your *malloc* is.
 
-##Building
+## Building
 
 TinyExpr is self-contained in two files: `tinyexpr.c` and `tinyexpr.h`. To use
 TinyExpr, simply add those two files to your project.
 
-##Short Example
+## Short Example
 
 Here is a minimal example to evaluate an expression at runtime.
 
@@ -40,7 +40,7 @@ Here is a minimal example to evaluate an expression at runtime.
 ```
 
 
-##Usage
+## Usage
 
 TinyExpr defines only four functions:
 
@@ -51,7 +51,7 @@ TinyExpr defines only four functions:
     void te_free(te_expr *expr);
 ```
 
-##te_interp
+## te_interp
 ```C
     double te_interp(const char *expression, int *error);
 ```
@@ -72,7 +72,7 @@ of the parse error on failure, and set `*error` to 0 on success.
     double c = te_interp("(5+5", &error); /* Returns NaN, error is set to 4. */
 ```
 
-##te_compile, te_eval, te_free
+## te_compile, te_eval, te_free
 ```C
     te_expr *te_compile(const char *expression, const te_variable *lookup, int lookup_len, int *error);
     double te_eval(const te_expr *n);
@@ -117,7 +117,7 @@ After you're finished, make sure to call `te_free()`.
 
 ```
 
-##Longer Example
+## Longer Example
 
 Here is a complete example that will evaluate an expression passed in from the command
 line. It also does error checking and binds the variables `x` and `y` to *3* and *4*, respectively.
@@ -178,7 +178,7 @@ This produces the output:
                 5.000000
 
 
-##Binding to Custom Functions
+## Binding to Custom Functions
 
 TinyExpr can also call to custom functions implemented in C. Here is a short example:
 
@@ -197,7 +197,7 @@ te_expr *n = te_compile("mysum(5, 6)", vars, 1, 0);
 ```
 
 
-##How it works
+## How it works
 
 `te_compile()` uses a simple recursive descent parser to compile your
 expression into a syntax tree. For example, the expression `"sin x + 1/4"`
@@ -216,7 +216,7 @@ and return the result of the expression.
 `te_free()` should always be called when you're done with the compiled expression.
 
 
-##Speed
+## Speed
 
 
 TinyExpr is pretty fast compared to C when the expression is short, when the
@@ -237,7 +237,7 @@ Here is some example performance numbers taken from the included
 
 
 
-##Grammar
+## Grammar
 
 TinyExpr parses the following grammar:
 
@@ -262,33 +262,39 @@ notation (e.g.  *1e3* for *1000*). A leading zero is not required (e.g. *.5*
 for *0.5*)
 
 
-##Functions supported
+## Functions supported
 
 TinyExpr supports addition (+), subtraction/negation (-), multiplication (\*),
 division (/), exponentiation (^) and modulus (%) with the normal operator
 precedence (the one exception being that exponentiation is evaluated
 left-to-right, but this can be changed - see below).
 
-In addition, the following C math functions are also supported:
+The following C math functions are also supported:
 
 - abs (calls to *fabs*), acos, asin, atan, atan2, ceil, cos, cosh, exp, floor, ln (calls to *log*), log (calls to *log10* by default, see below), log10, pow, sin, sinh, sqrt, tan, tanh
 
+The following functions are also built-in and provided by TinyExpr:
+
+- fac (factorials e.g. `fac 5` == 120)
+- ncr (combinations e.g. `ncr(6,2)` == 15)
+- npr (permutations e.g. `npr(6,2)` == 30)
+
 Also, the following constants are available:
 
 - `pi`, `e`
 
 
-##Compile-time options
+## Compile-time options
 
 
-By default, TinyExpr does exponentation from left to right. For example:
+By default, TinyExpr does exponentiation from left to right. For example:
 
 `a^b^c == (a^b)^c` and `-a^b == (-a)^b`
 
 This is by design. It's the way that spreadsheets do it (e.g. Excel, Google Sheets).
 
 
-If you would rather have exponentation work from right to left, you need to
+If you would rather have exponentiation work from right to left, you need to
 define `TE_POW_FROM_RIGHT` when compiling `tinyexpr.c`. There is a
 commented-out define near the top of that file. With this option enabled, the
 behaviour is:
@@ -300,7 +306,7 @@ That will match how many scripting languages do it (e.g. Python, Ruby).
 Also, if you'd like `log` to default to the natural log instead of `log10`,
 then you can define `TE_NAT_LOG`.
 
-##Hints
+## Hints
 
 - All functions/types start with the letters *te*.
 
diff --git a/test.c b/test.c
index b86a83f..c772950 100644
--- a/test.c
+++ b/test.c
@@ -213,6 +213,13 @@ void test_nans() {
         "1%0",
         "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;
@@ -234,6 +241,40 @@ void test_nans() {
 }
 
 
+void test_infs() {
+
+    const char *infs[] = {
+            "1/0",
+            "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;
+    for (i = 0; i < sizeof(infs) / sizeof(const char *); ++i) {
+        const char *expr = infs[i];
+
+        int err;
+        const double r = te_interp(expr, &err);
+        lequal(err, 0);
+        lok(r == r + 1);
+
+        te_expr *n = te_compile(expr, 0, 0, &err);
+        lok(n);
+        lequal(err, 0);
+        const double c = te_eval(n);
+        lok(c == c + 1);
+        te_free(n);
+    }
+}
+
+
 void test_variables() {
 
     double x, y, test;
@@ -587,17 +628,63 @@ 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},
+
+            {"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},
+    };
+
+
+    int i;
+    for (i = 0; i < sizeof(cases) / sizeof(test_case); ++i) {
+        const char *expr = cases[i].expr;
+        const double answer = cases[i].answer;
+
+        int err;
+        const double ev = te_interp(expr, &err);
+        lok(!err);
+        lfequal(ev, answer);
+
+        if (err) {
+            printf("FAILED: %s (%d)\n", expr, err);
+        }
+    }
+}
+
+
 int main(int argc, char *argv[])
 {
     lrun("Results", test_results);
     lrun("Syntax", test_syntax);
     lrun("NaNs", test_nans);
+    lrun("INFs", test_infs);
     lrun("Variables", test_variables);
     lrun("Functions", test_functions);
     lrun("Dynamic", test_dynamic);
     lrun("Closure", test_closure);
     lrun("Optimize", test_optimize);
     lrun("Pow", test_pow);
+    lrun("Combinatorics", test_combinatorics);
     lresults();
 
     return lfails != 0;
diff --git a/tinyexpr.c b/tinyexpr.c
index 5e50793..91a1848 100755
--- a/tinyexpr.c
+++ b/tinyexpr.c
@@ -39,11 +39,17 @@ For log = natural log uncomment the next line. */
 #include 
 #include 
 #include 
+#include 
 
 #ifndef NAN
 #define NAN (0.0/0.0)
 #endif
 
+#ifndef INFINITY
+#define INFINITY (1.0/0.0)
+#endif
+
+
 typedef double (*te_fun2)(double, double);
 
 enum {
@@ -113,6 +119,35 @@ void te_free(te_expr *n) {
 
 static double pi() {return 3.14159265358979323846;}
 static double e() {return 2.71828182845904523536;}
+static double fac(double a) {/* simplest version of fac */
+    if (a < 0.0)
+        return NAN;
+    if (a > UINT_MAX)
+        return INFINITY;
+    unsigned int ua = (unsigned int)(a);
+    unsigned long int result = 1, i;
+    for (i = 1; i <= ua; i++) {
+        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 */
@@ -126,6 +161,7 @@ static const te_variable functions[] = {
     {"cosh", cosh,    TE_FUNCTION1 | TE_FLAG_PURE, 0},
     {"e", e,          TE_FUNCTION0 | TE_FLAG_PURE, 0},
     {"exp", exp,      TE_FUNCTION1 | TE_FLAG_PURE, 0},
+    {"fac", fac,      TE_FUNCTION1 | TE_FLAG_PURE, 0},
     {"floor", floor,  TE_FUNCTION1 | TE_FLAG_PURE, 0},
     {"ln", log,       TE_FUNCTION1 | TE_FLAG_PURE, 0},
 #ifdef TE_NAT_LOG
@@ -134,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},
 
-#TinyExpr
+# TinyExpr
 
 TinyExpr is a very small recursive descent parser and evaluation engine for
 math expressions. It's handy when you want to add the ability to evaluation
@@ -12,7 +12,7 @@ math expressions at runtime without adding a bunch of cruft to you project.
 In addition to the standard math operators and precedence, TinyExpr also supports
 the standard C math functions and runtime binding of variables.
 
-##Features
+## Features
 
 - **ANSI C with no dependencies**.
 - Single source file and header file.
@@ -25,12 +25,12 @@ the standard C math functions and runtime binding of variables.
 - Easy to use and integrate with your code
 - Thread-safe, provided that your *malloc* is.
 
-##Building
+## Building
 
 TinyExpr is self-contained in two files: `tinyexpr.c` and `tinyexpr.h`. To use
 TinyExpr, simply add those two files to your project.
 
-##Short Example
+## Short Example
 
 Here is a minimal example to evaluate an expression at runtime.
 
@@ -40,7 +40,7 @@ Here is a minimal example to evaluate an expression at runtime.
 ```
 
 
-##Usage
+## Usage
 
 TinyExpr defines only four functions:
 
@@ -51,7 +51,7 @@ TinyExpr defines only four functions:
     void te_free(te_expr *expr);
 ```
 
-##te_interp
+## te_interp
 ```C
     double te_interp(const char *expression, int *error);
 ```
@@ -72,7 +72,7 @@ of the parse error on failure, and set `*error` to 0 on success.
     double c = te_interp("(5+5", &error); /* Returns NaN, error is set to 4. */
 ```
 
-##te_compile, te_eval, te_free
+## te_compile, te_eval, te_free
 ```C
     te_expr *te_compile(const char *expression, const te_variable *lookup, int lookup_len, int *error);
     double te_eval(const te_expr *n);
@@ -117,7 +117,7 @@ After you're finished, make sure to call `te_free()`.
 
 ```
 
-##Longer Example
+## Longer Example
 
 Here is a complete example that will evaluate an expression passed in from the command
 line. It also does error checking and binds the variables `x` and `y` to *3* and *4*, respectively.
@@ -178,7 +178,7 @@ This produces the output:
                 5.000000
 
 
-##Binding to Custom Functions
+## Binding to Custom Functions
 
 TinyExpr can also call to custom functions implemented in C. Here is a short example:
 
@@ -197,7 +197,7 @@ te_expr *n = te_compile("mysum(5, 6)", vars, 1, 0);
 ```
 
 
-##How it works
+## How it works
 
 `te_compile()` uses a simple recursive descent parser to compile your
 expression into a syntax tree. For example, the expression `"sin x + 1/4"`
@@ -216,7 +216,7 @@ and return the result of the expression.
 `te_free()` should always be called when you're done with the compiled expression.
 
 
-##Speed
+## Speed
 
 
 TinyExpr is pretty fast compared to C when the expression is short, when the
@@ -237,7 +237,7 @@ Here is some example performance numbers taken from the included
 
 
 
-##Grammar
+## Grammar
 
 TinyExpr parses the following grammar:
 
@@ -262,33 +262,39 @@ notation (e.g.  *1e3* for *1000*). A leading zero is not required (e.g. *.5*
 for *0.5*)
 
 
-##Functions supported
+## Functions supported
 
 TinyExpr supports addition (+), subtraction/negation (-), multiplication (\*),
 division (/), exponentiation (^) and modulus (%) with the normal operator
 precedence (the one exception being that exponentiation is evaluated
 left-to-right, but this can be changed - see below).
 
-In addition, the following C math functions are also supported:
+The following C math functions are also supported:
 
 - abs (calls to *fabs*), acos, asin, atan, atan2, ceil, cos, cosh, exp, floor, ln (calls to *log*), log (calls to *log10* by default, see below), log10, pow, sin, sinh, sqrt, tan, tanh
 
+The following functions are also built-in and provided by TinyExpr:
+
+- fac (factorials e.g. `fac 5` == 120)
+- ncr (combinations e.g. `ncr(6,2)` == 15)
+- npr (permutations e.g. `npr(6,2)` == 30)
+
 Also, the following constants are available:
 
 - `pi`, `e`
 
 
-##Compile-time options
+## Compile-time options
 
 
-By default, TinyExpr does exponentation from left to right. For example:
+By default, TinyExpr does exponentiation from left to right. For example:
 
 `a^b^c == (a^b)^c` and `-a^b == (-a)^b`
 
 This is by design. It's the way that spreadsheets do it (e.g. Excel, Google Sheets).
 
 
-If you would rather have exponentation work from right to left, you need to
+If you would rather have exponentiation work from right to left, you need to
 define `TE_POW_FROM_RIGHT` when compiling `tinyexpr.c`. There is a
 commented-out define near the top of that file. With this option enabled, the
 behaviour is:
@@ -300,7 +306,7 @@ That will match how many scripting languages do it (e.g. Python, Ruby).
 Also, if you'd like `log` to default to the natural log instead of `log10`,
 then you can define `TE_NAT_LOG`.
 
-##Hints
+## Hints
 
 - All functions/types start with the letters *te*.
 
diff --git a/test.c b/test.c
index b86a83f..c772950 100644
--- a/test.c
+++ b/test.c
@@ -213,6 +213,13 @@ void test_nans() {
         "1%0",
         "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;
@@ -234,6 +241,40 @@ void test_nans() {
 }
 
 
+void test_infs() {
+
+    const char *infs[] = {
+            "1/0",
+            "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;
+    for (i = 0; i < sizeof(infs) / sizeof(const char *); ++i) {
+        const char *expr = infs[i];
+
+        int err;
+        const double r = te_interp(expr, &err);
+        lequal(err, 0);
+        lok(r == r + 1);
+
+        te_expr *n = te_compile(expr, 0, 0, &err);
+        lok(n);
+        lequal(err, 0);
+        const double c = te_eval(n);
+        lok(c == c + 1);
+        te_free(n);
+    }
+}
+
+
 void test_variables() {
 
     double x, y, test;
@@ -587,17 +628,63 @@ 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},
+
+            {"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},
+    };
+
+
+    int i;
+    for (i = 0; i < sizeof(cases) / sizeof(test_case); ++i) {
+        const char *expr = cases[i].expr;
+        const double answer = cases[i].answer;
+
+        int err;
+        const double ev = te_interp(expr, &err);
+        lok(!err);
+        lfequal(ev, answer);
+
+        if (err) {
+            printf("FAILED: %s (%d)\n", expr, err);
+        }
+    }
+}
+
+
 int main(int argc, char *argv[])
 {
     lrun("Results", test_results);
     lrun("Syntax", test_syntax);
     lrun("NaNs", test_nans);
+    lrun("INFs", test_infs);
     lrun("Variables", test_variables);
     lrun("Functions", test_functions);
     lrun("Dynamic", test_dynamic);
     lrun("Closure", test_closure);
     lrun("Optimize", test_optimize);
     lrun("Pow", test_pow);
+    lrun("Combinatorics", test_combinatorics);
     lresults();
 
     return lfails != 0;
diff --git a/tinyexpr.c b/tinyexpr.c
index 5e50793..91a1848 100755
--- a/tinyexpr.c
+++ b/tinyexpr.c
@@ -39,11 +39,17 @@ For log = natural log uncomment the next line. */
 #include 
 #include 
 #include 
+#include 
 
 #ifndef NAN
 #define NAN (0.0/0.0)
 #endif
 
+#ifndef INFINITY
+#define INFINITY (1.0/0.0)
+#endif
+
+
 typedef double (*te_fun2)(double, double);
 
 enum {
@@ -113,6 +119,35 @@ void te_free(te_expr *n) {
 
 static double pi() {return 3.14159265358979323846;}
 static double e() {return 2.71828182845904523536;}
+static double fac(double a) {/* simplest version of fac */
+    if (a < 0.0)
+        return NAN;
+    if (a > UINT_MAX)
+        return INFINITY;
+    unsigned int ua = (unsigned int)(a);
+    unsigned long int result = 1, i;
+    for (i = 1; i <= ua; i++) {
+        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 */
@@ -126,6 +161,7 @@ static const te_variable functions[] = {
     {"cosh", cosh,    TE_FUNCTION1 | TE_FLAG_PURE, 0},
     {"e", e,          TE_FUNCTION0 | TE_FLAG_PURE, 0},
     {"exp", exp,      TE_FUNCTION1 | TE_FLAG_PURE, 0},
+    {"fac", fac,      TE_FUNCTION1 | TE_FLAG_PURE, 0},
     {"floor", floor,  TE_FUNCTION1 | TE_FLAG_PURE, 0},
     {"ln", log,       TE_FUNCTION1 | TE_FLAG_PURE, 0},
 #ifdef TE_NAT_LOG
@@ -134,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},