/* NUM.c
 *
 * Copyright (C) 1992-2002 Paul Boersma
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 * pb 2002/03/03
 * pb 2002/03/07 GPL
 */

#include "NUM.h"
#include "NUM2.h"
#include <stdlib.h>
#include "melder.h"
#define SIGN(x,s) ((s) < 0 ? -fabs (x) : fabs(x))
#define my  me ->
/*
	GSL is more accurate than the other routines, but makes
	"Sound: To Intensity..." 10 times slower...
*/
#define USE_GSL  0
#if USE_GSL
	#include "gsl_errno.h"
	#include "gsl_sf_bessel.h"
#endif

#ifdef _WIN32
	double pow_WIN32_CW11 (double base, double exponent) {
		return base == 0.0 ? exponent == 0.0 ? 1.0 : 0.0 : (pow) (base, exponent);
	}
#endif

#define NUMcomplexFastFourierTransform_MACRO(proc,type)  \
static void proc (type *data, long n, int isign) { \
	long mmax = 2, m, j = 1, i; \
	for (i = 1; i < n; i += 2) { \
		if (j > i) { \
			type dum; \
			dum = data [j], data [j] = data [i], data [i] = dum; \
			dum = data [j+1], data [j+1] = data [i+1], data [i+1] = dum; \
		} \
		m = n >> 1; \
		while (m >= 2 && j > m) { j -= m; m >>= 1; } \
		j += m; \
	} \
	while (n > mmax) { \
		long istep = 2 * mmax; \
		double theta = 2.0 * NUMpi / (isign * mmax); \
		double wr, wi, wtemp, wpr, wpi; \
		wtemp = sin (0.5 * theta); \
		wpr = -2.0 * wtemp * wtemp; \
		wpi = sin (theta); \
		wr = 1.0, wi = 0.0; \
		for (m = 1; m < mmax; m += 2) { \
			for (i = m; i <= n; i += istep) { \
				type tempr, tempi; \
				long i1 = i + 1, i2 = i + mmax, i3 = i2 + 1; \
				tempr = wr * data [i2] - wi * data [i3], tempi = wr * data [i3] + wi * data [i2]; \
				data [i2] = data [i] - tempr, data [i3] = data [i1] - tempi; \
				data [i] += tempr, data [i1] += tempi; \
			} \
			wtemp = wr, wr = wr * wpr - wi * wpi + wr, wi = wi * wpr + wtemp * wpi + wi; \
		} \
		mmax = istep; \
	} \
}
NUMcomplexFastFourierTransform_MACRO (NUMcomplexFastFourierTransform_f, float)
NUMcomplexFastFourierTransform_MACRO (NUMcomplexFastFourierTransform_d, double)

#define NUMforwardRealFastFourierTransform_MACRO(proc,complexFastFourierTransform, type)  \
void proc (type *data, long n) { \
	long i, nby4 = n >> 2, nplus3 = n + 3; \
	type h1r, h1i, h2r, h2i; \
	double theta = - NUMpi / (double) (n >> 1); \
	double wtemp = sin (0.5 * theta), wpr = -2.0 * wtemp * wtemp, wpi = sin (theta); \
	double wr = wpr + 1.0, wi = wpi; \
	complexFastFourierTransform (data, n, -1); \
	h1r = data [1]; \
	h1i = data [2]; \
	data [1] = h1r + h1i; \
	data [2] = h1r - h1i; \
	for (i = 2; i <= nby4; i ++) { \
		long i1 = i + i - 1, i2 = i1 + 1, i3 = nplus3 - i2, i4 = i3 + 1; \
		h1r = 0.5 * (data [i1] + data [i3]); \
		h1i = 0.5 * (data [i2] - data [i4]); \
		h2r = 0.5 * (data [i2] + data [i4]); \
		h2i = -0.5 * (data [i1] - data [i3]); \
		data [i1] = h1r + wr * h2r - wi * h2i; \
		data [i2] = h1i + wr * h2i + wi * h2r; \
		data [i3] = h1r - wr * h2r + wi * h2i; \
		data [i4] = wr * h2i + wi * h2r - h1i; \
		wtemp = wr; \
		wr += wtemp * wpr - wi * wpi; \
		wi += wi * wpr + wtemp * wpi; \
	} \
}
NUMforwardRealFastFourierTransform_MACRO (NUMforwardRealFastFourierTransform_f, NUMcomplexFastFourierTransform_f, float)
NUMforwardRealFastFourierTransform_MACRO (NUMforwardRealFastFourierTransform_d, NUMcomplexFastFourierTransform_d, double)

#define NUMreverseRealFastFourierTransform_MACRO(proc,complexFastFourierTransform, type)  \
void proc (type *data, long n) { \
	long i, nby4 = n >> 2, nplus3 = n + 3; \
	type h1r, h1i, h2r, h2i; \
	double theta = NUMpi / (double) (n >> 1); \
	double wtemp = sin (0.5 * theta), wpr = -2.0 * wtemp * wtemp, wpi = sin (theta); \
	double wr = wpr + 1.0, wi = wpi; \
	h1r = data [1]; \
	h1i = data [2]; \
	data [1] = 0.5 * (h1r + h1i); \
	data [2] = 0.5 * (h1r - h1i); \
	for (i = 2; i <= nby4; i ++) { \
		long i1 = i + i - 1, i2 = i1 + 1, i3 = nplus3 - i2, i4 = i3 + 1; \
		h1r = 0.5 * (data [i1] + data [i3]); \
		h1i = 0.5 * (data [i2] - data [i4]); \
		h2r = -0.5 * (data [i2] + data [i4]); \
		h2i = 0.5 * (data [i1] - data [i3]); \
		data [i1] = h1r + wr * h2r - wi * h2i; \
		data [i2] = h1i + wr * h2i + wi * h2r; \
		data [i3] = h1r - wr * h2r + wi * h2i; \
		data [i4] = wr * h2i + wi * h2r - h1i; \
		wtemp = wr; \
		wr += wtemp * wpr - wi * wpi; \
		wi += wi * wpr + wtemp * wpi; \
	} \
	complexFastFourierTransform (data, n, 1); \
}
NUMreverseRealFastFourierTransform_MACRO (NUMreverseRealFastFourierTransform_f, NUMcomplexFastFourierTransform_f, float)
NUMreverseRealFastFourierTransform_MACRO (NUMreverseRealFastFourierTransform_d, NUMcomplexFastFourierTransform_d, double)

void NUMrealft_f (float *data, long n, int isign) {
	if (isign == 1) NUMforwardRealFastFourierTransform_f (data, n);
	else NUMreverseRealFastFourierTransform_f (data, n);
}
void NUMrealft_d (double *data, long n, int isign) {
	if (isign == 1) NUMforwardRealFastFourierTransform_d (data, n);
	else NUMreverseRealFastFourierTransform_d (data, n);
}

void NUMmemcof (float data [], long n, long m, float *xms, float d []) {
	long k, j, i;
	float p = 0.0;
	float *wk1 = NUMfvector (1, n), *wk2 = NUMfvector (1, n), *wkm = NUMfvector (1, m);
	for (j = 1; j <= n; j ++) p += data [j] * data [j];
	*xms = p / n;
	wk1 [1] = data [1];
	wk2 [n - 1] = data [n];
	for (j = 2; j <= n - 1; j ++) { wk1 [j] = data [j]; wk2 [j - 1] = data [j]; }
	for (k = 1; k <= m; k ++) {
		float num = 0.0, denom = 0.0;
		for (j = 1; j <= n - k; j ++) {
			num += wk1 [j] * wk2 [j];
			denom += wk1[j] * wk1 [j] + wk2 [j] * wk2 [j];
		}
		Melder_assert (denom != 0.0);
		d [k] = 2.0 * num / denom;
		*xms *= 1.0 - d[k] * d [k];
		for (i = 1; i <= k - 1; i ++) d [i] = wkm [i] - d [k] * wkm [k - i];
		if (k == m) {
			NUMfvector_free (wkm, 1);
			NUMfvector_free (wk2, 1);
			NUMfvector_free (wk1, 1);
			return;
		}
		for (i = 1; i <= k; i ++) wkm [i] = d [i];
		for (j = 1; j <= n - k - 1; j ++)
			{ wk1 [j] -= wkm [k] * wk2 [j]; wk2 [j] = wk2 [j + 1] - wkm [k] * wk1 [j + 1]; }
	}
	Melder_fatal ("NUMmemcof: should never get here.");
}

void NUMfbtoa (double formant, double bandwidth, double dt, double *a1, double *a2) {
	*a1 = 2 * exp (- NUMpi * bandwidth * dt) * cos (2 * NUMpi * formant * dt);
	*a2 = exp (- 2 * NUMpi * bandwidth * dt);
}

void NUMfilterSecondOrderSection_a (float x [], long n, double a1, double a2) {
	long i;
	x [2] += a1 * x [1];
	for (i = 3; i <= n; i ++)
		x [i] += a1 * x [i - 1] - a2 * x [i - 2];
}

void NUMfilterSecondOrderSection_fb (float x [], long n, double dt, double formant, double bandwidth) {
	double a1, a2;
	NUMfbtoa (formant, bandwidth, dt, & a1, & a2);
	NUMfilterSecondOrderSection_a (x, n, a1, a2);
}

double NUMftopreemphasis (double f, double dt) {
	return exp (- 2.0 * NUMpi * f * dt);
}

void NUMpreemphasize_a (float x [], long n, double preemphasis) {
	long i;
	for (i = n; i >= 2; i --)
		x [i] -= preemphasis * x [i - 1];
}

void NUMdeemphasize_a (float x [], long n, double preemphasis) {
	long i;
	for (i = 2; i <= n; i ++)
		x [i] += preemphasis * x [i - 1];
}

void NUMpreemphasize_f (float x [], long n, double dt, double frequency) {
	NUMpreemphasize_a (x, n, NUMftopreemphasis (frequency, dt));
}

void NUMdeemphasize_f (float x [], long n, double dt, double frequency) {
	NUMdeemphasize_a (x, n, NUMftopreemphasis (frequency, dt));
}

void NUMautoscale (float x [], long n, double scale) {
	long i;
	double maximum = 0.0;
	for (i = 1; i <= n; i ++)
		if (fabs (x [i]) > maximum) maximum = fabs (x [i]);
	if (maximum > 0.0) {
		double factor = scale / maximum;
		for (i = 1; i <= n; i ++)
			x [i] *= factor;
	}
}

double NUMlnGamma (double x) {
	/*
		Press et al. advise to use a reflection formula for x < 1.0,
		but we are not so sure, as it gives 37 for all values below 10^-17.
		The series expansion simply goes on rising, giving e.g. 69 for 10^-30.
		The 69 is correct, since Gamma(10^30) = Gamma(1+10^-30)/10^-30
		which is approximately 10^30, and its logarithm is 30 ln (10).
		For such values, the procedure is accurate within 10^-10,
		as it is for values above 1 (according to Press et al.).
	*/
	if (x <= 0.0) {
		return NUMundefined;
	} else {
		double series = 1.0 + 76.18009173 / x - 86.50532033 / (x + 1.0) + 24.01409822 / (x + 2.0)
			- 1.231739516 / (x + 3.0) + 0.00120858003 / (x + 4.0) - 0.00000536382 / (x + 5.0);
		return - (x + 4.5) + (x - 0.5) * log (x + 4.5) + log (2.50662827465 * series);
	}
}

double NUMbeta (double z, double w) {
	if (z <= 0.0 || w <= 0.0) return NUMundefined;
	return exp (NUMlnGamma (z) + NUMlnGamma (w) - NUMlnGamma (z + w));
}

double NUMbinomialP (double p, double k, double n) {
	double binomialQ;
	if (p < 0.0 || p > 1.0 || n <= 0.0 || k < 0.0 || k > n) return NUMundefined;
	if (k == n) return 1.0;
	binomialQ = NUMincompleteBeta (k + 1, n - k, p);
	if (binomialQ == NUMundefined) return NUMundefined;
	return 1.0 - binomialQ;
}

