/*	Copyright (c) 1997 Lucent Technologies	*/
/*	  All Rights Reserved 	*/

/*	THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF LUCENT.	*/
/*	The copyright notice above does not evidence any   	*/
/*	actual or intended publication of such source code.	*/

/* sample funcadd */

#include "stdlib.h"	/* for atoi */
#include "math.h"
#include "funcadd.h"	/* includes "stdio1.h" */
#include "myalloc.h"

/* Sample function kth optionally prints its arguments using
 * a variant of stdio.h supplied by funcadd.h. */

#ifdef __cplusplus
/* The #ifdef __cplusplus lines allow this to be compiled
 * either by an ANSI C compiler or by a C++ compiler.
 */
extern "C" {
#endif

static void showvec( real *x, int n, AmplExports *ae);

static void showmat( real **A, int n, AmplExports *ae);

static void showmat2( real **A, int n, AmplExports *ae);

static void mx( real **A, real *x, real *y, int n);

static double dotprod( double *x, double *y, int n);

static int argmax( real *x, int n);

static real erf (real x);
static real gammp (real a, real x);
static void gser(real *gamser, real a, real x, real *gln);
static void gcf( real *gammcf, real a, real x, real *gln);

 static real
ginv(arglist *al)	/* generalized inverse of a single argument */
{
	real x = al->ra[0];
	x = x ? 1./x : 0.;
	if (al->derivs) {
		*al->derivs = -x*x;
		if (al->hes)
			*al->hes = -2.*x * *al->derivs;
		}
	return x;
	}

 static char *
sginv(arglist *al)	/* character-valued version of ginv */
{
	static char buf[32];
	AmplExports *ae = al->AE;	/* for sprintf */
	real x = al->ra[0];
	sprintf(buf, "x%.g", x ? 1./x : 0);
	return buf;
	}

 static real
myhypot(arglist *al)	/* myhypot(x,y) = sqrt(x*x + y*y) */
{
	real *d, *h, rv, x, x0, y, y0;

	x = x0 = al->ra[0];
	y = y0 = al->ra[1];

	if (x < 0.)
		x = -x;
	if (y < 0.)
		y = -y;
	rv = x;
	if (x < y) {
		rv = y;
		y = x;
		x = rv;
		}
	if (rv) {
		y /= x;
		rv = x * sqrt(1. + y*y);
		if (d = (real *)al->derivs) {
			d[0] = x0 / rv;
			d[1] = y0 / rv;
			if (h = al->hes) {
				h[0] =  d[1]*d[1] / rv;
				h[1] = -d[0]*d[1] / rv;
				h[2] =  d[0]*d[0] / rv;
				}
			}
		}
	else if (d = (real *)al->derivs) {
		d[0] = d[1] = 0;
		if (h = al->hes)
			h[0] = h[1] = h[2] = 0;
		}
	return rv;
	}

 static real
ncall(void)		/* returns it's invocation count */
{ static real x; return ++x; }

 static real
mean(register arglist *al)	/* mean of arbitrarily many arguments */
{
	real x, z;
	real *d, *de, *ra;
	int *at, i, j, n, nreal;
	char *se, *sym;
	AmplExports *ae = al->AE;

	n = al->n;
	if (n <= 0) return 0;
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	x = 0.0;
	nreal = 0;
	for(i=0; i<n; i++) {
		j = at[i];
		if (j >= 0) {
			x += ra[j];
			++nreal;
		} else {
			x += z = strtod(sym = al->sa[-(j+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"mean treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
	}
	if (d) {
		z = 1.0/n;
		for (i=0; i<nreal; i++) {d[i] = z;}
		/* The Hessian is == 0 and, if needed, has been */
		/* initialized to 0. */
	}
	return x/n;
}

 static real
mineig(register arglist *al)	/* minimum eigenvalue of a symmetric matrix */
{
	real z, lambdak, normx, normy, shift=0, cosine, ydotEkk, 
		lambda0=HUGE_VAL, hm, num;
	real *e0, *ek, *d, *h, *ra;
	int *at, i, ii, j, jj, k, kk, m, n, narg, cnt=0, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, **E=NULL;
	static real *lambda=NULL, *x=NULL, *y=NULL, *y0=NULL;

	/*
	srand48(0);
	*/

	narg = al->n;
	if (narg <= 0) return 0;
	n = (int)(sqrt(1+8*narg)-1)/2;
	if (n*(n+1)/2 != narg) {
	     fprintf(Stderr, 
	     "should be mineig({i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	/*
	if (A == NULL) {
	    CALLOC( A, n, real * );
	} else {
	    REALLOC( A, n, real * );
	}
	for (i=0; i<n; i++) REALLOC( A[i], n-i, real );
	*/
	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	/*
	REALLOC( A[0], n*(n+1)/2, real);
	for (i=0, k=0; i<n; k+=n-i, i++) {
	    A[i] = A[0]+k;
	}
	*/
	REALLOC( A[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A[i] = A[0]+k;

	if (E == NULL) {
	    REALLOC( E, n, real *);
	    E[0] = NULL;
	} else {
	    REALLOC( E, n, real *);
	}
	REALLOC( E[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) E[i] = E[0]+k;

	REALLOC( lambda, n, real );

	REALLOC( x, n, real );

	REALLOC( y, n, real );

	/*
	REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = drand48(); 
	*/
	REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = 0.5;

	k = 0;
	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		jj = at[k];
fprintf(Stderr, "i,j,k,jj: %3d %3d %3d %3d\n",i,j,k,jj); fflush(Stderr);
		if (jj >= 0) {
			A[i][j] = ra[jj];
		} else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"mineig treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
		k++;
	    }
	}
	/* showmat(A,n,ae); */

	for (i=0; i<n; i++) y[i]  = y0[i];
	do {
	    normy = sqrt(dotprod(y,y,n));
	    for (i=0; i<n; i++) { y[i] /= normy; }
	    cosine = dotprod(x,y,n);
	    for (i=0; i<n; i++) { x[i] = y[i]; }
	    mx( A, x, y, n);
	    i = argmax(y,n);
	    lambda[0] = y[i]/x[i];
	    cnt++;
	} while ( fabs(cosine) < 0.99999999999 );

	/*
	fprintf(Stderr, "lambda = %10.4f \n", lambda[0]);
	if (lambda[0] >= 0) {shift = lambda[0]/2 + y[0];} else {shift = y[0];}
	*/
	shift = y0[0];
	fprintf(Stderr, "shift = %f \n", shift);

	for (i=0; i<n; i++) A[i][i] -= shift;

	for (k=0; k<n; k++) {
	    for (i=0; i<n; i++) y[i] = y0[i];
	    /*
	    for (kk=0; kk<k; kk++) {
		ydotEkk = dotprod(y, E[kk], n);
		for (i=0; i<n; i++) y[i] -= ydotEkk*E[kk][i];
	    }
	    */
	    do {
		normy = sqrt(dotprod(y,y,n));
	        for (i=0; i<n; i++) { y[i] /= normy; }
	        /* showvec(y,n,ae); */
		cosine = dotprod(x,y,n);
	        for (i=0; i<n; i++) { x[i] = y[i]; }
		mx( A, x, y, n);
	        /* showvec(y,n,ae); */
		i = argmax(y,n);
	        lambdak = y[i]/x[i];
		cnt++;
	    } while ( fabs(cosine) < 0.9999999999 );

	    lambda[k] = lambdak + shift;
	    for (i=0; i<n; i++) E[k][i] = x[i];

	    if (lambda[k] < lambda0) {
		i0 = k;
		lambda0 = lambda[k];
	    }

	    fprintf(Stderr, "%9.4f %9.4f: ", lambda[k], lambdak);
	    for (i=0; i<n; i++) { fprintf(Stderr, "%7.4f ", E[k][i]); }
	    fprintf(Stderr, "\n");

	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    A[i][j] -= lambdak*x[i]*x[j];
		}
	    }

	    /* showmat(A,n,ae); */
	}
	fprintf(Stderr, "\n");


	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		A[i][j] = 0;
		for (k=0; k<n; k++) {
		    A[i][j] += lambda[k]*E[k][i]*E[k][j];
		}
	    }
	}
	showmat(A,n,ae);

	/*
	for (i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		fprintf(Stderr, "%2d %2d %10.4f \n", i, j,
		dotprod(E[i],E[j],n));
	    }
	}
	*/

	lambda0 = lambda[i0];
	e0 = E[i0];
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        d[m] = 2*e0[i]*e0[j];
			if (i==j) d[m] /= 2;
			m++;
		    }
		}
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    for (ii=i; ii<n; ii++) {
			for (jj=(ii==i)?j:ii; jj<n; jj++) {
			    hm = 0;
			    for (k=0; k<n; k++) {
				ek = E[k];
				if (k != i0) {
				    num = (e0[ i]*ek[ j] + e0[ j]*ek[ i])*
				          (e0[ii]*ek[jj] + e0[jj]*ek[ii]);
				    if ( i== j) num /= 2;
				    if (ii==jj) num /= 2;
				    hm += num/(lambda0-lambda[k]);
				}
			    }
			    h[m] = 2*hm;
			    m++;
			}
		    }
		}
	    }
	}
	return lambda0;
}

 static real
