/* Copyright (C) 1996 by Thomas Glen Smith.  All Rights Reserved. */
/* gammax APL2 V1.0.0 **************************************************
* Called by logamma and factorlx.                                      *
* Factorial (gamma) for complex numbers.                               *
* This is a translation from FORTRAN to C of the complex gamma         *
* function taken from the "Collected Algorithms from CACM," algorithm  *
* 421, by Hirondo Kuki.                                                *
***********************************************************************/
#define INCLUDES MATH+TRIGKEYS
#include "includes.h"
void gammax(carg,cans,lf0)
double *carg,*cans;
int lf0; /* 0 for logamma, 1 for gamma. */
{
	Factorlp;
	extern int aplerr;
	extern double fuzz;
	static double
		coef[7]={
			/* 1 */ +.641025641025641026e-2,
			/* 2 */ -.191752691752691753e-2,
			/* 3 */ +.841750841750841751e-3,
			/* 4 */ -.595238095238095238e-3,
			/* 5 */ +.793650793650793651e-3,
			/* 6 */ -.277777777777777778e-2,
			/* 7 */ +.833333333333333333e-1},
		f0=840.07385296052619e0, f1=20.001230821894200e0,
		g0=1680.1477059210524e0, g1=180.01477047052042e0,
		pi=3.14159265358979324e0, dpi=6.28318530717958648e0,
		hl2p=.918938533204672742e0, al2p=1.83787706640934548e0,
/*
	Constants eps and omega are machine dependent.  Eps is the basic
	round-off unit. For s/360 machines, it is chosen to be 16**-13.  For
	binary machine of n-bit accuracy set it to 2**(-n+1), and initialize
	de0 as 5.0 rather than as 2.0.  Omega is the largest number
	representable by the float point representation of the machine.  For
	S/360 machines, it is slightly less than 16**63.
*/
		eps=2.2e-16, omega=7.23700538e75;
	double a, al1, al2, b, delta, de0, de1, dn, h, h1, h2,
		t1, t2,
		u, uuu1, uuu2, uu1, uu2, u1, u2, vv1, vv2, v1, v2, w1, w2,
		y1, zz1, z1, z2;
	long int j,k,lf1,lf2,lf3,n;

	*cans = *(cans+1) = 0e0;
	z1 = *(carg);
	z2 = *(carg+1);
	if (((z2 >= 0e0) ? z2 : -z2)  < fuzz) {
		factorlp(carg,cans); /* Imaginary part is essentially 0. */
		return;
	}
	delta = fuzz;
	de0 = 2e0;
	de1 = 0e0;
/* Force sign of imaginary part of arg to non-negative. */
	lf1 = 0;
	if (!(z2 >= 0e0)) {
		lf1 = 1;
		z2 = -z2;
	}
/* 20 */
	lf2 = 0;
	if (!(z1 >= 0e0)) { /*IF(Z1.GE.0.0) GO TO 100*/
		/* case when real part of arg is negative */
		lf2 = 1;
		lf1--;
		t1 = al2p - pi * z2;
		t2 = pi*(.5 - z1);
		u = -dpi*z2;
		if (!(u >= -.1054e0)) { /* IF(U.GE.-0.1054D0) GO TO 40 */
			a = 0e0;
			/* If e**u .lt. 10**(-17), ignore it to save time and to avoid
			irrelevant underflow. */
			if (!(u <= -39.15))
				a = exp(u);
/* 30 */
			h1 = 1.0 - a;
		}
		else {
/* 40 */
			u2 = u * u;
			a  = -u * (f1 * u2 + f0);
			h1 = (a + a) / ((u2 + g1) * u2 + g0 + a);
			a  = 1.0 - h1;
		}
/* 50 */
		b  = z1 - (long)(z1 - .5); /* B=Z1-DINT(Z1-.5) */
		h2 = a * sin(dpi * b);
		b  = sin(pi * b);
		h1 = h1 + (b+b)*b*a;
		h  = ((h2 > 0.0) ? h2 : -h2) + h1 - dpi*a*delta;
		if (h <= 0e0) { /* IF (H.LE.0.0) GO TO 500 */
			aplerr = 130; return; }
		de0 += ((t1 > 0e0) ? t1 : -t1) + t2;
		de1 = pi + dpi*a/h;
		z1 = 1.0 - z1;
	}
/*
	Case when neither real part nor imaginary part of arg is
	negative.  Define threshold curve to be the broken lines
	connecting points 10f0*i, 10f4.142*i, 0.1f14.042*i, and
	0.1fomega*i.
*/
/* 100 */
	lf3 = 0;
	y1 = z1 - .5;
	w1 = 0e0;
	w2 = 0e0;
	k = 0;
	b = 14.142 - z2;
	if (b > 10.0) b = 10.0;
	else if (b < .1) b = .1; /* MAX(.1,MIN(10,14.142-Z2)) */
	b -= z1;
	for (;;) {
		if (b <= 0e0) break; /* IF(B.LE.0.0) GO TO 200 */
		/* Case when real part of arg is between 0 and threshold. */
		lf3 = 1;
		zz1 = z1;
		n = b + 1.0;
		dn = n;
		z1 += dn;
		a = z1*z1 + z2*z2;
		v1 = z1 / a;
		v2 = -z2 / a;
		/* Initialize u1+u2*i as the rightmost factor 1-1/(z+n) */
		u1 = 1.0 - v1;
		u2 = -v2;
		k = 6.0 - z2 * .6 - zz1;
		if (!(k <= 0)) { /* IF(K.LE.0) GO TO 120 */
			/* Forward assemble of factors (z+j-1)/(z+n) */
			n -= k;
			uu1 = (zz1*z1 + z2*z2) / a;
			uu2 = dn*z2/a;
			vv1 = 0e0;
			vv2 = 0e0;
			for (j = 1; j <= k; j++) { /* DO 110 J = 1,K */
				b  = u1*(uu1+vv1) - u2*(uu2+vv2);
				u2 = u1*(uu2+vv2) + u2*(uu1+vv1);
				u1 = b;
				vv1 += v1;
/* 110 */
				vv2 += v2;
			}
		}
/* 120 */
		if (!(n <= 1)) { /* IF(N.LE.1) GO TO 140 */
			/* Backward assembly of factors 1-j/(z+n) */
			vv1 = v1;
			vv2 = v2;
			for (j = 2; j <= n; j++) { /* DO 130 J = 2,N */
				vv1 += v1;
				vv2 += v2;
				b = u1*(1.0 - vv1) + u2*vv2;
				u2 = -u1*vv2 + u2*(1.0 - vv1);
/* 130 */
				u1 = b;
			}
		}
/* 140 */
		u = u1*u1 + u2*u2;
		if (u == 0e0) { /* IF (U.EQ.0.0) GO TO 500 */
			aplerr = 130; return; }
		if (!(lf0 == 0) && k <= 0) break; /* GO TO 200 */
/* 150 */
		al1 = log(u)*.5;
		if (!(lf0 != 0)) { /* IF (lf0.NE.0) GO TO 160 */
			w1 = al1;
			w2 = atan2(u2,u1); /* arctangent? */
			if (w2 < 0e0) w2 += dpi;
			if (k <= 0) break; /* IF(K.LE.0) GO TO 200 */
		}
/* 160 */
		a = zz1 + z2 - delta;
		if (a <= 0.0) { /* IF (A.LE.0.0) GO TO 500 */
			aplerr = 130; return; }
		de0 -= al1;
		de1 += (2.0 + 1.0/a);
		break;
	} /* end for(;;) */
/* Case when real part of arg is greater than threshold. */
/* 200 */
	a = z1*z1 + z2*z2;
	al1 = log(a)*.5;
	al2 = atan2(z2,z1);
	v1 = y1*al1 - z2*al2;
	v2 = y1*al2 + z2*al1;
/* Evaluate asymptotic terms.  Ignore this term, if abs val(arg) .gt. */
/* 10**9, to save time and to avoid irrelevant underflow. */
	vv1 = 0e0;
	vv2 = 0e0;
	if (!(a > 1e18)) { /* IF(A.GT.1E18) GO TO 220 */
		uu1 = z1 / a;
		uu2 = -z2 / a;
		uuu1 = uu1 * uu1 - uu2 * uu2;
		uuu2 = uu1 * uu2 * 2.0;
		vv1 = *coef;
		for (j = 2; j < 8; j++) { /* DO 210 J=2,7 */
			b = vv1*uuu1 - vv2*uuu2;
			vv2 = vv1*uuu2 + vv2*uuu1;
/* 210 */
			vv1 = b + *(coef + j - 1);
		}
		b = vv1*uu1 - vv2*uu2;
		vv2 = vv1*uu2 + vv2*uu1;
		vv1 = b;
	}
/* 220 */
	w1 = (((vv1 + hl2p) - w1) - z1) + v1;
	w2 = ((vv2 - w2) - z2) + v2;
	de0 += ((v1 > 0e0) ? v1 : -v1) + ((v2 > 0e0) ? v2 : -v2);
	if (k <= 0) de1 += al1;
	/* final assembly */
	for (;;) {
		if (!(lf2 != 0)) { /* IF(LF2.NE.0) GO TO 310 */
			if (lf0 == 0) break; /* IF(lf0.EQ.0) GO TO 400 */
			a = exp(w1);
			w1 = a * cos(w2);
			w2 = a * sin(w2);
			if (lf3 == 0) break; /* IF(LF3.EQ.0) GO TO 400 */
			b  = (w1*u1 + w2*u2) / u;
			w2 = (w2*u1 - w1*u2) / u;
			w1 = b;
			break; /* GO TO 400 */
		}
/* 310 */
		h = h1*h1 + h2*h2;
		if (h == 0e0) { /* IF(H.EQ.0.0) GO TO 500 */
			aplerr = 130; return; }
		if ((lf0 == 0) || !(h > 1e-2)) {
/* 320 */
			a = log(h)*.5;
			if (h <= 1e-2) de0 -= a;
			if (!(lf0 != 0)) {
				w1 = (t1 - a) - w1;
				w2 = (t2 - atan2(h2,h1)) - w2;
				break; /* GO TO 400 */
			}
		}
/* 330 */
		t1 -= w1;
		t2 -= w2;
		a = exp(t1);
		t1 = a * cos(t2);
		t2 = a * sin(t2);
		w1 = (t1*h1 - t1*h2)/h;
		w2 = (t2*h1 - t1*h2)/h;
		if (lf3 == 0) break; /* GO TO 400 */
		b  = w1*u1 - w2*u2;
		w2 = w1*u2 + w2*u1;
		w1 = b;
		break;
	}
/* 400 */
	if (lf1 != 0) w2 = -w2;
	/* truncation error of stirlings formula is up to 3*10**-17. */
	de1 = de0*eps + 3e-17 + de1*delta;
	*cans = w1;
	*(cans+1) = w2;
	return;
}