double NUMbinomialQ (double p, double k, double n) {
	if (p < 0.0 || p > 1.0 || n <= 0.0 || k < 0.0 || k > n) return NUMundefined;
	if (k == 0.0) return 1.0;
	return NUMincompleteBeta (k, n - k + 1, p);
}

#define ITMAX 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30
static double NUMridders3 (double (*func) (double, double, double), double constant, double df1, 
	double df2, double x1, double x2, double xacc)
{
	long iter, maxit = 100;
	double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew;

	fl = (*func)(x1, df1, df2) - constant;
	fh = (*func)(x2, df1, df2) - constant;
	if ((fl > 0 && fh < 0) || (fl < 0 && fh > 0))
	{
		xl = x1; xh = x2;
		ans = NUMundefined;
		for (iter=1; iter <= maxit; iter++)
		{
			xm = 0.5 * (xl + xh);
			fm = func (xm, df1, df2) - constant;
			s = sqrt (fm * fm - fl * fh);
			if (s == 0.0) return ans;
			xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s);
			if (fabs (xnew - ans) <= xacc) return ans;
			ans = xnew;
			fnew = func (ans, df1, df2) - constant;
			if (fnew == 0.0) return ans;
			if (SIGN (fm, fnew) != fm)
			{
				xl = xm; fl = fm;
				xh = ans; fh = fnew;
			}
			else if (SIGN (fl, fnew) != fl)
			{
				xh = ans; fh = fnew;
			}
			else if (SIGN (fh, fnew) != fh)
			{
				xl = ans; fl = fnew;
			}
			else
			{
				(void) Melder_error("NUMridders3: never get here.");
				return NUMundefined;
			}
			if (fabs (xh - xl) <= xacc) return ans;
		}
		(void) Melder_error("NUMridders3: exceeds maximum iterations.");
	}
	else
	{
		if (fl == 0.0) return x1;
		if (fh == 0.0) return x2;
		(void) Melder_error("NUMridders3: root must be bracketed.");
	}
	return NUMundefined;
}
#undef ITMAX
#undef EPS
#undef FPMIN

double NUMinvBinomialP (double p, double k, double n) {
	if (p < 0 || p > 1 || n <= 0 || k < 0 || k > n) return NUMundefined;
	if (k == n) return 1.0;
	return NUMridders3 (NUMbinomialP, p, k, n, 0.0, 1.0, 1.0e-10);
}

double NUMinvBinomialQ (double p, double k, double n) {
	if (p < 0 || p > 1 || n <= 0 || k < 0 || k > n) return NUMundefined;
	if (k == 0) return 0.0;
	return NUMridders3 (NUMbinomialQ, p, k, n, 0.0, 1.0, 1.0e-10);
}

static double NUMbessj0 (double x) {
	double ax, z, xx, y, ans1, ans2;
	if ((ax = fabs (x)) < 8.0) {
		y = x * x;
		ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7
				 + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
		ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718
				 + y * (59272.64853 + y * (267.8532712 + y * 1.0))));
		return ans1 / ans2;
	}
	z = 8.0 / ax;   /* <= 1.0 */
	y = z * z;
	xx = ax - 0.785398164;
	ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
			 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
	ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
			 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7)));
	return sqrt (0.636619772 / ax) * (cos (xx) * ans1 - z * sin (xx) * ans2);
}

static double NUMbessy0 (double x) {
	double z, xx, y, ans1, ans2;
	if (x < 8.0) {
		y = x * x;
		ans1 = -2957821389.0 + y * (7062834065.0 + y * (-512359803.6
				 + y * (10879881.29 + y * (-86327.92757 + y * 228.4622733))));
		ans2 = 40076544269.0 + y * (745249964.8 + y * (7189466.438
				 + y * (47447.26470 + y * (226.1030244 + y * 1.0))));
		return (ans1 / ans2) + 0.636619772 * NUMbessj0 (x) * log (x);
	}
	z = 8.0 / x;   /* <= 1.0 */
	y = z * z;
	xx = x - 0.785398164;
	ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
			 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
	ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
		 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 + y * (-0.934945152e-7))));
	return sqrt (0.636619772 / x) * (sin (xx) * ans1 + z * cos (xx) * ans2);
}

static double NUMbessj1 (double x) {
	double ax, z, xx, y, result, ans1, ans2;
	if ((ax = fabs (x)) < 8.0) {
		y = x * x;
		ans1 = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1
				 + y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
		ans2 = 144725228442.0 + y * (2300535178.0 + y * (18583304.74
				 + y * (99447.43394 + y * (376.9991397 + y * 1.0))));
		return ans1 / ans2;
	}
	z = 8.0 / ax;   /* <= 1.0 */
	y = z * z;
	xx = ax - 2.356194491;
	ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
			 + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
	ans2 = 0.04687499995 + y * (-0.2002690873e-3
			 + y * (0.8449199096e-5 + y * (-0.88228987e-6 + y * 0.105787412e-6)));
	result = sqrt (0.636619772 / ax) * (cos (xx) * ans1 - z * sin (xx) * ans2);
	if (x < 0.0) result = - result;
	return result;
}

static double NUMbessy1 (double x) {
	double z, xx, y, ans1, ans2;
	if (x < 8.0) {
		y = x * x;
		ans1 = x * (-0.4900604943e13 + y * (0.1275274390e13
				 + y * (-0.5153438139e11 + y * (0.7349264551e9
				 + y * (-0.4237922726e7 + y * 0.8511937935e4)))));
		ans2 = 0.2499580570e14 + y * (0.4244419664e12
				 + y * (0.3733650367e10 + y * (0.2245904002e8
				 + y * (0.1020426050e6 + y * (0.3549632885e3 + y)))));
		return (ans1 / ans2) + 0.636619772 * (NUMbessj1 (x) * log (x) - 1.0 / x);
	}
	z = 8.0 / x;   /* <= 1.0 */
	y = z * z;
	xx = x - 2.356194491;
	ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
			 + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
	ans2 = 0.04687499995 + y * (-0.2002690873e-3
			 + y * (0.8449199096e-5 + y * (-0.88228987e-6 + y * 0.105787412e-6)));
	return sqrt (0.636619772 / x) * (sin (xx) * ans1 + z * cos (xx) * ans2);
}

