/*
	Routines that need to be GPL-ed. The objective is to make this file empty.
	djmw 20020819
*/

#ifndef _NUM_h_
	#include "NUM.h"
#endif
#include "NUM2.h"
#include "SVD.h"
#include "NUMmachar.h"
#include "melder.h"

#define my me ->
#define SQR(a) ((temp=(a)) == 0.0 ? 0.0 : temp*temp)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
#define SIGN(x,s) ((s) < 0 ? -fabs (x) : fabs(x))

extern machar_Table NUMfpp;

/* only needed locally, MACRO_pythag is gpl */
#define MACRO_pythag(TYPE)\
{\
	TYPE absa, absb, temp;\
	absa = fabs(a);\
	absb = fabs(b);\
	if (absa > absb)\
	{\
		return absa * sqrt (1.0 + SQR (absb / absa));\
	}\
	else\
	{\
		return (absb == 0.0 ? 0.0 : absb * sqrt (1.0 + SQR (absa / absb)));\
	}\
}

static double pythag_d (double a, double b)
MACRO_pythag(double)

static float pythag_f (float a, float b)
MACRO_pythag(float)

#undef MACRO_pythag
#undef SQR

/* This routine is a corrected version 2.06 by djmw */
/*************************** sorting ************************************/

#define MACRO_NUMsort2(TYPE1, TYPE2)									\
{                                                                       \
	long l, j, ir, i;                                                   \
	TYPE1 rra; TYPE2 rrb;                                               \
	if (n < 2) return;                                                  \
	l = (n >> 1) + 1;                                                   \
	ir = n;                                                             \
	for (;;)                                                            \
	{                                                                   \
		if (l > 1)														\
		{                                                        		\
			rra = ra[--l]; rrb = rb[l];	                                \
		}																\
		else                                                            \
		{                                                               \
			rra = ra[ir]; rrb = rb[ir];                                 \
			ra [ir] = ra[1]; rb [ir] = rb[1];                           \
			if (--ir == 1)												\
			{															\
				ra[1] = rra;  rb[1] = rrb; return;						\
			}															\
		}                                                               \
		i = l;                                                          \
		j = l << 1; /* djmw  l + 1 */ 			                        \
		while (j <= ir)                                                 \
		{                                                               \
			if (j < ir && ra [j] < ra [j + 1]) j++;                     \
			if (rra < ra [j])                                           \
			{															\
				ra[i] = ra[j]; rb[i] = rb[j]; 							\
				j += (i = j);	/* djmw  i = j; j <<= 1; */				\
			}															\
			else j = ir + 1;                                            \
		}                                                               \
		ra[i] = rra; rb[i] = rrb;                                       \
	}                                                                   \
}
/*
void NUMsort2_fl (long n, float ra[], long rb[])
	MACRO_NUMsort2 (float, long)

void NUMsort2_ff (long n, float ra[], float rb[])
	MACRO_NUMsort2 (float, float)

void NUMsort2_dl (long n, double ra[], long rb[])
	MACRO_NUMsort2 (double, long)
	
void NUMsort2_dd (long n, double ra[], double rb[])
	MACRO_NUMsort2 (double, double)
	
void NUMsort2_ll (long n, long ra[], long rb[])
	MACRO_NUMsort2 (long, long)
*/	
#undef MACRO_NUMsort2


#define MACRO_NUMrank(TYPE)	\
{	\
	TYPE rank;	\
	long i, jt, j = 1;	\
	\
	while (j < n) {	\
		for (jt = j + 1; jt <= n && v[jt] == v[j]; jt++) {}	\
		rank = (j + jt - 1) * 0.5;	\
		for (i = j; i <= jt-1; i++) {v[i] = rank;}	\
		j = jt;	\
	} \
	if (j == n) v[n] = n;	\
}
/*
void NUMrank_f (long n, float v[])
	MACRO_NUMrank(float)
	
void NUMrank_d (long n, double v[])
	MACRO_NUMrank(double)
*/
#undef MACRO_NUMrank

#define MACRO_NUMrank(TYPE,TS)			\
{	\
	long i, j, nr = re - rb + 1, *index = NULL;	\
	TYPE *v = NULL;	\
	\
	v = NUM##TS##vector (1, nr);	\
	if (v == NULL) return 0;	\
	index = NUMlvector (1, nr);	\
	if (index == NULL) goto end;	\
\
	for (j = cb; j <= ce; j++)	\
	{	\
		for (i = 1; i <= nr; i++) v[i] = m[rb + i - 1][j]; \
		for (i = 1; i <= nr; i++) index[i] = i;	\
		NUMsort2_##TS##l (nr, v, index);	\
		NUMrank_##TS (nr, v);	\
		for (i = 1; i <= nr; i++) m[rb + index[i] - 1][j] = v[i];	\
	}	\
end:	\
	NUM##TS##vector_free (v, 1);	\
	NUMlvector_free (index, 1);	\
	return ! Melder_hasError();	\
}
/*
int NUMrankColumns_f (float **m, long rb, long re, long cb, long ce)
	MACRO_NUMrank(float,f)

int NUMrankColumns_d (double **m, long rb, long re, long cb, long ce)
	MACRO_NUMrank(double,d)
*/
#undef MACRO_NUMrank

#define M 7
#define NSTACK 50

