Commit 2652e10e authored by Morten Welinder's avatar Morten Welinder

LAMBERTW: New function.

This computes W_0 or W_-1.

This probably isn't perfect near the branch point.
parent ea4c921a
......@@ -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
......
......@@ -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 },
......
......@@ -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"/>
......
......@@ -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:
......
......@@ -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. */
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment