//rlp_math.cpp, Copyright (c) 2004-2007 R.Lackner
//
//    This file is part of RLPlot.
//
//    RLPlot 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.
//
//    RLPlot 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 RLPlot; if not, write to the Free Software
//    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
#include "rlplot.h"
#include <math.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <time.h>

#define SWAP(a,b) {double temp=(a);(a)=(b);(b)=temp;}
#define _PREC 1.0e-12

extern Default defs;

static char *MRQ_error = 0L;
static double sqrt2pi = sqrt(_PI*2.0);

//---------------------------------------------------------------------------
//utilitity functions for memory allocation
double **dmatrix(int nrl, int nrh, int ncl, int nch)
{
	int i;
	double **m;

	m = (double **)malloc(nrh * sizeof(double*));
	//Allocate rows and set pointers to them
	for(i = 0; i < nrh; i++) {
		m[i] = (double *)malloc(nrh * sizeof(double));
		}
	return m;
}
void free_dmatrix(double **m, int nrl, int nrh, int ncl, int)
{
	int i;

	for(i = 0; i < nrh; i++) free(m[i]);
	free(m);
}

//---------------------------------------------------------------------------
//The routine gaussj solves linear equations by Gauss-Jordan elimination
bool gaussj(double **a, int n, double **b, int m)
{
	int *indxc, *indxr, *ipiv;
	int i, icol, irow, j, k, l, ll;
	double big, dum, pivinv;

	indxc = (int*)malloc(n*sizeof(int*));
	indxr = (int*)malloc(n*sizeof(int*));
	ipiv = (int*)malloc(n*sizeof(int*));
	for (j = 0; j < n; j++) ipiv[j] = 0;
	for (i = 0; i < n; i++) {				//This is the main loop over the
		big = 0.0;							//    columns to be reduced
		for(j = 0; j < n; j ++)				//This is the outer loop of the search
			if(ipiv[j] != 1)				//    for a pivot element
				for(k = 0; 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) {
						MRQ_error = "Singular Matrix (1)";
						free(ipiv);		free(indxr);	free(indxc);
						return false;
						}
				}
		++(ipiv[icol]);
		//We now have the pivot element, so we interchange rows, if needed,
		// to put the pivot element on the diagonal.
		if(irow != icol) {
			for(l = 0; l < n; l++) SWAP(a[irow][l], a[icol][l])
			for(l = 0; l < m; l++) SWAP(b[irow][l], b[icol][l])
			}
		indxr[i] = irow;		indxc[i] = icol;
		if(a[icol][icol] == 0.0) {
			MRQ_error = "Singular Matrix (2)";
			free(ipiv);		free(indxr);	free(indxc);
			return false;
			}
		pivinv = 1.0/a[icol][icol];
		a[icol][icol] = 1.0;
		for(l = 0; l < n; l++) a[icol][l] *= pivinv;
		for(l = 0; l < m; l++) b[icol][l] *= pivinv;
		for(ll = 0; ll <  n; ll++)
			if(ll != icol) { 							//Next, we reduce the rows
				dum = a[ll][icol];
				a[ll][icol] = 0.0;
				for(l = 0; l < n; l++) a[ll][l] -= a[icol][l]*dum;
				for(l = 0; l < m; l++) b[ll][l] -= b[icol][l]*dum;
				}
		}											// This is the end of the main loop
	for (l = n; l > 0; l--) {						//   over columns of the reduction.
		if(indxr[l] != indxc[l]) 					//   Unscramble the solution
			for(k = 0; k < n; k++) SWAP (a[k][indxr[l]], a[k][indxc[l]]);
		}											//And we are done.
	free(ipiv);		free(indxr);	free(indxc);
	return true;
}

//---------------------------------------------------------------------------
//The routine mrqcof is called by mrqmin to evaluate the linearized fitting
// matrix alpha and vector beta
void mrqcof(double x[], double y[], double z[], int ndata, double **a, int ma,
	int lista[], int mfit, double **alpha, double beta[], double *chisq,
	void (*funcs)(double, double, double **, double *, double *, int))
{
	int k, j, i;
	double ymod, wt, dy;
	double *dyda;

	dyda = (double*)malloc(ma*sizeof(double));
	for(j = 0; j < mfit; j++) {					//Initialize (symmetric) alpha, beta
		for(k = 0; k <= j; k++) alpha[j][k] = 0.0;
		beta[j] = 0.0;
		}
	*chisq = 0.0;
	for (i = 0; i < ndata; i++) {		 		//Summation loop over all data
		(*funcs)(x[i], z ? z[i] : 0.0, a, &ymod, dyda, ma);
		if(ymod != 0.0) dy = y[i]-ymod;			//functions = 0.0 if out of range
		else dy = 0.0;
		for(j = 0; j < mfit; j++) {
			wt = dyda[lista[j]];
			for (k = 0; k <= j; k++){
				alpha[j][k] += wt*dyda[lista[k]];
				}
			beta[j] += dy*wt;
			}
		(*chisq) += dy*dy; 							//And find X^2 if function o.k.
		}
	for(j = 0; j < mfit; j++)						//Fill the symmetric side
		for(k = 0; k <= j; k++) alpha[k][j]=alpha[j][k];
	free(dyda);
}

//---------------------------------------------------------------------------
//The routine mrqmin performs one iteration of Marquart's method for nonlinear
// parameter estimation
bool mrqmin(double *x, double *y, double *z, int ndata, double **a, int ma,
	int *lista, int mfit, double **covar, double **alpha, double *chisq,
	void (*funcs)(double, double, double **, double *, double *, int), double *alamda)
{
	int k, kk, j, ihit;
	static double *da, *atry, *beta, ochisq;
	static double **oneda, **atryref;

	if (*alamda < 0.0) {								//Initialization
		MRQ_error = 0L;
		oneda = dmatrix(1, mfit, 1, 1);
		atry = (double *)malloc(ma * sizeof(double));
		atryref = (double**)malloc(ma * sizeof(double*));
		for(j=0; j < ma; atryref[j++] = &atry[j]);
		da = (double*)malloc(ma *sizeof(double));
		beta = (double*)malloc(ma *sizeof(double));
		kk = mfit+1;
		for(j = 0; j < ma; j++) { 						//Does lista contain a proper
			ihit = 0;									//   permutation of the
			for(k = 0; k < mfit; k++)					//   coefficients ?
				if(lista[k] == j) ihit++;
			if(ihit == 0)
				lista[kk++] = j;
			else if (ihit >1) ErrorBox("Bad LISTA permutations in MRQMIN-1");
			}
		if(kk != ma+1) ErrorBox("Bad LISTA permutations in MRQMIN-2");
		*alamda = 0.001;
		mrqcof(x, y, z, ndata, a, ma, lista, mfit, alpha, beta, chisq, funcs);
		ochisq=(*chisq);
		}
	for (j = 0; j < mfit; j++) {						//Alter linearized fitting matrix
		for(k = 0; k < mfit; k++) covar[j][k] = alpha[j][k];	// by augmenting
		covar[j][j] = alpha[j][j]*(1.0+(*alamda));		// diagaonal elements
		oneda[j][0] = beta[j];
		}
	if (!gaussj(covar, mfit, oneda, 1)) return false;	//Matrix solution ?
	for(j = 0; j < mfit; j++) da[j] = oneda[j][0];
	if(*alamda == 0.0) {								//Once converged evaluate
														//  covariance matrix with
		free(beta);										//  alamda = 0.
		free(da);
		free(atry);
		free(atryref);
		free_dmatrix(oneda, 1, mfit, 1, 1);
		return true;
		}
	for(j = 0; j < ma; j++) atry[j] = *a[j];
	for(j = 0; j < mfit; j++)							//Did the trial succeed ?
		atry[lista[j]] = *a[lista[j]] + da[j];
	mrqcof(x, y, z, ndata, atryref, ma, lista, mfit, covar, da, chisq, funcs);
	if(*chisq < ochisq) {								//Success, accept the new solution
		*alamda *= 0.1;
		ochisq=(*chisq);
		for(j = 0; j < mfit; j++) {
			for(k = 0; k < mfit; k++) alpha[j][k] = covar[j][k];
			beta[j] = da[j];
			*a[lista[j]] = atry[lista[j]];
			}
		}
	else {												//Failure, increase almda and
		*alamda *= 10.0;								//    return.
		*chisq = ochisq;
		}
	return true;
}

bool Check_MRQerror()
{
	bool bRet;

	if(bRet = MRQ_error != 0L) ErrorBox(MRQ_error);
	MRQ_error = 0L;
	return bRet;
}

//---------------------------------------------------------------------------
//Use heap sort to sort elements of an float array
//W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1988/1989)
//Numerical Recipes in C, Cambridge University Press, ISBN 0-521-35465-X
// p. 245
void SortArray(int n, double *vals)
{
	int l, j, ir, i;
	double rra, *ra = vals-1;

	if(n < 2 || !vals) return;
	l=(n >> 1) + 1;				ir = n;
	for( ; ; ) {
		if(l > 1) rra = ra[--l];
		else {
			rra = ra[ir];		ra[ir] = ra[1];
			if(--ir == 1) {
				ra[1] = rra;	return;
				}
			}
		i = l;					j = l << 1;
		while (j <= ir) {
			if (j < ir && ra[j] < ra[j+1]) ++j;
			if (rra < ra[j]) {
				ra[i] = ra[j];	j += (i=j);
				}
			else j = ir + 1;
			}
		ra[i] = rra;
		}
}

//sorts array v1 making the corresponding rearrangement of v2
void SortArray2(int n, double *v1, double *v2)
{
	int l, j, ir, i;
	double rra, rrb, *ra = v1-1, *rb = v2-1;

	if(n < 2 || !v1 || !v2) 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;
		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);
				}
			else j = ir + 1;
			}
		ra[i] = rra;			rb[i] = rrb;
		}
}

//Use heap sort to sort elements of an xy array
void SortFpArray(int n, lfPOINT *vals)
{
	int l, j, ir, i;
	lfPOINT rra, *ra = vals-1;

	if(n < 2) return;
	l=(n >> 1) + 1;					ir = n;
	for( ; ; ) {
		if(l > 1) {
			rra.fx = ra[--l].fx; rra.fy = ra[l].fy;
			}
		else {
			rra.fx = ra[ir].fx;		rra.fy = ra[ir].fy;
			ra[ir].fx = ra[1].fx;	ra[ir].fy = ra[1].fy;	
			if(--ir == 1) {
				ra[1].fx = rra.fx;	ra[1].fy = rra.fy;
				return;
				}
			}
		i = l;					j = l << 1;
		while (j <= ir) {
			if (j < ir && ra[j].fx < ra[j+1].fx) ++j;
			if (rra.fx < ra[j].fx) {
				ra[i].fx = ra[j].fx;	ra[i].fy = ra[j].fy;
				j += (i=j);
				}
			else j = ir + 1;
			}
		ra[i].fx = rra.fx;				ra[i].fy = rra.fy;
		}
}

//randomize array
double *randarr(double *v0, int n, long *seed)
{
	double r, *v, *v_tmp;
	int i, j, l;

	if(!(v = (double*)malloc(n *sizeof(double)))) return 0L;
	if(!(v_tmp = (double*)memdup(v0, n *sizeof(double),0))) return 0L;
	for(l = n, i = 0; i < n; ) {
		r = ran2(seed);			j = (int)(r *((double)l));
		if(j < l) {
			v[i++] = v_tmp[j];
			if(j < l)memcpy(v_tmp+j, v_tmp+j+1, (l-j)*sizeof(double));
			l--;
			}
		}
	return v;
}

//resample array
double *resample(double *v0, int n, long *seed)
{
	double r, *v;
	int i, j;

	if(!(v = (double*)malloc(n *sizeof(double)))) return 0L;
	for(i = 0; i < n; ) {
		r = ran2(seed);			j = (int)(r *((double)n));
		if(j < n) v[i++] = v0[j];
		}
	return v;
}

//---------------------------------------------------------------------------
// Cubic Spline Interpolation
// Ref.: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Rcipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 96 ff.
void spline(lfPOINT *v, int n, double *y2)
{
	int i, k;
	double p, qn, sig, un, *u;

	u = (double *)malloc(n * sizeof(double));
	y2[0] = u[0] = 0.0;
	for(i = 1; i < (n-1); i++) {
		sig = (v[i].fx-v[i-1].fx)/(v[i+1].fx-v[i-1].fx);
		p = sig*y2[i-1]+2.0;			y2[i]=(sig-1.0)/p;
		u[i]=(v[i+1].fy-v[i].fy)/(v[i+1].fx-v[i].fx)-(v[i].fy-v[i-1].fy)/(v[i].fx-v[i-1].fx);
		u[i]=(6.0*u[i]/(v[i+1].fx-v[i-1].fx)-sig*u[i-1])/p;
		}
	qn = un = 0.0;
	y2[n-1] = (un - qn * u[n-2])/(qn*y2[n-2]+1.0);
	for(k = n-2; k >= 0; k--) {
		y2[k] = y2[k]*y2[k+1]+u[k];
		}
	free(u);
}

