[gnumeric] POCHHAMMER: Extend domain.



commit fe63c11c2cd407d24c2c7d25bcfa4abfa8f2f691
Author: Morten Welinder <terra gnome org>
Date:   Sat Nov 30 16:03:43 2013 -0500

    POCHHAMMER: Extend domain.

 ChangeLog                   |    5 +++
 NEWS                        |    3 ++
 plugins/fn-math/functions.c |    4 +--
 samples/amath.gnumeric      |  Bin 43761 -> 44307 bytes
 src/sf-gamma.c              |   74 ++++++++++++++++++------------------------
 src/sf-gamma.h              |    2 +-
 tools/ChangeLog             |    5 +++
 tools/process-amath.pl      |    2 +-
 8 files changed, 48 insertions(+), 47 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 3b449cf..fb9cbb4 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-11-30  Morten Welinder  <terra gnome org>
+
+       * src/sf-gamma.c (pochhammer): Drop give_log arguments.  Extend to
+       negative values.
+
 2013-11-28  Morten Welinder <terra gnome org>
 
        * configure.ac: Post-release bump.
diff --git a/NEWS b/NEWS
index 67c9f37..57d35ed 100644
--- a/NEWS
+++ b/NEWS
@@ -3,6 +3,9 @@ Gnumeric 1.12.10
 Andreas:
        * Fix handling of dashes on ODF import. [#719509]
 
+Morten:
+       * Extend POCHHAMMER to negative values.
+
 --------------------------------------------------------------------------
 Gnumeric 1.12.9
 
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index 8d31129..6d3b327 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -1257,7 +1257,6 @@ static GnmFuncHelp const help_pochhammer[] = {
         { GNM_FUNC_HELP_NAME, F_("POCHHAMMER:the value of GAMMA(@{x}+ {n})/GAMMA(@{x})")},
         { GNM_FUNC_HELP_ARG, F_("x:number")},
         { GNM_FUNC_HELP_ARG, F_("n:number")},
-       { GNM_FUNC_HELP_ARG, F_("give_log:if true, log of the result will be returned instead") },
         { GNM_FUNC_HELP_EXAMPLES, "=POCHHAMMER(1,5)" },
         { GNM_FUNC_HELP_EXAMPLES, "=POCHHAMMER(6,0.5)" },
         { GNM_FUNC_HELP_SEEALSO, "GAMMA"},
@@ -1269,9 +1268,8 @@ gnumeric_pochhammer (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 {
        gnm_float x = value_get_as_float (argv[0]);
        gnm_float n = value_get_as_float (argv[1]);
-       gboolean give_log = argv[2] ? value_get_as_checked_bool (argv[2]) : FALSE;
 
-       return value_new_float (pochhammer (x, n, give_log));
+       return value_new_float (pochhammer (x, n));
 }
 
 /***************************************************************************/
diff --git a/samples/amath.gnumeric b/samples/amath.gnumeric
index e00f855..c48fe93 100644
Binary files a/samples/amath.gnumeric and b/samples/amath.gnumeric differ
diff --git a/src/sf-gamma.c b/src/sf-gamma.c
index 1c759d5..f606083 100644
--- a/src/sf-gamma.c
+++ b/src/sf-gamma.c
@@ -657,18 +657,18 @@ static gnm_float
 pochhammer_naive (gnm_float x, int n)
 {
        void *state = gnm_quad_start ();
-       GnmQuad qp, qx, qone;
+       GnmQuad qp, qx;
        gnm_float r;
 
-       gnm_quad_init (&qone, 1);
-       qp = qone;
+       qp = gnm_quad_one;
        gnm_quad_init (&qx, x);
        while (n-- > 0) {
                gnm_quad_mul (&qp, &qp, &qx);
-               gnm_quad_add (&qx, &qx, &qone);
+               gnm_quad_add (&qx, &qx, &gnm_quad_one);
        }
        r = gnm_quad_value (&qp);
        gnm_quad_end (state);
+
        return r;
 }
 
@@ -681,7 +681,7 @@ pochhammer_naive (gnm_float x, int n)
  */
 
 gnm_float
-pochhammer (gnm_float x, gnm_float n, gboolean give_log)
+pochhammer (gnm_float x, gnm_float n)
 {
        gnm_float rn, rx, lr;
        GnmQuad m1, m2;
@@ -690,12 +690,8 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
        if (gnm_isnan (x) || gnm_isnan (n))
                return gnm_nan;
 
-       /* This isn't a fundamental restriction, but one we impose.  */
-       if (x <= 0 || x + n <= 0)
-               return gnm_nan;
-
        if (n == 0)
-               return give_log ? 0 : 1;
+               return 1;
 
        rx = gnm_floor (x + 0.5);
        rn = gnm_floor (n + 0.5);
@@ -705,10 +701,8 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
         * We don't want to use this if x is also an integer
         * (but we might do so below if x is insanely large).
         */
-       if (n == rn && x != rx && n >= 0 && n < 40) {
-               gnm_float r = pochhammer_naive (x, (int)n);
-               return give_log ? gnm_log (r) : r;
-       }
+       if (n == rn && x != rx && n >= 0 && n < 40)
+               return pochhammer_naive (x, (int)n);
 
        if (!qfactf (x + n - 1, &m1, &e1) &&
            !qfactf (x - 1, &m2, &e2)) {
@@ -721,9 +715,19 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
                r = gnm_quad_value (&qr);
                gnm_quad_end (state);
 
-               return give_log
-                       ? gnm_log (r) + M_LN2gnum * de
-                       : gnm_ldexp (r, de);
+               return gnm_ldexp (r, de);
+       }
+
+       if (x == rx && x <= 0) {
+               if (n != rn)
+                       return 0;
+               if (x == 0)
+                       return (n > 0)
+                               ? 0
+                               : ((gnm_fmod (-n, 2) == 0 ? +1 : -1) /
+                                  gnm_fact (-n));
+               if (n > -x)
+                       return gnm_nan;
        }
 
        /*
@@ -731,23 +735,14 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
         * insanely big, possibly both.
         */
 
-       if (gnm_abs (x) < 1) {
-               /* n is big. */
-               if (give_log)
-                       goto via_log;
-               else
-                       return gnm_pinf;
-       }
+       if (gnm_abs (x) < 1)
+               return gnm_pinf;
 
-       if (n < 0) {
-               gnm_float r = pochhammer (x + n, -n, give_log);
-               return give_log ? -r : 1 / r;
-       }
+       if (n < 0)
+               return 1 / pochhammer (x + n, -n);
 
-       if (n == rn && n >= 0 && n < 100) {
-               gnm_float r = pochhammer_naive (x, (int)n);
-               return give_log ? gnm_log (r) : r;
-       }
+       if (n == rn && n >= 0 && n < 100)
+               return pochhammer_naive (x, (int)n);
 
        if (gnm_abs (n) < 1) {
                /* x is big.  */
@@ -757,17 +752,16 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
                pochhammer_small_n (x, n, &qr);
                r = gnm_quad_value (&qr);
                gnm_quad_end (state);
-               return give_log ? gnm_log (r) : r;
+               return r;
        }
 
        /* Panic mode.  */
        g_printerr ("x=%.20g  n=%.20g\n", x, n);
-via_log:
        lr = ((x - 0.5) * gnm_log1p (n / x) +
              n * gnm_log (x + n) -
              n +
              (lgammacor (x + n) - lgammacor (x)));
-       return give_log ? lr : gnm_exp (lr);
+       return gnm_exp (lr);
 }
 
 /* ------------------------------------------------------------------------- */
@@ -981,7 +975,7 @@ combin (gnm_float n, gnm_float k)
                return c;
        }
 
-       if (k < 30) {
+       if (k < 100) {
                void *state = gnm_quad_start ();
                GnmQuad p, a, b;
                gnm_float c;
@@ -1001,11 +995,7 @@ combin (gnm_float n, gnm_float k)
                return c;
        }
 
-       {
-               gnm_float lp = pochhammer (n - k + 1, k, TRUE) -
-                       gnm_lgamma (k + 1);
-               return gnm_floor (0.5 + gnm_exp (lp));
-       }
+       return pochhammer (n - k + 1, k) / gnm_fact (k);
 }
 
 gnm_float
@@ -1014,7 +1004,7 @@ permut (gnm_float n, gnm_float k)
        if (k < 0 || k > n || n != gnm_floor (n) || k != gnm_floor (k))
                return gnm_nan;
 
-       return pochhammer (n - k + 1, k, FALSE);
+       return pochhammer (n - k + 1, k);
 }
 
 /* ------------------------------------------------------------------------- */
diff --git a/src/sf-gamma.h b/src/sf-gamma.h
index 21bc872..b3e2023 100644
--- a/src/sf-gamma.h
+++ b/src/sf-gamma.h
@@ -14,7 +14,7 @@ gnm_float gnm_lbeta (gnm_float a, gnm_float b);
 gnm_float gnm_beta (gnm_float a, gnm_float b);
 gnm_float gnm_lbeta3 (gnm_float a, gnm_float b, int *sign);
 
-gnm_float pochhammer (gnm_float x, gnm_float n, gboolean give_log);
+gnm_float pochhammer (gnm_float x, gnm_float n);
 gnm_float combin (gnm_float n, gnm_float k);
 gnm_float permut (gnm_float n, gnm_float k);
 
diff --git a/tools/ChangeLog b/tools/ChangeLog
index 57265a8..93a082d 100644
--- a/tools/ChangeLog
+++ b/tools/ChangeLog
@@ -1,3 +1,8 @@
+2013-11-30  Morten Welinder  <terra gnome org>
+
+       * process-amath.pl (def_expr_handler): Allow negative arguments
+       for POCHHAMMER.
+
 2013-11-28  Morten Welinder <terra gnome org>
 
        * Release 1.12.9
diff --git a/tools/process-amath.pl b/tools/process-amath.pl
index a65a325..b7ccd63 100755
--- a/tools/process-amath.pl
+++ b/tools/process-amath.pl
@@ -140,7 +140,6 @@ my %expr_handlers =
     ('beta' => \&non_negative_handler,
      'gammaln' => \&non_negative_handler,
      'factdouble' => \&non_negative_handler,
-     'pochhammer' => \&positive_handler, # We shouldn't need this
      'combin' => \&non_negative_handler,
      'r.dcauchy' => sub { &reorder_handler ("3,1,2", @_); },
      'r.pcauchy' => sub { &reorder_handler ("3,1,2", @_); },
@@ -226,6 +225,7 @@ sub output_test {
     my ($gfunc,$expr,$res) = @_;
 
     my $gfunc0 = ($gfunc eq $last_func) ? '' : $gfunc;
+    $res = "=$res" if $res =~ m{/};
 
     my $N = $test_row++;
     print 
"\"$gfunc0\",\"=$expr\",\"$res\",\"=IF(B$N=C$N,\"\"\"\",IF(C$N=0,-LOG10(ABS(B$N)),-LOG10(ABS((B$N-C$N)/C$N))))\"\n";


[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]