Implement factorial powered by the Gamma function.

src/operator_functions.cpp

@@ -7,6 +7,22 @@ FUNCTION_MPFR2(divide, div)
 FUNCTION_MPFR2(pow, pow);
 
 FUNCTION_MPFR(negate, neg);
+FUNCTION(factorial) {
+    number* left = (number*) libab_unwrap_param(params, 0);   
+    mpfr_t plus_one;
+    mpfr_t gamma;
+
+    mpfr_init2(plus_one, PRECISION);
+    mpfr_init2(gamma, PRECISION);
+    mpfr_add_ui(plus_one, left->value, 1, MPFR_RNDN);
+    mpfr_gamma(gamma, plus_one, MPFR_RNDN);
+
+    number* to_return = new number(std::move(gamma));
+    ref return_value = create_value<number>(ab, to_return);
+    libab_ref_copy(return_value, into);
+
+    return LIBAB_SUCCESS;
+}
 
 FUNCTION_COMPARE(lt, <);
 FUNCTION_COMPARE(lte, <=);