//---------------------------------------------------------------------------
// The Gamma Function: return the ln(G(xx)) for xx > 0
// Ref: B.W. Brown, J. Lovato, K. Russel (1994)
//    DCDFLIB.C, Library of C Routinesfor Cumulative Distribution Functions,
//    Inverses, and other Parameters.

double devlpl(double a[], int n, double x)
{
	double term;
	int i;

	for(term = a[n-1], i= n-2; i>=0; i--) term = a[i] + term * x;
	return term;
}


double gammln(double x)
{
	static double coef[] = {0.83333333333333023564e-1,-0.27777777768818808e-2, 
	0.79365006754279e-3, -0.594997310889e-3, 0.8065880899e-3};
static double scoefd[] = {0.62003838007126989331e2, 0.9822521104713994894e1,
	-0.8906016659497461257e1, 0.1000000000000000000e1};
static double scoefn[] = {0.62003838007127258804e2, 0.36036772530024836321e2,
	0.20782472531792126786e2, 0.6338067999387272343e1,0.215994312846059073e1,
	0.3980671310203570498e0, 0.1093115956710439502e0,0.92381945590275995e-2,
	0.29737866448101651e-2};
	double offset, prod, xx;
	int i,n;

    if(x < 6.0) {
		prod = 1.0e0;	    xx = x;
		while(xx > 3.0) {
			xx -= 1.0;			prod *= xx;
			}
		if(x <= 2.0) while(xx < 2.0) {
			prod /= xx;			xx += 1.0;
			}
		// compute rational approximation to gamma(x)
		return log(devlpl(scoefn, 9, xx-2.0) / devlpl(scoefd, 4, xx-2.0) * prod);
		}
	else {
		offset = 0.91893853320467274178;	// hln2pi
		// if necessary make x at least 12 and carry correction in offset
		if(n = 13.0 >= x ? (int)(12.0 - x) : 0) xx = x;
		else {
			for(i=1, prod = 1.0; i<= n; i++) prod *= (x+(double)(i-1));
			offset -= log(prod);			xx = x+(double)n;
			}
		// compute power series
		return devlpl(coef, 5, 1.0/(xx*xx)) / xx + (offset+(xx-0.5)*log(xx)-xx);
		}
}

//---------------------------------------------------------------------------
// Special Functions
// Ref.: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Rcipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 166 ff.

//The Factorial Function: return n!
double factrl(int n)
{
	static int ntop = 4;
	static double a[33]={1.0, 1.0, 2.0, 6.0, 24.0};
	int j;

	if(n < 0) return 0.0;		//error: no factorial for negative numbers
	if(n > 32) return exp(gammln(n+1.0));
	while(ntop < n) {			//fill in table up to desired value
		j = ntop++;		a[ntop]=a[j] * ntop;
		}
	return a[n];
}

//returns the incomplete gamma function evaluated by its series representation
void gser(double *gamser, double a, double x, double *gln)
{
	int n;
	double sum, del, ap;

	*gln = gammln(a);
	if(x <= 0) {
		*gamser = 0.0;			return;
		}
	else {
		ap = a;					del = sum = 1.0/a;
		for(n = 1; n <= 100; n++) {
			ap += 1.0;			del *= x/ap;		sum += del;
			if(fabs(del) <= fabs(sum) * _PREC) {
				*gamser = sum * exp(-x + a * log(x)-(*gln));
				return;
				}
			}
		// maximum number of iterations exceeded
		*gamser = sum * exp(-x + a * log(x)-(*gln));
		}

}

//returns the incomplete gamma function evaluated by its continued fraction representation
void gcf(double *gammcf, double a, double x, double *gln)
{
	int n;
	double gold=0.0, g, fac=1.0, b1=1.0, b0=0.0, anf, ana, an, a1, a0=1.0;

	*gln=gammln(a);		a1=x;
	for(n=1; n <= 100; 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) <= _PREC) {
				*gammcf = exp(-x + a * log(x) -(*gln)) * g;
				return;
				}
			gold = g;
			}
		}
	// maximum number of iterations exceeded
	*gammcf = exp(-x + a * log(x) -(*gln)) * gold;
}

//returns the incomplete gamma function P(a,x)
double gammp(double a, double x)
{
	double gamser, gammcf, gln;

	if(x < 0.0 || a <= 0.0) return 0.0;
	if(x < (a+1.0)) {
		gser(&gamser, a, x, &gln);			return gamser;
		}
	else {
		gcf(&gammcf, a, x, &gln);			return 1.0-gammcf;
		}
	return 0.0;
}

//returns the complementary incomplete gamma function Q(a,x)
double gammq(double a, double x)
{
	double gamser, gammcf, gln;

	if(x < 0.0 || a <= 0.0) return 0.0;
	if(x < (a+1.0)) {
		gser(&gamser, a, x, &gln);			return 1.0-gamser;
		}
	else {
		gcf(&gammcf, a, x, &gln);			return gammcf;
		}
	return 0.0;
}

//continued fraction for incomplete beta function, used by betai()
double betacf(double a, double b, double x)
{
	double qap, qam, qab, em, tem, d, bz, bm = 1.0, bp, bpp, az = 1.0, am = 1.0, ap, app, aold;
	int m;

	qab = a+b;		qap = a+1.0;		qam = a-1.0;	bz = 1.0-qab*x/qap;
	for(m = 1; m <= 100; m++) {
		em = (double)m;			tem = em+em;
		d = em*(b-em)*x/((qam+tem)*(a+tem));
		ap = az + d * am;		bp = bz + d *bm;
		d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
		app = ap + d * az;		bpp = bp + d * bz;
		aold = az;				am = ap/bpp;
		bm = bp/bpp;			az = app/bpp;
		bz = 1.0;
		if(fabs(az-aold) <= (_PREC * fabs(az))) return az;	//success: return
		}
	return az;												//fail: iterations exceeded
}

//The incomplete beta function Ix(a,b) for 0 <= x <= 1
double betai(double a, double b, double x)
{
	double bt;

	if(x < 0.0 || x > 1.0) return 0.0;		//range !
	if(x == 0.0 || x == 1.0) bt = 0.0;
	else
		bt = exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
	if(x < (a+1.0)/(a+b+2.0)) return bt * betacf(a, b, x)/a;
	else return 1.0 - bt * betacf(b, a, 1.0 - x)/b;
}

//The following relations are obviously based on:
//  Abramowitz, M. & Stegun I.A. (1964): Hanbook of Mathematical Functions.
//    Applied Mathematics Series, vol. 55 (Washington: National Bureau
//    of Standards).

//the binomial coefficient
double bincof(double n, double k)
{
	if(n<0 || k<0 || k > n) return 0.0;
	return exp(gammln(n+1.0) - gammln(k+1.0) - gammln(n-k+1.0));
}

//the cumulative binomial distribution
double binomdistf(double k, double n, double p)
{
	if(k > n || n < 0.0 || p < 0.0 || p >1.0) return 0.0;
	return betai(n-k, k+1, p);
}

//the beta function
double betaf(double z, double w)
{
	return exp(gammln(z)+gammln(w)-gammln(z+w));
}

//the error function: not all compilers have a built in erf()
double errf(double x)
{
	return x < 0.0 ? -gammp(0.5, x*x) : gammp(0.5, x*x);
}

//the complementary error function
double  errfc(double x)
{
//	return x < 0.0 ? 2.0 - gammq(0.5, x*x) : gammq(0.5, x*x);
	return x < 0.0 ? 1.0 + gammp(0.5, x*x) : gammq(0.5, x*x);
}

//cumulative normal distribution
double norm_dist(double x, double m, double s)
{
	return 0.5 + errf((x - m)/(s * _SQRT2))/2.0;
}

//normal distribution
double norm_freq(double x, double m, double s)
{
	double ex;

	ex = (x-m)/s;	ex = exp(-0.5*ex*ex);
	return ex/(s*sqrt2pi);
}

//cumulative exponential distribution
double exp_dist(double x, double l, double s)
{
	if(x >= 0.0 && l > 0.0) return 1.0-exp(-x*l);
	else return 0.0;
}

//inverse exponential distribution
double exp_inv(double p, double l, double s)
{
	if(p >= 1.0) return HUGE_VAL;
	if(l <= 0.0) return 0.0;
	return -log(1.0-p)/l;
}

//exponential distribution
double exp_freq(double x, double l, double s)
{
	if(x >= 0.0 && l > 0.0) return l*exp(-x*l);
	else return fabs(l);
}

//cumulative lognormal distribution
double lognorm_dist(double x, double m, double s)
{
	return 0.5 + errf((log(x) - m)/(s * _SQRT2))/2.0;
}

//lognormal distribution
double lognorm_freq(double x, double m, double s)
{
	double tmp;

	if(x > 0.0 && m > 0.0 && s > 0.0) {
		tmp = (log(x)-m)/s;
		return exp(-0.5*tmp*tmp)/(x*s*sqrt2pi);
		}
	return 0.0;
}

//chi square distribution
double chi_dist(double x, double df, double)
{
	if(x <= 0.0) return 1.0;
	return gammq(df/2.0, x/2.0);
}

double chi_freq(double x, double df)
{
	if(x < 0.0 || df <= 0.0) return 0.0;
	if(x < 1.0e-32) x = 1.0e-32;
//formula by Wikipedia 
//	return exp(log(2.0)*(1.0-df/2.0)+log(x)*(df-1.0)+x*x/-2.0-gammln(df/2.0));
//formula by StatSoft's STATISTICA documentation
	return exp(-x/2.0+log(x)*(df/2.0-1.0)-log(2.0)*df/2.0-gammln(df/2.0));
}

//t-distribution
double t_dist(double t, double df, double)
{
	return betai(df/2.0, 0.5, (df/(df+t*t)));
}

double t_freq(double t, double df)
{
	double a, b, c, d;
 
	a = gammln((df+1.0)/2.0);		b = log(sqrt(df * _PI));
	c = gammln(df/2.0);				d = log(1.0+t*t/df) * (df+1)/2.0;
	return exp(a-b-c-d);
}

//poisson distribution
double pois_dist(double x, double m, double)
{
	return gammq(x+1.0, m);
}

//f-distribution
double f_dist(double f, double df1, double df2)
{
	return f > 0.0 ? betai(df2/2.0, df1/2.0, df2/(df2+df1*f)): 1.0;
}

double f_freq(double x, double df1, double df2)
{
	double a, b, c, d;

	a = gammln((df1+df2)/2.0);		b = gammln(df1/2.0) + gammln(df2/2.0);
	c = log(df1/df2) * df1/2.0 + log(x) * (df1/2.0-1.0);
	d = log(1+(df1/df2)*x) * (-(df1+df2)/2.0);
	return exp(a-b+c+d);
}

//---------------------------------------------------------------------------
// The Weibull distribution
//---------------------------------------------------------------------------
double weib_dist(double x, double shape, double scale)
{
	double dn=1.0, sum, term, tmp;

	if(shape <= 0.0 || scale <= 0.0) return HUGE_VAL;
	if(x <= 0.0) return 0.0;
	term = -pow(x/scale, shape);		tmp = fabs(term);
	if(tmp < 2.22e-16) return tmp;
	if (tmp > 0.697) return -exp(term)+1.0;
	x = sum = term;
	do {				//do taylor series
		dn += 1.0 ;		term *= x/dn;		sum += term;
		}while (fabs(term) > fabs(sum) * 2.22e-16) ;
	return -sum;
}

double weib_freq(double x, double shape, double scale)
{
	double tmp1, tmp2;

	if (shape <= 0.0 || scale <= 0.0) return HUGE_VAL;
	if (x < 0) return 0.0;
	if(x > -HUGE_VAL && x < HUGE_VAL) {
		if(x == 0.0 && shape < 1.0) return HUGE_VAL;
		tmp1 = pow(x / scale, shape - 1.0);
		tmp2 = tmp1 * (x / scale);
		return shape * tmp1 * exp(-tmp2) / scale;
		}
	return HUGE_VAL;
}

//---------------------------------------------------------------------------
// The Cauchy (Lorentz) distribution
//---------------------------------------------------------------------------
double cauch_dist(double x, double loc, double scale)
{
	double y;

	if(scale < 0.0) return HUGE_VAL;
	x = (x - loc) / scale;
	if(x > -HUGE_VAL && x < HUGE_VAL) {
		if (fabs(x) > 1.0) {
			y = atan(1.0/x)/_PI;		return (x > 0) ? 1.0-y : -y;
			} 
		else return 0.5 + atan(x)/_PI;
		}
	return HUGE_VAL;
}