linear(arglist *al)	/* generalized inverse of a single argument */
{
	real x = al->ra[0];
	if (al->derivs) {
		*al->derivs = 1;
		if (al->hes)
			*al->hes = 0;
		}
	return x;
}

static real R_real, R_imag;
static real lx, ly, lz, qx, qy, qz, zmin, zmax;

static void
R(arglist *al) /* inverse of 3D spectral density on a spherical shell */
{
	int j, n=1000; 
/*	int j, n=20; */

	real pi = 4*atan(1);
	real lx, ly, lz, qx, qy, qz, zmin, zmax, z;
	real z_perp, 
	     lq, 	     /* component of l in direction of q */
	     lox, loy, loz,  /* component of l orthogonal to q */
	     l_mag, q_mag;

	lx = al->ra[0];
	ly = al->ra[1];
	lz = al->ra[2];
	qx = al->ra[3];
	qy = al->ra[4];
	qz = al->ra[5];
	zmin = al->ra[6];
	zmax = al->ra[7];

	q_mag = sqrt(qx*qx + qy*qy + qz*qz);
	qx /= q_mag;
	qy /= q_mag;
	qz /= q_mag;

	R_real=0;
	R_imag=0;

	for (j=0; j<n; j++) {
	    z = zmin + (zmax-zmin)*(j+0.5)/n;
	    z_perp = sqrt(1-z*z);
	    lq = lx*qx + ly*qy + lz*qz;
	    lox = lx - lq*qx;
	    loy = ly - lq*qy;
	    loz = lz - lq*qz;
	    l_mag = sqrt( lox*lox + loy*loy + loz*loz );
	    R_real += j0(2*pi*z_perp*l_mag)*cos(-2*pi*z*lq);
	    R_imag += j0(2*pi*z_perp*l_mag)*sin(-2*pi*z*lq);
	}
	R_real *= 2*pi*(zmax-zmin)/n;
	R_imag *= 2*pi*(zmax-zmin)/n;

	/*
	if (al->derivs) {
		*al->derivs = 1;
		if (al->hes)
			*al->hes = 0;
		}
	*/
}

 static real