#define MACRO_INDEXX(TYPE) \
{	\
	long i, indxt, ir = n, temp, j, k, l = 1; \
	long istack[51 /*NSTACK+1*/];	\
	int jstack = 0; 	\
	TYPE a;	\
	\
	if (n < 1) return;	\
	for (j = 1; j <= n; j++) indx[j] = j;	\
	if (n == 1) return;	\
	else if (n == 2)	\
	{	\
		if (arr[1] > arr[2]) { indx[1] = 2; indx[2] = 1; }	\
		return;	\
	}	\
	for (j = 1; j <= NSTACK; j++) istack[j] = 0; \
	for (;;)	\
	{	\
		if (ir-l < M)	\
		{	\
			for (j=l+1; j <= ir; j++)	\
			{	\
				indxt = indx[j]; a = arr[indxt];	\
				for (i=j-1; i >= l; i--)	\
				{	\
					if (arr[indx[i]] < a) break;	\
					indx[i+1] = indx[i];	\
				}	\
				indx[i+1] = indxt;	\
			}	\
			if (jstack == 0) break;	\
			ir = istack[jstack--];	\
			l = istack[jstack--];	\
		}	\
		else	\
		{	\
			k= (l + ir) >> 1;	\
			SWAP (indx[k], indx[l+1]);	\
			if (arr[indx[l]] > arr[indx[ir]])	\
			{	\
				SWAP (indx[l],indx[ir])	\
			}	\
			if (arr[indx[l+1]] > arr[indx[ir]])	\
			{	\
				SWAP (indx[l+1], indx[ir])	\
			}	\
			if (arr[indx[l]] > arr[indx[l+1]])	\
			{	\
				SWAP (indx[l],indx[l+1])	\
			}	\
			i = l+1; j = ir; indxt = indx[l+1]; a = arr[indxt];	\
			for(;;)	\
			{	\
				do i++; while (arr[indx[i]] < a);	\
				do j--; while (arr[indx[j]] > a);	\
				if (j < i) break;	\
				SWAP (indx[i], indx[j])	\
			}	\
			indx[l+1] = indx[j]; indx[j] = indxt; jstack += 2;	\
			if (jstack > NSTACK) return; \
			if (ir-i+1 >= j-l)	\
			{	\
				istack[jstack] = ir; istack[jstack-1] = i; ir = j-1;	\
			}	\
			else	\
			{	\
				istack[jstack] = j-1; istack[jstack-1] = l; l = i;	\
			}	\
		}	\
	}	\
	return;	\
}
/*
void NUMindexx_f (const float arr[], long n, long indx[])
	MACRO_INDEXX(float)

void NUMindexx_d (const double arr[], long n, long indx[])
	MACRO_INDEXX(double)
*/
#undef MACRO_INDEXX

/*
int NUMindexx_s (char *arr[], long n, long indx[])
{
	long i, indxt, ir = n, temp, j, k, l = 1, *istack;
	int jstack = 0; char *a;

	if (n < 1) return 0;
	for (j=1; j <= n; j++) indx[j] = j;
	if (n == 1) return 1;
	else if (n == 2)
	{
		if (NUMstrcmp (arr[1], arr[2]) > 0) { indx[1] = 2; indx[2] = 1; }
		return 1;
	}
	if ((istack = NUMlvector (1, NSTACK)) == NULL) return 0;
	for (;;)
	{
		if (ir-l < M)
		{
			for (j=l+1; j <= ir; j++)
			{
				indxt = indx[j]; a = (char *) arr[indxt];
				for (i=j-1; i >= l; i--)
				{
					if (NUMstrcmp (arr[indx[i]], a) < 0) break;
					indx[i+1] = indx[i];
				}
				indx[i+1] = indxt;
			}
			if (jstack == 0) break;
			ir = istack[jstack--];
			l = istack[jstack--];
		}
		else
		{
			k= (l + ir) >> 1;
			SWAP (indx[k], indx[l+1]);
			if (NUMstrcmp (arr[indx[l]], arr[indx[ir]]) > 0)
			{
				SWAP (indx[l],indx[ir])
			}
			if (NUMstrcmp (arr[indx[l+1]], arr[indx[ir]]) > 0)
			{
				SWAP (indx[l+1], indx[ir])
			}
			if (NUMstrcmp (arr[indx[l]], arr[indx[l+1]]) > 0)
			{
				SWAP (indx[l],indx[l+1])
			}
			i = l+1; j = ir; indxt = indx[l+1]; a = (char *) arr[indxt];
			for(;;)
			{
				do i++; while (NUMstrcmp (arr[indx[i]], a) < 0);
				do j--; while (NUMstrcmp (arr[indx[j]], a) > 0);
				if (j < i) break;
				SWAP (indx[i], indx[j])
			}
			indx[l+1] = indx[j]; indx[j] = indxt; jstack += 2;
			if (jstack > NSTACK) { NUMlvector_free (istack, 1); return 0; }
			if (ir-i+1 >= j-l)
			{
				istack[jstack] = ir; istack[jstack-1] = i; ir = j-1;
			}
			else
			{
				istack[jstack] = j-1; istack[jstack-1] = l; l = i;
			}
		}
	}
	NUMlvector_free (istack, 1);
	return 1;
}

#undef M
#undef NSTACK
*/