double NUMbesselY (long n, double x) {
	long j;
	double by, bym, byp, tox;
	Melder_assert (x > 0);
	if (n == 0) return NUMbessy0 (x);
	if (n == 1) return NUMbessy1 (x);
	Melder_assert (n >= 2);
	tox = 2.0 / x;
	by = NUMbessy1 (x);
	bym = NUMbessy0 (x);
	for (j = 1; j < n; j ++) { byp = j * tox * by - bym; bym = by; by = byp; }
	return by;
}

#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10
double NUMbesselJ (long n, double x) {
	int jsum;
	long j, m;
	double ax, bj, bjm, bjp, sum, tox, result;
	if (n == 0) return NUMbessj0 (x);
	if (n == 1) return NUMbessj1 (x);
	Melder_assert (n >= 2);
	ax = fabs (x);
	if (ax == 0.0) return 0.0;
	tox = 2.0 / ax;
	if (ax > n) {
		bjm = NUMbessj0 (ax);
		bj = NUMbessj1 (ax);
		for (j = 1; j < n; j ++) { bjp = j * tox * bj - bjm; bjm = bj; bj = bjp; }
		return bj;
	}
	m = 2 * ((n + (long) sqrt (ACC * n)) / 2);
	jsum = 0;
	bjp = result = sum = 0.0;
	bj = 1.0;
	for (j = m; j > 0; j --) {
		bjm = j * tox * bj - bjp; bjp = bj; bj = bjm;
		if (fabs (bj) > BIGNO) { bj *= BIGNI; bjp *= BIGNI; result *= BIGNI; sum *= BIGNI; }
		if (jsum) sum += bj;
		jsum = ! jsum;
		if (j == n) result = bjp;
	}
	sum = 2.0 * sum - bj;
	result /= sum;
	return x < 0.0 && (n & 1) ? - result : result;
}
#undef ACC
#undef BIGNO
#undef BIGNI

/* Modified Bessel function I0. Abramowicz & Stegun, p. 378.*/
static double NUMbessel_i0 (double x) {
	double t;
	if (x < 0.0) return NUMbessel_i0 (- x);
	if (x < 3.75) {
		/* Formula 9.8.1. Accuracy 1.6e-7. */
		t = x / 3.75;
		t *= t;
		return 1.0 + t * (3.5156229 + t * (3.0899424 + t * (1.2067492
			+ t * (0.2659732 + t * (0.0360768 + t * 0.0045813)))));
	}
	/*
		otherwise: x >= 3.75
	*/
	/* Formula 9.8.2. Accuracy of the polynomial factor 1.9e-7. */
	t = 3.75 / x;   /* <= 1.0 */
	return exp (x) / sqrt (x) * (0.39894228 + t * (0.01328592
		+ t * (0.00225319 + t * (-0.00157565 + t * (0.00916281
		+ t * (-0.02057706 + t * (0.02635537 + t * (-0.01647633
		+ t * 0.00392377))))))));
}

/* Modified Bessel function I1. Abramowicz & Stegun, p. 378. */
static double NUMbessel_i1 (double x) {
	double t;
	if (x < 0.0) return - NUMbessel_i1 (- x);
	if (x < 3.75) {
		/* Formula 9.8.3. Accuracy of the polynomial factor 8e-9. */
		t = x / 3.75;
		t *= t;
		return x * (0.5 + t * (0.87890594 + t * (0.51498869 + t * (0.15084934
			+ t * (0.02658733 + t * (0.00301532 + t * 0.00032411))))));
	}
	/*
		otherwise: x >= 3.75
	*/
	/* Formula 9.8.4. Accuracy of the polynomial factor 2.2e-7. */
	t = 3.75 / x;   /* <= 1.0 */
	return exp (x) / sqrt (x) * (0.39894228 + t * (-0.03988024
		+ t * (-0.00362018 + t * (0.00163801 + t * (-0.01031555
		+ t * (0.02282967 + t * (-0.02895312 + t * (0.01787654
		+ t * (-0.00420059)))))))));
}

double NUMbesselI (long n, double x) {
#if USE_GSL
	gsl_sf_result result;
	int status = gsl_sf_bessel_In_e (n, x, & result);
	return status == GSL_SUCCESS ? result. val : NUMundefined;
#else
	#define ACC 40.0
	#define BIGNO 1.0e10
	#define BIGNI 1.0e-10
	long j;
	double bi, bim, bip, tox, result;
	if (n == 0) return NUMbessel_i0 (x);
	if (n == 1) return NUMbessel_i1 (x);
	Melder_assert (n >= 2);
	if (x == 0.0) return 0.0;
	tox = 2.0 / fabs (x);
	bip = result = 0.0;
	bi = 1.0;
	for (j = 2 * (n + (long) sqrt (ACC * n)); j > 0; j --) {
		bim = bip + j * tox * bi; bip = bi; bi = bim;
		if (fabs (bi) > BIGNO) { result *= BIGNI; bi *= BIGNI; bip *= BIGNI; }
		if (j == n) result = bip;
	}
	result *= NUMbessel_i0 (x) / bi;
	return x < 0.0 && (n & 1) ? - result : result;
	#undef ACC
	#undef BIGNO
	#undef BIGNI
#endif
}

/* Modified Bessel function K0. Abramowicz & Stegun, p. 379. */
static double NUMbessel_k0 (double x) {
	double t;
	if (x <= 0.0) return NUMundefined;   /* Positive infinity. */
	if (x <= 2.0) {
		/* Formula 9.8.5. Accuracy 1e-8. */
		double x2 = 0.5 * x, t = x2 * x2;
		return - log (x2) * NUMbessel_i0 (x) + (-0.57721566 + t * (0.42278420
			+ t * (0.23069756 + t * (0.03488590 + t * (0.00262698
			+ t * (0.00010750 + t * 0.00000740))))));
	}
	/*
		otherwise: 2 < x < positive infinity
	*/
	/* Formula 9.8.6. Accuracy of the polynomial factor 1.9e-7. */
	t = 2.0 / x;   /* < 1.0 */
	return exp (- x) / sqrt (x) * (1.25331414 + t * (-0.07832358
		+ t * (0.02189568 + t * (-0.01062446 + t * (0.00587872
		+ t * (-0.00251540 + t * 0.00053208))))));
}