C(arglist *al) /* inverse of 3D spectral density on a spherical shell */
{
	if (lx != al->ra[0] ||
	    ly != al->ra[1] ||
	    lz != al->ra[2] ||
	    qx != al->ra[3] ||
	    qy != al->ra[4] ||
	    qz != al->ra[5] ||
	    zmin != al->ra[6] ||
	    zmax != al->ra[7] ) {
		R(al);
	    }

	return R_real;
}

 static real
S(arglist *al) /* inverse of 3D spectral density on a spherical shell */
{
	if (lx != al->ra[0] ||
	    ly != al->ra[1] ||
	    lz != al->ra[2] ||
	    qx != al->ra[3] ||
	    qy != al->ra[4] ||
	    qz != al->ra[5] ||
	    zmin != al->ra[6] ||
	    zmax != al->ra[7] ) {
		R(al);
	    }

	return R_imag;
}

 static real
kth_eig(register arglist *al)	/* minimum kth of a symmetric matrix */
{
	real z, lambdak, normx, normy, shift=0, cosine, ydotEkk, 
		lambda0=HUGE_VAL, hm, num;
	real *e0, *ek, *d, *h, *ra;
	int *at, i, ii, j, jj, k, kk, m, n, N, cnt=0, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, **E=NULL;
	static real *lambda=NULL, *x=NULL, *y=NULL, *y0=NULL;

	/*
	srand48(0);
	*/

	N = al->n - 1;
	if (N <= 0) return 0;
	n = (int)(sqrt(1+8*N)-1)/2;
	if (n*(n+1)/2 != N) {
	     fprintf(Stderr, 
	     "should be kth_eig(k, {i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	i0 = (int)al->ra[0];

	/*
	if (A == NULL) {
	    CALLOC( A, n, real * );
	} else {
	    REALLOC( A, n, real * );
	}
	for (i=0; i<n; i++) REALLOC( A[i], n-i, real );
	*/
	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	/*
	REALLOC( A[0], n*(n+1)/2, real);
	for (i=0, k=0; i<n; k+=n-i, i++) {
	    A[i] = A[0]+k;
	}
	*/
	REALLOC( A[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A[i] = A[0]+k;

	if (E == NULL) {
	    REALLOC( E, n, real *);
	    E[0] = NULL;
	} else {
	    REALLOC( E, n, real *);
	}
	REALLOC( E[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) E[i] = E[0]+k;

	REALLOC( lambda, n, real );

	REALLOC( x, n, real );

	REALLOC( y, n, real );

	if (y0 == NULL) {
	    /*
	    REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = drand48(); 
	    */
	    REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = 0.5; 
	} else {
	    REALLOC( y0, n, real ); 
	}

	k = 1;
	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		jj = at[k];
		if (jj >= 0) {
			A[i][j] = ra[jj];
		} else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_eig treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
		k++;
	    }
	}
	/* showmat(A,n,ae); */

	for (i=0; i<n; i++) y[i]  = y0[i];
	do {
	    normy = sqrt(dotprod(y,y,n));
	    for (i=0; i<n; i++) { y[i] /= normy; }
	    cosine = dotprod(x,y,n);
	    for (i=0; i<n; i++) { x[i] = y[i]; }
	    mx( A, x, y, n);
	    i = argmax(y,n);
	    lambda[0] = y[i]/x[i];
	    cnt++;
	} while ( fabs(cosine) < 0.99999999999 );

	shift = y0[0];

	for (i=0; i<n; i++) A[i][i] -= shift;

	for (k=0; k<n; k++) {
	    for (i=0; i<n; i++) y[i] = y0[i];
	    do {
		normy = sqrt(dotprod(y,y,n));
	        for (i=0; i<n; i++) { y[i] /= normy; }
	        /* showvec(y,n,ae); */
		cosine = dotprod(x,y,n);
	        for (i=0; i<n; i++) { x[i] = y[i]; }
		mx( A, x, y, n);
	        /* showvec(y,n,ae); */
		i = argmax(y,n);
	        lambdak = y[i]/x[i];
		cnt++;
	    } while ( fabs(cosine) < 0.9999999999 );

	    lambda[k] = lambdak + shift;
	    for (i=0; i<n; i++) E[k][i] = x[i];

	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    A[i][j] -= lambdak*x[i]*x[j];
		}
	    }

	    /* showmat(A,n,ae); */
	}

	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		A[i][j] = 0;
		for (k=0; k<n; k++) {
		    A[i][j] += lambda[k]*E[k][i]*E[k][j];
		}
	    }
	}
	/* showmat(A,n,ae); */

	i0--;	/* convert from 1..n to 0..n-1 */
	lambda0 = lambda[i0];
	e0 = E[i0];
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        d[m] = 2*e0[i]*e0[j];
			if (i==j) d[m] /= 2;
			m++;
		    }
		}
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    for (ii=i; ii<n; ii++) {
			for (jj=(ii==i)?j:ii; jj<n; jj++) {
			    hm = 0;
			    for (k=0; k<n; k++) {
				ek = E[k];
				if (k != i0) {
				    num = (e0[ i]*ek[ j] + e0[ j]*ek[ i])*
				          (e0[ii]*ek[jj] + e0[jj]*ek[ii]);
				    if ( i== j) num /= 2;
				    if (ii==jj) num /= 2;
				    hm += num/(lambda0-lambda[k]);
				}
			    }
			    h[m] = 2*hm;
			    m++;
			}
		    }
		}
	    }
	}
	return lambda0;
}

 /************** needs to be coded ***************************/

 static real
