/* GNUPLOT - complexfun.c */

/*
 * FILE CONTENTS
 *
 * -	BSD 2-clause license covering material in this file
 *
 * -	f_Sign
 *		Sign(z) = z/|z| for complex z
 *
 * -	f_LambertW lambert_initial LambertW
 *		Lambert W function for complex numbers
 *
 * -	f_lnGamma lnGamma
 *		log Gamma for complex argument z
 *	- 14 term Lanczos approximation
 *
 * -	f_Igamma Igamma Igamma_GL Igamma_negative_z Igamma_Poincare
 *		lower incomplete gamma function P(a, z)
 * 	- adapted from previous gnuplot real-valued function igamma(a,x)
 *	- REFERENCE ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3
 *	  B. L. Shea "Chi-Squared and Incomplete Gamma Integral"
 *	- Poincar expansion for large z (not currently used) based on
 *	  Gil et al (2016) ACM TOMS 43:3 Article 26
 *	- Coefficients for Gauss-Legendre quadrature used for a > 100
 *	  Press et al, Numerical Recipes (3rd Ed.) Section 6.2
 *
 * -	Riemann_zeta f_zeta
 *		Riemann zeta function zeta(s) for general complex argument
 *	polynomial series using	algorithm 3 from Borwein [2000] MR1777614
 */

/*[
 * Copyright Ethan A Merritt 2019-2021
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * Redistributions of source code must retain the above copyright notice, this
 * list of conditions and the following disclaimer.  Redistributions in binary
 * form must reproduce the above copyright notice, this list of conditions and
 * the following disclaimer in the documentation and/or other materials provided
 * with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
]*/

#include "syscfg.h"
#include "gp_types.h"
#include "eval.h"
#include "stdfn.h"
#include "util.h"	/* for int_error() */

#ifdef HAVE_COMPLEX_FUNCS

#include 	/* C99 _Complex */
#include "complexfun.h"

#ifdef HAVE_FENV_H
#include 
#endif

/*
 * Various complex functions like cexp may set errno on underflow
 * We would prefer to return 0.0 rather than NaN
 */
#ifdef HAVE_FENV_H
#define initialize_underflow( who ) \
    if (errno) \
	int_error(NO_CARET, "%s: error present on entry (errno %d %s)", who, errno, strerror(errno)); \
    else feclearexcept(FE_ALL_EXCEPT);
#else
#define initialize_underflow( who ) \
    if (errno) \
	int_error(NO_CARET, "%s: error present on entry (errno %d %s)", who, errno, strerror(errno));
#endif

#ifdef HAVE_FENV_H
#define handle_underflow( who, var ) \
    if (errno) { \
	if (fetestexcept(FE_UNDERFLOW)) { \
	    var = 0.0; \
	    errno = 0; \
	} else { \
	    fprintf(stderr,"%s: errno = %d\n", who, errno); \
	} \
    }
#else
#define handle_underflow( who, var ) int_error(NO_CARET, "%s: errno = %d", who, errno);
#endif

/* internal prototypes */
static complex double lnGamma( complex double z );
static complex double Igamma( complex double a, complex double z );
static double complex Igamma_GL( double complex a, double complex z );
static double complex Igamma_negative_z( double a, double complex z );

#undef IGAMMA_POINCARE
#ifdef IGAMMA_POINCARE
static double complex Igamma_Poincare( double a, double complex z );
#endif

/* wrapper for Igamma so that when it replaces igamma
 * there is still something for old callers who want to call
 * it with real arguments rather than complex.
 */
double
igamma( double a, double z )
{
    return creal( Igamma( (complex double)a, (complex double)z ) );
}

/*
 * Complex Sign function
 * Sign(z) = z/|z| for z non-zero
 */
