[gnumeric] LAMBERTW: New function.

[Morten Welinder <mortenw src gnome org>]

commit 2652e10ea76f1d92904abf9671c743d3547935bc
Author: Morten Welinder <terra gnome org>
Date:   Fri Nov 16 13:47:51 2018 -0500

    LAMBERTW: New function.
    
    This computes W_0 or W_-1.
    
    This probably isn't perfect near the branch point.

 NEWS                          |  1 +
 plugins/fn-math/functions.c   | 29 ++++++++++++++
 plugins/fn-math/plugin.xml.in |  1 +
 src/mathfunc.c                | 88 +++++++++++++++++++++++++++++++++++++++++++
 src/mathfunc.h                |  1 +
 5 files changed, 120 insertions(+)
---
diff --git a/NEWS b/NEWS
index 807776e83..4b4f4289b 100644
--- a/NEWS
+++ b/NEWS
@@ -9,6 +9,7 @@ Morten:
        * Fix conditional style crash.
        * Fix applix locale problem.  [#362]
        * Fix applix encoding and escape problems.  [#363]
+       * New LAMBERTW function.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.43
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index b89626742..3b4d6f0cd 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -1253,6 +1253,32 @@ gnumeric_int (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 
 /***************************************************************************/
 
+static GnmFuncHelp const help_lambertw[] = {
+       { GNM_FUNC_HELP_NAME, F_("LAMBERTW:the Lambert W function")},
+       { GNM_FUNC_HELP_ARG, F_("x:number")},
+       { GNM_FUNC_HELP_ARG, F_("k:branch")},
+       { GNM_FUNC_HELP_NOTE, F_("@k defaults to 0, the principal branch.") },
+       { GNM_FUNC_HELP_NOTE, F_("@k must be either 0 or -1.") },
+       { GNM_FUNC_HELP_EXAMPLES, "=LAMBERTW(3)" },
+       { GNM_FUNC_HELP_EXAMPLES, "=LAMBERTW(3,-1)" },
+       { GNM_FUNC_HELP_SEEALSO, "EXP"},
+       { GNM_FUNC_HELP_END}
+};
+
+static GnmValue *
+gnumeric_lambertw (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
+{
+       gnm_float x = value_get_as_float (argv[0]);
+       gnm_float k = argv[1] ? value_get_as_float (argv[1]) : 0;
+
+       if (k != 0 && k != -1)
+               return value_new_error_NUM (ei->pos);
+
+       return value_new_float (gnm_lambert_w (x, (int)k));
+}
+
+/***************************************************************************/
+
 static GnmFuncHelp const help_log[] = {
        { GNM_FUNC_HELP_NAME, F_("LOG:logarithm of @{x} with base @{base}")},
        { GNM_FUNC_HELP_ARG, F_("x:positive number")},
@@ -3574,6 +3600,9 @@ GnmFuncDescriptor const math_functions[] = {
          gnumeric_int, NULL,
          GNM_FUNC_SIMPLE + GNM_FUNC_AUTO_FIRST,
          GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
+       { "lambertw", "f|f", help_lambertw,
+         gnumeric_lambertw, NULL,
+         GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
        { "lcm", NULL, help_lcm,
          NULL, gnumeric_lcm,
          GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
diff --git a/plugins/fn-math/plugin.xml.in b/plugins/fn-math/plugin.xml.in
index 3778e63e1..7dfcb9e20 100644
--- a/plugins/fn-math/plugin.xml.in
+++ b/plugins/fn-math/plugin.xml.in
@@ -59,6 +59,7 @@
                                <function name="hypot"/>
                                <function name="igamma"/>
                                <function name="int"/>
+                               <function name="lambertw"/>
                                <function name="lcm"/>
                                <function name="linsolve"/>
                                <function name="ln"/>
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 0c4e6735b..f8546def1 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -5019,6 +5019,94 @@ gnm_agm (gnm_float a, gnm_float b)
        return a / scale;
 }
 
+/**
+ * gnm_lambert_w:
+ * @x: a number
+ * @k: branch, either 0 or -1
+ *
+ * Returns: The arithmetic-geometric mean of @a and @b.
+ */
+gnm_float
+gnm_lambert_w (gnm_float x, int k)
+{
+       gnm_float w;
+       static const gnm_float one_over_e = 1 / M_Egnum;
+       static const gnm_float sqrt_one_over_e = gnm_sqrt (1 / M_Egnum);
+       static const gboolean debug = FALSE;
+       gnm_float wmin, wmax;
+       int i, imax = 20;
+
+       if (gnm_isnan (x) || x < -one_over_e)
+               return gnm_nan;
+       else if (x == -one_over_e)
+               return -1;
+
+       if (k == 0) {
+               if (x == gnm_pinf)
+                       return gnm_pinf;
+               if (x < 0)
+                       w = 1.5 * (gnm_sqrt (x + one_over_e) - sqrt_one_over_e);
+               else if (x < 10)
+                       w = gnm_sqrt (x) / 1.7;
+               else {
+                       gnm_float l1 = gnm_log (x);
+                       gnm_float l2 = gnm_log (l1);
+                       w = l1 - l2;
+               }
+               wmin = -1;
+               wmax = gnm_pinf;
+       } else if (k == -1) {
+               if (x >= 0)
+                       return (x == 0) ? gnm_ninf : gnm_nan;
+               if (x < -0.1)
+                       w = -1 - 3 * gnm_sqrt (x + one_over_e);
+               else {
+                       gnm_float l1 = gnm_log (-x);
+                       gnm_float l2 = gnm_log (-l1);
+                       w = l1 - l2;
+               }
+
+               wmin = gnm_ninf;
+               wmax = -1;
+       } else
+               return gnm_nan;
+
+       if (debug) g_printerr ("x = %.20g    w=%.20g\n", x, w);
+       for (i = 0; i < imax; i++) {
+               gnm_float ew = gnm_exp (w);
+               gnm_float wew = w * ew;
+               gnm_float d1 = ew * (w + 1);
+               gnm_float d2 = ew * (w + 2);
+               gnm_float dw;
+               gnm_float wold = w;
+
+               dw = (-2 * ((wew - x) * d1) / (2 * d1 * d1 - (wew - x) * d2));
+               w += dw;
+
+               if (w <= wmin || w >= wmax) {
+                       // We overshot
+                       gnm_float l = (w < wmin ? wmin : wmax);
+                       g_printerr (" (%2d w = %.20g)\n", i, w);
+                       dw = (l - wold) * 15 / 16;
+                       w = wold + dw;
+               }
+
+               if (debug) {
+                       g_printerr ("  %2d w = %.20g\n", i, w);
+                       if (i == imax - 1) {
+                               g_printerr ("  wew = %.20g\n", wew);
+                               g_printerr ("  d1  = %.20g\n", d1);
+                               g_printerr ("  d2  = %.20g\n", d2);
+                       }
+               }
+
+               if (gnm_abs (dw) <= 2 * GNM_EPSILON * gnm_abs (w))
+                       break;
+       }
+
+       return w;
+}
+
 
 /**
  * pow1p:
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 32a68ab8a..997c8ba7d 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -38,6 +38,7 @@ gnm_float gnm_owent (gnm_float h, gnm_float a);
 gnm_float gnm_logcf (gnm_float x, gnm_float i, gnm_float d);
 gnm_float expmx2h (gnm_float x);
 gnm_float gnm_agm(gnm_float a, gnm_float b);
+gnm_float gnm_lambert_w(gnm_float x, int k);
 
 /* "d": density.  */
 /* "p": distribution function.  */