kth_diag(register arglist *al)	/* kth diagonal element of D in LDL^T */
				/* must be called in order: 1,2,...,n */
{
	real z, diagval;
	real *d, *h, *ra;
	int *at, i, ii, j, jj, k, m, n, N;
	static int k0=-1;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, Ident, **A0=NULL;
	static int convex_flag=1;

	N = al->n;
	if (N <= 0) return 0;
	n = (int)(sqrt(1+8*N)-1)/2;
	if (n*(n+1)/2 != N) {
	     fprintf(Stderr, 
	     "should be kth_diag({i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	k0++;  if (k0==n) k0 = 0;

/*
printf("IN KTH_DIAG \n");
printf("k0 = %3d \n", k0);
*/

	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	REALLOC( A[0], 2*n*n, real);
	for (i=0, k=0; i<n; k+=2*n, i++) A[i] = A[0]+k;

	if (A0 == NULL) {
	    REALLOC( A0, n, real *);
	    A0[0] = NULL;
	} else {
	    REALLOC( A0, n, real *);
	}
	REALLOC( A0[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A0[i] = A0[0]+k;

	if (k0 == 0) {
	    k = 0;
	    for(i=0; i<n; i++) {
	        for (j=i; j<n; j++) {
		    jj = at[k];
		    if (jj >= 0) {
			A[i][j] = ra[jj];
		    } else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_diag treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		    }
		    A[j][i] = A[i][j];
		    k++;
	        }
	        for(j=n; j<2*n; j++) {
		    A[i][j] = 0;
	        }
		A[i][n+i] = 1;
	    }
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    A0[i][j] = A[i][j];
		}
	    }
	    /*
	    showmat2(A,n,ae); 
	    */
	}

/*
printf("Gradient \n");
*/
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        if (i<k0 && j<k0) {
			    d[m] = A[i][k0]*A[j][k0];
			    if (j!=i) d[m] *= 2;
		        } else 
		        if (i<k0 && j==k0) {
			    d[m] = -2*A[i][k0];
		        } else 
		        if (i==k0 && j==k0) {
			    d[m] = 1;
		        } else {
			    d[m] = 0;
		        }
/*
printf("%2d %2d %10f \n", i,j,d[m]);
*/
			m++;
		    }
		}
	}
/*
printf("Hessian \n");
*/
	if (h) {
	    m = 0;
	    for (ii=0; ii<n; ii++) {
		for (jj=ii; jj<n; jj++) {
		    for (i=0; i<=ii; i++) {
			for (j=i; j<=(i==ii?jj:n-1); j++) {
			    if (i<k0 && j<k0 && ii<k0 && jj<k0) {
			        if ( i!=j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0]

			                   - A[jj][k0]*A[ii][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+jj]*A[ii][k0];
				} else 
				if ( i!=j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0];
				} else 
				if ( i==j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0];
				} else 
				if ( i==j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0];
				}
			    } else 
			    if (i<k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = A[i][k0]*A[j][n+ii]
				     + A[j][k0]*A[ii][n+i];
				if (i!=j) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+i]
				     + A[jj][k0]*A[i][n+ii];
				if (ii!=jj) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+i ] - A[i ][n+ii];
			    } else {
				h[m] = 0;
			    }
/*
printf("%2d %2d %2d %2d %10f \n", i,j,ii,jj,h[m]);
*/
			    m++;
			}
		    }
		}
	    }
	}

	for (i=0; i<k0; i++) {
	    for (j=k0; j<n+k0; j++) {
		A[k0][j] -= A[k0][i]*A[i][j];
	    }
	    A[k0][i] = 0;
	}
	diagval = A[k0][k0];
/*
printf("val=%10f \n", diagval);
*/
	for (j=k0; j<=n+k0; j++) {
	    A[k0][j] /= diagval;
	}
	for (i=0; i<k0; i++) {
	    for (j=k0+1; j<=n+k0; j++) {
		A[i][j] -= A[i][k0]*A[k0][j];
	    }
	    A[i][k0] = 0;
	}

	/*
	fprintf(stdout,"k = %3d \n", k0);
	showmat2(A,n,ae);
	fprintf(stdout,"\n");
	fflush(stdout);
	*/

	if (k0==n-1) {
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    Ident = 0;
		    for (k=0; k<n; k++) {
			Ident += A0[i][k]*A[k][n+j];
		    }
		    if ( (i==j && fabs(Ident - 1) > 0.001)
		       ||(i!=j && fabs(Ident) > 0.001) ) {
			   fprintf(stdout,"error in Gaussian elimination \n");
		       }
		}
	    }
	}

	if (diagval < -1.0e-12) {convex_flag = 0;}
	if (k0 == n-1) {
	    if (convex_flag == 0) {
		fprintf(stdout,"matrix X not positive semidefinite \n");
	    }
	    convex_flag = 1;
	}

	return diagval;
}

 static real
