[gnumeric] ILOG: new function.

[Morten Welinder <mortenw src gnome org>]

commit 0bada92137f14ebb4ece553a595c3f331613db26
Author: Morten Welinder <terra gnome org>
Date:   Tue Oct 12 20:41:39 2021 -0400

    ILOG: new function.
    
    This is the integer log function.  Implemented for bases 2 and 10 so far.

 NEWS                          |  1 +
 plugins/fn-math/functions.c   | 34 ++++++++++++++++++++++
 plugins/fn-math/plugin.xml.in |  1 +
 src/mathfunc.c                | 67 +++++++++++++++++++++++++++++++++++++++++++
 src/mathfunc.h                |  1 +
 5 files changed, 104 insertions(+)
---
diff --git a/NEWS b/NEWS
index 8ab05ab52..774038ad4 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,7 @@ Morten:
        * Fix --with-long-double configuration.
        * Fix applix import/locale problem.
        * Make sure CEIL doesn't return -0.  [#603]
+       * Add ILOG function.  [#591]
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.50
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index 1726d6846..20f6d2ce8 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -1339,6 +1339,37 @@ gnumeric_log (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 
 /***************************************************************************/
 
+static GnmFuncHelp const help_ilog[] = {
+       { GNM_FUNC_HELP_NAME, F_("ILOG:integer logarithm of @{x} with base @{base}")},
+       { GNM_FUNC_HELP_ARG, F_("x:positive number")},
+       { GNM_FUNC_HELP_ARG, F_("base:base of the logarithm, defaults to 10")},
+       { GNM_FUNC_HELP_NOTE, F_("@{base} must be positive and not equal to 1.") },
+       { GNM_FUNC_HELP_NOTE, F_("If @{x} \xe2\x89\xa4 0, LOG returns #NUM! error.") },
+       { GNM_FUNC_HELP_NOTE, F_("This function returns the logarithm of @{x} using @{base} rounded down to 
nearest integer.  Unlike FLOOR(LOG(@{x},@{base})), this function is not subject error of representation of 
the intermediate result.") },
+       { GNM_FUNC_HELP_NOTE, F_("This function is not implemented for all possible value.  #VALUE! will be 
returned for such arguments.") },
+       { GNM_FUNC_HELP_EXAMPLES, "=ILOG(2^32,2)" },
+       { GNM_FUNC_HELP_EXAMPLES, "=LOG(10^15)" },
+       { GNM_FUNC_HELP_SEEALSO, "LOG"},
+       { GNM_FUNC_HELP_END}
+};
+
+static GnmValue *
+gnumeric_ilog (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
+{
+       gnm_float x = value_get_as_float (argv[0]);
+       gnm_float base = argv[1] ? value_get_as_float (argv[1]) : 10;
+
+       if (base == 1. || base <= 0.)
+               return value_new_error_NUM (ei->pos);
+
+       if (x <= 0.0)
+               return value_new_error_NUM (ei->pos);
+
+       return value_new_float (gnm_ilog (x, base));
+}
+
+/***************************************************************************/
+
 static GnmFuncHelp const help_ln[] = {
        { GNM_FUNC_HELP_NAME, F_("LN:the natural logarithm of @{x}")},
        { GNM_FUNC_HELP_ARG, F_("x:positive number")},
@@ -3619,6 +3650,9 @@ GnmFuncDescriptor const math_functions[] = {
        { "igamma",    "ff|bbb",  help_igamma,
          gnumeric_igamma, NULL,
          GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
+       { "ilog",     "f|f",  help_ilog,
+         gnumeric_ilog, NULL,
+         GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
        { "int",     "f",     help_int,
          gnumeric_int, NULL,
          GNM_FUNC_SIMPLE + GNM_FUNC_AUTO_FIRST,
diff --git a/plugins/fn-math/plugin.xml.in b/plugins/fn-math/plugin.xml.in
index b5bd4bd8b..70b221097 100644
--- a/plugins/fn-math/plugin.xml.in
+++ b/plugins/fn-math/plugin.xml.in
@@ -59,6 +59,7 @@
                                <function name="gd"/>
                                <function name="hypot"/>
                                <function name="igamma"/>
+                               <function name="ilog"/>
                                <function name="int"/>
                                <function name="lambertw"/>
                                <function name="lcm"/>
diff --git a/src/mathfunc.c b/src/mathfunc.c
index f51b72261..c2222ec14 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -6041,3 +6041,70 @@ gnm_owent (gnm_float h, gnm_float a)
 }
 
 /* ------------------------------------------------------------------------- */
+
+gnm_float
+gnm_ilog (gnm_float x, gnm_float b)
+{
+       if (gnm_isnan (x) || x < 0 ||
+           gnm_isnan (b) || b == 1 || b <= 0 || b == gnm_pinf)
+               return gnm_nan;
+
+       if (x == 0)
+               return b < 1 ? gnm_pinf : gnm_ninf;
+
+       if (x == gnm_pinf)
+               return b < 1 ? gnm_ninf : gnm_pinf;
+
+       if (b == 2) {
+               int e;
+               gnm_float m = gnm_frexp (x, &e);
+               (void)m;
+               return e - 1;
+       }
+
+       if (b == 10) {
+               if (x >= 1 && x <= 1e22) {
+                       // This code relies on 10^i being exact
+                       int l10 = (int)(gnm_log10 (x));
+                       if (gnm_pow10 (l10) > x)
+                               l10++;
+                       return l10;
+               } else {
+                       void *state = gnm_quad_start ();
+                       GnmQuad qx, q10, qlog10, qfudge;
+
+                       gnm_quad_init (&q10, 10);
+                       gnm_quad_log (&qlog10, &q10);
+
+                       gnm_quad_init (&qx, x);
+                       gnm_quad_log (&qx, &qx);
+                       gnm_quad_div (&qx, &qx, &qlog10);
+
+                       // This looks bad, but actually isn't because the
+                       // true logarithm cannot be too close to an integer
+                       // while still being less.
+                       //
+                       // Let eps = 1ulp for 10^i (roughly 10^i * GNM_EPSILON)
+                       //
+                       // log10(10^i-eps) =
+                       //    i + log10(1-eps/10^i) =
+                       //    1 - eps/10^i/log(10) + O((eps/10^i)^2)
+                       // As long as we add something smaller than
+                       // eps/10^i/log(10) (roughly GNM_EPSILON/3), we
+                       // should be fine.  *should*
+                       // Verification needed.
+                       gnm_quad_init (&qfudge, GNM_EPSILON / 32);
+                       gnm_quad_add (&qx, &qx, &qfudge);
+                       gnm_quad_floor (&qx, &qx);
+
+                       gnm_quad_end (state);
+
+                       return gnm_quad_value (&qx);
+               }
+       }
+
+       // Not implemented.
+       return gnm_nan;
+}
+
+/* ------------------------------------------------------------------------- */
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 997c8ba7d..cbdf32cf2 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -39,6 +39,7 @@ 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);
+gnm_float gnm_ilog (gnm_float x, gnm_float b);
 
 /* "d": density.  */
 /* "p": distribution function.  */