/* Modified Bessel function K1. Abramowicz & Stegun, p. 379. */
static double NUMbessel_k1 (double x) {
	double t;
	if (x <= 0.0) return NUMundefined;   /* Positive infinity. */
	if (x <= 2.0) {
		/* Formula 9.8.7. Accuracy  of the polynomial factor 8e-9. */
		double x2 = 0.5 * x, t = x2 * x2;
		return log (x2) * NUMbessel_i1 (x) + (1.0 / x) * (1.0 + t * (0.15443144
			+ t * (-0.67278579 + t * (-0.18156897 + t * (-0.01919402
			+ t * (-0.00110404 + t * (-0.00004686)))))));
	}
	/*
		otherwise: 2 < x < positive infinity
	*/
	/* Formula 9.8.8. Accuracy of the polynomial factor 2.2e-7. */
	t = 2.0 / x;   /* < 1.0 */
	return exp (- x) / sqrt (x) * (1.25331414 + t * (0.23498619
			 + t * (-0.03655620 + t * (0.01504268 + t * (-0.00780353
			 + t * (0.00325614 + t * (-0.00068245)))))));
}

double NUMbesselK (long n, double x) {
	long j;
	double bk, bkm, bkp, tox;
	Melder_assert (x > 0);
	if (n == 0) return NUMbessel_k0 (x);
	if (n == 1) return NUMbessel_k1 (x);
	Melder_assert (n >= 2);
	tox = 2.0/x;
	bkm = NUMbessel_k0 (x);
	bk = NUMbessel_k1 (x);
	for (j = 1; j < n; j ++) { bkp = bkm + j * tox * bk; bkm = bk; bk = bkp; }
	return bk;
}

double NUMsigmoid (double x)
	{ return x > 0.0 ? 1 / (1 + exp (- x)) : 1 - 1 / (1 + exp (x)); }

double NUMerfcc (double x) {
	double t = 1 / (1 + 0.5 * fabs (x));
	double result = t * exp (- x * x - 1.26551223 + t * (1.00002368 + t * (0.37409196 +
		t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 +
		t * (1.48851587 + t * (-0.82215223 + t * 0.17087277)))))))));
	return x >= 0.0 ? result : 2 - result;
}

double NUMgaussP (double z) {
	return 1 - 0.5 * NUMerfcc (NUMsqrt1_2 * z);
}

double NUMgaussQ (double z) {
	return 0.5 * NUMerfcc (NUMsqrt1_2 * z);
}

static double NUMgammaSeries (double a, double x) {
	long n, ITMAX = 100;
	double EPS = 3e-7, ap = a, sum = 1.0 / a, del = sum;
	double lnGamma = NUMlnGamma (a);
	if (x < 0.0) return NUMundefined;
	if (x == 0.0) return 0.0;
	del = sum = 1.0 / a;
	for (n = 1; n <= ITMAX; n ++) {
		++ ap;
		del *= x / ap;
		sum += del;
		if (fabs (del) < fabs (sum) * EPS)
			return sum * exp (- x + a * log (x) - lnGamma);
	}
	return NUMundefined;
}

static double NUMgammaContinuedFraction (double a, double x) {
	int i, ITMAX = 100;
	double an ,b, c, d, del, h, EPS = 3e-7, FPMIN = 1e-30;
	double lnGamma = NUMlnGamma (a);
	b = x + 1.0 - a;
	c = 1.0 / FPMIN;
	d = 1.0 / b;
	h = d;
	for (i = 1; i <= ITMAX; i ++) {
		an = -i * (i-a); b += 2.0; d = an * d + b;
		if (fabs (d) < FPMIN) d = FPMIN;
		c = b + an / c;
		if (fabs (c) < FPMIN) c = FPMIN;
		d = 1.0 / d; del = d * c; h *= del;
		if (fabs (del - 1.0) < EPS) break;
	}
	if (i > ITMAX) return NUMundefined;
	return exp (- x + a * log (x) - lnGamma) * h;
}

double NUMincompleteGammaP (double a, double x) {
	double q;
	if (a <= 0.0 || x < 0.0) return NUMundefined;
	if (x < a + 1.0)
		return NUMgammaSeries (a, x);
	q = NUMgammaContinuedFraction (a, x);
	return 1 - q;
}

double NUMincompleteGammaQ (double a, double x) {
	if (a <= 0.0 || x < 0.0) return NUMundefined;
	if (x < a + 1.0) {
		double p = NUMgammaSeries (a, x);
		if (p == NUMundefined) return NUMundefined;
		return 1 - p;
	}
	return NUMgammaContinuedFraction (a, x);
}

double NUMchiSquareP (double chiSquare, long degreesOfFreedom) {
	if (chiSquare < 0 || degreesOfFreedom <= 0) return NUMundefined;
	return NUMincompleteGammaP (0.5 * degreesOfFreedom, 0.5 * chiSquare);
}

double NUMchiSquareQ (double chiSquare, long degreesOfFreedom) {
	if (chiSquare < 0 || degreesOfFreedom <= 0) return NUMundefined;
	return NUMincompleteGammaQ (0.5 * degreesOfFreedom, 0.5 * chiSquare);
}

double NUMcombinations (long n, long k) {
	double result = 1.0;
	long i;
	if (k > n / 2) k = n - k;
	for (i = 1; i <= k; i ++) result *= n - i + 1;
	for (i = 2; i <= k; i ++) result /= i;
	return result;
}

#define BRENT_MAXIMUM_NUMBER_OF_ITERATIONS 100