kth_diag_eps(register arglist *al)	/* kth diagonal element of D in LDL^T */
				/* must be called in order: 1,2,...,n */
{
	real z, diagval;
	real *d, *h, *ra;
	int *at, i, ii, j, jj, k, m, n, N;
	static int k0=-1;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, Ident, **A0=NULL;

	N = al->n;
	if (N <= 0) return 0;
	n = (int)(sqrt(1+8*N)-1)/2;
	if (n*(n+1)/2 != N) {
	     fprintf(Stderr, 
	     "should be kth_diag({i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	k0++;  if (k0==n) k0 = 0;

/*
printf("IN KTH_DIAG \n");
printf("k0 = %3d \n", k0);
*/

	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	REALLOC( A[0], 2*n*n, real);
	for (i=0, k=0; i<n; k+=2*n, i++) A[i] = A[0]+k;

	if (A0 == NULL) {
	    REALLOC( A0, n, real *);
	    A0[0] = NULL;
	} else {
	    REALLOC( A0, n, real *);
	}
	REALLOC( A0[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A0[i] = A0[0]+k;

	if (k0 == 0) {
	    k = 0;
	    for(i=0; i<n; i++) {
	        for (j=i; j<n; j++) {
		    jj = at[k];
		    if (jj >= 0) {
			A[i][j] = ra[jj];
		    } else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_diag treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		    }
		    A[j][i] = A[i][j];
		    k++;
	        }
	        for(j=n; j<2*n; j++) {
		    A[i][j] = 0;
	        }
		A[i][n+i] = 1;

		A[i][i] += 1.0e-6;
	    }
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    A0[i][j] = A[i][j];
		}
	    }
	    /*
	    showmat2(A,n,ae); 
	    */
	}

/*
printf("Gradient \n");
*/
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        if (i<k0 && j<k0) {
			    d[m] = A[i][k0]*A[j][k0];
			    if (j!=i) d[m] *= 2;
		        } else 
		        if (i<k0 && j==k0) {
			    d[m] = -2*A[i][k0];
		        } else 
		        if (i==k0 && j==k0) {
			    d[m] = 1;
		        } else {
			    d[m] = 0;
		        }
/*
printf("%2d %2d %10f \n", i,j,d[m]);
*/
			m++;
		    }
		}
	}
/*
printf("Hessian \n");
*/
	if (h) {
	    m = 0;
	    for (ii=0; ii<n; ii++) {
		for (jj=ii; jj<n; jj++) {
		    for (i=0; i<=ii; i++) {
			for (j=i; j<=(i==ii?jj:n-1); j++) {
			    if (i<k0 && j<k0 && ii<k0 && jj<k0) {
			        if ( i!=j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0]

			                   - A[jj][k0]*A[ii][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+jj]*A[ii][k0];
				} else 
				if ( i!=j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0];
				} else 
				if ( i==j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0];
				} else 
				if ( i==j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0];
				}
			    } else 
			    if (i<k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = A[i][k0]*A[j][n+ii]
				     + A[j][k0]*A[ii][n+i];
				if (i!=j) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+i]
				     + A[jj][k0]*A[i][n+ii];
				if (ii!=jj) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+i ] - A[i ][n+ii];
			    } else {
				h[m] = 0;
			    }
/*
printf("%2d %2d %2d %2d %10f \n", i,j,ii,jj,h[m]);
*/
			    m++;
			}
		    }
		}
	    }
	}

	for (i=0; i<k0; i++) {
	    for (j=k0; j<n+k0; j++) {
		A[k0][j] -= A[k0][i]*A[i][j];
	    }
	    A[k0][i] = 0;
	}
	diagval = A[k0][k0];
/*
printf("val=%10f \n", diagval);
*/
	for (j=k0; j<=n+k0; j++) {
	    A[k0][j] /= diagval;
	}
	for (i=0; i<k0; i++) {
	    for (j=k0+1; j<=n+k0; j++) {
		A[i][j] -= A[i][k0]*A[k0][j];
	    }
	    A[i][k0] = 0;
	}

	/*
	fprintf(stdout,"k = %3d \n", k0);
	showmat2(A,n,ae);
	fprintf(stdout,"\n");
	fflush(stdout);
	*/

	if (k0==n-1) {
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    Ident = 0;
		    for (k=0; k<n; k++) {
			Ident += A0[i][k]*A[k][n+j];
		    }
		    if ( (i==j && fabs(Ident - 1) > 0.001)
		       ||(i!=j && fabs(Ident) > 0.001) ) {
			   fprintf(stdout,"error in Gaussian elimination \n");
		       }
		}
	    }
	}

	return diagval;
}

 static real