void
f_Sign(union argument *arg)
{
    struct value result;
    struct value a;
    complex double z;

    pop(&a);	/* Complex argument z */
    if (a.type == INTGR) {
	push(Gcomplex(&result, sgn(a.v.int_val), 0.0));
    } else if (a.type == CMPLX) {
	z = a.v.cmplx_val.real + I*a.v.cmplx_val.imag;
	if (z != 0.0)
	    z = z/cabs(z);
	push(Gcomplex(&result, creal(z), cimag(z)));
    } else
	int_error(NO_CARET, "z must be numeric");
}


/*
 * Lambert W function for complex numbers
 *
 * W(z) is a multi-valued function with the defining property
 * 
 *     z = W(z) exp(W(z))   for complex z
 *
 * LambertW( z, k ) is the kth branch of W
 *
 * This implementation guided by C++ code by Istvn Mez 
 * See also
 *   R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth,
 *   On the Lambert W function, Adv. Comput. Math. 5 (1996), no. 4, 329359.
 *   DOI 10.1007/BF02124750.
 */

/* Internal Prototypes */
complex double lambert_initial(complex double z, int k);
complex double LambertW(complex double z, int k);

void
f_LambertW(union argument *arg)
{
    struct value result;
    struct value a;

    struct cmplx z;	/* gnuplot complex parameter z */
    int k;		/* gnuplot integer parameter k */

    complex double w;	/* C99 _Complex representation */

    pop(&a);		/* Integer argument k */
    if (a.type != INTGR)
	int_error(NO_CARET, "k must be integer");
    k = a.v.int_val;
    pop(&a);		/* Complex argument z */
    if (a.type != CMPLX)
	int_error(NO_CARET, "z must be real or complex");
    z = a.v.cmplx_val;

    w = z.real + I*z.imag;
    w = LambertW( w, k );

    push(Gcomplex(&result, creal(w), cimag(w)));
}

/*
 * First and second derivatives for z * e^z
 * dzexpz( z )  = first derivative of ze^z = e^z + ze^z
 * ddzexpz( z ) = second derivative of ze^z = e^z + e^z + ze^z
 */
#define dzexpz(z)  (cexp(z) + z * cexp(z))
#define ddzexpz(z) (2. * cexp(z) + z * cexp(z))

/*
 * The hard part is choosing a starting point
 * since Halley's method does not have a large radius of convergence
 * EAM: The domain windows in which special case starting points are used
 *      as found in the Mez code produced glitches in my tests.
 *      I adjusted them empirically but I have no justification for
 *      the specific window sizes or thresholds.
 */
complex double
lambert_initial( complex double z, int k )
{
    complex double e = 2.71828182845904523536;
    complex double branch = 2 * M_PI * I * k;
    complex double ip;
    double close;

    double case1_window = 1.2;	/* see note above, was 1.0 */
    double case2_window = 0.9;	/* see note above, was 1.0 */
    double case3_window = 0.5;	/* see note above, was 0.5 */

    /* Initial term of Eq (4.20) from Corless et al */
    ip = clog(z) + branch - clog(clog(z) + branch);

    /* Close to a branch point use (4.22) from Corless et al */
    close = cabs(z - (-1/e));
    if (close  0 || close < case2_window)
		ip = -1. + p - (1./3.) * p*p + (11./72.) * p*p*p;
	}
#if (0)
	/* This treatment empirically causes more glitches than it removes */
	if (k == 1 && cimag(z) < 0.0) {
	    if (creal(z) > 0 && close < case2_window) {
		ip = -1. - p - (1./3.) * p*p - (11./72.) * p*p*p;
		if (cimag(z) > -0.1)
		    ip += (-43./540.) * p*p*p*p;
	    }
	}
#endif
	if (k == -1 && cimag(z) > 0.) {
	    if (close < case2_window)
		ip = -1. - p - (1./3.) * p*p - (11./72.) * p*p*p;
	}
    }

    /* Pad approximant for W(0,a) */
    if (k == 0 && cabs(z - 0.5)