/*
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);   \
        a[k][l]=h+s*(g-h*tau);

int NUMjacobi (float **a, long n, float d[], float **v, long *nrot)
{
    long j, iq, ip, i;
	int status = 0;
    float tresh, theta, tau, t, sm, s, h, g, c, *b, *z;
	
    if (((b = NUMfvector (1, n)) == NULL) ||
		((z = NUMfvector (1, n)) == NULL)) goto end;
		 
    for (ip = 1; ip <= n; ip++)
    {
        for (iq = 1; iq <= n; iq++) v[ip][iq] = 0.0;
        v[ip][ip] = 1.0;
    }
    for (ip = 1; ip <= n; ip++) b[ip] = d[ip] = a[ip][ip];
    *nrot = 0;
    for (i = 1; i <= 50; i++)
    {
        sm=0.0;
        for (ip=1; ip <= n-1; ip++)
		{
			for (iq=ip+1; iq <= n; iq++) sm += fabs (a[ip][iq]);
		}
        if (sm == 0.0)
		{
			eigenSort_f (d, v, n); status = 1; goto end;
		}
        if (i < 4) tresh = 0.2 * sm / ( n * n);
        else tresh = 0.0;
        for (ip=1; ip <= n-1; ip++)
        {
            for (iq=ip+1; iq <= n; iq++)
            {
                g = 100.0 * fabs (a[ip][iq]);
                if (i > 4 && ( fabs (d[ip]) + g) == fabs (d[ip])
                    && (fabs (d[iq]) + g) == fabs (d[iq])) a[ip][iq] = 0.0;
                else if (fabs (a[ip][iq]) > tresh)
                {
                    h = d[iq] - d[ip];
                    if (fabs (h) + g == fabs (h)) t = a[ip][iq] / h;
                    else
                    {
                        theta = 0.5 * h / a[ip][iq];
                        t = 1.0 / ( fabs (theta) + sqrt (1.0 + theta * theta));
                        if (theta < 0.0) t = -t;
                    }
                    c = 1.0 / sqrt (1 + t * t); s = t * c;
                    tau = s / ( 1.0 + c); h = t * a[ip][iq];
                    z[ip] -= h; z[iq] += h;
                    d[ip] -= h; d[iq] += h;
                    a[ip][iq] = 0.0;
                    for (j=1; j <= ip-1; j++)
                    {
                        ROTATE (a, j, ip, j, iq)
                    }
                    for (j=ip+1; j <= iq-1; j++)
                    {
                        ROTATE (a, ip, j, j, iq)
                    }
                    for (j=iq+1; j <= n; j++)
                    {
                        ROTATE (a, ip, j, iq, j)
                    }
                    for (j=1; j <= n; j++)
                    {
                        ROTATE (v, j, ip, j, iq)
                    }
                    ++(*nrot);
                }
            }
        }
        for (ip=1; ip <= n; ip++)
		{
			b[ip] += z[ip]; d[ip] = b[ip]; z[ip] = 0.0;
		}
    }
    *nrot *= -1;
    (void) Melder_error ("NUMjacobi: too many rotations");
end:
    NUMfvector_free (b, 1);
    NUMfvector_free (z, 1);
    return status;
}

#undef ROTATE
*/

/* Voor eigen van symmetrische nxn matrix;
NUMtred2 gevolgd door NUMtqli
Transition_eigen, Matrix_eigen, Eigen_initFromSymmetricMatrix, SSCP_to_PCA,
NUMeigensystem
*/
void NUMtred2_f (float **a, long n, float d[], float e[])
{
	long l, k, j, i;
	float scale, hh, h, g, f;

	for (i=n; i >= 2; i--)
	{
		l = i - 1; h = scale = 0.0;
		if (l > 1)
		{
			for (k=1; k <= l; k++) scale += fabs (a[i][k]);
			if (scale == 0.0) e[i] = a[i][l];
			else
			{
				for (k=1; k <= l; k++)
				{
					a[i][k] /= scale;
					h += a[i][k] * a[i][k];
				}
				f = a[i][l]; g = ( f >= 0.0 ? -sqrt(h) : sqrt(h));
				e[i] = scale * g; h -= f * g; a[i][l] = f - g; f = 0.0;
				for (j=1;j<=l;j++)
				{
					a[j][i] = a[i][j] / h;
					g = 0.0;
					for (k=1; k <= j; k++) g += a[j][k] * a[i][k];
					for (k=j+1; k <= l; k++) g += a[k][j] * a[i][k];
					e[j] = g/h;
					f += e[j] * a[i][j];
				}
				hh=f/(h+h);
				for (j=1; j <= l; j++)
				{
					f = a[i][j]; e[j] = g = e[j] - hh * f;
					for (k=1; k <= j; k++) a[j][k] -= ( f * e[k] + g * a[i][k]);
				}
			}
		} else
			e[i] = a[i][l];
		d[i] = h;
	}
	d[1] = e[1] = 0.0;
	/* Contents of this loop can be omitted if eigenvectors not
			wanted except for statement d[i]=a[i][i]; */
	for (i=1; i <= n; i++)
	{
		l = i - 1;
		if (d[i])
		{
			for (j=1; j <= l; j++)
			{
				g = 0.0;
				for (k=1; k <= l; k++) g += a[i][k] * a[k][j];
				for (k=1; k <= l; k++) a[k][j] -= g * a[k][i];
			}
		}
		d[i] = a[i][i]; a[i][i] = 1.0;
		for (j=1; j <= l; j++) a[j][i] = a[i][j] = 0.0;
	}
}

void NUMtred2_d (double **a, long n, double d[], double e[])
{
	long l, k, j, i;
	double scale, hh, h, g, f;

	for (i=n; i >= 2; i--)
	{
		l = i - 1; h = scale = 0.0;
		if (l > 1)
		{
			for (k=1; k <= l; k++)
			{
				scale += fabs (a[i][k]);
			}
			if (scale == 0.0)
			{
				e[i] = a[i][l];
			}
			else
			{
				for (k=1; k <= l; k++)
				{
					a[i][k] /= scale;
					h += a[i][k] * a[i][k];
				}
				f = a[i][l]; g = ( f >= 0.0 ? -sqrt(h) : sqrt(h));
				e[i] = scale * g; h -= f * g; a[i][l] = f - g; f = 0.0;
				for (j=1;j<=l;j++)
				{
					a[j][i] = a[i][j] / h;
					g = 0.0;
					for (k=1; k <= j; k++) g += a[j][k] * a[i][k];
					for (k=j+1; k <= l; k++) g += a[k][j] * a[i][k];
					e[j] = g/h;
					f += e[j] * a[i][j];
				}
				hh=f/(h+h);
				for (j=1; j <= l; j++)
				{
					f = a[i][j]; e[j] = g = e[j] - hh * f;
					for (k=1; k <= j; k++) a[j][k] -= ( f * e[k] + g * a[i][k]);
				}
			}
		} else
			e[i] = a[i][l];
		d[i] = h;
	}
	d[1] = e[1] = 0.0;
	/* Contents of this loop can be omitted if eigenvectors not
			wanted except for statement d[i]=a[i][i]; */
	for (i=1; i <= n; i++)
	{
		l = i - 1;
		if (d[i])
		{
			for (j=1; j <= l; j++)
			{
				g = 0.0;
				for (k=1; k <= l; k++) g += a[i][k] * a[k][j];
				for (k=1; k <= l; k++) a[k][j] -= g * a[k][i];
			}
		}
		d[i] = a[i][i]; a[i][i] = 1.0;
		for (j=1; j <= l; j++) a[j][i] = a[i][j] = 0.0;
	}
}