double cauch_freq(double x, double loc, double scale)
{
	double y;

	if(scale < 0.0) return HUGE_VAL;
	if(x > -HUGE_VAL && x < HUGE_VAL) {
		y = (x - loc) / scale;
		return 1.0 / (_PI * scale * (1.0 + y*y));
		}
	return HUGE_VAL;
}

//---------------------------------------------------------------------------
// The Logistic distribution
//---------------------------------------------------------------------------
double logis_dist(double x, double loc, double scale)
{
	if(scale < 0.0) return HUGE_VAL;
	x = exp(-(x - loc) / scale);
	if(x > -HUGE_VAL && x < HUGE_VAL) {
		return 1.0/(1.0 + x);
		}
	return HUGE_VAL;
}

double logis_freq(double x, double loc, double scale)
{
	double e, f;

	x = fabs((x - loc) / scale);
	if(x > -HUGE_VAL && x < HUGE_VAL) {
		e = exp(-x);     f = 1.0 + e;	
		return  e / (scale * f*f);
		}
	return HUGE_VAL;
}

//---------------------------------------------------------------------------
// Shapiro-Wilk W test and its significance level
// Algorithm AS 394, 1995, Appl. Statist. 44(4), 547-551
//
static int do_swilk(double (*func)(double, double, double), double p1, double p2, 
	double *x, int n, int n1, int n2, double *a, double *w, double *pw)
{

//initialized data
const static double z90 = 1.2816;		//tinv(0.2, inf)
const static double z95 = 1.6449;		//tinv(0.1, inf)
const static double z99 = 2.3263;		//tinv(.05, inf)
const static double zm = 1.7509;		//(z90 + z95 + z99)/3
const static double zss = 0.56268;
const static double bf1 = 0.8378;
const static double xx90 = 0.556;
const static double xx95 = 0.622;
const static double sqrth = 0.70711;	//sqrt(0.5)
const static double smal = 1.0e-19;		//small value
const static double pi6 = 1.909859;
const static double stqr = 1.047198;	//pi / 3

//polynomial coefficients
static double g[2] = {-2.273, 0.459};
static double c1[6] = {0.0, 0.221157, -0.147981, -2.07119, 4.434685, -2.706056};
static double c2[6] = {0.0, 0.042981, -0.293762, -1.752461, 5.682633, -3.582633};
static double c3[4] = {0.544, -0.39978, 0.025054, -6.714e-4};
static double c4[4] = {1.3822, -0.77857, 0.062767, -0.0020322};
static double c5[4] = {-1.5861, -0.31082, -0.083751, 0.0038915};
static double c6[3] = {-0.4803, -0.082676, 0.0030302};
static double c7[2] = {0.164, 0.533};
static double c8[2] = {0.1736, 0.315};
static double c9[2] = {0.256, -0.00635};

	//local variables
	int i, j, ncens, i1, nn2;
	double zbar, ssassx, summ2, ssumm2, gamma, delta, range;
	double a1, a2, an, bf, ld, m, s, sa, xi, sx, xx, y, w1;
	double fac, asa, an25, ssa, z90f, sax, zfm, z95f, zsd, z99f, rsn, ssx, xsx;

	//parameter adjustment
	--a;

	*pw = 1.0;
	if(*w >= 0.0) *w = 1.0;
	an = (double)(n);			nn2 = n>>1;
	if(n2 < nn2) return 3;
	if(n < 3) return 1;
	// calculate coefficients a[]
	if(true) {
		if(n == 3) a[1] = sqrth;
		else {
			for(i = 1, summ2 = 0.0, an25 = an + 0.25; i <= n2; ++i) {
				a[i] = distinv(func, p1, p2, (i-0.375)/an25, 0);
				summ2 += (a[i] * a[i]);
				}
			summ2 *= 2.0;			ssumm2 = sqrt(summ2);
			rsn = 1.0 / sqrt(an);	a1 = devlpl(c1, 6, rsn) -a[1]/ssumm2;
			//normalize a[]
			if(n > 5) {
				i1 = 3;
				a2 = -a[2] / ssumm2 + devlpl(c2, 6, rsn);
				fac = sqrt((summ2 - 2.0*a[1]*a[1] - 2.0*a[2]*a[2])
					/ (1.0 - 2.0*a1*a1 - 2.0*a2*a2));
				a[2] = a2;
				}
			else {
				i1 = 2;
				fac = sqrt((summ2 -2.0*a[1]*a[1]) / (1.0 - 2.0*a1*a1));
				}
			a[1] = a1;
			for(i = i1; i <= nn2; ++i) a[i] /= -fac;
			} 
		}
	if(n1 < 3) return 1;
	ncens = n - n1;
	if(ncens < 0 || (ncens > 0 && n < 20)) return 4;
	delta = (double)ncens / an;
	if(delta > 0.8) return 5;
	//if w input as negative, calculate significance level of -w
	if(*w < 0.0) { 
		w1 = 1.0 + *w;
		goto sw_prob;
		}
	//check for zero range
	if((range = x[n1-1] -x[0]) < smal) return 6;
	//check for sort order
	xx = x[0]/range;	sx = xx;	sa = -a[1];		j = n -1;
	for(i = 1; i < n1; --j) {
		xi = x[i] / range;			sx += xi;			++i;
		if(i != j) sa += i > j ? a[i < j ? i : j] : -a[i < j ? i : j];
		xx = xi;
		}
	//calculate w statistic as squared correlation between data and coefficients
	sa /= n1;		sx /= n1;		ssa = ssx = sax = 0.0;		j = n -1;
	for(i = 0; i < n1; ++i, --j) {
		if(i > j) asa = a[1+j] - sa;
		else if(i < j) asa = -a[1+i] - sa;
		else asa = -sa;
		xsx = x[i] / range - sx;		ssa += asa * asa;
		ssx += xsx * xsx;				sax += asa * xsx;
		}
	ssassx = sqrt(ssa * ssx);
	w1 = (ssassx - sax) * (ssassx + sax) / (ssa * ssx);
sw_prob:
	*w = 1.0 - w1;			//reduce rounding errors
	if(n == 3) {
		*pw = pi6 * (asin(sqrt(*w)) - stqr);
		return 0;
		}
	y = log(w1);
	xx = log(an);
	if(n <= 11) {
		gamma = devlpl(g, 2, an);
		if(y >= gamma) {
			*pw = smal;		return 0;
			}
		y = -log(gamma - y);		m = devlpl(c3, 4, an);
		s = exp(devlpl(c4, 4, an));
		}
	else {					//n >= 12
		m = devlpl(c5, 4, xx);		s = exp(devlpl(c6, 3, xx));
		}
	//Censoring by proportion  NCENS/N
	if(ncens > 0) {
		ld = -log(delta);			bf = 1.0 + xx * bf1;
		z90f = z90 + bf * pow(devlpl(c7, 2, pow(xx90, xx)), ld);
		z95f = z95 + bf * pow(devlpl(c8, 2, pow(xx95, xx)), ld);
		z99f = z99 + bf * pow(devlpl(c9, 2, xx), ld);
		//Regress z90f ... z99f on normal deviates z90 ... z99
		//   to get pseudo-mean and pseudo-sd of z as the slope and intercept 
		zfm = (z90f + z95f + z99f)/3.0;
		zsd = (z90 * (z90f - zfm) + z95 * (z95f - zfm) + z99 * (z99f - zfm)) / zss;
		zbar = zfm - zsd * zm;		m += zbar * s;		s *= zsd;
		}
	*pw = 1.0 - norm_dist(y, m, s);
	return 0;
}

void swilk1(int n, double *v0, double (*func)(double, double, double), double p1, double p2, 
	bool bsorted, double *w, double *p)
{
	double *v, *a;

	if(!n || !w || !p) return;			*w = *p = 1.0;
	if(!(a = (double*)malloc(n *sizeof(double)))) return;
	if(!bsorted && (v = (double*)memdup(v0, n*sizeof(double), 0)))SortArray(n, v);
	else if(bsorted) v = v0;
	else return;
	if(do_swilk(func, p1, p2, v, n, n, n>>1, a, w, p)){
		//an error occured
		*w = *p = -1.0;
		}
	free(a);	if(v != v0) free(v);
}

//Kolmogorov-Smirnov's test and distribution of D
// (1) Miller L. (1956) Journal of the American Statistical Association.  51: 111-121
// (2) Mises R. (1964) Mathematical Theory of Probability and Statistics (New York: Academic Press)
//     Chapters IX(C) and IX(E)
// (3) Press W.H., Flannery B.P.,Teukolsky S.A., Vetterling W.T. (1988/1989)
//     Numerical Recipes in C, Cambridge University Press, ISBN 0-521-35465-X, pp. 490 ff.
//
double ks_dist(int n, double d)
{
	double j, jn, sum, las, q, r, s, dn = (double)n;

	las = floor(dn - dn * d);
	for (j = sum = 0.0; j <= las; j += 1.0) {
		jn = j / dn;							q = gammln(dn+1) - gammln(j+1) - gammln(dn-j+1.0);
		r = (dn - j) * log( 1 - d - jn );		s = (j - 1.0) * log( d + jn );
		sum += exp(q + r + s);
		}
	return(d*sum);
}

void KolSmir(int n, double *v0, double (*func)(double, double, double), double p1, double p2, 
	bool bsorted, double *d, double *p)
{
	int i;
	double *v, *dev, *x, ff, dt, dt1, dt2;
	double dn = (double)n, f0 = 0.0;

	if(!n || !d || !p) return;			*d = *p = 0.0;
	if(!(dev = (double*)malloc(n*sizeof(double)))) return;
	if(!(x = (double*)malloc(n*sizeof(double)))){
		free(dev);						return;
		}
	if(!bsorted && (v = (double*)memdup(v0, n*sizeof(double), 0)))SortArray(n, v);
	else if(bsorted) v = v0;
	else return;
	for(i = 0, *d = 0.0; i < n; i++) {
		x[i] = (double)(i+1)/dn;		ff = (*func)(v[i], p1, p2); 
		dt1 = fabs(f0-ff);				dt2 = fabs(dev[i] = (f0 = x[i])-ff);
		dt = dt1 > dt2 ? dt1 : dt2;		if(dt > *d) *d = dt;
		}
	free(dev);	free(x);
	*p = ks_dist(n, *d);
	if(v != v0) free(v);
}

//---------------------------------------------------------------------------
// Inverse of statitistical functions:
// funcd supplies the function value fn and the derivative df of the function sf at x
void funcd(double x, double *fn, double *df, double (*sf)(double, double, double), 
		   double df1, double df2, double p)
{
	double y1, y2;

	*fn = (sf)(x, df1, df2);
	if(sf == norm_dist) *df = norm_freq(x, df1,df2);
	else if(sf == chi_dist) *df = -chi_freq(x, df1);
	else if(sf == t_dist) *df = -2.0 * t_freq(x, df1);
	else if(sf == f_dist) *df = -1.0 * f_freq(x, df1, df2);
	else if(sf == lognorm_dist) *df = lognorm_freq(x, df1, df2);
	else if(sf == weib_dist) *df = weib_freq(x, df1, df2);
	else if(sf == cauch_dist) *df = cauch_freq(x, df1, df2);
	else if(sf == logis_dist) *df = logis_freq(x, df1, df2);
	else {		//numerical differentiation
		y1 = (sf)(x * 0.995, df1, df2);		y2 = (sf)(x * 1.005, df1, df2);
		*df = (y2-y1)*100.0/x;
		}
	*fn = *fn - p;
}

//distinv does actual Newton-Raphson root finding
double distinv(double (*sf)(double, double, double), double df1, double df2, double p, double x0)
{
	int i, j;
	double df, df0, adf, dx, f, rtn;

	for(j = 0, rtn = dx = x0; j < 200; j++) {
		for(i = 0, df0 = 0.0; i < 20; i++) {
			funcd(rtn, &f, &df, sf, df1, df2, p);
			if((adf=fabs(df)) > 1.0e-12 || df0 > adf) break;
			rtn += (dx = dx/2.0);				df0 = adf;
			if(i >= 19) return HUGE_VAL;
			}
		dx = f/df*(0.01*(double)(100-j));		rtn -= dx;
		if(fabs(dx) < _PREC && j > 3)return rtn; 
		}
	return HUGE_VAL;
}

