[gnumeric] INFACT, IMGAMMA: New functions.

[Morten Welinder <mortenw src gnome org>]

commit f0a7d3ae3d98b6ddb0e519649f8e9feffe6d9a87
Author: Morten Welinder <terra gnome org>
Date:   Mon Dec 2 23:06:55 2013 -0500

    INFACT, IMGAMMA: New functions.
    
    Because we can.  And because the complex gamma functions plays
    a part in the complex incomplete gamma function which, in turn,
    plays a part in bessel function series.

 NEWS                             |    1 +
 plugins/fn-complex/functions.c   |   58 +++++++++++++++++++++++++
 plugins/fn-complex/plugin.xml.in |    2 +
 src/complex.c                    |   88 ++++++++++++++++++++++++++++++++++++++
 src/complex.h                    |    3 +
 5 files changed, 152 insertions(+), 0 deletions(-)
---
diff --git a/NEWS b/NEWS
index 0b9eff5..00ff9c8 100644
--- a/NEWS
+++ b/NEWS
@@ -7,6 +7,7 @@ Morten:
        * Extend POCHHAMMER to negative values.
        * Improve accuracy of BETA and BETALN.
        * Improve accuracy of BESSELJ and BESSELY.
+       * New functions: IMGAMMA, IMFACT.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.9
diff --git a/plugins/fn-complex/functions.c b/plugins/fn-complex/functions.c
index 8da2f45..93373a1 100644
--- a/plugins/fn-complex/functions.c
+++ b/plugins/fn-complex/functions.c
@@ -1094,6 +1094,56 @@ gnumeric_imsqrt (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
 
 /***************************************************************************/
 
+static GnmFuncHelp const help_imfact[] = {
+        { GNM_FUNC_HELP_NAME, F_("IMFACT:the factorial of the complex number @{z}") },
+        { GNM_FUNC_HELP_ARG, F_("z:a complex number") },
+        { GNM_FUNC_HELP_NOTE, F_("If @{z} is not a valid complex number, #VALUE! is returned.") },
+        { GNM_FUNC_HELP_EXAMPLES, "=IMFACT(\"1+j\")" },
+        { GNM_FUNC_HELP_SEEALSO, "IMGAMMA" },
+        { GNM_FUNC_HELP_END}
+};
+
+
+static GnmValue *
+gnumeric_imfact (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
+{
+       complex_t c, res;
+       char imunit;
+
+       if (value_get_as_complex (argv[0], &c, &imunit))
+               return value_new_error_NUM (ei->pos);
+
+       complex_fact (&res, &c);
+       return value_new_complex (&res, imunit);
+}
+
+/***************************************************************************/
+
+static GnmFuncHelp const help_imgamma[] = {
+        { GNM_FUNC_HELP_NAME, F_("IMGAMMA:the gamma function of the complex number @{z}") },
+        { GNM_FUNC_HELP_ARG, F_("z:a complex number") },
+        { GNM_FUNC_HELP_NOTE, F_("If @{z} is not a valid complex number, #VALUE! is returned.") },
+        { GNM_FUNC_HELP_EXAMPLES, "=IMGAMMA(\"1+j\")" },
+        { GNM_FUNC_HELP_SEEALSO, "IMGAMMA" },
+        { GNM_FUNC_HELP_END}
+};
+
+
+static GnmValue *
+gnumeric_imgamma (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
+{
+       complex_t c, res;
+       char imunit;
+
+       if (value_get_as_complex (argv[0], &c, &imunit))
+               return value_new_error_NUM (ei->pos);
+
+       complex_gamma (&res, &c);
+       return value_new_complex (&res, imunit);
+}
+
+/***************************************************************************/
+
 static GnmFuncHelp const help_imsub[] = {
         { GNM_FUNC_HELP_NAME, F_("IMSUB:the difference of two complex numbers") },
         { GNM_FUNC_HELP_ARG, F_("z1:a complex number") },
@@ -1355,5 +1405,13 @@ GnmFuncDescriptor const complex_functions[] = {
         { "imarccoth",    "S",    help_imarccoth,
          gnumeric_imarccoth, NULL, NULL, NULL,
          GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
+
+        { "imfact",       "S",    help_imfact,
+         gnumeric_imfact, NULL, NULL, NULL,
+         GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
+        { "imgamma",      "S",    help_imgamma,
+         gnumeric_imgamma, NULL, NULL, NULL,
+         GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
+
         {NULL}
 };
diff --git a/plugins/fn-complex/plugin.xml.in b/plugins/fn-complex/plugin.xml.in
index 094b034..7d7b7e2 100644
--- a/plugins/fn-complex/plugin.xml.in
+++ b/plugins/fn-complex/plugin.xml.in
@@ -16,6 +16,8 @@
                                 <function name="imaginary"/>
                                 <function name="imargument"/>
                                 <function name="imconjugate"/>
+                                <function name="imfact"/>
+                                <function name="imgamma"/>
                                 <function name="iminv"/>
                                 <function name="imneg"/>
                                 <function name="imcos"/>
diff --git a/src/complex.c b/src/complex.c
index 739324e..fd707ce 100644
--- a/src/complex.c
+++ b/src/complex.c
@@ -14,6 +14,7 @@
 #include <stdlib.h>
 #include <errno.h>
 #include <mathfunc.h>
+#include <sf-gamma.h>
 
 /* ------------------------------------------------------------------------- */
 
@@ -177,3 +178,90 @@ complex_invalid_p (complex_t const *src)
 }
 
 /* ------------------------------------------------------------------------- */
+
+void
+complex_gamma (complex_t *dst, complex_t const *src)
+{
+       if (complex_real_p (src)) {
+               complex_init (dst, gnm_gamma (src->re), 0);
+       } else {
+               complex_t z = *src, f;
+
+               complex_fact (&f, src);
+               complex_div (dst, &f, &z);
+       }
+}
+
+/* ------------------------------------------------------------------------- */
+
+void
+complex_fact (complex_t *dst, complex_t const *src)
+{
+       if (complex_real_p (src)) {
+               complex_init (dst, gnm_fact (src->re), 0);
+       } else if (src->re < 0.5) {
+               /* Fact(z) = -pi / (sin(pi*z) * Gamma(-z)) */
+               complex_t a, b, mz;
+
+               complex_init (&mz, -src->re, -src->im);
+               complex_gamma (&a, &mz);
+
+               complex_init (&b,
+                             M_PIgnum * gnm_fmod (src->re, 2),
+                             M_PIgnum * src->im);
+               complex_sin (&b, &b);
+
+               complex_mul (&a, &a, &b);
+
+               complex_init (&b, -M_PIgnum, 0);
+
+               complex_div (dst, &b, &a);
+       } else {
+               static const gnm_float c[] = {
+                       GNM_const (0.99999999999999709182),
+                       GNM_const (57.156235665862923517),
+                       GNM_const (-59.597960355475491248),
+                       GNM_const (14.136097974741747174),
+                       GNM_const (-0.49191381609762019978),
+                       GNM_const (0.33994649984811888699e-4),
+                       GNM_const (0.46523628927048575665e-4),
+                       GNM_const (-0.98374475304879564677e-4),
+                       GNM_const (0.15808870322491248884e-3),
+                       GNM_const (-0.21026444172410488319e-3),
+                       GNM_const (0.21743961811521264320e-3),
+                       GNM_const (-0.16431810653676389022e-3),
+                       GNM_const (0.84418223983852743293e-4),
+                       GNM_const (-0.26190838401581408670e-4),
+                       GNM_const (0.36899182659531622704e-5)
+               };
+               const gnm_float g = GNM_const(607.0) / 128;
+               const gnm_const sqrt2pi =
+                       gnm_const (2.506628274631000502415765284811045253006986740609938316629923);
+               complex_t zph, zpgh, s, f;
+               int i;
+
+               complex_init (&zph, src->re + 0.5, src->im);
+               complex_init (&zpgh, src->re + g + 0.5, src->im);
+               complex_init (&s, 0, 0);
+
+               for (i = G_N_ELEMENTS(c) - 1; i >= 1; i--) {
+                       complex_t d, q;
+                       complex_init (&d, src->re + i, src->im);
+                       complex_init (&q, c[i], 0);
+                       complex_div (&q, &q, &d);
+                       complex_add (&s, &s, &q);
+               }
+               s.re += c[0];
+
+               complex_init (&f, sqrt2pi, 0);
+               complex_mul (&s, &s, &f);
+               complex_exp (&f, &zpgh);
+               complex_div (&s, &s, &f);
+               complex_pow (&f, &zpgh, &zph);
+               complex_mul (&s, &s, &f);
+
+               *dst = s;
+       }
+}
+
+/* ------------------------------------------------------------------------- */
diff --git a/src/complex.h b/src/complex.h
index 7efafc4..f06dea4 100644
--- a/src/complex.h
+++ b/src/complex.h
@@ -66,6 +66,9 @@ int complex_from_string (complex_t *dst, char const *src, char *imunit);
 
 int complex_invalid_p (complex_t const *src);
 
+void complex_gamma (complex_t *dst, complex_t const *src);
+void complex_fact (complex_t *dst, complex_t const *src);
+
 /* ------------------------------------------------------------------------- */
 
 G_END_DECLS