#define MACRO_NUMtqli(T,TD,TYPE,EPS) \
{\
	long m, l, iter, i, k; \
	TYPE s, r, p, g, f, dd, c, b; \
\
	if (NUMfpp == NULL) NUMmachar (); \
\
	for (i = 2; i <= n; i++) \
	{ \
		e[i-1] = e[i]; \
	} \
	e[n] = 0.0; \
	for (l = 1; l <= n; l++) \
	{\
		iter = 0;\
		do\
		{\
			for (m = l; m <= n - 1; m++)\
			{\
				dd = fabs (d[m]) + fabs (d[m+1]);\
				/* !!!! let op aanpassen fabs(e[m])/dd < eps */\
				if (((TYPE)(fabs (e[m]) / dd)) < EPS) break;\
			}\
			if (m != l)\
			{\
				if (iter++ == 30) return Melder_error\
					("NUMtqli: too many iterations. Maybe your matrix is singular?");\
				g = (d[l+1] - d[l]) / ( 2.0 * e[l]);\
				r = pythag_##T (g, 1.0);\
				g = d[m] - d[l] + e[l] / (g + SIGN (r, g));\
				s = c = 1.0;\
				p = 0.0;\
				for (i = m - 1; i >= l; i--)\
				{\
					f = s * e[i];\
					b = c * e[i];\
					e[i+1] = (r = pythag_##T (f, g));\
					if (r == 0.0)\
					{\
						d[i+1] -= p;\
						e[m] = 0.0;\
						break;\
					}\
					s = f / r;\
					c = g / r;\
					g = d[i+1] - p;\
					r = (d[i] - g) * s + 2.0 * c * b;\
					d[i+1] = g + (p = s * r);\
					g = c * r - b;\
					if (z != NULL)\
					{\
						for (k = 1; k <= n; k++)\
						{\
							f = z[k][i+1];\
							z[k][i+1] = s * z[k][i] + c * f;\
							z[k][i] = c * z[k][i] - s * f;\
						}\
					}\
				}\
				if (r == 0.0 && i >= l) continue;\
				d[l] -= p;\
				e[l] = g;\
				e[m] = 0.0;\
			}\
		} while (m != l);\
	}\
	return 1;\
}

int NUMtqli_f (float d[], float e[], long n, float **z)
MACRO_NUMtqli(f,_f,float,5.96e-8)

int NUMtqli_d (double d[], double e[], long n, double **z)
MACRO_NUMtqli(d,_d,double,NUMfpp->eps)

#undef MACRO_NUMtqli

int NUMgaussj_d (double **a, long n, double **b, long m)
{
	int status = 0; long *indxc = NULL, *indxr = NULL, *ipiv = NULL;
	long i, icol, irow, j, k, l, ll;
	double big,dum,pivinv,temp;

	if (((indxc = NUMlvector (1, n)) == NULL) ||
		((indxr = NUMlvector (1, n)) == NULL) ||
		((ipiv = NUMlvector (1, n)) == NULL)) goto end;
	for (j = 1; j <= n; j++) ipiv[j]=0;
	for (i = 1; i <= n; i++) {
		big = 0.0;
		for (j = 1; j <= n; j++)
			if (ipiv[j] != 1)
				for (k = 1; k <= n; k++) {
					if (ipiv[k] == 0) {
						if (fabs(a[j][k]) >= big) {
							big=fabs(a[j][k]);
							irow=j;
							icol=k;
						}
					}
					else if (ipiv[k] > 1)
					{
						Melder_error("NUMgaussj: Singular Matrix-1");
						goto end;
					}
				}
		++(ipiv[icol]);
		if (irow != icol) {
			for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
			for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
		}
		indxr[i]=irow;
		indxc[i]=icol;
		if (a[icol][icol] == 0.0)
		{
			Melder_error("NUMgaussj: Singular Matrix-2"); goto end;
		}
		pivinv=1.0/a[icol][icol];
		a[icol][icol]=1.0;
		for (l=1;l<=n;l++) a[icol][l] *= pivinv;
		for (l=1;l<=m;l++) b[icol][l] *= pivinv;
		for (ll=1;ll<=n;ll++)
			if (ll != icol) {
				dum=a[ll][icol];
				a[ll][icol]=0.0;
				for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
				for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
			}
	}
	for (l=n;l>=1;l--) {
		if (indxr[l] != indxc[l])
			for (k=1;k<=n;k++)
				SWAP(a[k][indxr[l]],a[k][indxc[l]]);
	}
	status = 1;
end:
	NUMlvector_free (ipiv, 1);
	NUMlvector_free (indxr, 1);
	NUMlvector_free (indxc, 1);
	return status;
}

int NUMgaussj (float **a, long n)
{
	long *indxc = NULL, *indxr = NULL, *ipiv = NULL;
	long i, icol, irow, j, k, l, ll, status = 0;
	float big, dum, pivinv, temp;

	if (((indxc = NUMlvector (1, n)) == NULL) ||
		((indxr = NUMlvector (1, n)) == NULL) ||
		((ipiv = NUMlvector (1, n)) == NULL)) goto end;
	for (j=1; j <= n; j++) ipiv[j] = 0;
	for (i=1; i <= n; i++)
	{
		big = 0.0;
		for (j=1; j <= n; j++)
			if (ipiv[j] != 1)
				for (k=1; k <= n; k++)
				{
					if (ipiv[k] == 0)
					{
						if (fabs(a[j][k]) >= big)
						{
							big = fabs(a[j][k]); irow = j; icol = k;
						}
					}
					else if (ipiv[k] > 1)
					{
						Melder_error("NUMgaussj: Singular Matrix-1");
						goto end;
					}
				}
		++(ipiv[icol]);
		if (irow != icol) for (l=1; l <= n; l++) SWAP(a[irow][l],a[icol][l])
		indxr[i] = irow; indxc[i] = icol;
		if (a[icol][icol] == 0.0)
		{
			Melder_error("NUMgaussj: Singular Matrix-2"); goto end;
		}
		pivinv = 1.0 / a[icol][icol]; a[icol][icol] = 1.0;
		for (l=1; l <= n; l++) a[icol][l] *= pivinv;
		for (ll=1; ll <= n; ll++)
			if (ll != icol)
			{
				dum = a[ll][icol]; a[ll][icol] = 0.0;
				for (l=1; l <= n; l++) a[ll][l] -= a[icol][l]*dum;
			}
	}
	for (l=n; l >= 1;l--)
	{
		if (indxr[l] != indxc[l])
		{
			for (k=1; k <= n; k++) SWAP(a[k][indxr[l]],a[k][indxc[l]]);
		}
	}
	status = 1;
end:
	NUMlvector_free (ipiv, 1);
	NUMlvector_free (indxr, 1);
	NUMlvector_free (indxc, 1);
	return status;
}

int NUMsvdcmp(double **a, long m, long n, double w[], double **v)
{
	long flag, i, its, j, jj, k, l, nm;
	double anorm, c, f, g, h, s, scale, x, y, z, *rv1 = NULL, tmp;

	if ((rv1 = NUMdvector (1, n)) == NULL) return 0;
	if (! NUMfpp) NUMmachar ();
	
	/*
		Special case matrix == 'column vector'
	*/
	g = scale = anorm = 0.0;
	
	if (n == 1)
	{
		for (k = 1; k <= m; k++)
		{
			anorm += a[k][1] * a[k][1];
		}
		anorm = sqrt (anorm);	
		if (anorm == 0)
		{
			NUMdvector_free (rv1, 1);
			return Melder_error ("NUMsvdcmp: matrix too singular (norm=0).");
		}
		for (k = 1; k <= m; k++)
		{
			a[k][1] /= anorm;
		}
		v[1][1] = 1;
		w[1] = anorm;
		return 1;
	}
	
	for (i = 1; i <= n; i++)
	{
		l = i + 1;
		rv1[i] = scale * g;
		g = s = scale = 0.0;
		if (i <= m)
		{
			for (k = i; k <= m; k++) scale += fabs (a[k][i]);
			if (scale)
			{
				for (k = i; k <= m; k++)
				{
					a[k][i] /= scale;
					s += a[k][i] * a[k][i];
				}
				f = a[i][i];
				g = -SIGN (sqrt(s), f);
				h = f * g - s;
				a[i][i] = f - g;
				for (j = l; j <= n; j++)
				{
					for (s = 0.0, k = i; k <= m; k++) s += a[k][i] * a[k][j];
					f = s / h;
					for (k = i; k <= m; k++) a[k][j] += f * a[k][i];
				}
				for (k = i; k <= m; k++) a[k][i] *= scale;
			}
		}
		w[i] = scale * g;
		g = s = scale = 0.0;
		if (i <= m && i != n)
		{
			for (k = l; k <= n; k++) scale += fabs (a[i][k]);
			if (scale)
			{
				for (k = l; k <= n; k++)
				{
					a[i][k] /= scale;
					s += a[i][k] * a[i][k];
				}
				f = a[i][l];
				g = -SIGN (sqrt(s), f);
				h = f * g - s;
				a[i][l] = f - g;
				for (k = l; k <= n; k++) rv1[k] = a[i][k] / h;
				for (j = l; j <= m; j++)
				{
					for (s = 0.0, k = l; k <= n; k++) s += a[j][k] * a[i][k];
					for (k = l; k <= n; k++) a[j][k] += s * rv1[k];
				}
				for (k = l; k <= n; k++) a[i][k] *= scale;
			}
		}
		anorm = anorm > (tmp = (fabs (w[i]) + fabs (rv1[i]))) ? anorm : tmp;
	}
	
	if (anorm == 0)
	{
		NUMdvector_free (rv1, 1);
		return Melder_error ("NUMsvdcmp: matrix too singular (norm=0).");
	}
	
	for (i = n; i >= 1; i--)
	{
		if (i < n)
		{
			if (g)
			{
				for (j = l; j <= n; j++) v[j][i] = (a[i][j] / a[i][l]) / g;
				for (j = l; j <= n; j++)
				{
					for (s = 0.0, k = l; k <= n; k++) s += a[i][k] * v[k][j];
					for (k = l; k <= n; k++) v[k][j] += s * v[k][i];
				}
			}
			for (j = l; j <= n; j++) v[i][j] = v[j][i] = 0.0;
		}
		v[i][i] = 1.0;
		g = rv1[i];
		l=i;
	}
	for (i = m < n ? m : n; i >= 1; i--)
	{
		l = i + 1;
		g = w[i];
		for (j = l; j <= n; j++) a[i][j] = 0.0;
		if (g)
		{
			g = 1.0 / g;
			for (j = l; j <= n; j++)
			{
				for (s = 0.0, k = l; k <= m; k++) s += a[k][i] * a[k][j];
				f = ( s / a[i][i]) * g;
				for (k = i; k <= m; k++) a[k][j] += f * a[k][i];
			}
			for (j = i; j <= m; j++) a[j][i] *= g;
		}
		else
		{
			for (j = i; j <= m; j++) a[j][i] = 0.0;
		}
		++a[i][i];
	}
	for (k = n; k >= 1; k--)
	{
		for (its = 1; its <= 30; its++)
		{
			flag = 1;
			for (l = k; l >= 1; l--)
			{
				nm = l - 1;
/* The following tests ''fails'' on machines whose floating point 
   registers have more than 64 bits and when the comparants are always 
   in a register: 
   if ((fabs (rv1[l]) + anorm) == anorm) { flag=0; break; }
   if ((fabs (w[nm]) + anorm) == anorm) break;
   We have to make the following substitutions:
 */   
				if (fabs (rv1[l])/ anorm < NUMfpp -> eps)
				{
					flag=0; break;
				}
				if (nm == 0) 
				{
					NUMdvector_free (rv1, 1);
					return Melder_error ("NUMsvdcmp: very ill-conditioned.");
				}
				if (fabs (w[nm]) / anorm < NUMfpp -> eps) break;
			}
			if(flag)
			{
				c = 0.0;
				s = 1.0;
				for (i = l; i <= k; i++)
				{
					f = s * rv1[i];
					rv1[i] = c * rv1[i];
/*					if ((fabs (f) + anorm) == anorm) break; */
					if (fabs (f) / anorm < NUMfpp -> eps) break;
					g = w[i];
					h = pythag_d (f, g);
					w[i] = h;
					h = 1.0 / h;
					c = g * h;
					s = -f * h;
					for (j = 1; j <= m; j++)
					{
						y = a[j][nm];
						z = a[j][i];
						a[j][nm] = y * c + z * s;
						a[j][i]  = z * c - y * s;
					}
				}
			}
			z = w[k];
			if (l == k)
			{
				if (z < 0.0)
				{
					w[k] = -z; 
					for (j = 1; j <= n; j++) v[j][k] = -v[j][k];
				}
				break;
			}
			if (its == 30)
			{ 
				NUMdvector_free (rv1, 1);
				return Melder_error ("NUMsvdcmp: no convergence in "
					"30 iterations");
			}
			x = w[l];
			nm = k - 1;
			y = w[nm];
			g = rv1[nm];
			h = rv1[k];
			f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
			g = pythag_d (f, 1.0);
			f = ((x - z) * (x + z) + h * ((y / (f + SIGN (g, f))) - h)) / x;
			c = s = 1.0;
			for (j = l; j <= nm; j++)
			{
				i = j + 1;
				g = rv1[i];
				y = w[i];
				h = s * g;
				g = c * g;
				z = pythag_d (f, h);
				rv1[j] = z;
				c = f / z;
				s = h / z;
				f = x * c + g * s;
				g = g * c - x * s;
				h = y * s; 
				y *= c;
				for (jj = 1; jj <= n; jj++)
				{
					x = v[jj][j];
					z = v[jj][i];
					v[jj][j] = x * c + z * s;
					v[jj][i] = z * c - x * s;
				}
				z = pythag_d (f, h);
				w[j] = z;
				if (z)
				{
					z = 1.0 / z;
					c = f * z;
					s = h * z;
				}
				f = c * g + s * y;
				x = c * y - s * g;
				for (jj = 1; jj <= m; jj++)
				{
					y = a[jj][j];
					z = a[jj][i];
					a[jj][j] = y * c + z * s;
					a[jj][i] = z * c - y * s;
				}
			}
			rv1[l] = 0.0;
			rv1[k] = f;
			w[k] = x;
		}
	}
	NUMdvector_free (rv1, 1);
	return 1;
}

int NUMsvbksb (double **u, double w[], double **v, long m, long n,
	double b[], double x[])
{
	long jj,j,i; double s, *tmp;

	if ((tmp = NUMdvector (1, n)) == NULL) return 0;
	for (j=1; j <= n; j++)
	{
		s = 0.0;
		if (w[j])
		{
			for (i=1; i <= m; i++) s += u[i][j] * b[i];
			s /= w[j];
		}
		tmp[j] = s;
	}
	for (j=1; j <= n; j++)
	{
		s = 0.0;
		for (jj=1; jj <= n; jj++) s += v[j][jj] * tmp[jj];
		x[j] = s;
	}
	NUMdvector_free (tmp, 1);
	return 1;
}

#define TINY 1.0e-20;
int NUMludcmp_f(float **a, long n, long *indx, float *d)
{
	long i,imax,j,k;
	float big,dum,sum,temp;
	float *vv;

	if ((vv = NUMfvector (1, n)) == NULL)  return 0;
	*d = 1.0;
	for (i=1; i <= n; i++)
	{
		big=0.0;
		for (j=1; j <= n; j++) if ((temp = fabs (a[i][j])) > big) big = temp;
		if (big == 0.0)
		{
			NUMfvector_free (vv, 1);
			return Melder_error ("NUMludcmp: Singular matrix.");
		}
		vv[i] = 1.0/big;
	}
	for (j=1; j <= n; j++)
	{
		for (i=1; i < j; i++)
		{
			sum = a[i][j];
			for (k=1; k < i; k++) sum -= a[i][k] * a[k][j];
			a[i][j] = sum;
		}
		big = 0.0;
		for (i=j; i <= n; i++)
		{
			sum = a[i][j];
			for (k=1; k < j; k++) sum -= a[i][k] * a[k][j];
			a[i][j] = sum;
			if ( (dum=vv[i]*fabs(sum)) >= big) { big=dum; imax=i; }
		}
		if (j != imax)
		{
			for (k=1; k <= n; k++)
			{
				dum = a[imax][k]; a[imax][k] = a[j][k]; a[j][k] = dum;
			}
			*d = -(*d);
			vv[imax] = vv[j];
		}
		indx[j] = imax;
		if (a[j][j] == 0.0) a[j][j] = TINY;
		if (j != n) {
			dum = 1.0 / a[j][j];
			for (i=j+1; i <= n; i++) a[i][j] *= dum;
		}
	}
	NUMfvector_free (vv, 1);
	return 1;
}

int NUMludcmp (double **a, long n, long *indx, double *d)
{
	long i, imax, j, k;
	double big, dum, sum, temp, *vv;

	if ((vv = NUMdvector (1, n)) == NULL)  return 0;
	*d = 1.0;
	for (i=1; i <= n; i++)
	{
		big=0.0;
		for (j=1; j <= n; j++)
		{
			if ((temp = fabs (a[i][j])) > big) big = temp;
		}
		if (big == 0.0)
		{
			NUMdvector_free (vv, 1);
			return Melder_error ("NUMludcmp: Singular matrix.");
		}
		vv[i] = 1.0 / big;
	}
	for (j=1; j <= n; j++)
	{
		for (i=1; i < j; i++)
		{
			sum = a[i][j];
			for (k=1; k < i; k++) sum -= a[i][k] * a[k][j];
			a[i][j] = sum;
		}
		big = 0.0;
		for (i=j; i <= n; i++)
		{
			sum = a[i][j];
			for (k=1; k < j; k++) sum -= a[i][k] * a[k][j];
			a[i][j] = sum;
			if ((dum = vv[i] * fabs(sum)) >= big)
			{
				big = dum; imax = i;
			}
		}
		if (j != imax)
		{
			for (k=1; k <= n; k++)
			{
				dum = a[imax][k]; a[imax][k] = a[j][k]; a[j][k] = dum;
			}
			*d = -(*d);
			vv[imax] = vv[j];
		}
		indx[j] = imax;
		if (a[j][j] == 0.0) a[j][j] = TINY;
		if (j != n)
		{
			dum = 1.0 / a[j][j];
			for (i=j+1; i <= n; i++) a[i][j] *= dum;
		}
	}
	NUMdvector_free (vv, 1);
	return 1;
}

/*
int NUMcholeskyDecomposition (double **a, long n, double d[])
{
	long i, j, k; 
	double sum;

	for (i = 1; i <= n; i++)
	{
		for (j = i; j <= n; j++)
		{
			sum = a[i][j];
			for (k = i - 1; k >= 1; k--)
			{
				sum -= a[i][k] * a[j][k];
			}
			if (i == j)
			{
				if (sum <= 0.0) return Melder_error
					("NUMcholeskyDecomposition: Matrix not positive definite.");
				d[i] = sqrt (sum);
			}
			else a[j][i] = sum / d[i];
		}
	}
	return 1;
}
*/
#define ITMAX 100

int NUMrtsafe (void (*funcd) (void *, double, double *, double *), void *data, 
	double x1, double x2, double xacc, double *rts)
{
	long j; double df, dx, dxold, f, fh, fl, temp, xh, xl;
	*rts = 0;
	(*funcd) (data, x1, &fl, &df);
	(*funcd) (data, x2, &fh, &df);
	if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
		return Melder_error("NUMrtsafe: root must be bracketed.");
	if (fl == 0.0) { *rts = x1; return 1; }
	if (fh == 0.0) { *rts = x2; return 1; }
	if (fl < 0.0)
	{
		xl = x1; xh = x2;
	}
	else
	{
		xh = x1; xl = x2;
	}
	*rts = 0.5 * (x1 + x2);
	dxold = fabs (x2 - x1);
	dx = dxold;
	(*funcd) (data, *rts, &f, &df);
	for (j=1; j <= ITMAX; j++)
	{
		if ((((*rts - xh) * df - f) * ((*rts - xl) * df - f) >= 0.0)
			|| (fabs (2.0 * f) > fabs (dxold * df)))
		{
			dxold = dx; dx = 0.5 * (xh - xl); *rts = xl + dx;
			if (xl == *rts) return 1;
		}
		else
		{
			dxold = dx; dx = f / df; temp = *rts; *rts -= dx;
			if (temp == *rts) return 1;
		}
		if (fabs(dx) < xacc) return 1;
		(*funcd) (data, *rts, &f, &df);
		if (f < 0.0) xl = *rts; else xh = *rts;
	}
	/* return the value found */
	return Melder_error ("NUMrtsafe: maximum number of iterations exceeded.");
}


#define EPS 3.0e-7
#define FPMIN 1.0e-30
/*
static void NUMgser(double *gamser, double a, double x, double *gln)
{
	long n;
	double sum,del,ap;

	*gln = NUMlnGamma(a);
	if (x == 0.0) { *gamser=0.0; return; }
	else
	{
		ap = a; del = sum = 1.0 / a;
		for (n=1; n <= ITMAX; n++)
		{
			++ap;
			del *= x / ap;
			sum += del;
			if (fabs(del) < fabs(sum)*EPS)
			{
				*gamser = sum * exp (-x + a*log(x) - (*gln));
				return;
			}
		}
		Melder_warning ("NUMgser: a too large, ITMAX too small.");
		return;
	}
}
*/
double NUMbetaContinuedFraction(double a, double b, double x)
{
	long m, m2;
	double aa, c, d, del, h, qab, qam, qap;

	qab = a + b; qap = a + 1.0; qam = a - 1.0;
	c = 1.0; d = 1.0 - qab * x / qap;
	if (fabs (d) < FPMIN) d = FPMIN;
	d = 1.0 / d; h = d;
	for (m=1; m <= ITMAX; m++)
	{
		m2 = 2 * m;
		aa = m * (b-m) * x / ( (qam+m2) * (a+m2));
		d = 1.0 + aa * d;
		if (fabs(d) < FPMIN) d = FPMIN;
		c = 1.0 + aa / c;
		if (fabs(c) < FPMIN) c = FPMIN;
		d = 1.0 / d; h *= d * c;
		aa = -(a+m) * (qab+m) * x / ( (a+m2) * (qap+m2));
		d = 1.0 + aa * d;
		if (fabs(d) < FPMIN) d = FPMIN;
		c = 1.0 + aa / c;
		if (fabs(c) < FPMIN) c=FPMIN;
		d = 1.0 / d; del = d * c; h *= del;
		if (fabs(del-1.0) < EPS) break;
	}
	if (m > ITMAX)
	{
		Melder_warning("NUMbetaContinuedFraction: a or b too big, or "
			"ITMAX too small.");
		return 0;
	}
	return h;
}

double NUMridders2dw (double (*func) (double, long), double constant,
	long df, double x1, double x2, double xacc)
{
	long iter, maxit = 60;
	double ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew;

	fl = (*func)(x1, df) - constant;
	fh = (*func)(x2, df) - 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, df) - 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, df) - 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("NUMridders2: never get here.");
				return NUMundefined;
			}
			if (fabs (xh - xl) <= xacc) return ans;
		}
		(void) Melder_error("NUMridders2: exceed maximum iterations.");
		return NUMundefined;
	}
	else
	{
		if (fl == 0.0) return x1;
		if (fh == 0.0) return x2;
		(void) Melder_error("NUMridders2: root must be bracketed.");
		return NUMundefined;
	}
	return 0.0;
}