double NUM_minimize
	(double ax, double bx, double cx, double (*f) (double x, void *closure), void *closure,
	 double absoluteTolerance, double fractionalTolerance, double *xmin)
{
	int iter;
	double a, b, d, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm;
	/* There may be a compiler warning "d used before set." Can occur in case of zero tolerance. */
	double e = 0.0, goldenSection2 = 1.0 - (0.5 * sqrt (5.0) - 0.5);

	a = ax < cx ? ax : cx;
	b = ax > cx ? ax : cx;
	x = w = v = bx;
	fw = fv = fx = f (x, closure);
	for (iter = 1; iter <= BRENT_MAXIMUM_NUMBER_OF_ITERATIONS; iter ++) {
		xm = 0.5 * (a + b);
		tol1 = absoluteTolerance + fractionalTolerance * fabs (x);
		tol2 = 2.0 * tol1;
		if (fabs (x - xm) <= (tol2 - 0.5 * (b - a))) { *xmin = x; return fx; }
		if (fabs (e) > tol1) {
			r = (x - w) * (fx - fv);
			q = (x - v) * (fx - fw);
			p = (x - v) * q - (x - w) * r;
			q = 2.0 * (q - r);
			if (q > 0.0) p = -p;
			q = fabs (q);
			etemp = e;
			e = d;
			if (fabs (p) >= fabs (0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x)) {
				e = x >= xm ? a - x : b - x;
				d = goldenSection2 * e;
			} else {
				d = p / q;
				u = x + d;
				if (u - a < tol2 || b - u < tol2) d = xm > x ? tol1 : - tol1;
			}
		} else {
			d = goldenSection2 * (e = (x >= xm ? a - x : b - x));
		}
		u = fabs (d) >= tol1 ? x + d : x + ( d > 0.0 ? tol1 : - tol1 );
		fu = f (u, closure);
		if (fu <= fx) {
			if (u >= x) a = x; else b = x;
			v = w; w = x; x = u;
			fv = fw; fw = fx; fx = fu;
		} else {
			if (u < x) a = u; else b = u;
			if (fu <= fw || w == x) { v = w; w = u; fv = fw; fw = fu; }
			else if (fu <= fv || v == x || v == w) { v = u; fv = fu; }
		}
	}
	Melder_warning ("(NUM_minimize:) Too many iterations.");
	*xmin = x;
	return fx;
}

#define NUM_interpolate_simple_cases \
	if (nx < 1) return NUMundefined; \
	if (x > nx) return y [nx]; \
	if (x < 1) return y [1]; \
	if (x == midleft) return y [midleft]; \
	/* 1 < x < nx && x not integer: interpolate. */ \
	if (maxDepth > midright - 1) maxDepth = midright - 1; \
	if (maxDepth > nx - midleft) maxDepth = nx - midleft; \
	if (maxDepth <= 0) return y [(long) floor (x + 0.5)]; \
	if (maxDepth == 1) return y [midleft] + (x - midleft) * (y [midright] - y [midleft]); \
	if (maxDepth == 2) { \
		double yl = y [midleft], yr = y [midright]; \
		double dyl = 0.5 * (yr - y [midleft - 1]), dyr = 0.5 * (y [midright + 1] - yl); \
		double fil = x - midleft, fir = midright - x; \
		return yl * fir + yr * fil - fil * fir * (0.5 * (dyr - dyl) + (fil - 0.5) * (dyl + dyr - 2 * (yr - yl))); \
	}

#if defined (sgi) || defined (sun) || defined (HPUX) || __POWERPC__
#define NUM_interpolate_sinc_MACRO(proc,type) \
double proc (type y [], long nx, double x, long maxDepth) { \
	long ix, midleft = floor (x), midright = midleft + 1, left, right; \
	double result = 0.0, a, halfsina, aa, daa, cosaa, sinaa, cosdaa, sindaa; \
	NUM_interpolate_simple_cases \
	left = midright - maxDepth, right = midleft + maxDepth; \
	a = NUMpi * (x - midleft); \
	halfsina = 0.5 * sin (a); \
	aa = a / (x - left + 1); cosaa = cos (aa); sinaa = sin (aa); \
	daa = NUMpi / (x - left + 1); cosdaa = cos (daa); sindaa = sin (daa); \
	for (ix = midleft; ix >= left; ix --) { \
		double d = halfsina / a * (1.0 + cosaa), help; \
		result += y [ix] * d; \
		a += NUMpi; \
		help = cosaa * cosdaa - sinaa * sindaa; \
		sinaa = cosaa * sindaa + sinaa * cosdaa; \
		cosaa = help; \
		halfsina = - halfsina; \
	} \
	a = NUMpi * (midright - x); \
	halfsina = 0.5 * sin (a); \
	aa = a / (right - x + 1); cosaa = cos (aa); sinaa = sin (aa); \
	daa = NUMpi / (right - x + 1); cosdaa = cos (daa); sindaa = sin (daa); \
	for (ix = midright; ix <= right; ix ++) { \
		double d = halfsina / a * (1.0 + cosaa), help; \
		result += y [ix] * d; \
		a += NUMpi;	 \
		help = cosaa * cosdaa - sinaa * sindaa; \
		sinaa = cosaa * sindaa + sinaa * cosdaa; \
		cosaa = help; \
		halfsina = - halfsina; \
	} \
	return result; \
}
#else
#define NUM_interpolate_sinc_MACRO(proc,type) \
double proc (type y [], long nx, double x, long maxDepth) { \
	long ix, midleft = floor (x), midright = midleft + 1, left, right; \
	double result = 0.0, a, halfsina, aa, daa; \
	NUM_interpolate_simple_cases \
	left = midright - maxDepth, right = midleft + maxDepth; \
	a = NUMpi * (x - midleft); \
	halfsina = 0.5 * sin (a); \
	aa = a / (x - left + 1); \
	daa = NUMpi / (x - left + 1); \
	for (ix = midleft; ix >= left; ix --) { \
		double d = halfsina / a * (1.0 + cos (aa)); \
		result += y [ix] * d; \
		a += NUMpi;	 \
		aa += daa;	 \
		halfsina = - halfsina; \
	} \
	a = NUMpi * (midright - x); \
	halfsina = 0.5 * sin (a); \
	aa = a / (right - x + 1); \
	daa = NUMpi / (right - x + 1); \
	for (ix = midright; ix <= right; ix ++) { \
		double d = halfsina / a * (1.0 + cos (aa)); \
		result += y [ix] * d; \
		a += NUMpi;	 \
		aa += daa; \
		halfsina = - halfsina; \
	} \
	return result; \
}
#endif
NUM_interpolate_sinc_MACRO (NUM_interpolate_sinc_f, float)
NUM_interpolate_sinc_MACRO (NUM_interpolate_sinc_d, double)