kth_diag2(register arglist *al)	/* kth diagonal element of D in LDL^T */
				/* must be called in order: 1,2,...,n */
{
	real z, diagval;
	real *d, *h, *ra;
	int *at, i, ii, j, jj, k, m, n, N;
	static int k0=-1;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, Ident, **A0=NULL;

	N = al->n - 1;
	if (N <= 0) return 0;
	n = (int)(sqrt(N));
	if (n*n != N) {
	     fprintf(Stderr, 
	     "should be kth_diag2(k, {i in 1..n, j in 1..n} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	k = (int)al->ra[0];
	k--;	/* convert from 1..n to 0..n-1 */
	if (k != k0+1 && k != 0) {
	     fprintf(Stderr, 
	     "k in kth_diag2 not in increasing order \n");
	     fflush(Stderr);
	}
	k0 = k;

	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	REALLOC( A[0], 2*n*n, real);
	for (i=0, k=0; i<n; k+=2*n, i++) A[i] = A[0]+k;

	if (A0 == NULL) {
	    REALLOC( A0, n, real *);
	    A0[0] = NULL;
	} else {
	    REALLOC( A0, n, real *);
	}
	REALLOC( A0[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A0[i] = A0[0]+k;

	if (k0 == 0) {
	    k = 1;
	    for(i=0; i<n; i++) {
	        for (j=0; j<n; j++) {
		    jj = at[k];
		    if (jj >= 0) {
			A[i][j] = ra[jj];
		    } else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_diag2 treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		    }
		    k++;
	        }
	        for(j=n; j<2*n; j++) {
		    A[i][j] = 0;
	        }
		A[i][n+i] = 1;
	    }
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    A0[i][j] = A[i][j];
		}
	    }
	    showmat2(A,n,ae);
	}

	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=0; j<n; j++) {
		        if (i<k0 && j<k0) {
			    d[m] = A[i][k0]*A[j][k0];
		        } else 
		        if (i<k0 && j==k0) {
			    d[m] = -2*A[i][k0];
		        } else 
		        if (i==k0 && j<k0) {
			    d[m] = -2*A[j][k0];
		        } else 
		        if (i==k0 && j==k0) {
			    d[m] = 1;
		        } else {
			    d[m] = 0;
		        }
			m++;
		    }
		}
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    for (ii=i; ii<n; ii++) {
			for (jj=(ii==i)?j:0; jj<n; jj++) {
			    if (i<k0 && j<k0 && ii<k0 && jj<k0) {
			        h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				       - A[i ][k0]*A[j ][n+ii]*A[jj][k0];
			    } else 
			    if (i<k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = A[i][k0]*A[j][n+ii]
				     + A[j][k0]*A[ii][n+i];
			    } else 
			    if (i<k0 && j<k0 && ii==k0 && jj<k0) {
				h[m] = A[i][k0]*A[j][n+jj]
				     + A[j][k0]*A[jj][n+i];
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+i]
				     + A[jj][k0]*A[i][n+ii];
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+i ] - A[i ][n+ii];
			    } else 
			    if (i<k0 && j==k0 && ii==k0 && jj<k0) {
				h[m] = - A[jj][n+i ] - A[i ][n+jj];
			    } else 
			    if (i==k0 && j<k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+j]
				     + A[jj][k0]*A[j][n+ii];
			    } else 
			    if (i==k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+j ] - A[j ][n+ii];
			    } else 
			    if (i==k0 && j<k0 && ii==k0 && jj<k0) {
				h[m] = - A[jj][n+j ] - A[j ][n+jj];
			    } else {
				h[m] = 0;
			    }
			    m++;
			}
		    }
		}
	    }
	}

	for (i=0; i<k0; i++) {
	    for (j=k0; j<n+k0; j++) {
		A[k0][j] -= A[k0][i]*A[i][j];
	    }
	    A[k0][i] = 0;
	}
	diagval = A[k0][k0];
	for (j=k0; j<=n+k0; j++) {
	    A[k0][j] /= diagval;
	}
	for (i=0; i<k0; i++) {
	    for (j=k0+1; j<=n+k0; j++) {
		A[i][j] -= A[i][k0]*A[k0][j];
	    }
	    A[i][k0] = 0;
	}

	/*
	fprintf(stdout,"k = %3d \n", k0);
	showmat2(A,n,ae);
	fprintf(stdout,"\n");
	fflush(stdout);
	*/

	if (k0==n-1) {
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    Ident = 0;
		    for (k=0; k<n; k++) {
			Ident += A0[i][k]*A[k][n+j];
		    }
		    if ( (i==j && fabs(Ident - 1) > 0.00001)
		       ||(i!=j && fabs(Ident) > 0.00001) ) {
			   fprintf(stdout,"error in Gaussian elimination \n");
		       }
		}
	    }
	}

	return diagval;
}

 static real
myerf(arglist *al)	/* 1/sqrt(2 pi) int_-infty^x e^{-t^2/2} dt */
{
	int j;
	double h = 1.0/128.0, x0 = -10;
	double h2 = h/2;
	double myerf=0;
	double pi = 4*atan(1);
	double rtp = sqrt(2*pi);
	double rt  = sqrt(2);
	double deriv;

	real x = al->ra[0];

	myerf = (erf(x/rt) + 1)/2;

	if (al->derivs) {
		deriv = exp(-x*x/2)/rtp;
		*al->derivs = deriv;
		if (al->hes)
			*al->hes = -x * deriv;
		}
	return myerf;
	}

 static real