double NUMridders3dw (double (*func) (double, long, long), double constant,
	long df1, long 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.");
		return NUMundefined;
	}
	else
	{
		if (fl == 0.0) return x1;
		if (fh == 0.0) return x2;
		(void) Melder_error("NUMridders3: root must be bracketed.");
		return NUMundefined;
	}
	return 0.0;
}

#undef EPS
#undef FPMIN

double NUMincompleteBeta (double a, double b, double x)
{
	double bt = 0;

	if (a <= 0 || b <= 0 || x < 0 || x > 1) return NUMundefined;
	
	if (x > 0 && x < 1) bt = exp (NUMlnGamma (a + b) - NUMlnGamma (a) -
		NUMlnGamma (b) + a * log (x) + b * log (1 - x));
	
	return (x < (a + 1) / (a + b + 2)) ? bt *
		NUMbetaContinuedFraction (a, b, x) / a :
		1 - bt * NUMbetaContinuedFraction (b, a, 1 - x) / b;
}

void NUMhunt_f (float xx[], long n, float x, long *jlo)
{
	long jm, jhi, inc;
	int ascnd;

	ascnd = (xx[n] >= xx[1]);
	if (*jlo <= 0 || *jlo > n) { *jlo = 0; jhi = n + 1; }
	else
	{
		inc = 1;
		if (x >= xx[*jlo] == ascnd)
		{
			if (*jlo == n) return;
			jhi = *jlo + 1;
			while (x >= xx[jhi] == ascnd)
			{
				*jlo = jhi;
				inc += inc;
				jhi = *jlo + inc;
				if (jhi > n) { jhi = n + 1; break; }
			}
		}
		else
		{
			if (*jlo == 1) { *jlo = 0; return; }
			jhi=(*jlo)--;
			while (x < xx[*jlo] == ascnd)
			{
				jhi = *jlo;
				inc *= 2;
				if (inc >= jhi) { *jlo = 0; break; }
				else *jlo = jhi - inc;
			}
		}
	}
	while (jhi - *jlo != 1)
	{
		jm = (jhi + *jlo) / 2;
		if (x >= xx[jm] == ascnd) *jlo = jm;
		else jhi = jm;
	}
	if (x == xx[n]) *jlo = n - 1;
	if (x == xx[1]) *jlo = 1;
}