/********** Improving extrema **********/

struct improve_params {
	int depth;
	float *y_f;
	double *y_d;
	long ixmax;
	int isMaximum;
};

static double improve_evaluate_f (double x, void *closure) {
	struct improve_params *me = closure;
	double y = NUM_interpolate_sinc_f (my y_f, my ixmax, x, my depth);
	return my isMaximum ? - y : y;
}
static double improve_evaluate_d (double x, void *closure) {
	struct improve_params *me = closure;
	double y = NUM_interpolate_sinc_d (my y_d, my ixmax, x, my depth);
	return my isMaximum ? - y : y;
}

double NUMimproveExtremum_f (float *y, long nx, long ixmid, int interpolation, double *ixmid_real, int isMaximum) {
	struct improve_params params;
	if (ixmid <= 1) { *ixmid_real = 1; return y [1]; }
	if (ixmid >= nx) { *ixmid_real = nx; return y [nx]; }
	if (interpolation <= NUM_PEAK_INTERPOLATE_NONE) { *ixmid_real = ixmid; return y [ixmid]; }
	if (interpolation == NUM_PEAK_INTERPOLATE_PARABOLIC) {
		double dy = 0.5 * (y [ixmid + 1] - y [ixmid - 1]);
		double d2y = 2 * y [ixmid] - y [ixmid - 1] - y [ixmid + 1];
		*ixmid_real = ixmid + dy / d2y;
		return y [ixmid] + 0.5 * dy * dy / d2y;
	}
	/* Sinc interpolation. */
	params. y_f = y;
	params. depth = interpolation == NUM_PEAK_INTERPOLATE_CUBIC ? 2 :
		interpolation == NUM_PEAK_INTERPOLATE_SINC70 ? 70 : 700;
	params. ixmax = nx;
	params. isMaximum = isMaximum;
	return isMaximum ?
		- NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_f, & params, 1e-10, 1e-7, ixmid_real) :
		 NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_f, & params, 1e-10, 1e-7, ixmid_real);
}

double NUMimproveMaximum_f (float *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_f (y, nx, ixmid, interpolation, ixmid_real, 1); }
double NUMimproveMinimum_f (float *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_f (y, nx, ixmid, interpolation, ixmid_real, 0); }

double NUMimproveExtremum_d (double *y, long nx, long ixmid, int interpolation, double *ixmid_real, int isMaximum) {
	struct improve_params params;
	if (ixmid <= 1) { *ixmid_real = 1; return y [1]; }
	if (ixmid >= nx) { *ixmid_real = nx; return y [nx]; }
	if (interpolation <= NUM_PEAK_INTERPOLATE_NONE) { *ixmid_real = ixmid; return y [ixmid]; }
	if (interpolation == NUM_PEAK_INTERPOLATE_PARABOLIC) {
		double dy = 0.5 * (y [ixmid + 1] - y [ixmid - 1]);
		double d2y = 2 * y [ixmid] - y [ixmid - 1] - y [ixmid + 1];
		*ixmid_real = ixmid + dy / d2y;
		return y [ixmid] + 0.5 * dy * dy / d2y;
	}
	/* Sinc interpolation. */
	params. y_d = y;
	params. depth = interpolation == NUM_PEAK_INTERPOLATE_SINC70 ? 70 : 700;
	params. ixmax = nx;
	params. isMaximum = isMaximum;
	return isMaximum ?
		- NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_d, & params, 1e-10, 1e-11, ixmid_real) :
		NUM_minimize (ixmid - 1, ixmid, ixmid + 1, improve_evaluate_d, & params, 1e-10, 1e-11, ixmid_real);
}

double NUMimproveMaximum_d (double *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_d (y, nx, ixmid, interpolation, ixmid_real, 1); }
double NUMimproveMinimum_d (double *y, long nx, long ixmid, int interpolation, double *ixmid_real)
	{ return NUMimproveExtremum_d (y, nx, ixmid, interpolation, ixmid_real, 0); }

/********** Viterbi **********/

int NUM_viterbi (
	long numberOfFrames, long maxnCandidates,
	long (*getNumberOfCandidates) (long iframe, void *closure),
	double (*getLocalCost) (long iframe, long icand, void *closure),
	double (*getTransitionCost) (long iframe, long icand1, long icand2, void *closure),
	void (*putResult) (long iframe, long place, void *closure),
	void *closure)
{
	double **delta = NULL, maximum;
	long **psi = NULL, *numberOfCandidates = NULL, iframe, icand, icand1, icand2, place;
	if (! (delta = NUMdmatrix (1, numberOfFrames, 1, maxnCandidates))) goto end;
	if (! (psi = NUMlmatrix (1, numberOfFrames, 1, maxnCandidates))) goto end;
	if (! (numberOfCandidates = NUMlvector (1, numberOfFrames))) goto end;
	for (iframe = 1; iframe <= numberOfFrames; iframe ++) {
		numberOfCandidates [iframe] = getNumberOfCandidates (iframe, closure);
		for (icand = 1; icand <= numberOfCandidates [iframe]; icand ++)
			delta [iframe] [icand] = - getLocalCost (iframe, icand, closure);
	}
	for (iframe = 2; iframe <= numberOfFrames; iframe ++) {
		for (icand2 = 1; icand2 <= numberOfCandidates [iframe]; icand2 ++) {
			maximum = -1e300;
			place = 0;
			for (icand1 = 1; icand1 <= numberOfCandidates [iframe - 1]; icand1 ++) {
				double value = delta [iframe - 1] [icand1] + delta [iframe] [icand2]
					- getTransitionCost (iframe, icand1, icand2, closure);
				if (value > maximum) { maximum = value; place = icand1; }
			}
			delta [iframe] [icand2] = maximum;
			psi [iframe] [icand2] = place;
		}
	}
	/*
	 * Find the end of the most probable path.
	 */
	maximum = delta [numberOfFrames] [place = 1];
	for (icand = 2; icand <= numberOfCandidates [numberOfFrames]; icand ++)
		if (delta [numberOfFrames] [icand] > maximum)
			maximum = delta [numberOfFrames] [place = icand];
	/*
	 * Backtrack.
	 */
	for (iframe = numberOfFrames; iframe >= 1; iframe --) {
		putResult (iframe, place, closure);
		place = psi [iframe] [place];
	}
end:
	NUMdmatrix_free (delta, 1, 1);
	NUMlmatrix_free (psi, 1, 1);
	NUMlvector_free (numberOfCandidates, 1);
	if (Melder_hasError ()) return 0;
	return 1;
}