//---------------------------------------------------------------------------
//some statistical basics
//do quartiles, median of data
void d_quartile(int n, double *v, double *q1, double *q2, double *q3)
{
	int n2, n3;
	double f1, f2;

	if(!v || n<2) return;
	SortArray(n, v);			n2 = n >> 1;
	if(q1) {
		n3 = n2 >> 1;
		switch(n%4) {
		case 3:		n3 ++;		f1 = 2.0;		f2 = 2.0;		break;
		case 2:		n3 ++;		f1 = 3.0;		f2 = 1.0;		break;
		case 1:		n3 ++;		f1 = 4.0;		f2 = 0.0;		break;
		default:	f1 = 1.0;	f2 = 3.0;						break;
			}
		*q1 = (f1*v[n3-1] + f2*v[n3])/4.0;
		}
	if(q2) {
		if(n & 1) *q2 = v[n2];
		else *q2 = (v[n2-1] + v[n2])/2.0;
		}
	if(q3) {
		n3 = n2 >> 1;
		switch(n%4) {
		case 3:		n3++;		f1 = 2.0;		f2 = 2.0;	break;
		case 2:		f1 = 3.0;	f2 = 1.0;					break;
		case 1:		f1 = 4.0;	f2 = 0.0;					break;
		default:	f1 = 1.0;	f2 = 3.0;					break;
			}
		n3 += n2;
		*q3 = (f2*v[n3-1] + f1*v[n3])/4.0;
		}
}

// statistical basics partly based on
// Davies, J. and Gogh, B. (2000), GSL-1.7 - The GNU scientific library
//
//do variance
double d_variance(int n, double *v, double *mean, double *ss)
{
	int i;
	double d, m, va, e;

	for(i = 0, m = 0.0, d = 1.0; i < n; i++, d += 1.0) {
		m += (v[i] - m)/d;
		}
	if (mean) *mean = m;
	for(i = 0, va = 0.0, d = 1.0; i < n; i++, d += 1.0) {
		e = v[i] - m;		va += (e * e - va)/d;
		}
	if (ss) *ss = va * (double)n;
	return va * ((double)n/((double)(n-1)));
}

//do arithmethic mean
double d_amean(int n, double *v)
{
	int i;
	double d, mean;

	for(i = 0, mean = 0.0, d = 1.0; i < n; i++, d += 1.0) {
		mean += (v[i] - mean)/d;
		}
	return mean;
}


//do geometric mean
double d_gmean(int n, double *v)
{
	int i;
	double sum;

	for(i = 0, sum = 0.0; i < n; i++) {
		if(v[i] <= 0.0) return 0.0;
		sum += log(v[i]);
		}
	return exp(sum/n);
}

//do harmonic mean
double d_hmean(int n, double *v)
{
	int i;

	double sum;

	for(i = 0, sum = 0.0; i < n; i++) {
		if(v[i] == 0.0) return 0.0;
		sum += 1.0/(v[i]);
		}
	return (n/sum);
}

//kurtosis
double d_kurt(int n, double *v)
{
	double sum, avg, sd, tmp, dn = n;
	int i;

	for(i = 0, sum = 0.0; i < n; i++) sum += v[i];
	for(i = 0, avg = sum/dn, sum = 0.0; i < n; i++) sum += (tmp = v[i]-avg) * tmp;
	for(i = 0, sd = sqrt(sum/(dn-1.0)), sum=0.0; i < n; i++) sum += ((tmp = (v[i]-avg)/sd)*tmp*tmp*tmp);
	sum *= ((dn*(dn+1.0))/((dn-1.0)*(dn-2.0)*(dn-3.0)));
	tmp = (3.0 * (dn-1.0) * (dn-1.0))/((dn-2.0)*(dn-3.0));
	return sum - tmp;
}

//skewness
double d_skew(int n, double *v)
{
	double sum, avg, sd, tmp, dn = n;
	int i;

	for(i = 0, avg = 0.0; i < n; i++) avg += ((v[i]-avg)/((double)(i+1)));
	for(i = 0, sum = 0.0; i < n; i++) sum += (tmp = v[i]-avg) * tmp;
	for(i = 0, sd = sqrt(sum/(dn-1.0)), sum=0.0; i < n; i++) sum += ((tmp = (v[i]-avg)/sd)*tmp*tmp);
	return sum * dn/((dn-1.0)*(dn-2.0));
}

//---------------------------------------------------------------------------
// Create a frequency distribution by counting the elements which may be 
// assigned to a bin
double d_classes(DataObj *d, double start, double step, double *v, int nv, char *range)
{
	int i, j, r, c, nc, *f;
	AccRange *ar;

	if(!range || !nv || !v || step <= 0.0 || !(ar = new AccRange(range))) return 0.0;
	if(!(nc = ar->CountItems()) || !ar->GetFirst(&c, &r) || !(f=(int*)calloc(nc, sizeof(int)))) {
		delete ar;				return 0.0;
		}
	for(i = 0; i < nv; i++) {
		j = (int)(floor((v[i] - start)/step));
		if(j < 0) j = 0;		if(j >= nc) j = (nc-1);
		f[j]++;
		}
	for( ; nc > 0 && !(f[nc-1]); nc--);
	for(i = 0; ar->GetNext(&c, &r) && i < nc; i++) {
		d->SetValue(r, c, (double)f[i]);
		}
	free(f);					return ((double)nv);
}