void NUMhunt_d (double xx[], long n, double x, long *jlo)
{
	long jm, jhi, inc;
	int ascnd;

	ascnd = (xx[n] >= xx[1]);
	if (*jlo <= 0 || *jlo > n) { *jlo = 0; jhi = n + 1; }
	else
	{
		inc = 1;
		if (x >= xx[*jlo] == ascnd)
		{
			if (*jlo == n) return;
			jhi = *jlo + 1;
			while (x >= xx[jhi] == ascnd)
			{
				*jlo = jhi;
				inc += inc;
				jhi = *jlo + inc;
				if (jhi > n) { jhi = n + 1; break; }
			}
		}
		else
		{
			if (*jlo == 1) { *jlo = 0; return; }
			jhi=(*jlo)--;
			while (x < xx[*jlo] == ascnd)
			{
				jhi = *jlo;
				inc *= 2;
				if (inc >= jhi) { *jlo = 0; break; }
				else *jlo = jhi - inc;
			}
		}
	}
	while (jhi - *jlo != 1)
	{
		jm = (jhi + *jlo) / 2;
		if (x >= xx[jm] == ascnd) *jlo = jm;
		else jhi = jm;
	}
	if (x == xx[n]) *jlo = n - 1;
	if (x == xx[1]) *jlo = 1;
}

#undef SWAP
#undef SIGN

/* End of file NUM_nongpl.c */