[genius] Thu Feb 26 14:30:40 2015 Jiri (George) Lebl <jirka 5z com>

[George Lebl <jirka src gnome org>]

commit 26b89bb72eca553bab9abe4edd4877a2fa44b2c7
Author: Jiri (George) Lebl <jiri lebl gmail com>
Date:   Thu Feb 26 14:31:06 2015 -0600

    Thu Feb 26 14:30:40 2015  Jiri (George) Lebl <jirka 5z com>
    
        * lib/functions/elementary.gel, src/funclib.c: Use MPFR version of
          erf for real values, and move the hacky computation using power
          series as gel snippet inside funclib.c.  Should be replaced by
          a proper computation

 ChangeLog                    |    7 +++++
 lib/functions/elementary.gel |   56 ++++++++++++++++++++++-------------------
 lib/library-strings.c        |    3 +-
 src/funclib.c                |   56 ++++++++++++++++++++++++++++++------------
 4 files changed, 78 insertions(+), 44 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 13510aa..5b20f2d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+Thu Feb 26 14:30:40 2015  Jiri (George) Lebl <jirka 5z com>
+
+       * lib/functions/elementary.gel, src/funclib.c: Use MPFR version of
+         erf for real values, and move the hacky computation using power
+         series as gel snippet inside funclib.c.  Should be replaced by
+         a proper computation
+
 Mon Feb 16 14:26:02 2015  Jiri (George) Lebl <jirka 5z com>
 
        * examples/fourier-series-plotting.gel: Replace the hand coded
diff --git a/lib/functions/elementary.gel b/lib/functions/elementary.gel
index 9c3a315..176ad3e 100644
--- a/lib/functions/elementary.gel
+++ b/lib/functions/elementary.gel
@@ -233,34 +233,38 @@ function log(x,b...) = (
        ln(x)/ln(base)
 );
 
+# This is still used for complex values in the hacky computation in funclib.c
+# although for real values we use MPFR
 parameter ErrorFunctionTolerance = 10.0^(-10);
-SetHelp ("ErrorFunctionTolerance", "parameters", "Tolerance of the ErrorFunction")
+SetHelp ("ErrorFunctionTolerance", "parameters", "Tolerance of the ErrorFunction (used for complex values 
only)")
 