//---------------------------------------------------------------------------
// Pearsons linear correlation
// (1) W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Rcipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 503 ff.
// (2) B. Gough (2000), linear.c, gsl-1.7 the GNU scientific library
double d_pearson(double *x, double *y, int n, char *dest, DataObj *data, double *ra)
{
	int j, r, c;
	double yt, xt, t, df, res[4];
	double syy=0.0, sxy=0.0, sxx=0.0, ay=0.0, ax=0.0;
	AccRange *rD;


	for(j = 0;	j < n; j++) {				// find means
		ax += (x[j] - ax) / (j+1);			ay += (y[j] - ay) / (j+1);
		}
	for(j = 0; j < n; j++) {				// correlation
		xt = x[j] - ax;						yt = y[j] - ay;
		sxx += (xt*xt-sxx) / (j+1);			syy += (yt*yt-syy) / (j+1);
		sxy += (xt*yt-sxy) / (j+1);
		}
	res[0] = sxy/sqrt(sxx*syy);				//pearsons r
	if(dest || ra) {
		res[1] = 0.5 * log((1.0+res[0]+_PREC)/(1.0-res[0]+_PREC));	//Fishers z-transform
		df = n-2;
		t = res[0]*sqrt(df/((1.0-res[0]+_PREC)*(1.0+res[0]+_PREC)));	//Student's t
		res[2] = betai(0.5*df, 0.5, df/(df+t*t));					//probability
		res[3] = n;
		}
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(j = 0; j < 4 && rD->GetNext(&c, &r); j++) {
			data->SetValue(r, c, res[j]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if (ra){
		memcpy(ra, res, 4 * sizeof(double));
		}
	return res[0];
}

//---------------------------------------------------------------------------
// Given an array w, rank returns the rank of v1 in v
// if v1 is not found in v 0 is returned
double d_rank(int n, double *v, double v1)
{
	double *sv;
	int i, j;

	if(!n || !v) return 0.0;		if(n < 2) return 1.0;
	if(!(sv = (double*)memdup(v, n * sizeof(double), 0))) return 0.0;
	SortArray(n, sv);
	for(i = j = 0; i < n; i++) {
		if(v1 == sv[i]) {
			for( ;(i+j)<n; j++) if(sv[i+j] > v1) break;
			free(sv);				return (double)i + 1.0 + (((double)j-1.0)/2.0);
			}
		}
	free(sv);						return 0.0;
}

//---------------------------------------------------------------------------
// Spearman rank-order correlation
// Ref.: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Recipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 507 ff.

//Given a sorted array w, crank replaces the elements by their rank
void crank(int n, double *w0, double *s)
{
	int j=1, ji, jt;
	double t, rank, *w = w0-1;

	*s = 0.0;
	while (j < n) {
		if(w[j+1] != w[j]) {
			w[j] = j;		++j;
			}
		else {
			for(jt = j+1; jt <= n; jt++) if(w[jt] != w[j]) break;
			rank = 0.5 * (j+jt-1);
			for(ji = j; ji <= (jt-1); ji++) w[ji] = rank;
			t = jt -j;		*s += t*t*t -t;				j = jt;
			}
		}
	if(j == n) w[n] = n;
}

//the actual rank correlation
double d_spearman(double *sx, double *sy, int n, char *dest, DataObj *data, double *ra)
{
	int j, r, c;
	double *x, *y, vard, t, sg, sf, fac, en3n, en, df, aved, tmp;
	double res[6];
	AccRange *rD;

	if(!(x = (double*)memdup(sx, n*sizeof(double), 0)) 
		|| !(y = (double*)memdup(sy, n*sizeof(double), 0)))return 0.0;
	SortArray2(n, x, y);			crank(n, x, &sf);
	SortArray2(n, y, x);			crank(n, y, &sg);
	for(j = 0, res[0] = 0.0; j < n; j++) res[0] += ((tmp = (x[j]-y[j]))*tmp);
	en = n;						en3n = en*en*en -en;
	aved = en3n/6.0 - (sf+sg)/12.0;
	fac = (1.0-sf/en3n)*(1.0-sg/en3n);
	vard = ((en-1.0)*en*en*((tmp = (en+1.0))*tmp)/36.0)*fac;
	res[1] = (res[0]-aved)/sqrt(vard);
	res[2] = errfc(fabs(res[1])/_SQRT2);
	res[3] = (1.0-(6.0/en3n)*(res[0]+0.5*(sf+sg)))/fac;
	t = res[3]*sqrt((en-2.0)/((res[3]+1.0)*(1.0-res[3])));
	df = en-2.0;	res[5] = (double)n;
    res[4] = betai(0.5*df, 0.5, df/(df+t*t));
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(j = 0; j < 6 && rD->GetNext(&c, &r); j++) {
			data->SetValue(r, c, res[j]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if(ra) {
		memcpy(ra, res, 6 * sizeof(double));
		}
	free(x);						free(y);
	return res[3];
}

//---------------------------------------------------------------------------
// Kendal's non-parametric correlation
// Ref.: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (1989), 
//    Numerical Recipies in C. The Art of Scientific Computing, 
//    Cambridge University Press, ISBN 0-521-35465, pp. 510 ff.

double d_kendall(double *x, double *y, int n, char *dest, DataObj *data, double *ra)
{
	int j, k, n1, n2, is, r, c;
	double aa, a1, a2, sv, res[4];
	AccRange *rD;

	for (j = n1 = n2 = is = 0; j < (n-1); j++) {
		for(k = j+1; k < n; k++) {
			a1 = x[j] - x[k];		a2 = y[j] - y[k];		aa = a1*a2;
			if(aa != 0.0) {
				n1++;				n2++;
				if (aa > 0.0) is++;
				else is--;
				}
			else {
				if(a1 != 0.0) n1++;	if(a2 != 0.0) n2++;
				}
			}
		}
	res[0] = ((double)is)/(sqrt((double)n1) * sqrt((double)n2));
	sv = (4.0 * ((double)n) + 10.0)/(9.0*((double)n)*((double)(n-1)));
	res[1] = res[0]/sqrt(sv);	res[2] = errfc(fabs(res[1])/_SQRT2);
	res[3] = n;			
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(j = 0; j < 4 && rD->GetNext(&c, &r); j++) {
			data->SetValue(r, c, res[j]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if (ra){
		memcpy(ra, res, 4 * sizeof(double));
		}
	return res[0];
}


//linear regression
double d_regression(double *x, double *y, int n, char *dest, DataObj *data, double *ra)
{
	double sx, sy, dx, dy, sxy, sxx, syy, sdy, df;
	double res[10];		// slope, intercept, mean x, mean y, SE of slope, 
						//   variance(x), variance(y), variance(fit), F of regression, significance
	int i, j, r, c;
	AccRange *rD;

	if(n < 2) return 0.0;
	for(i = 0, 	sx = sy = 0.0; i < n; i++) {
		sx += x[i];			sy += y[i];
		}
	res[2] = sx /n;			res[3] = sy/n;
	sxy = sxx = syy = 0.0;
	for(i = 0; i < n; i++) {
		dx = x[i]-res[2];	dy = y[i]-res[3];
		sxx += (dx*dx);		syy += (dy*dy);		sxy += (dx*dy);
		}
	res[0] = sxy / sxx;		res[1] = res[3] - res[0] * res[2];
	for(i = 0, sdy = 0.0; i < n; i++) {
		dy = y[i] - (res[1] + x[i] *res[0]);
		sdy += (dy * dy);
		}
	sdy = sdy/(n-2);		res[4] = sqrt(sdy/sxx);		df = (n-2);
	res[5] = sxx/(n-1);		res[6] = syy/(n-1);			res[7] = sdy;
	res[8] = sxy/sdy*sxy/sxx;
	res[9] = betai(df/2.0, 0.5, df/(df+res[8]));
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(j = 0; j < 10 && rD->GetNext(&c, &r); j++) {
			data->SetValue(r, c, res[j]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if (ra)	memcpy(ra, res, 10 * sizeof(double));
	return n;
}

//covariance
double d_covar(double *x, double *y, int n, char *dest, DataObj *data)
{
	int i;
	double sx, sy, dx, dy, sxy;

	if(n < 2) return 0.0;
	for(i = 0, 	sx = sy = 0.0; i < n; i++) {
		sx += x[i];			sy += y[i];
		}
	sx /= n;		sy /= n;		sxy = 0.0;
	for(i = 0; i < n; i++) {
		dx = x[i]-sx;		dy = y[i]-sy;
		sxy += (dx*dy - sxy) / (i+1);
		}
	return sxy;
}

//Mann-Whitney U Test
double d_utest(double *x, double *y, int n1, int n2, char *dest, DataObj *data, double *ra)
{
	double *da, *ta, u1, u2, su, su1, ts, dn1 = n1, dn2 = n2;
	double res[9];
	AccRange *rD;
	int i, j, n, r, c;

	if(!x || !y || n1 < 2 || n2 < 2) return 0.0;
	da = (double*)malloc((n = (n1+n2)) * sizeof(double));
	ta = (double*)malloc(n * sizeof(double));
	if(!da || !ta) {
		if(da) free(da);	if(ta) free(ta); return 0.0;
		}
	for(i = 0; i < n1; i++) {
		da[i] = x[i];		ta[i] = 1.0;
		}
	for(j = 0; j < n2; j++) {
		da[i] = y[j];		ta[i++] = 2.0;
		}
	SortArray2(n, da, ta);	crank(n, da, &ts);
	for(i = 0, res[0] = res[1] = 0.0; i < n; i++) {
		if(ta[i] == 1.0) res[0] += da[i];
		else res[1] += da[i];
		}
	free(da);										free(ta);
	u1 = (dn1*dn2 + (dn1*(dn1+1))/2.0) - res[0];	u2 = (dn1*dn2 + ((dn2+1)*dn2)/2.0) - res[1];
	su = sqrt((dn1*dn2*(dn1+dn2+1))/12.0);			res[2] = u2 > u1 ? u2 : u1;
	su1 = ((dn1*dn2)/((dn1+dn2)*(dn1+dn2-1))) * (((dn1+dn2)*(dn1+dn2)*(dn1+dn2)-(dn1+dn2)-ts)/12.0);
	su1 = sqrt(su1);
	res[3] = (res[2] - (n1*n2)/2.0)/su;			res[6] = errfc(res[3]/_SQRT2);
	res[4] = n1;								res[5] = n2;
	res[7] = (res[2] - (n1*n2)/2.0)/su1;		res[8] = errfc(res[7]/_SQRT2);
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(i = 0; i < 9 && rD->GetNext(&c, &r); i++) {
			data->SetValue(r, c, res[i]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if (ra)	memcpy(ra, res, 9 * sizeof(double));
	return res[8];
}

//t-test
double d_ttest(double *x, double *y, int n1, int n2, char *dest, DataObj *data, double *results)
{
	int i, r, c;
	double sx, sy, mx, my, d, df, p;
	double res[9];			// mean1, SD1, n1, mean2, SD2, n2, p if variances equal,
	AccRange *rD;			//    corrected df, corrected p

	d_variance(n1, x, &mx, &sx);		d_variance(n2, y, &my, &sy);
	d = ((sx+sy)/(n1+n2-2)) * ((double)(n1+n2)/(double)(n1*n2));
	d = (mx-my)/sqrt(d);	//Student's t

	//Welch's correction for differences in variance
	df = (sx/(double)n1)*(sx/(double)n1)/(double)(n1+1)+(sy/(double)n2)*(sy/(double)n2)/(double)(n2+1);
	df = (sx/(double)n1+sy/(double)n2)*(sx/(double)n1+sy/(double)n2)/df;
	df -= 2.0;		df = floor(df);

//	an alternative formula for correction
//	p = (sx/(double)n1)*(sx/(double)n1)/(double)(n1-1) + (sy/(double)n2)*(sy/(double)n2)/(double)(n2-1);
//	df = (sx/(double)n1 + sy/(double)n2) * (sx/(double)n1 + sy/(double)n2) / p;

	p = betai(df/2.0, 0.5, (df/(df+d*d)));
	if((dest) && (data) && (rD = new AccRange(dest))) {
		res[0] = mx;	res[1] = sqrt(sx/(double)(n1-1));	res[2] = n1;
		res[3] = my;	res[4] = sqrt(sy/(double)(n2-1));	res[5] = n2;
		res[7] = df;	df = (n1-1) + (n2-1);	res[6] = betai(df/2.0, 0.5, (df/(df+d*d)));
		res[8] = p;
		rD->GetFirst(&c, &r);
		for(i = 0; i < 9 && rD->GetNext(&c, &r); i++) {
			data->SetValue(r, c, res[i]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if(results) {
		results[0] = mx;	results[1] = sqrt(sx/(double)(n1-1));	results[2] = n1;
		results[3] = my;	results[4] = sqrt(sy/(double)(n2-1));	results[5] = n2;
		results[7] = df;	df = (n1-1) + (n2-1);	results[6] = betai(df/2.0, 0.5, (df/(df+d*d)));
		results[8] = p;		results[9] = d;
		}
	return p;
}

//t-test for paired samples
double d_ttest2(double *x, double *y, int n, char *dest, DataObj *data, double *ra)
{
	double sx, sy, mx, my, df, cov, sd, t, p;
	int i, r, c;
	double res[6];			// mean1, SD1, mean2, SD2, n, p 
	AccRange *rD;

	d_variance(n, x, &mx, &sx);		d_variance(n, y, &my, &sy);
	sx = d_variance(n, x, &mx);		sy = d_variance(n, y, &my);
	cov = d_covar(x, y, n, 0L, 0L) * ((double)n/(double)(n-1));
	sd = sqrt((sx+sy-2*cov)/n);
	t = (mx-my)/sd;					df = (n-1);
	p = betai(0.5*df, 0.5, df/(df+t*t));
	res[0] = mx;	res[1] = sqrt(sx);	res[5] = p;
	res[2] = my;	res[3] = sqrt(sy);	res[4] = n;
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(i = 0; i < 6 && rD->GetNext(&c, &r); i++) {
			data->SetValue(r, c, res[i]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if (ra)	memcpy(ra, res, 6 * sizeof(double));
	return p;
}

//f-test
double d_ftest(double *x, double *y, int n1, int n2, char *dest, DataObj *data, double *ra)
{
	int i, r, c;
	double sx, sy, mx, my, d, df1, df2, p;
	double res[6];			// mean1, SD1, n1, mean2, SD2, n2
	AccRange *rD;

	for(i=0, sx = 0.0; i < n1; sx += x[i], i++);				mx = sx/n1;
	for(i=0, sy = 0.0; i < n2; sy += y[i], i++);				my = sy/n2;
	for(i=0, sx = 0.0; i < n1; sx += ((d=(x[i]-mx))*d), i++);	sx /= (n1-1);
	for(i=0, sy = 0.0; i < n2; sy += ((d=(y[i]-my))*d), i++);	sy /= (n2-1);
	if(sx > sy) {
		d = sx/sy;		df1 = n1-1;		df2 = n2-1;
		}
	else {
		d = sy/sx;		df1 = n2-1;		df2 = n1-1;
		}
	p = 2.0 * betai(df2/2.0, df1/2.0, df2/(df2+df1*d));
	if(p > 1.0) p = 2.0-p;
	res[0] = mx;	res[1] = sqrt(sx);	res[2] = n1;
	res[3] = my;	res[4] = sqrt(sy);	res[5] = n2;
	if((dest) && (data) && (rD = new AccRange(dest))) {
		rD->GetFirst(&c, &r);
		for(i = 0; i < 6 && rD->GetNext(&c, &r); i++) {
			data->SetValue(r, c, res[i]);
			}
		data->Command(CMD_UPDATE, 0L, 0L);
		delete rD;
		}
	if (ra)	memcpy(ra, res, 6 * sizeof(double));
	return p;
}
//---------------------------------------------------------------------------
// Simple one way anova
//---------------------------------------------------------------------------
bool do_anova1(int n, int *nv, double **vals, double **res_tab, double *gm, double **means, double **ss)
{
	int i, j, ntot;
	double tmp, *csums, *css, ssa, ssw, sst, mtot, d;

	if(!(csums = (double*)calloc(n+1, sizeof(double)))
		|| !(css = (double*)calloc(n+1, sizeof(double)))) return false;

	for(i = ntot = 0, mtot = 0.0, d = 1.0; i< n; i++){
		for(j = 0, csums[i] = 0.0, tmp = 1.0; j < nv[i]; j++, d+=1.0, tmp +=1.0) {
			mtot += (vals[i][j] - mtot)/d;		
			csums[i] += (vals[i][j] -csums[i])/tmp;
			}
		ntot += nv[i];
		}
	for(i = 0; i < n; i++) {
		for(j = 0, css[i] = 0.0; j < nv[i]; j++) {
			tmp = vals[i][j] - csums[i];	css[i] += (tmp*tmp);
			}
		}
	for(i = 0, ssa = ssw = sst = 0.0;  i < n; i++) {
		tmp =(csums[i] - mtot);		ssa += (tmp*tmp) * ((double)nv[i]);
		ssw += css[i];
		}
	sst = ssa + ssw;
	res_tab[0][0] = n - 1;				res_tab[1][0] = ntot - n;
	res_tab[2][0] = ntot -1;			res_tab[0][1] = ssa;
	res_tab[1][1] = ssw;				res_tab[2][1] = sst;
	res_tab[0][2] = ssa/res_tab[0][0];	res_tab[1][2] = ssw/res_tab[1][0];
	res_tab[0][3] = res_tab[0][2]/res_tab[1][2];
	res_tab[0][4] = f_dist(res_tab[0][3], res_tab[0][0], res_tab[1][0]);
	if(gm) *gm = mtot;
	if(means) *means = csums;			else free(csums);
	if(ss) *ss = css;					else free(css);
	return true;
}

//---------------------------------------------------------------------------
// Bartlett's Test for homogeneity of variances
// RR Sokal & FJ Rohlf: Biometry, 3rd ed., pp. 398 ff.
//---------------------------------------------------------------------------
bool bartlett(int n, int *nc, double *ss, double *chi2)
{
	int i, sdf, df;
	double mss, mlss, *lnss, cf;

	if(!n || !nc || !ss || !chi2) return false;
	if(!(lnss = (double*)malloc(n * sizeof(double))))return false;
	for(i = sdf = 0, mss = mlss = cf = 0.0; i < n; i++) {
		sdf += (df = nc[i]-1);				lnss[i] = log(ss[i]);
		mss += (ss[i] * ((double)df));		mlss += (lnss[i] * ((double)df)); 
		cf += (1.0/((double)df));
		}
	*chi2 = ((double)sdf) * log(mss/((double)sdf)) - mlss;
	cf -= (1.0/((double)sdf));				cf = 1.0 + cf/(3.0 * ((double)(n-1)));
	*chi2 /= cf;
	// P = chi_dist(*chi2, n-1, 0);
	free(lnss);		return true;
}
//---------------------------------------------------------------------------
// Leven's Test for homogeneity of variances
//---------------------------------------------------------------------------
bool levene(int type, int n, int *nv, double *means, double **vals, double *F, double *P)
{
	int i, j;
	bool bRet = false;
	double cm, **res_tab, **cols;
	
	if(!n || !nv || !means || !vals) return false;
	//setup matrix for results
	if((res_tab = (double**)calloc(3, sizeof(double*)))
		&& (res_tab[0] = (double*) malloc(5*sizeof(double)))
		&& (res_tab[1] = (double*) malloc(5*sizeof(double)))
		&& (res_tab[2] = (double*) malloc(5*sizeof(double)))
		&& (cols = (double**)calloc(n+1, sizeof(double*)))) bRet = true;
	//allocate mem for data
	for(i = 0; bRet && i<n; i++) {
		if(!(cols[i]=(double*)malloc((nv[i]+1)*sizeof(double)))) bRet = false;
		}
	//data are absolute differences to mean ...
	for(i = 0, cm = 0.0; bRet && i < n; i++) {
		switch(type) {
			case 1:			//use means
				cm = means[i];								break;
			case 2:			//use medians
				d_quartile(nv[i], vals[i], 0L, &cm, 0L);	break;
			}
		for(j = 0; j < nv[i]; j++) {
			cols[i][j] = vals[i][j] > cm ? vals[i][j] - cm : cm - vals[i][j];
			}
		}
	//Levene's test statistic is based on ANOVA of the differences
	if(bRet && (bRet = do_anova1(n, nv, cols, res_tab, 0L, 0L, 0L))){
		if(F) *F = res_tab[0][3];				if(P) *P = res_tab[0][4];
		}
	//clean up
	if(bRet) {
		for(i = 0; i < n; i++) if(cols[i]) free(cols[i]);
		for(i = 0; i < 3; i++) if(res_tab[i]) free(res_tab[i]);
		free(cols);								free(res_tab);	
		}
	return bRet;
}

//---------------------------------------------------------------------------
// Modules from the R-project
//
//---------------------------------------------------------------------------
#define M_1_SQRT_2PI	0.398942280401432677939946059934	/* 1/sqrt(2pi) */
/*
 *  Copyright (C) 1998       Ross Ihaka
 *  Copyright (C) 2000--2005 The R Development Core Team
 *
 *  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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 *  DESCRIPTION
 *    Computes the probability that the maximum of rr studentized
 *    ranges, each based on cc means and with df degrees of freedom
 *    for the standard error, is less than q.
 *    The algorithm is based on that of the reference.
 *
 *  REFERENCE
 *    Copenhaver, Margaret Diponzio & Holland, Burt S.
 *    Multiple comparisons of simple effects in
 *    the two-way analysis of variance with fixed effects.
 *    Journal of Statistical Computation and Simulation,
 *    Vol.30, pp.1-15, 1988.
 */

double wprob(double w, double rr, double cc)
{
/*  wprob() :

	This function calculates probability integral of Hartley's
	form of the range.

	w     = value of range
	rr    = no. of rows or groups
	cc    = no. of columns or treatments
	ir    = error flag = 1 if pr_w probability > 1
	pr_w = returned probability integral from (0, w)

	program will not terminate if ir is raised.

	bb = upper limit of legendre integration
	iMax = maximum acceptable value of integral
	nleg = order of legendre quadrature
	ihalf = int ((nleg + 1) / 2)
	wlar = value of range above which wincr1 intervals are used to
	       calculate second part of integral,
	       else wincr2 intervals are used.
	C1, C2, C3 = values which are used as cutoffs for terminating
           or modifying a calculation.
	xleg = legendre 12-point nodes
	aleg = legendre 12-point coefficients
 */
#define nleg	12
#define ihalf	6

    /* looks like this is suboptimal for double precision.
       (see how C1-C3 are used) <MM> */
    /* const double iMax  = 1.; not used if = 1*/
    const static double C1 = -30.0, C2 = -50.0, C3 = 60.;
    const static double bb = 8.0, wlar = 3.0, wincr1 = 2.0, wincr2 = 3.;
    const static double xleg[ihalf] = {	0.981560634246719250690549090149,
	0.904117256370474856678465866119,	0.769902674194304687036893833213,
	0.587317954286617447296702418941,	0.367831498998180193752691536644,
	0.125233408511468915472441369464};
    const static double aleg[ihalf] = {	0.047175336386511827194615961485,
	0.106939325995318430960254718194,	0.160078328543346226334652529543,
	0.203167426723065921749064455810,	0.233492536538354808760849898925,
	0.249147045813402785000562436043};
    double a, ac, pr_w, b, binc, blb, bub, c, cc1, einsum, elsum,
		pminus, pplus, qexpo, qsqz, rinsum, wi, wincr, xx;
    int j, jj;

    qsqz = w * 0.5;

    // if w >= 16 then the integral lower bound (occurs for c=20)
    // is 0.99999999999995 so return a value of 1
	if (qsqz >= bb)	return 1.0;

	// find (f(w/2) - 1) ^ cc
    // (first term in integral of hartley's form). 
	pr_w = 2.0 * norm_dist(qsqz, 0.0, 1.0) -1.0;
    // if pr_w ^ cc < 2e-22 then set pr_w = 0 
    if (pr_w >= exp(C2 / cc)) pr_w = pow(pr_w, cc);
    else pr_w = 0.0;
    // if w is large then the second component of the
    // integral is small, so fewer intervals are needed.
    if (w > wlar) wincr = wincr1;
    else wincr = wincr2;

    /* find the integral of second term of hartley's form */
    /* for the integral of the range for equal-length */
    /* intervals using legendre quadrature.  limits of */
    /* integration are from (w/2, 8).  two or three */
    /* equal-length intervals are used. */
    /* blb and bub are lower and upper limits of integration. */
    blb = qsqz;			    binc = (bb - qsqz) / wincr;
    bub = blb + binc;	    einsum = 0.0;

    // integrate over each interval
    cc1 = cc - 1.0;
    for (wi = 1; wi <= wincr; wi++) {
		elsum = 0.0;		a = 0.5 * (bub + blb);
		// legendre quadrature with order = nleg
		b = 0.5 * (bub - blb);
		for (jj = 1; jj <= nleg; jj++) {
			if (ihalf < jj) {
				j = (nleg - jj) + 1;		xx = xleg[j-1];
				}
			else {
				j = jj;						xx = -xleg[j-1];
				}
			c = b * xx;					    ac = a + c;
			// if exp(-qexpo/2) < 9e-14, then doesn't contribute to integral
			if ((qexpo = ac * ac) > C3) break;
			pplus = 2.0 * norm_dist(ac, 0.0, 1.0);    pminus= 2.0 * norm_dist(ac, w, 1.0);
			// if rinsum ^ (cc-1) < 9e-14, then doesn't contribute to integral
			rinsum = (pplus * 0.5) - (pminus * 0.5);
			if (rinsum >= exp(C1 / cc1)) {
				rinsum = (aleg[j-1] * exp(-(0.5 * qexpo))) * pow(rinsum, cc1);
				elsum += rinsum;
				}
			}
		elsum *= (((2.0 * b) * cc) * M_1_SQRT_2PI);
		einsum += elsum;		blb = bub;			bub += binc;
		}
	// if pr_w ^ rr < 9e-14, then return 0 */
	pr_w = einsum + pr_w;
	if (pr_w <= exp(C1 / rr))return 0.;
    pr_w = pow(pr_w, rr);
 	return pr_w < 1.0 ? pr_w : 1.0;
}

double ptukey(double q, double rr, double cc, double df, int lower_tail, int log_p)
{
/* 	q = value of studentized range
	rr = no. of rows or groups
	cc = no. of columns or treatments
	df = degrees of freedom of error term
	ir[0] = error flag = 1 if wprob probability > 1
	ir[1] = error flag = 1 if qprob probability > 1

	All references in wprob to Abramowitz and Stegun
	are from the following reference:
		Abramowitz, Milton and Stegun, Irene A.
		Handbook of Mathematical Functions.
		New York:  Dover publications, Inc. (1970).
	All constants taken from this text are given to 25 significant digits.

	nlegq = order of legendre quadrature
	ihalfq = int ((nlegq + 1) / 2)
	eps = max. allowable value of integral
	eps1 & eps2 = values below which there is no contribution to integral.

	d.f. <= dhaf:	integral is divided into ulen1 length intervals.  else
	d.f. <= dquar:	integral is divided into ulen2 length intervals.  else
	d.f. <= deigh:	integral is divided into ulen3 length intervals.  else
	d.f. <= dlarg:	integral is divided into ulen4 length intervals.

	d.f. > dlarg:	the range is used to calculate integral.

	xlegq = legendre 16-point nodes
	alegq = legendre 16-point coefficients

	The coefficients and nodes for the legendre quadrature used in
	qprob and wprob were calculated using the algorithms found in:
		Stroud, A. H. and Secrest, D.,	Gaussian Quadrature Formulas.
		Englewood Cliffs, New Jersey:  Prentice-Hall, Inc, 1966.

	All values matched the tables (provided in same reference)
	to 30 significant digits.

	f(x) = .5 + erf(x / sqrt(2)) / 2      for x > 0
	f(x) = erfc( -x / sqrt(2)) / 2	      for x < 0
	where f(x) is standard normal c. d. f.

	if degrees of freedom large, approximate integral with range distribution.
 */
#define nlegq	16
#define ihalfq	8

/*  const double eps = 1.0; not used if = 1 */
    const static double eps1 = -30.0, eps2 = 1.0e-14;
    const static double dhaf  = 100.0, dquar = 800.0, deigh = 5000.0, dlarg = 25000.0;
    const static double ulen1 = 1.0, ulen2 = 0.5, ulen3 = 0.25, ulen4 = 0.125;
    const static double xlegq[ihalfq] = { 0.989400934991649932596154173450,
	0.944575023073232576077988415535, 0.865631202387831743880467897712,
	0.755404408355003033895101194847, 0.617876244402643748446671764049,
	0.458016777657227386342419442984, 0.281603550779258913230460501460,
	0.950125098376374401853193354250e-1};
    const static double alegq[ihalfq] = {0.271524594117540948517805724560e-1,
	0.622535239386478928628438369944e-1, 0.951585116824927848099251076022e-1,
	0.124628971255533872052476282192, 0.149595988816576732081501730547,
	0.169156519395002538189312079030, 0.182603415044923588866763667969,
	0.189450610455068496285396723208};
    double ans, f2, f21, f2lf, ff4, otsum, qsqz, rotsum, t1, twa1, ulen, wprb;
    int i, j, jj;

    if (df > dlarg)	return wprob(q, rr, cc);
    f2 = df * 0.5;								// calculate leading constant
    f2lf = ((f2 * log(df)) - (df * log(2.0))) - gammln(f2);
    f21 = f2 - 1.0;
    // integral is divided into unit, half-unit, quarter-unit, or eighth-unit length intervals 
	//    depending on the value of the degrees of freedom.
    ff4 = df * 0.25;
    if	    (df <= dhaf)	ulen = ulen1;
    else if (df <= dquar)	ulen = ulen2;
    else if (df <= deigh)	ulen = ulen3;
    else					ulen = ulen4;
    f2lf += log(ulen);
    for (i = 1, ans = 0.0; i <= 50; i++) {		// integrate over each subinterval
		otsum = 0.0;
		// legendre quadrature with order = nlegq, nodes (stored in xlegq) are symmetric around zero.
		twa1 = (2 * i - 1) * ulen;
		for (jj = 1; jj <= nlegq; jj++) {
			if (ihalfq < jj) {
				j = jj - ihalfq - 1;
				t1 = (f2lf + (f21 * log(twa1 + (xlegq[j] * ulen)))) - (((xlegq[j] * ulen) + twa1) * ff4);
				} 
			else {
				j = jj - 1;
				t1 = (f2lf + (f21 * log(twa1 - (xlegq[j] * ulen)))) + (((xlegq[j] * ulen) - twa1) * ff4);
				}
			if (t1 >= eps1) {			// if exp(t1) < 9e-14, then doesn't contribute to integral 
				if (ihalfq < jj) qsqz = q * sqrt(((xlegq[j] * ulen) + twa1) * 0.5);
				else qsqz = q * sqrt(((-(xlegq[j] * ulen)) + twa1) * 0.5);
				wprb = wprob(qsqz, rr, cc);		// call wprob to find integral of range portion
				rotsum = (wprb * alegq[j]) * exp(t1);			otsum += rotsum;
				}
			}									// end legendre integral for interval i
		// If integral for interval i < 1e-14, then stop. However, in order to avoid small area 
		//    under left tail, at least  1 / ulen  intervals are calculated.
		if (i * ulen >= 1.0 && otsum <= eps2) break;
		ans += otsum;							//end of interval i 
		}
	return ans > 1.0 ? 1.0 : ans;
 }

 /*
 *  Copyright (C) 1998 	     Ross Ihaka
 *  Copyright (C) 2000--2005 The R Development Core Team
 *  based in part on AS70 (C) 1974 Royal Statistical Society
 *
 *  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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 *  DESCRIPTION
 *	Computes the quantiles of the maximum of rr studentized
 *	ranges, each based on cc means and with df degrees of freedom
 *	for the standard error, is less than q.
 *	The algorithm is based on that of the reference.
 *
 *  REFERENCE
 *	Copenhaver, Margaret Diponzio & Holland, Burt S., Multiple comparisons of simple
 *	effects in the two-way analysis of variance with fixed effects.
 *	Journal of Statistical Computation and Simulation, Vol.30, pp.1-15, 1988.
 */

/* qinv() :
 *	this function finds percentage point of the studentized range
 *	which is used as initial estimate for the secant method.
 *	function is adapted from portion of algorithm as 70
 *	from applied statistics (1974) ,vol. 23, no. 1
 *	by odeh, r. e. and evans, j. o.
 *	  p = percentage point
 *	  c = no. of columns or treatments
 *	  v = degrees of freedom
 *	  qinv = returned initial estimate
 *	vmax is cutoff above which degrees of freedom
 *	is treated as infinity.
 */

static double qinv(double p, double c, double v)
{
    const static double p0 = 0.322232421088, q0 = 0.993484626060e-01;
    const static double p1 = -1.0, q1 = 0.588581570495;
    const static double p2 = -0.342242088547, q2 = 0.531103462366;
    const static double p3 = -0.204231210125, q3 = 0.103537752850;
    const static double p4 = -0.453642210148e-04, q4 = 0.38560700634e-02;
    const static double c1 = 0.8832, c2 = 0.2368, c3 = 1.214, c4 = 1.208, c5 = 1.4142;
    const static double vmax = 120.0;
    double ps, q, t, yi;

    ps = 0.5 - 0.5 * p;
    yi = sqrt (log (1.0 / (ps * ps)));
    t = yi + (((( yi * p4 + p3) * yi + p2) * yi + p1) * yi + p0)
	   / (((( yi * q4 + q3) * yi + q2) * yi + q1) * yi + q0);
    if (v < vmax) t += (t * t * t + t) / v / 4.0;
    q = c1 - c2 * t;
    if (v < vmax) q += -c3 / v + c4 * t / v;
    return t * (q * log (c - 1.0) + c5);
}

/*
 *  Copenhaver, Margaret Diponzio & Holland, Burt S.
 *  Multiple comparisons of simple effects in
 *  the two-way analysis of variance with fixed effects.
 *  Journal of Statistical Computation and Simulation,
 *  Vol.30, pp.1-15, 1988.
 *
 *  Uses the secant method to find critical values.
 *
 *  p = confidence level (1 - alpha)
 *  rr = no. of rows or groups
 *  cc = no. of columns or treatments
 *  df = degrees of freedom of error term
 *
 */
double qtukey(double p, double rr, double cc, double df, int lower_tail, int log_p)
{
    const int maxiter = 50;
    double ans = HUGE_VAL, valx0, valx1, x0, x1;
    int iter;

    // df must be > 1 ; there must be at least two values 
	if(p >= 1.0 || df < 2 || rr < 1 || cc < 2) return HUGE_VAL;
	if(p < 0.0) p = 0.0;
    x0 = qinv(p, cc, df);									// Initial value
    valx0 = ptukey(x0, rr, cc, df, true, false) - p;		// Find prob(value < x0)
    // Find the second iterate and prob(value < x1). If the first iterate has probability value 
    // exceeding p then second iterate is 1 less than first iterate; otherwise it is 1 greater.
	x1 = valx0 > 0.0 ? (x1 = x0 > 1.0 ? x0-1.0 : 0.0) : (x0 + 1.0);
    valx1 = ptukey(x1, rr, cc, df, true, false) - p;
    for(iter=1; iter < maxiter ; iter++) {					// Iterate
		ans = x1 - ((valx1 * (x1 - x0)) / (valx1 - valx0));
		valx0 = valx1;		x0 = x1;
		if (ans < 0.0) {									// New iterate must be >= 0
			ans = 0.0;			valx1 = -p;
			}
		valx1 = ptukey(ans, rr, cc, df, true, false) - p;	//  Find prob(value < new iterate)
		x1 = ans;
		if (fabs(x1 - x0) < _PREC)	return ans;				// Convergence ?
		}
    //The process did not converge in 'maxiter' iterations 
    return ans;
}
//---------------------------------------------------------------------------
// END Modules from the R-project


//---------------------------------------------------------------------------
// Calendar, Date- and Time functions
// The following characters are used as format specifiers in a format string,
//    all other characters are either ignored or copyied to the output
//
//    Y   four digits year               y    two digits year
//    X   month's full name              x    three character month name
//    Z   two digits day of month        z    same as Z but no leading zero
//    V   two digit month number         v    number of month
//    W   single letter month
//    D   full name of day               d    three characters for day name
//    E   two digits weekday             e    one or two digits weekday
//    F   single character day name
//    H   two digits for hours           h    hours with no leading zero
//    M   two digits for minutes         m    minutes with no leading zero
//    S   two digits for seconds         s    seconds with no leading zero
//    T   two digits seconds, two dec.   t    same as T but no leading zero
//    U   full precision seconds

static char *dt_month[] = {"January", "February", "March", "April", "May", "June",
	"July", "August", "September", "October", "November", "December"};

static char *dt_months[] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug",
	"Sep", "Oct", "Nov", "Dec"};

static int dt_monthl[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

static char *dt_day[] = {"Sunday", "Monday", "Tuesday", "Wednesday", "Thursday",
	"Friday", "Saturday"};

static char *dt_days[] = {"Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"};

static bool leapyear(int year) {
	return ((year % 4 == 0 && year % 100 != 0) || year % 400 == 0);
}

int year2aday(int y)
{
	int aday, y1;

	y1 = y - 1900;
	aday = y1 * 365;
	aday += ((y1-1) >> 2 );
	aday -= (y1 / 100);
	aday += ((y/400)-4);
	return aday;
}

static void set_dow(rlp_datetime *dt)
{
	dt->dow = (dt->aday %7)+1;
}

void add_date(rlp_datetime *base, rlp_datetime *inc)
{
	int i, dom;

	if(base) {
		if(base->month < 1) base->month = 1;
		if(inc) {
			base->seconds += inc->seconds;
			if(base->seconds >= 60.0) {
				base->minutes++;		base->seconds -= 60.0;
				}
			base->minutes += inc->minutes;
			if(base->minutes >= 60) {
				base->hours++;			base->minutes -= 60;
				}
			base->hours += inc->hours;
			if(base->hours >= 24) {
				base->dom++;			base->hours -= 24;
				}
			base->year += inc->year;	base->dom += inc->dom;
			base->month += inc->month;
			}
		dom = dt_monthl[base->month-1];
		if(leapyear(base->year) && base->month == 2) dom = 29;
		if(base->dom > dom) {
			base->month++;			base->dom -= dom;
			}
		if(base->month > 12) {
			base->year++;			base->month -= 12;
			}
		base->aday = year2aday(base->year);
		for(i = base->doy = 0; i < (base->month-1); i++) {
			dom = dt_monthl[i];
			if(i == 1 && leapyear(base->year)) dom = 29;
			base->doy += dom;
			}
		base->doy += base->dom;
		base->aday += base->doy;	set_dow(base);
		}
}

static int parse_date (rlp_datetime *dt, char *src, char *fmt)
{
	int i, j, k;
	char tmp_str[10];

	if(!src || !src[0] || !fmt || !fmt[0]) return 0;
	if(*src == '\'') src++;
	for(i = j = 0; fmt[i] && src[j]; i++) {
		switch (fmt[i]) {
		case 'Y':		case 'y':			// year is numeric
			if(j && src[j] == '-' || src[j] == '/' || src[j] == '.') j++;
#ifdef USE_WIN_SECURE
			if(sscanf_s(src+j, "%d", &dt->year)) {
#else
			if(sscanf(src+j, "%d", &dt->year)) {
#endif
				if(dt->year < 0) return 0;
				while(isdigit(src[j])) j++;
				if(dt->year<60) dt->year += 2000;
				else if(dt->year <99) dt->year += 1900;
				}
			else return 0;
			break;
		case 'X':		case 'x':			// month can be text
			if(j && src[j] == '-' || src[j] == '/' || src[j] == '.') j++;
			tmp_str[0] = toupper(src[j]);
			tmp_str[1] = tolower(src[j+1]);
			tmp_str[2] = tolower(src[j+2]);
			tmp_str[3] = 0;
			for(k = dt->month = 0; k < 12; k++) {
				if(0 == strcmp(tmp_str,dt_months[k])) {
					dt->month = k+1;			break;
					}
				}
			if(dt->month) while(isalpha(src[j])) j++;
			else return 0;
			break;
		case 'V':		case 'v':			//    or numeric
			if(j && src[j] == '-' || src[j] == '/' || src[j] == '.') j++;
#ifdef USE_WIN_SECURE
			if(sscanf_s(src+j, "%d", &dt->month)) {
#else
			if(sscanf(src+j, "%d", &dt->month)) {
#endif
				if(dt->month <= 0 || dt->month > 12) return 0;
				j++;				if(isdigit(src[j])) j++;
				}
			else return 0;
			break;
		case 'Z':		case 'z':			// day of month is numeric
			if(j && src[j] == '-' || src[j] == '/' || src[j] == '.') j++;
#ifdef USE_WIN_SECURE
			if(sscanf_s(src+j, "%d", &dt->dom)) {
#else
			if(sscanf(src+j, "%d", &dt->dom)) {
#endif
				if(dt->dom <= 0 || dt->dom > 31) return 0;
				j++;				if(isdigit(src[j])) j++;
				}
			else return 0;
			break;
		case 'H':		case 'h':			// hours are numeric
#ifdef USE_WIN_SECURE
			if(sscanf_s(src+j, "%2d", &dt->hours)) {
#else
			if(sscanf(src+j, "%2d", &dt->hours)) {
#endif
				if(dt->hours < 0 || dt->hours > 23) return 0;
				j++;				if(isdigit(src[j])) j++;
				}
			else return 0;
			break;
		case 'M':		case 'm':			// minutes are numeric
			if(j && src[j] == ' ' || src[j] == ':') j++;
#ifdef USE_WIN_SECURE
			if(sscanf_s(src+j, "%2d", &dt->minutes)) {
#else
			if(sscanf(src+j, "%2d", &dt->minutes)) {
#endif
				if(dt->minutes < 0 || dt->minutes >= 60) return 0;
				j++;				if(isdigit(src[j])) j++;
				}
			else return 0;
			break;
		case 'S':		case 's':			// seconds are numeric
		case 'T':		case 't':
			if(j && src[j] == ' ' || src[j] == ':') j++;
#ifdef USE_WIN_SECURE
			if(sscanf_s(src+j, "%lf", &dt->seconds)) {
#else
			if(sscanf(src+j, "%lf", &dt->seconds)) {
#endif
				if(dt->seconds < 0.0 || dt->seconds >= 60.0) return 0;
				while(isdigit(src[j]) || src[j] == '.') j++;
				}
			else return 0;
			dt->seconds += 1.0e-12;
			break;
		default:
			if(fmt[i] && fmt[i] == src[j]) j++;
			}
		}
	if(dt->year && dt->month && dt->dom) {
		for(dt->doy = 0, i = dt->month-2; i >= 0; i--) {
			if(i == 1) dt->doy += leapyear(dt->year) ? 29 : 28;
			else dt->doy += dt_monthl[i]; 
			}
		dt->doy += dt->dom;
		if(dt->year >= 1900) dt->aday = year2aday(dt->year);
		dt->aday += dt->doy;
		}
	return j;
}

char *date2text(rlp_datetime *dt, char *fmt)
{
	static char res[80];
	int i, pos;
	double secs;

	res[0] = 0;
	if(!fmt || !fmt[0] || !dt) return res;
	set_dow(dt);
	secs = dt->seconds;
	if (secs > 59.4999) secs = 59.4999;
	for(pos = i = 0; fmt[i] && pos < 70; i++) {
#ifdef USE_WIN_SECURE
		switch(fmt[i]) {
		case 'Y':
			if(dt->year) pos+=sprintf_s(res+pos, 80-pos, "%4d", dt->year);
			else pos += sprintf_s(res+pos, 80-pos, "####");			break;
		case 'y':
			if(dt->year) pos+=sprintf_s(res+pos, 80-pos, "%02d", (dt->year %100));
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'Z':
			if(dt->dom) pos+=sprintf_s(res+pos, 80-pos, "%02d", dt->dom);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'z':
			if(dt->dom) pos+=sprintf_s(res+pos, 80-pos, "%d", dt->dom);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'X':
			if(dt->month >0 && dt->month < 13) pos+=sprintf_s(res+pos, 80-pos, "%s", dt_month[dt->month-1]);
			else pos += sprintf_s(res+pos, 80-pos, "###");			break;
		case 'x':
			if(dt->month >0 && dt->month < 13) pos+=sprintf_s(res+pos, 80-pos, "%s", dt_months[dt->month-1]);
			else pos += sprintf_s(res+pos, 80-pos, "###");			break;
		case 'V':
			if(dt->month >0 && dt->month < 13) pos+=sprintf_s(res+pos, 80-pos, "%02d", dt->month);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'v':
			if(dt->month >0 && dt->month < 13) pos+=sprintf_s(res+pos, 80-pos, "%d", dt->month);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'W':
			if(dt->month >0 && dt->month < 13) pos+=sprintf_s(res+pos, 80-pos, "%c", dt_month[dt->month-1][0]);
			else pos += sprintf_s(res+pos, 80-pos, "#");			break;
		case 'D':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf_s(res+pos, 80-pos, "%s", dt_day[dt->dow-1]);
			else pos += sprintf_s(res+pos, 80-pos, "###");			break;
		case 'd':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf_s(res+pos, 80-pos, "%s", dt_days[dt->dow-1]);
			else pos += sprintf_s(res+pos, 80-pos, "###");			break;
		case 'E':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf_s(res+pos, 80-pos, "%02d", dt->dow);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'e':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf_s(res+pos, 80-pos, "%d", dt->dow);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'F':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf_s(res+pos, 80-pos, "%c", dt_day[dt->dow-1][0]);
			else pos += sprintf_s(res+pos, 80-pos, "#");			break;
		case 'H':
			if(dt->hours >=0 && dt->hours < 24) pos+=sprintf_s(res+pos, 80-pos, "%02d", dt->hours);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'h':
			if(dt->hours >=0 && dt->hours < 24) pos+=sprintf_s(res+pos, 80-pos, "%d", dt->hours);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'M':
			if(dt->minutes >=0 && dt->minutes < 60) pos+=sprintf_s(res+pos, 80-pos, "%02d", dt->minutes);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'm':
			if(dt->minutes >=0 && dt->minutes < 60) pos+=sprintf_s(res+pos, 80-pos, "%d", dt->minutes);
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'S':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf_s(res+pos, 80-pos, "%02d", iround(secs));
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 's':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf_s(res+pos, 80-pos, "%d", iround(secs));
			else pos += sprintf_s(res+pos, 80-pos, "##");			break;
		case 'T':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf_s(res+pos, 80-pos, "%02.2lf", dt->seconds);
			else pos += sprintf_s(res+pos, 80-pos, "##.##");		break;
		case 't':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf_s(res+pos, 80-pos, "%.2lf", dt->seconds);
			else pos += sprintf_s(res+pos, 80-pos, "##.##");		break;
		default:
			pos += sprintf_s(res+pos, 80-pos, "%c", fmt[i]);		break;
			}
#else
		switch(fmt[i]) {
		case 'Y':
			if(dt->year) pos+=sprintf(res+pos, "%4d", dt->year);
			else pos += sprintf(res+pos, "####");		break;
		case 'y':
			if(dt->year) pos+=sprintf(res+pos, "%02d", (dt->year %100));
			else pos += sprintf(res+pos, "##");			break;
		case 'Z':
			if(dt->dom) pos+=sprintf(res+pos, "%02d", dt->dom);
			else pos += sprintf(res+pos, "##");			break;
		case 'z':
			if(dt->dom) pos+=sprintf(res+pos, "%d", dt->dom);
			else pos += sprintf(res+pos, "##");			break;
		case 'X':
			if(dt->month >0 && dt->month < 13) pos+=sprintf(res+pos, "%s", dt_month[dt->month-1]);
			else pos += sprintf(res+pos, "###");		break;
		case 'x':
			if(dt->month >0 && dt->month < 13) pos+=sprintf(res+pos, "%s", dt_months[dt->month-1]);
			else pos += sprintf(res+pos, "###");		break;
		case 'V':
			if(dt->month >0 && dt->month < 13) pos+=sprintf(res+pos, "%02d", dt->month);
			else pos += sprintf(res+pos, "##");			break;
		case 'v':
			if(dt->month >0 && dt->month < 13) pos+=sprintf(res+pos, "%d", dt->month);
			else pos += sprintf(res+pos, "##");			break;
		case 'W':
			if(dt->month >0 && dt->month < 13) pos+=sprintf(res+pos, "%c", dt_month[dt->month-1][0]);
			else pos += sprintf(res+pos, "#");			break;
		case 'D':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf(res+pos, "%s", dt_day[dt->dow-1]);
			else pos += sprintf(res+pos, "###");		break;
		case 'd':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf(res+pos, "%s", dt_days[dt->dow-1]);
			else pos += sprintf(res+pos, "###");		break;
		case 'E':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf(res+pos, "%02d", dt->dow);
			else pos += sprintf(res+pos, "##");			break;
		case 'e':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf(res+pos, "%d", dt->dow);
			else pos += sprintf(res+pos, "##");			break;
		case 'F':
			if(dt->dow >0 && dt->dow < 8) pos+=sprintf(res+pos, "%c", dt_day[dt->dow-1][0]);
			else pos += sprintf(res+pos, "#");			break;
		case 'H':
			if(dt->hours >=0 && dt->hours < 24) pos+=sprintf(res+pos, "%02d", dt->hours);
			else pos += sprintf(res+pos, "##");			break;
		case 'h':
			if(dt->hours >=0 && dt->hours < 24) pos+=sprintf(res+pos, "%d", dt->hours);
			else pos += sprintf(res+pos, "##");			break;
		case 'M':
			if(dt->minutes >=0 && dt->minutes < 60) pos+=sprintf(res+pos, "%02d", dt->minutes);
			else pos += sprintf(res+pos, "##");			break;
		case 'm':
			if(dt->minutes >=0 && dt->minutes < 60) pos+=sprintf(res+pos, "%d", dt->minutes);
			else pos += sprintf(res+pos, "##");			break;
		case 'S':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf(res+pos, "%02d", iround(secs));
			else pos += sprintf(res+pos, "##");			break;
		case 's':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf(res+pos, "%d", iround(secs));
			else pos += sprintf(res+pos, "##");			break;
		case 'T':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf(res+pos, "%02.2lf", dt->seconds);
			else pos += sprintf(res+pos, "##.##");		break;
		case 't':
			if(dt->seconds >=0.0 && dt->seconds < 60.0) pos+=sprintf(res+pos, "%.2lf", dt->seconds);
			else pos += sprintf(res+pos, "##.##");		break;
		default:
			pos += sprintf(res+pos, "%c", fmt[i]);		break;
			}
#endif
	}
	res[pos] = 0;
	return res;
}

double date2value(rlp_datetime *dt)
{
	double res;

	if(!dt) return 0.0;

	res = dt->seconds/60.0 + (double)dt->minutes;
	res = res/60.0 + (double)dt->hours;
	res = res/24.0 + (double)dt->aday;
	return res;
}

void parse_datevalue(rlp_datetime *dt, double dv)
{
	int i, j, d;

	if(!dt || dv < 0.0) return;
	if(dv > 1.0) {
		dt->aday = (int)floor(dv);
		dt->year = (int)(dv/365.2425);
		d = (int)floor(dv);
		do {
			dt->doy = d - 365*dt->year;
			dt->doy -= ((dt->year-1)>>2);
			dt->doy += ((dt->year)/100);
			dt->doy -= ((dt->year+300)/400);
			if(dt->doy < 1) dt->year--;
			}while(dt->doy < 1);
		dt->year += 1900;
		for(i = dt->month = 0, d = dt->doy; i < 12 && d > 0; i++) {
			if(i == 1 && d > (j = (leapyear(dt->year)) ? 29 : 28)) d -= j;
			else if(i != 1 && d > dt_monthl[i]) d -= dt_monthl[i];
			else break;
			}
		dt->month = i+1;				dt->dom = d;
		}
	dv -= floor(dv);				dv *= 24.0;
	dt->hours = (int)floor(dv);		dv -= floor(dv);
	dv *= 60.0;						dt->minutes = (int)floor(dv); 
	dv -= floor(dv);				dt->seconds = dv *60.0 + 1.0e-12; 
	if(dt->seconds > 59.9999) {
		dt->seconds = 0.0;			dt->minutes++;
		if(dt->minutes == 60) {
			dt->hours++;			dt->minutes = 0;
			}
		}
}

static char *dt_popfmt[] = {"Z.V.Y H:M:S", "Z/V/Y H:M:S", "Z-V-Y H:M:S", "Z.X.Y H:M:S",
	"Y.V.Z H:M:S", "Y-X-Z H:M:S", "H:M:S", 0L};

bool date_value(char *desc, char *fmt, double *value)
{
	int i;
	rlp_datetime dt;

	dt.year = dt.aday = dt.doy = dt.month = dt.dom = dt.dow = dt.hours = dt.minutes = 0;
	dt.seconds = 0.0;
	if(!value || !desc || !desc[0]) return false;
	if(fmt && fmt[0]) {
		if(parse_date(&dt, desc, fmt)) {
			*value = date2value(&dt);	return true;
			}
		}
	else {
		if(parse_date(&dt, desc, defs.fmt_datetime)) {
			*value = date2value(&dt);	return true;
			}
		}
	for(i=0; dt_popfmt[i]; i++) {
		if(parse_date(&dt, desc, dt_popfmt[i])) {
			*value = date2value(&dt);	return true;
			}
		}
	return false;
}

char *value_date(double dv, char *fmt)
{
	rlp_datetime dt;

	parse_datevalue(&dt, dv);
	return date2text(&dt, fmt ? fmt : defs.fmt_date);
}

double now_today()
{
	double res = 0.0;
	time_t ti = time(0L);
#ifdef USE_WIN_SECURE
	char dtbuff[80];

	ctime_s(dtbuff, 80, &ti);
	date_value(dtbuff+4, "x z H:M:S Y", &res);
#else
	date_value(ctime(&ti)+4, "x z H:M:S Y", &res);
#endif
	return res;
}

void split_date(double dv, int *y, int *mo, int *dom, int *dow, int *doy, int *h, int *m, double *s)
{
	rlp_datetime dt;

	parse_datevalue(&dt, dv);
	set_dow(&dt);
	if(y) *y = dt.year;				if(mo) *mo = dt.month;
	if(dom) *dom = dt.dom;			if(dow) *dow = dt.dow;
	if(doy) *doy = dt.doy;			if(h) *h = dt.hours;
	if(m) *m = dt.minutes;			if(s) *s = dt.seconds;
}

//---------------------------------------------------------------------------
// Use the Delauney triangulation to create a 3D mesh of dispersed data
//
Triangle* Triangulate1(char *xr, char *yr, char *zr, DataObj *data)
{
	AccRange *rX, *rY, *rZ;
	int i, j, n, rx, cx, ry, cy, rz, cz;
	double zMin;
	fPOINT3D *da;
	fRECT lim;
	Triangle *trl, *trn;
	Triangulate *tria;

	rX = rY = rZ = 0L;				trl = trn  = 0L;
	if((rX = new AccRange(xr)) && (rY = new AccRange(yr)) && (rZ = new AccRange(zr))
		&& rX->GetFirst(&cx, &rx) && rY->GetFirst(&cy, &ry) && rZ->GetFirst(&cz, &rz)
		&& (n = rX->CountItems()) && (da = (fPOINT3D*)malloc(n * sizeof(fPOINT3D)))
		&& (trl = new Triangle()) && (trn = new Triangle())) {
		//get minima and maxima
		for(i = j = 0; i < n; i++) {
			if(rX->GetNext(&cx, &rx) && rY->GetNext(&cy, &ry) && rZ->GetNext(&cz, &rz)) {
				data->GetValue(rx, cx, &da[j].fx);	data->GetValue(ry, cy, &da[j].fy);
				data->GetValue(rz, cz, &da[j].fz);	j++;
				}
			}
		if(!j) {
			free(da); delete rX;	delete rY;	delete rZ;	return trl;
			}
		for(i = 0, j = n; i < n; i++) {
			if(i) {
				if(da[i].fx < lim.Xmin) lim.Xmin = da[i].fx;	if(da[i].fx > lim.Xmax) lim.Xmax = da[i].fx;
				if(da[i].fy < lim.Ymin) lim.Ymin = da[i].fy;	if(da[i].fy > lim.Ymax) lim.Ymax = da[i].fy;
				if(da[i].fz < zMin) zMin = da[i].fz;
				}
			else {
				lim.Xmax = lim.Xmin = da[i].fx;		lim.Ymax = lim.Ymin = da[i].fy;		zMin = da[i].fz;
				}
			}
		//setup two super triangles
		trl->pt[0].fz = trl->pt[1].fz = trl->pt[2].fz = zMin;
		trn->pt[0].fz = trn->pt[1].fz = trn->pt[2].fz = zMin;
		trl->pt[0].fx = trn->pt[0].fx = trl->pt[2].fx = lim.Xmin;
		trl->pt[0].fy = trn->pt[0].fy = trn->pt[1].fy = lim.Ymin;
		trl->pt[1].fx = trn->pt[2].fx = trn->pt[1].fx = lim.Xmax;
		trl->pt[1].fy = trn->pt[2].fy = trl->pt[2].fy = lim.Ymax;
		trl->SetRect();			trn->SetRect();
		trl->next = trn;		trn->next = 0L;
		//do triangulation
		tria = new Triangulate(trl);
		for(i = 0; i < n; i++) {
			tria->AddVertex(&da[i]);
			}
		free(da);
		}
	if(rX) delete rX;	if(rY) delete rY;	if(rZ) delete rZ;
	trl = tria->trl;	delete tria;		return trl;
}

Ribbon *SurfTria(GraphObj *parent, DataObj *data, char *xr, char *yr, char *zr)
{
	Triangle *trl, *trc, *trn;
	int i, j, n, npl;
	double tmp;
	Plane3D **planes;

	trl = Triangulate1(xr, zr, yr, data);
	for(i = 0, trc = trl; trc; i++) trc = trc->next;
	if((n = i) && (planes = (Plane3D**)malloc(n*sizeof(Plane3D*)))) 
		for(i = npl = 0, trc = trl; trc && i < n; i++) {
		for(j = 0; j < 4; j++) {	//swap y and z values;
			tmp = trc->pt[j].fz;	trc->pt[j].fz = trc->pt[j].fy;	trc->pt[j].fy = tmp;
			}
		planes[npl++] = new Plane3D(0L, data, trc->pt, 4);
		trn = trc->next;	delete trc;		trc = trn;
		}
	if(npl) return new Ribbon(parent, data, (GraphObj**)planes, npl);
	return 0L;
}