stderf(arglist *al)	/* 2/sqrt(pi) int_0^x e^{-t^2} dt */
{
	int j;
	double h = 1.0/128.0, x0 = -10;
	double h2 = h/2;
	double stderf=0;
	double pi = 4*atan(1);
	double rp = sqrt(pi);
	double rt  = sqrt(2);
	double deriv;

	real x = al->ra[0];

	stderf = erf(x);

	if (al->derivs) {
		deriv = 2*exp(-x*x)/rp;
		*al->derivs = deriv;
		if (al->hes)
			*al->hes = -2*x * deriv;
		}
	return stderf;
	}

static void showvec( real *x, int n, AmplExports *ae)
{
	int i, j;

	for(i=0; i<n; i++) {
	    fprintf(stderr, "%10.2e ", x[i]);
	}
	fprintf(Stderr, "\n");
	fflush(Stderr);
}

static void showmat( real **A, int n, AmplExports *ae)
{
	int i, j;

	for(i=0; i<n; i++) {
	    for (j=0; j<n; j++) {
		if (j>=i) {
		    fprintf(stderr, "%10.4f ", A[i][j]);
		} else {
		    fprintf(stderr, "           ");
		}
	    }
	    fprintf(Stderr, "\n");
	    fflush(Stderr);
	}
}

static void showmat2( real **A, int n, AmplExports *ae)
{
	int i, j;

	for(i=0; i<n; i++) {
	    for (j=0; j<2*n; j++) {
		    fprintf(stdout, "%10.4f ", A[i][j]);
	    }
	    fprintf(stdout, "\n");
	    fflush(stdout);
	}
}

static void mx( real **A, real *x, real *y, int n)
{
	int i,j;

	for (i=0; i<n; i++) y[i] = 0;
	for (i=0; i<n; i++) {
	    y[i] += A[i][i]*x[i];
	    for (j=i+1; j<n; j++) {
		y[i] += A[i][j]*x[j];
		y[j] += A[i][j]*x[i];
	    }
	}
}

static real dotprod(         /* inner product between n-vectors x and y */
        real  *x,
        real  *y,
        int     n
)
{
        int i;
        real dotprod=0.0e0;

        for (i=0; i<n; i++) dotprod += x[i]*y[i];

        return dotprod;
}

static int argmax(         /* inner product between n-vectors x and y */
        real  *x,
        int     n
)
{
        int i, ii = -1;
        real max = -1.0e0;

        for (i=0; i<n; i++) if (fabs(x[i]) > max) {ii=i; max=fabs(x[i]);}

        return ii;
}

 static real
mymin(register arglist *al)	/* minimum eigenvalue of a symmetric matrix */
{
	real z, x0=HUGE_VAL;
	real *d, *h, *ra;
	int *at, i, j, jj, k, m, n, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real *x=NULL;

	n = al->n;
	if (n <= 0) return 0;
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	REALLOC( x, n, real );

	for(i=0; i<n; i++) {
		jj = at[i];
		if (jj >= 0) {
			x[i] = ra[jj];
		} else {
			x[i] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"mymin treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
	}

	for (k=0; k<n; k++) {
	    if (x[k] < x0) {
		i0 = k;
		x0 = x[k];
	    }
	}
/*
fprintf(stdout, "min var = %d, val = %10.5f \n", i0, x0); fflush(stdout);
*/

	if (d) {
		for (i=0; i<n; i++) { d[i] = 0; }
		d[i0] = 1;
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    h[m] = 0;
		    m++;
		}
	    }
	}

	return x0;
}

 static real
nonneg(register arglist *al)	/* minimum eigenvalue of a symmetric matrix */
{
	real z, x0=HUGE_VAL;
	real *d, *h, *ra;
	int *at, i, j, jj, k, m, n, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real *x=NULL;

	n = al->n - 1;
	if (n <= 0) return 0;
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	REALLOC( x, n+1, real );

	i0 = (int)al->ra[0];

	for(i=1; i<=n; i++) {
		jj = at[i];
		if (jj >= 0) {
			x[i] = ra[jj];
		} else {
			x[i] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"nonneg treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
	}

	if (d) {
		for (i=0; i<n; i++) { d[i] = 0; }
		d[i0-1] = 1;
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    h[m] = 0;
		    m++;
		}
	    }
	}

	return x[i0];
}

 static char *
kth(register arglist *al)	/* kth(k,a1,a2,...,an) return ak */
{
	int i, j, k, n;
	char *comma;
	static char buf[32];
	AmplExports *ae = al->AE;

	k = al->at[0] ? atoi(al->sa[0]) : (int)al->ra[0];
	n = al->n;
	if (k < 0) {
		fprintf(Stderr, "kth(");
		comma = "";
		for(i = 0; i < n; i++, comma = ", ")
			if ((j = al->at[i]) >= 0)
				fprintf(Stderr, "%s%.g", comma, al->ra[j]);
			else
				fprintf(Stderr, "%s%s", comma, al->sa[-(j+1)]);
		fprintf(Stderr, ")\n");
		fflush(Stderr);
		k = -k;
		}
	if (n <= 1 || k <= 0 || k >= n)
		return "";
	if ((j = al->at[k]) >= 0) {
		sprintf(buf, "%.g", al->ra[j]);
		return buf;
		}
	return al->sa[-(j+1)];
	}

 void