-SetHelp("ErrorFunction","functions","The error function, 2/sqrt(pi) * int_0^x e^(-t^2) dt")
-function ErrorFunction (x) = (
-       if IsMatrix (x) then
-               return ApplyOverMatrix (x, ErrorFunction)
-       else if not IsValue (x) then
-               (error("ErrorFunction: argument not a value");bailout);
-       twosqrtpi = 2/sqrt(pi);
-       a = 1;
-       s = 0;
-       n = 0;
-       f = 1;
-       xx = x;
-       xsq = x^2;
-       do (
-               t = xx * a * twosqrtpi;
-               s = s + t;
-               increment n;
-               f = f * n;
-               a = ((-1)^n) / (((2*n)+1) * f);
-               xx = xx * xsq
-       ) while (|t| > ErrorFunctionTolerance);
-       s
-);
-erf = ErrorFunction
-SetHelpAlias ("ErrorFunction", "erf");
+# This is actually done for complex values inside funclib.c
+# This should as some point be replaced by a proper version of erf
+#SetHelp("ErrorFunction","functions","The error function, 2/sqrt(pi) * int_0^x e^(-t^2) dt")
+#function ErrorFunction (x) = (
+#      if IsMatrix (x) then
+#              return ApplyOverMatrix (x, ErrorFunction)
+#      else if not IsValue (x) then
+#              (error("ErrorFunction: argument not a value");bailout);
+#      twosqrtpi = 2/sqrt(pi);
+#      a = 1;
+#      s = 0;
+#      n = 0;
+#      f = 1;
+#      xx = x;
+#      xsq = x^2;
+#      do (
+#              t = xx * a * twosqrtpi;
+#              s = s + t;
+#              increment n;
+#              f = f * n;
+#              a = ((-1)^n) / (((2*n)+1) * f);
+#              xx = xx * xsq
+#      ) while (|t| > ErrorFunctionTolerance);
+#      s
+#);
+#erf = ErrorFunction
+#SetHelpAlias ("ErrorFunction", "erf");
 
 #FIXME: Should probably be in a separate source file
 SetHelp("NewtonsMethodPoly","polynomial","Attempt to find a root of a polynomial using Newton's method, 
returning after two successive values are within epsilon or after maxn tries (then returns null)")
diff --git a/lib/library-strings.c b/lib/library-strings.c
index 35bb600..d2addcb 100644
--- a/lib/library-strings.c
+++ b/lib/library-strings.c
@@ -9,7 +9,7 @@ char *fake = N_("Tolerance for continuity of functions and for calculating the l
 char *fake = N_("How many iterations to try to find the limit for derivative");
 char *fake = N_("How many successive steps to be within tolerance for calculation of derivative");
 char *fake = N_("Tolerance for calculating the derivatives of functions");
-char *fake = N_("Tolerance of the ErrorFunction");
+char *fake = N_("Tolerance of the ErrorFunction (used for complex values only)");
 char *fake = N_("Tolerance of the GaussDistribution function");
 char *fake = N_("The function used for numerical integration in NumericalIntegral");
 char *fake = N_("Steps to perform in NumericalIntegral");
@@ -214,7 +214,6 @@ char *fake = N_("Compute two-sided derivative using three-point formula");
 char *fake = N_("argument (angle) of complex number");
 char *fake = N_("Dirichlet kernel of order n");
 char *fake = N_("Returns 1 if and only if all elements are zero");
-char *fake = N_("The error function, 2/sqrt(pi) * int_0^x e^(-t^2) dt");
 char *fake = N_("Fejer kernel of order n");
 char *fake = N_("Returns 1 if and only if all elements are equal");
 char *fake = N_("Principal branch of the Lambert W function for real values greater than or equal to -1/e");
diff --git a/src/funclib.c b/src/funclib.c
index 4132b3f..4a419cb 100644
--- a/src/funclib.c
+++ b/src/funclib.c
@@ -42,6 +42,7 @@
 
 /* FIXME:static GelEFunc *_internal_ln_function = NULL; */
 static GelEFunc *_internal_exp_function = NULL;
+static GelEFunc *_internal_erf_function = NULL;
 
 static GelEFunc *conj_function = NULL;
 static GelEFunc *sin_function = NULL;
@@ -65,7 +66,7 @@ static GelEFunc *Numerator_function = NULL;
 static GelEFunc *Denominator_function = NULL;
 static GelEFunc *Re_function = NULL;
 static GelEFunc *Im_function = NULL;
-/*static GelEFunc *ErrorFunction_function = NULL;*/
+static GelEFunc *ErrorFunction_function = NULL;
 static GelEFunc *RiemannZeta_function = NULL;
 static GelEFunc *GammaFunction_function = NULL;
 static GelEFunc *BesselJ0_function = NULL;
@@ -1675,7 +1676,6 @@ GoldenRatio_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
        return gel_makenum (golden_ratio_cache);
 }
 
-/*  FIXME: I have bad GEL implementation that handles complex values
 static GelETree *
 ErrorFunction_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
 {
@@ -1695,24 +1695,48 @@ ErrorFunction_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
        if G_UNLIKELY ( ! check_argument_number (a, 0, "ErrorFunction"))
                return NULL;
        if G_UNLIKELY (mpw_is_complex (a[0]->val.value)) {
-               gel_errorout (_("%s: Not implemented (yet) for complex values"),
-                             "ErrorFunction");
-               return NULL;
-       }
+               /* FIXME: this is evil! */
+               if G_UNLIKELY (_internal_erf_function == NULL) {
+                       _internal_erf_function = d_makeufunc(d_intern("<internal>exp"),
+                                                            gel_parseexp
+                                                            ("twosqrtpi = 2/sqrt(pi); "
+                                                             "a = 1; "
+                                                             "s = 0; "
+                                                             "n = 0; "
+                                                             "f = 1; "
+                                                             "xx = x; "
+                                                             "xsq = x^2; "
+                                                             "do ( "
+                                                             " t = xx * a * twosqrtpi; "
+                                                             " s = s + t; "
+                                                             " increment n; "
+                                                             " f = f * n; "
+                                                             " a = ((-1)^n) / (((2*n)+1) * f); "
+                                                             " xx = xx * xsq "
+                                                             ") while (|t| > ErrorFunctionTolerance); "
+                                                             "s ",
+                                                             NULL, FALSE, FALSE,
+                                                             NULL, NULL),
+                                                            g_slist_append(NULL,d_intern("x")),1,
+                                                            NULL);
+               }
+               return gel_funccall(ctx,_internal_erf_function,a,1);
 
-       MPW_MPF_REAL (num, a[0]->val.value, tmp);
+               return NULL;
+       } else {
+               MPW_MPF_REAL (num, a[0]->val.value, tmp);
 
-       mpfr_init (ret);
-       mpfr_erf (ret, num, GMP_RNDN);
+               mpfr_init (ret);
+               mpfr_erf (ret, num, GMP_RNDN);
 
-       MPW_MPF_KILL (num, tmp);
+               MPW_MPF_KILL (num, tmp);
 
-       mpw_init (retw);
-       mpw_set_mpf_use (retw, ret);
+               mpw_init (retw);
+               mpw_set_mpf_use (retw, ret);
 
-       return gel_makenum (retw);
+               return gel_makenum (retw);
+       }
 }
-*/
 
 static GelETree *
 RiemannZeta_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
@@ -7061,9 +7085,9 @@ gel_funclib_addall(void)
        FUNC (CatalanConstant, 0, "", "constants",
              N_("Catalan's Constant (0.915...)"));
 
-       /*FUNC (ErrorFunction, 1, "x", "functions", N_("The error function, 2/sqrt(2) * int_0^x e^(-t^2) dt 
(only real values implemented)"));
+       FUNC (ErrorFunction, 1, "x", "functions", N_("The error function, 2/sqrt(2) * int_0^x e^(-t^2) dt"));
        ErrorFunction_function = f;
-       ALIAS (erf, 1, ErrorFunction);*/
+       ALIAS (erf, 1, ErrorFunction);
        FUNC (RiemannZeta, 1, "x", "functions", N_("The Riemann zeta function (only real values 
implemented)"));
        f->no_mod_all_args = 1;
        RiemannZeta_function = f;