add Riemann zeta function and Gauss error function

4 files changed, 53 insertions, 0 deletions

diff --git a/src/mp-equation.c b/src/mp-equation.c
index 4778298..df9d872 100644
--- a/src/mp-equation.c
+++ b/src/mp-equation.c
@@ -174,6 +174,7 @@ function_is_defined(ParserState *state, const char *name)
         strcmp(lower_name, "sin⁻¹") == 0 || strcmp(lower_name, "cos⁻¹") == 0 || strcmp(lower_name, "tan⁻¹") == 0 ||
         strcmp(lower_name, "sinh") == 0 || strcmp(lower_name, "cosh") == 0 || strcmp(lower_name, "tanh") == 0 ||
         strcmp(lower_name, "sinh⁻¹") == 0 || strcmp(lower_name, "cosh⁻¹") == 0 || strcmp(lower_name, "tanh⁻¹") == 0 ||
+        strcmp(lower_name, "erf") == 0 || strcmp(lower_name, "zeta") == 0 ||
         strcmp(lower_name, "asinh") == 0 || strcmp(lower_name, "acosh") == 0 || strcmp(lower_name, "atanh") == 0 ||
         strcmp(lower_name, "ones") == 0 ||
         strcmp(lower_name, "twos") == 0) {
@@ -261,6 +262,10 @@ get_function(ParserState *state, const char *name, const MPNumber *x, MPNumber *
         mp_acosh(x, z);
     else if (strcmp(lower_name, "tanh⁻¹") == 0 || strcmp(lower_name, "atanh") == 0)
         mp_atanh(x, z);
+    else if (strcmp(lower_name, "erf") == 0)
+        mp_erf(x, z);
+    else if (strcmp(lower_name, "zeta") == 0)
+        mp_zeta(x, z);
     else if (strcmp(lower_name, "ones") == 0)
         mp_ones_complement(x, state->options->wordlen, z);
     else if (strcmp(lower_name, "twos") == 0)
diff --git a/src/mp.c b/src/mp.c
index a25280b..1286f56 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -583,6 +583,40 @@ mp_xpowy_integer(const MPNumber *x, int64_t n, MPNumber *z)
     mpc_pow_si(z->num, x->num, (long) n, MPC_RNDNN);
 }
 
+void
+mp_erf(const MPNumber *x, MPNumber *z)
+{
+    if (mp_is_complex(x))
+    {   /* Translators: Error displayed when error function (erf) value is undefined */
+        mperr(_("The error function is only defined for real numbers"));
+        mp_set_from_integer(0, z);
+        return;
+    }
+
+    mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+    mpfr_erf(mpc_realref(z->num), mpc_realref(x->num), MPFR_RNDN);
+}
+
+void
+mp_zeta(const MPNumber *x, MPNumber *z)
+{
+    MPNumber one = mp_new();
+
+    mp_set_from_integer(1, &one);
+    if (mp_is_complex(x) || mp_compare(x, &one) == 0)
+    {   /* Translators: Error displayed when zeta function value is undefined */
+        mperr(_("The Riemann zeta function is only defined for real numbers ≠1"));
+        mp_set_from_integer(0, z);
+        mp_clear(&one);
+        return;
+    }
+
+    mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+    mpfr_zeta(mpc_realref(z->num), mpc_realref(x->num), MPFR_RNDN);
+
+    mp_clear(&one);
+}
+
 GList*
 mp_factorize(const MPNumber *x)
 {
diff --git a/src/mp.h b/src/mp.h
index 1748c03..f57ff16 100644
--- a/src/mp.h
+++ b/src/mp.h
@@ -298,6 +298,12 @@ void   mp_acosh(const MPNumber *x, MPNumber *z);
 /* Sets z = tanh⁻¹ x */
 void   mp_atanh(const MPNumber *x, MPNumber *z);
 
+/* Sets z to the value of the error function of x */
+void   mp_erf(const MPNumber *x, MPNumber *z);
+
+/* Sets z to the value of the Riemann Zeta function of x */
+void   mp_zeta(const MPNumber *x, MPNumber *z);
+
 /* Returns true if x is cannot be represented in a binary word of length 'wordlen' */
 bool   mp_is_overflow(const MPNumber *x, int wordlen);
 
diff --git a/src/test-mp-equation.c b/src/test-mp-equation.c
index becf742..861c445 100644
--- a/src/test-mp-equation.c
+++ b/src/test-mp-equation.c
@@ -590,6 +590,14 @@ test_equations(void)
     options.angle_units = MP_GRADIANS;
     test("sin 100", "1", 0);
 
+    test("zeta 0", "−0.5", 0);
+    test("zeta 2", "1.644934067", 0);
+    test("zeta −2", "0", 0);
+
+    test("erf 0", "0", 0);
+    test("erf 1", "0.842700793", 0);
+    test("erf −1", "−0.842700793", 0);
+
     /* Complex numbers */
     options.angle_units = MP_DEGREES;
     test("i", "i", 0);