/******************/

struct parm2 {
	int ntrack;
	long ncomb;
	long **indices;
	double (*getLocalCost) (long iframe, long icand, int itrack, void *closure);
	double (*getTransitionCost) (long iframe, long icand1, long icand2, int itrack, void *closure);
	void (*putResult) (long iframe, long place, int itrack, void *closure);
	void *closure;
};

static long getNumberOfCandidates_n (long iframe, void *closure) {
	struct parm2 *me = closure;
	(void) iframe;
	return my ncomb;
}
static double getLocalCost_n (long iframe, long jcand, void *closure) {
	struct parm2 *me = closure;
	int itrack;
	double localCost = 0.0;
	for (itrack = 1; itrack <= my ntrack; itrack ++)
		localCost += my getLocalCost (iframe, my indices [jcand] [itrack], itrack, my closure);
	return localCost;
}
static double getTransitionCost_n (long iframe, long jcand1, long jcand2, void *closure) {
	struct parm2 *me = closure;
	int itrack;
	double transitionCost = 0.0;
	for (itrack = 1; itrack <= my ntrack; itrack ++)
		transitionCost += my getTransitionCost (iframe,
			my indices [jcand1] [itrack], my indices [jcand2] [itrack], itrack, my closure);
	return transitionCost;
}
static void putResult_n (long iframe, long jplace, void *closure) {
	struct parm2 *me = closure;
	int itrack;
	for (itrack = 1; itrack <= my ntrack; itrack ++)
		my putResult (iframe, my indices [jplace] [itrack], itrack, my closure);
}

int NUM_viterbi_multi (
	long nframe, long ncand, int ntrack,
	double (*getLocalCost) (long iframe, long icand, int itrack, void *closure),
	double (*getTransitionCost) (long iframe, long icand1, long icand2, int itrack, void *closure),
	void (*putResult) (long iframe, long place, int itrack, void *closure),
	void *closure)
{
	double ncomb;
	struct parm2 parm;
	int itrack, jtrack;
	long *icand = NULL, jcomb;
	parm.indices = NULL;

	if (ntrack > ncand) return Melder_error ("(NUMviterbin:) "
		"Number of tracks (%d) should not exceed number of candidates (%ld).", ntrack, ncand);
	ncomb = NUMcombinations (ncand, ntrack);
	if (ncomb > 10000000) return Melder_error ("(NUMviterbin:) "
		"Unrealistically high number of combinations (%ld).", ncomb);
	parm. ntrack = ntrack;
	parm. ncomb = ncomb;

	/*
	 * For ncand == 5 and ntrack == 3, parm.indices is going to contain:
	 *   1 2 3
	 *   1 2 4
	 *   1 2 5
	 *   1 3 4
	 *   1 3 5
	 *   1 4 5
	 *   2 3 4
	 *   2 3 5
	 *   2 4 5
	 *   3 4 5
	 */
	if (! (parm. indices = NUMlmatrix (1, parm. ncomb, 1, ntrack)) ||
	    ! (icand = NUMlvector (1, ntrack))) goto end;
	for (itrack = 1; itrack <= ntrack; itrack ++)
		icand [itrack] = itrack;   /* Start out with "1 2 3". */
	jcomb = 0;
	for (;;) {
		jcomb ++;
		for (itrack = 1; itrack <= ntrack; itrack ++)
			parm. indices [jcomb] [itrack] = icand [itrack];
		for (itrack = ntrack; itrack >= 1; itrack --) {
			if (++ icand [itrack] <= ncand - (ntrack - itrack)) {
				for (jtrack = itrack + 1; jtrack <= ntrack; jtrack ++)
					icand [jtrack] = icand [itrack] + jtrack - itrack;
				break;
			}
		}
		if (itrack == 0) break;
	}
	Melder_assert (jcomb == ncomb);
	parm. getLocalCost = getLocalCost;
	parm. getTransitionCost = getTransitionCost;
	parm. putResult = putResult;
	parm. closure = closure;
	if (! NUM_viterbi (nframe, ncomb, getNumberOfCandidates_n, getLocalCost_n, getTransitionCost_n,
		putResult_n, & parm)) goto end;
end:
	NUMlmatrix_free (parm. indices, 1, 1);
	NUMlvector_free (icand, 1);
	if (Melder_hasError ()) return 0;
	return 1;
}

int NUMrotationsPointInPolygon (double x0, double y0, long n, float x [], float y []) {
	long nup = 0, i;
	int upold = y [n] > y0, upnew;
	for (i = 1; i <= n; i ++) if ((upnew = y [i] > y0) != upold) {
		long j = i == 1 ? n : i - 1;
		if (x0 < x [i] + (x [j] - x [i]) * (y0 - y [i]) / (y [j] - y [i]))
			if (upnew) nup ++; else nup --;
		upold = upnew;
	}
	return nup;
}

/* End of file NUM.c */