funcadd(AmplExports *ae){
	/* Insert calls on addfunc here... */

/* Arg 3, called type, must satisfy 0 <= type <= 6:
 * type&1 == 0:	0,2,4,6	==> force all arguments to be numeric.
 * type&1 == 1:	1,3,5	==> pass both symbolic and numeric arguments.
 * type&6 == 0:	0,1	==> the function is real valued.
 * type&6 == 2:	2,3	==> the function is char * valued; static storage
			    suffices: AMPL copies the return value.
 * type&6 == 4:	4,5	==> the function is random (real valued).
 * type&6 == 6: 6	==> random, real valued, pass nargs real args,
 *				0 <= nargs <= 2.
 *
 * Arg 4, called nargs, is interpretted as follows:
 *	>=  0 ==> the function has exactly nargs arguments
 *	<= -1 ==> the function has >= -(nargs+1) arguments.
 */

	addfunc("ginv",     (rfunc)ginv,     0,  1, 0);
	addfunc("sginv",    (rfunc)sginv,    2,  1, 0);
	addfunc("hypot",    (rfunc)myhypot,  0,  2, 0);
	addfunc("ncall",    (rfunc)ncall,    0,  0, 0);
	addfunc("rncall",   (rfunc)ncall,    4,  0, 0); /*could change 4 to 6*/
	addfunc("mean0",    (rfunc)mean,     0, -1, 0);
	addfunc("mean",     (rfunc)mean,     1, -1, 0);
	addfunc("mineig",   (rfunc)mineig,   1, -1, 0);
	addfunc("kth_eig",  (rfunc)kth_eig,  1, -1, 0);
	addfunc("linear",   (rfunc)linear,   1, -1, 0);
	addfunc("C",   	    (rfunc)C,        0,  8, 0);
	addfunc("S",   	    (rfunc)S,        0,  8, 0);
	addfunc("kth_diag", (rfunc)kth_diag, 1, -1, 0);
	addfunc("kth_diag_eps",(rfunc)kth_diag_eps,1, -1, 0);
	addfunc("kth_diag2",(rfunc)kth_diag2,1, -1, 0);
	addfunc("myerf",    (rfunc)myerf,    0,  1, 0);
	addfunc("stderf",   (rfunc)stderf,   0,  1, 0);
	addfunc("mymin",    (rfunc)mymin,    1, -1, 0);
	addfunc("nonneg",   (rfunc)nonneg,   1, -1, 0);
	addfunc("kth",      (rfunc)kth,      3, -2, 0);
	}

static double erf (double x)
{
    return x < 0.0 ? -gammp(0.5, x*x) : gammp(0.5, x*x);
}

static double gammln(double xx)
{
    double x, tmp, ser;
    static double cof[6]={76.18009173, -86.50532033, 24.01409822,
			  -1.231739516, 0.120858003e-2, -0.536382e-5};
    int j;
    x = xx-1.0;
    tmp = x+5.5;
    tmp -= (x+0.5)*log(tmp);
    ser=1.0;
    for (j=0; j<=5;j++) {
	x += 1.0;
	ser += cof[j]/x;
    }
    return -tmp+log(2.50662827465*ser);
}

static double gammp (double a, double x)
{
    double gamser, gammcf, gln;

    if (x < 0.0 || a <= 0.0) {
	/* fprintf(stderr,"invalid arguments in gammp\n"); exit(0); */
    }
    if (x < (a+1.0)) {
	gser(&gamser,a,x,&gln);
	return gamser;
    } else {
	gcf(&gammcf,a,x,&gln);
	return 1.0-gammcf;
    }
}

#define ITMAX 100
#define EPS 3.0e-7

static void gser(double *gamser, double a, double x, double *gln)
{
    int n;
    double sum, del, ap;
    *gln=gammln(a);
    if (x <= 0.0) {
	/* if (x < 0.0) {fprintf(stderr,"x < 0 in gser\n"); exit(0);} */
	*gamser=0.0;
	return;
    } else {
	ap = a;
	del = sum = 1.0/a;
	for (n=1; n<= ITMAX; n++) {
	    ap += 1.0;
	    del *= x/ap;
	    sum += del;
	    if (fabs(del) < fabs(sum)*EPS) {
		*gamser=sum*exp(-x+a*log(x)-(*gln));
		return;
	    }
	}
	/*
	fprintf(stderr,"a too large or ITMAX too small in gser\n"); exit(0);
	*/
    }
}

static void gcf( double *gammcf, double a, double x, double *gln)
{
    int n;
    double gold=0.0,g,fac=1.0, b1=1.0;
    double b0=0.0, anf, ana, an, a1, a0=1.0;
    *gln = gammln(a);
    a1 = x;
    for (n=1; n<=ITMAX; n++) {
	an= (double) n;
	ana = an - a;
	a0 = (a1+a0*ana)*fac;
	b0 = (b1+b0*ana)*fac;
	anf = an*fac;
	a1 = x*a0+anf*a1;
	b1 = x*b0+anf*b1;
	if (a1) {
	    fac = 1.0/a1;
	    g = b1*fac;
	    if (fabs((g-gold)/g) < EPS) {
		*gammcf = exp(-x+a*log(x)-(*gln))*g;
		return;
	    }
	    gold = g;
	}
    }
    /*
    fprintf(stderr,"a too large or ITMAX too small in gcf\n"); exit(0);
    */
}

#ifdef __cplusplus
	}
#endif