/* To define printing from R use Rprintf, otherwise printf */
//#define LINKTOR
#ifdef LINKTOR
//#include <C:\\Program Files/R/rw2001/include/R_ext/Print.h>
void Rprintf(char*, ...);
void _Rprintf(char*, ...);
#define PRINT Rprintf
#else
#define PRINT printf
#endif
#include <math.h>
#include <malloc.h>
#include <stdio.h>

	/* Standard defines */
#ifdef __unix__
#define EXPORT void
#define INT long
#define IFMT %ld
#else
#define EXPORT __declspec(dllexport)void
#define INT int
#define IFMT %d
#endif
#define MAKEARRAY(size,type) ((type *)malloc((size)*sizeof(type)))
#define MIN(a,b) ((a)<(b)? (a):(b))
#define MAX(a,b) ((a)<(b)? (b):(a))
	/* TRUE represented as 1 in Splus */
#define FALSE 0
#define TRUE 1

/*====================================================================== */
/* Three functions for dealing with symmetric positive definite overlapping
block diagonal matrices.  Note that it would be possible to write versions of
these that were a few times faster for specific block sizes and overlaps.
This should be done if the code starts getting a lot of serious use on
very large data sets.  It could be done by inserting print statements (either
in this C code or in the R code on files development*.R with the C option
appearing to be more reliable) to write information
to assist with the writing of code.  The setting up of the matrices could
also be done in C if this happens. */

EXPORT BLOCKldrFACTORIZE(INT *nblocks,INT *blocksize,INT *overlap,
	double *diagelts,double *offdiag,INT *errorcode)
{
/*	Find LDR factorization (in place) of block symmetric matrix with
	nblocks blocks of size blocksize with specified overlap.
	Note that diagonal elements in result are stored as reciprocals.
	Errorcode returned as 1 if integer parameters are unsatisfactory. */
INT i,j,k,iblock,delta_vector,delta_offdiag,Ndiag,Noffdiag,blockoffsetZ,
    	blockoffdiagoffsetZ,blockoffset,blockoffdiagoffset,oi,oj;
INT* offset;
double* temp;
double sum;

offset = MAKEARRAY(*blocksize,INT);
temp = MAKEARRAY(*blocksize,double);

/* Compute the following:
delta_vector:	Increment in vector index between adjacent blocks
delta_offdiag:	Increment in index to off-diagonal elements between adjacent blocks
offset:		Offsets for rows in set of offdiagonal elements
Ndiag:		Number of rows = number of columns in matrices
Noffdiag:	Number of offdiagonal elements stored (includes some structural
	zeros for the first block so that its structure is the same as other blocks)
blockoffsetZ:	Vector index offset for use with last block
blockoffdiagoffsetZ:	Offset to index for off-diagonal elements for use \
		with last block */

*errorcode=0;
if(*nblocks < 2) {*errorcode=1;
	PRINT("BLOCKldrFACTORIZE: nblocks must be at least 2\n");
	return;}
if(*overlap < 1) {*errorcode=1;
	PRINT("BLOCKldrFACTORIZE: overlap must be at least 1\n");
	return;}
if(2 * *overlap > *blocksize) {*errorcode=1;
	PRINT("BLOCKldrFACTORIZE: overlap must not exceed 1/2 blocksize\n");
	return;}

delta_vector = *blocksize - *overlap;
delta_offdiag = delta_vector * (*blocksize + *overlap-1)/2;
Ndiag = *nblocks * delta_vector + *overlap;
Noffdiag=*nblocks*delta_offdiag+(*overlap*(*overlap-1))/2+*overlap*delta_vector;
blockoffsetZ = Ndiag - *blocksize;
blockoffdiagoffsetZ = (*nblocks-1) * delta_offdiag;

offset[0]=delta_vector;
for (i=1; i<*blocksize; i++) {
	k=i-1;
	if(i< *overlap)k=k+delta_vector;
	offset[i]=offset[i-1]+k;
}

for (i=0; i< *blocksize; i++){
	oi = offset[i];
	sum = diagelts[i];
	for (j=0; j<i; j++){
		temp[j] = offdiag[oi+j];
		oj = offset[j];
		for (k=0;k<j; k++){temp[j] -= temp[k]*offdiag[oj+k];}
		offdiag[oi+j] = temp[j]*diagelts[j];
		sum -= offdiag[oi+j]*temp[j];
	}
	diagelts[i] = 1/sum;
}
blockoffset = delta_vector;
blockoffdiagoffset = delta_offdiag;
for (iblock=1; iblock< *nblocks; iblock++){
   for (i= *overlap; i< *blocksize; i++){
	oi = blockoffdiagoffset+offset[i];
	sum = diagelts[blockoffset+i];
	for (j=0; j<i; j++){
		temp[j] = offdiag[oi+j];
		oj = blockoffdiagoffset+offset[j];
		for (k=0; k<j; k++){temp[j] -= temp[k]*offdiag[oj+k];}
		offdiag[oi+j] = temp[j]*diagelts[blockoffset+j];
		sum -= offdiag[oi+j]*temp[j];
	}
	diagelts[blockoffset+i] = 1/sum;
   }
   blockoffset += delta_vector;
   blockoffdiagoffset += delta_offdiag;
}

return;
}
/*====================================================================== */
EXPORT BLOCKldrSOLVE(INT *nblocks,INT *blocksize,INT *overlap,
	double *diagelts,double *offdiag,double *v,INT *errorcode)
{
/*	Solve matrix equation with LHS matrix given by LDR factorization
	of block symmetric matrix with nblocks blocks of size blocksize with
	specified overlap.  RHS vector is v and will be replaced by solution.
	Note that diagonal elements in result are stored as reciprocals.
	Errorcode returned as 1 if integer parameters are unsatisfactory. */
INT i,j,k,iblock,delta_vector,delta_offdiag,Ndiag,Noffdiag,blockoffsetZ,
    	blockoffdiagoffsetZ,blockoffset,blockoffdiagoffset,oi;
INT* offset;
offset = MAKEARRAY(*blocksize,INT);

/* Compute the following:
delta_vector:	Increment in vector index between adjacent blocks
delta_offdiag:	Increment in index to off-diagonal elements between adjacent blocks
offset:		Offsets for rows in set of offdiagonal elements
Ndiag:		Number of rows = number of columns in matrices
Noffdiag:	Number of offdiagonal elements stored (includes some structural
	zeros for the first block so that its structure is the same as other blocks)
blockoffsetZ:	Vector index offset for use with last block
blockoffdiagoffsetZ:	Offset to index for off-diagonal elements for use \
		with last block */

*errorcode=0;
if(*nblocks < 2) {*errorcode=1;
	PRINT("BLOCKldrSOLVE: nblocks must be at least 2\n");
	return;}
if(*overlap < 1) {*errorcode=1;
	PRINT("BLOCKldrSOLVE: overlap must be at least 1\n");
	return;}
if(2 * *overlap > *blocksize) {*errorcode=1;
	PRINT("BLOCKldrSOLVE: overlap must not exceed 1/2 blocksize\n");
	return;}

delta_vector = *blocksize - *overlap;
delta_offdiag = delta_vector * (*blocksize + *overlap-1)/2;
Ndiag = *nblocks * delta_vector + *overlap;
Noffdiag=*nblocks*delta_offdiag+(*overlap*(*overlap-1))/2+*overlap*delta_vector;
blockoffsetZ = Ndiag - *blocksize;
blockoffdiagoffsetZ = (*nblocks-1) * delta_offdiag;

offset[0]=delta_vector;
for (i=1; i<*blocksize; i++) {
	k=i-1;
	if(i< *overlap)k=k+delta_vector;
	offset[i]=offset[i-1]+k;
}

for (i=0; i< *blocksize; i++){
	oi = offset[i];
	for (j=0; j<i; j++){v[i] -= v[j]*offdiag[oi+j];}
}
blockoffset = delta_vector;
blockoffdiagoffset = delta_offdiag;
for (iblock=1; iblock < *nblocks; iblock++ ){
   for (i = *overlap; i< *blocksize; i++){
	oi = blockoffdiagoffset+offset[i];
	for (j=0; j<i; j++){
		v[blockoffset+i] -= v[blockoffset+j]*offdiag[oi+j];
	}
   }
   blockoffset += delta_vector;
   blockoffdiagoffset += delta_offdiag;
}

/* Multiply by inverse of diagonal matrix */
for (i=0; i<Ndiag; i++){v[i] *= diagelts[i];}

/* Compute LU inverse v */
blockoffset = blockoffsetZ;
blockoffdiagoffset = blockoffdiagoffsetZ;
for (iblock= *nblocks-1; iblock>0; iblock -= 1){
   for (i=*blocksize-1; i>=*overlap; i -= 1){
	oi = blockoffdiagoffset+offset[i];
	for (j=0; j< i; j++){
		v[blockoffset+j] -= v[blockoffset+i]*offdiag[oi+j];
	}
   }
   blockoffset -= delta_vector;
   blockoffdiagoffset -= delta_offdiag;
}
for (i=*blocksize-1; i>=0; i -= 1){
	oi = offset[i];
	for (j=0; j<i; j++){v[j] -= v[i]*offdiag[oi+j];}
}

return;
}
/*====================================================================== */
EXPORT BLOCKldrINVERSE(INT *nblocks,INT *blocksize,INT *overlap,
	double *diagelts,double *offdiag,
	double *invdiag,double *invoff,INT *errorcode)
{
/*	Find portion of inverse matrix from LDR factorization of block
	symmetric matrix with nblocks blocks of size blocksize with
	specified overlap.  Note that diagonal elements in input are stored
	as reciprocals.  Result returned in invoff and invdiag.
	Errorcode returned as 1 if integer parameters are unsatisfactory. */
INT i,j,k,iblock,delta_vector,delta_offdiag,Ndiag,Noffdiag,blockoffsetZ,
    	blockoffdiagoffsetZ,blockoffset,blockoffdiagoffset,oi,ok;
INT* offset;
double sum;
offset = MAKEARRAY(*blocksize,INT);

/* Compute the following:
delta_vector:	Increment in vector index between adjacent blocks
delta_offdiag:	Increment in index to off-diagonal elements between adjacent blocks
offset:		Offsets for rows in set of offdiagonal elements
Ndiag:		Number of rows = number of columns in matrices
Noffdiag:	Number of offdiagonal elements stored (includes some structural
	zeros for the first block so that its structure is the same as other blocks)
blockoffsetZ:	Vector index offset for use with last block
blockoffdiagoffsetZ:	Offset to index for off-diagonal elements for use \
		with last block */

*errorcode=0;
if(*nblocks < 2) {*errorcode=1;
	PRINT("BLOCKldrINVERSE: nblocks must be at least 2\n");
	return;}
if(*overlap < 1) {*errorcode=1;
	PRINT("BLOCKldrINVERSE: overlap must be at least 1\n");
	return;}
if(2 * *overlap > *blocksize) {*errorcode=1;
	PRINT("BLOCKldrINVERSE: overlap must not exceed 1/2 blocksize\n");
	return;}

delta_vector = *blocksize - *overlap;
delta_offdiag = delta_vector * (*blocksize + *overlap-1)/2;
Ndiag = *nblocks * delta_vector + *overlap;
Noffdiag=*nblocks*delta_offdiag+(*overlap*(*overlap-1))/2+*overlap*delta_vector;
blockoffsetZ = Ndiag - *blocksize;
blockoffdiagoffsetZ = (*nblocks-1) * delta_offdiag;

offset[0]=delta_vector;
for (i=1; i<*blocksize; i++) {
	k=i-1;
	if(i< *overlap)k=k+delta_vector;
	offset[i]=offset[i-1]+k;
}

/* In reverse order of rows, compute non-structural-zero offdiagonal
	elements of L inverse. */
blockoffset = blockoffsetZ;
blockoffdiagoffset = blockoffdiagoffsetZ;
for (iblock=*nblocks; iblock>1; iblock -= 1){
   for (i= *blocksize-1; i>= *overlap; i -= 1){
	oi = blockoffdiagoffset+offset[i];
	for (j=i-1; j >= 0; j -= 1){
		sum = offdiag[oi+j];
		for (k=j+1; k<=i-1; k++){
		   sum += invoff[oi+k]*offdiag[blockoffdiagoffset+offset[k]+j];}
		invoff[oi+j] = -sum;
	}
   }
   blockoffset -= delta_vector;
   blockoffdiagoffset -= delta_offdiag;
}
for (i=*blocksize-1; i>=0; i -= 1){
	oi = offset[i];
	for (j=i-1; j>=0; j -= 1){
		sum = offdiag[oi+j];
		for (k=j+1; k<=i-1; k++){
			sum += invoff[oi+k]*offdiag[offset[k]+j];}
		invoff[oi+j] = -sum;
	}
}

/* Now pre-multiply by D inverse */
for (i=0; i< Ndiag; i++){invdiag[i] = diagelts[i];}
for (i=1; i< *blocksize; i++){
	for (j=0; j<=i-1; j++){invoff[offset[i]+j] *= diagelts[i];}
}
blockoffset = delta_vector;
blockoffdiagoffset = delta_offdiag;
for (iblock=2; iblock <= *nblocks; iblock++){
   for (i=*overlap; i< *blocksize; i++){
	for (j=0; j<=i-1; j++){
	   invoff[blockoffdiagoffset+offset[i]+j] *= diagelts[blockoffset+i];}
   }
   blockoffset += delta_vector;
   blockoffdiagoffset += delta_offdiag;
}

/* Now add (I - U)Z in place */
blockoffset = blockoffsetZ;
blockoffdiagoffset = blockoffdiagoffsetZ;
for (j= *blocksize-2; j>=0; j -= 1){
	for (i= *blocksize-2; i>j; i -= 1){
		oi = blockoffdiagoffset+offset[i];
		sum = 0.;
		for (k=i+1; k< *blocksize; k++){
			ok = blockoffdiagoffset+offset[k];
			sum += invoff[ok+j]*offdiag[ok+i];}
		invoff[oi+j] -= sum;
	}
	sum = 0.;
	for (k=j+1; k< *blocksize; k++){
		ok = blockoffdiagoffset+offset[k];
		sum += invoff[ok+j]*offdiag[ok+j];}
	invdiag[blockoffset+j] -= sum;
}

for (iblock=*nblocks-1; iblock > 0; iblock -= 1){
   blockoffset -= delta_vector;
   blockoffdiagoffset -= delta_offdiag;
   for (j=delta_vector-1; j>=0; j -= 1){
	for (i= *blocksize-1; i>= delta_vector; i -= 1){
		oi = blockoffdiagoffset+offset[i];
	/* Compute using elements in of Z in i'th column, including diagonal */
		sum = invdiag[blockoffset+i]*offdiag[oi+j];
		for (k=i+1; k< *blocksize; k++){
			ok = blockoffdiagoffset+offset[k];
			sum += invoff[ok+i]*offdiag[ok+j];}
		for (k=j+1; k<i; k++){
			ok = blockoffdiagoffset+offset[k];
			sum += invoff[oi+k]*offdiag[ok+j];}
		invoff[oi+j] = -sum;
	}
	for (i=delta_vector-1; i>j; i -= 1){
		oi = blockoffdiagoffset+offset[i];
		sum = 0.;
		for (k=i+1; k< *blocksize; k++){
			ok = blockoffdiagoffset+offset[k];
			sum += invoff[ok+j]*offdiag[ok+i];}
		invoff[oi+j] -= sum;
	}
	sum = 0.;
	for (k=j+1; k< *blocksize; k++){
		ok = blockoffdiagoffset+offset[k];
		sum += invoff[ok+j]*offdiag[ok+j];}
	invdiag[blockoffset+j] -= sum;
   }
}

return;
}
/*====================================================================== */
EXPORT blockvariog(double *block,double *time,double *y,
	double *v,INT *count,INT *n,INT *maxlag)
{
/*	To compute sample variogram into v and number of
		contributions into count.
	Data is given as vectors of length n in block then time order:
		block	: a set of block identifiers
		time	: the times (as integers) of measurements
		y	: a set of measurements
*/
INT lagm1,i,j,thismaxlag,lag;
double temp,thisblock;
thismaxlag = *maxlag;

for (lagm1=0; lagm1<thismaxlag; lagm1++){
	/* Use v as an accumulator for sum of squared deviations. */
	v[lagm1]=0.;
	count[lagm1]=0;
	}
for (i=0; i< *n; i++){
	thisblock  = block[i];
	for (j=i+1; j< *n; j++) {
		if(block[j] != thisblock) break;
		lag = (INT)(time[j]-time[i]);
		if(lag > thismaxlag) break;
		if(lag > 0){
			count[lag-1] ++;
			temp = y[i]-y[j];
			v[lag-1] += temp*temp;
			}
	}
}
for (lagm1=0; lagm1< thismaxlag; lagm1++){
	if(count[lagm1] > 0){
		/* variogram is half of average squared deviation. */
		v[lagm1]=v[lagm1]*.5/count[lagm1];
		}
	}
return;
}
/*====================================================================== */
EXPORT ewma(double *x,double *lambda,double *result,INT *n)
{
/*	To compute EWMA of x with parameter lambda into result. */
INT i;
double omlambda;
omlambda=1.-*lambda;
result[0]=x[0];
for (i=1; i<*n; i++){
	result[i]=omlambda*result[i-1]+*lambda*x[i];
	}
return;
}
/*====================================================================== */
EXPORT KF(INT *n,double *y,double *VARy,double *VARdf,
	double *result,double *VARr)
{
/*	Computes Kalman Filter of y into result
	 and its estimation variance into VARr.*/
/*	Number of data points is n */
/*	Data points in vector y */
/*	Error variances of data in VARy */
/*	Variances of changes in Brownian motion trend in VARdf */
INT i;
double old,recip;

result[0]=y[0];
VARr[0]=VARy[0];
for (i=1; i<*n; i++){
	old=VARr[i-1] + VARdf[i-1];
	recip=1./(old+VARy[i]);
	VARr[i]=old*VARy[i]*recip;
	result[i]=(VARy[i]*result[i-1]+old*y[i])*recip;
	}
return;
}
/*====================================================================== */
EXPORT KFLS(INT *n,INT *nsmooth,double *y,double *VARy,double *VARdf,
        double *result,double *VARresult)
{
/*      Computes Kalman Filter plus (limited) fixed lag Smoothing */
/*      Number of data points is n */
/*      Maximum amount of smoothing is nsmooth (at least 2) */
/*      Data points in vector y */
/*      Error variances of data in VARy */
/*      Variances of changes in Brownian motion trend in VARdf */
/*      Output point estimate in result */
/*      Output estimation variances in result */
INT i,j;
double old,recip,Ai,Vip1;
double *backest,*backvar;

backest=MAKEARRAY(*nsmooth,double);
backvar=MAKEARRAY(*nsmooth,double);

/*      Kalman Filter for entire data set       */
result[0]=y[0];
VARresult[0]=VARy[0];
for (i=1; i<*n; i++){
        old=VARresult[i-1] + VARdf[i-1];
        recip=1./(old+VARy[i]);
        VARresult[i]=old*VARy[i]*recip;
        result[i]=(VARy[i]*result[i-1]+old*y[i])*recip;
        }
/*      Fixed interval smoothing algorithm for short data segments      */
for (i=0; i< *n- *nsmooth; i++){
        for (j=0; j< *nsmooth; j++){
                backest[j]=result[i+j];
                backvar[j]=VARresult[i+j];
                }
        for (j= *nsmooth-2; j>=0; j=j-1){
                Vip1=backvar[j]+VARdf[i+j];
                Ai=backvar[j]/Vip1;
                backvar[j] -= Ai*Ai*(Vip1-backvar[j+1]);
                backest[j] += Ai*(backest[j+1]-backest[j]);
                }
        VARresult[i]=backvar[0];
        result[i]=backest[0];
        }
/*      Fixed interval smoothing algorithm for final data segment       */
for (i=*n-2; i>= *n- *nsmooth; i=i-1){
        Vip1=VARresult[i]+VARdf[i];
        Ai=VARresult[i]/Vip1;
        VARresult[i] -= Ai*Ai*(Vip1-VARresult[i+1]);
        result[i] += Ai*(result[i+1]-result[i]);
        }
free(backest);
free(backvar);
return;
}
/*====================================================================== */
EXPORT KFLSuniform(INT *n,INT *nsmooth,double *y,double *VARy,double *VARdf,
        double *result,double *VARresult)
{
/*      Computes Kalman Filter plus fixed lag smoothing (uniform version)*/
/*      Number of data points is n */
/*      Maximum amount of smoothing is nsmooth (at least 2) */
/*      Data points in vector y */
/*      Error variance of data in VARy */
/*      Variance of changes in Brownian motion trend in VARdf */
/*      Output point estimate in result */
/*      Output estimation variances in result */
INT i,j;
double old,recip,Ai,Vip1;
double *backest,*backvar;

backest=MAKEARRAY(*nsmooth,double);
backvar=MAKEARRAY(*nsmooth,double);

/*      Kalman Filter for entire data set       */
result[0]=y[0];
VARresult[0]= *VARy;
for (i=1; i<*n; i++){
        old=VARresult[i-1] + *VARdf;
        recip=1./(old+ *VARy);
        VARresult[i]=old* *VARy*recip;
        result[i]=(*VARy*result[i-1]+old*y[i])*recip;
        }
/*      Fixed interval smoothing algorithm for short data segments      */
for (i=0; i< *n- *nsmooth; i++){
        for (j=0; j< *nsmooth; j++){
                backest[j]=result[i+j];
                backvar[j]=VARresult[i+j];
                }
        for (j= *nsmooth-2; j>=0; j=j-1){
                Vip1=backvar[j]+ *VARdf;
                Ai=backvar[j]/Vip1;
                backvar[j] -= Ai*Ai*(Vip1-backvar[j+1]);
                backest[j] += Ai*(backest[j+1]-backest[j]);
                }
        VARresult[i]=backvar[0];
        result[i]=backest[0];
        }
/*      Fixed interval smoothing algorithm for final data segment       */
for (i=*n-2; i>= *n- *nsmooth; i=i-1){
        Vip1=VARresult[i]+ *VARdf;
        Ai=VARresult[i]/Vip1;
        VARresult[i] -= Ai*Ai*(Vip1-VARresult[i+1]);
        result[i] += Ai*(result[i+1]-result[i]);
        }
free(backest);
free(backvar);
return;
}
/*====================================================================== */
EXPORT KFuniform(INT *n,double *y,double *VARyS,double *VARdfS,
	double *result,double *VARresult)
{
/*	Computes Kalman Filter of y into result
		and its estimation variance into VARresult for the case
		of equal precision, equally spaced data.*/
/*	Number of data points is n */
/*	Data points in vector y */
/*	Error variance (scalar) of data in VARyS */
/*	Variance of changes in Brownian motion trend (scalar) in VARdfS */
INT i;
double old,recip,VARy,VARdf;
VARy= *VARyS;
VARdf= *VARdfS;

/*	Kalman Filter	*/
result[0]=y[0];
VARresult[0]=VARy;
for (i=1; i<*n; i++){
	old=VARresult[i-1] + VARdf;
	recip=1./(old+VARy);
	VARresult[i]=old*VARy*recip;
	result[i]=(VARy*result[i-1]+old*y[i])*recip;
	}
return;
}
/*====================================================================== */
EXPORT KS(INT *n,double *y,double *VARy,double *VARdf,
	double *result,double *VARresult)
{
/*	Computes Kalman Smoother (i.e. Kalman Filter followed by fixed
		interval smoothing algorithm) of y into result
		and its estimation variance into VARresult.*/
/*	Number of data points is n */
/*	Data points in vector y */
/*	Error variances of data in VARy */
/*	Variances of changes in Brownian motion trend in VARdf */
INT i;
double old,recip,Ai,Vip1;

/*	Kalman Filter	*/
result[0]=y[0];
VARresult[0]=VARy[0];
for (i=1; i<*n; i++){
	old=VARresult[i-1] + VARdf[i-1];
	recip=1./(old+VARy[i]);
	VARresult[i]=old*VARy[i]*recip;
	result[i]=(VARy[i]*result[i-1]+old*y[i])*recip;
	}
/*	Now: fixed interval smoothing algorithm	*/
for (i=*n-2; i>=0; i=i-1){
	Vip1=VARresult[i]+VARdf[i];
	Ai=VARresult[i]/Vip1;
	VARresult[i] -= Ai*Ai*(Vip1-VARresult[i+1]);
	result[i] += Ai*(result[i+1]-result[i]);
	}
return;
}
/*====================================================================== */
EXPORT KSuniform(INT *n,double *y,double *VARyS,double *VARdfS,
	double *result,double *VARresult)
{
/*	Computes Kalman Smoother (i.e. Kalman Filter followed by fixed
		interval smoothing algorithm) of y into result
		and its estimation variance into VARresult for the case
		of equal precision, equally spaced data.*/
/*	Number of data points is n */
/*	Data points in vector y */
/*	Error variance (scalar) of data in VARyS */
/*	Variance of changes in Brownian motion trend (scalar) in VARdfS */
INT i;
double old,recip,Ai,Vip1,VARy,VARdf;
VARy= *VARyS;
VARdf= *VARdfS;

/*	Kalman Filter	*/
result[0]=y[0];
VARresult[0]=VARy;
for (i=1; i<*n; i++){
	old=VARresult[i-1] + VARdf;
	recip=1./(old+VARy);
	VARresult[i]=old*VARy*recip;
	result[i]=(VARy*result[i-1]+old*y[i])*recip;
	}
/*	Now: fixed interval smoothing algorithm	*/
for (i=*n-2; i>=0; i=i-1){
	Vip1=VARresult[i]+VARdf;
	Ai=VARresult[i]/Vip1;
	VARresult[i] -= Ai*Ai*(Vip1-VARresult[i+1]);
	result[i] += Ai*(result[i+1]-result[i]);
	}
return;
}
/*====================================================================== */
EXPORT ldrTSM(double *x,INT *errorcode)
{
/*	Find LDR factorization of tridiagonal symmetric matrix (in place)
	Note that D is stored as reciprocals.
	Errorcode returned as 1 if matrix not positive definite. */
INT i,n;
double *below,*diag,temp,thisdiag;
n=(INT)x[0];
diag=&x[1];
below=diag+n-1;
*errorcode=0;

if(diag[0] <= 0) {*errorcode=1;
	PRINT("ldrTSM: Matrix not positive definite\n");
	return;}
diag[0]=1/diag[0];
for (i=1; i<n; i++){
	temp=below[i]*diag[i-1];
	thisdiag=diag[i]-below[i]*temp;
	if(thisdiag <= 0) {*errorcode=1;
		PRINT("ldrTSM: Matrix not positive definite\n");
		return;}
	diag[i]=1/thisdiag;
	below[i]=temp;
	}
return;
}
/*====================================================================== */
EXPORT robustblockvariog(double *block,double *time,double *y,
	double *v,INT *count,INT *n,INT *maxlag)
{
/*	To compute robust estimate of variogram into v and number of
		contributions into count.
	Data is given as vectors of length n in block then time order:
		block	: a set of block identifiers
		time	: the times (as integers) of measurements
		y	: a set of measurements
*/
INT lagm1,i,j,thismaxlag,lag;
double temp,thisblock;
thismaxlag = *maxlag;

for (lagm1=0; lagm1<thismaxlag; lagm1++){
	/* Use v as an accumulator for sum of absolute deviations. */
	v[lagm1]=0.;
	count[lagm1]=0;
	}
for (i=0; i< *n; i++){
	thisblock  = block[i];
	for (j=i+1; j< *n; j++) {
		if(block[j] != thisblock) break;
		lag = (INT)(time[j]-time[i]);
		if(lag > thismaxlag) break;
		if(lag > 0){
			count[lag-1] ++;
			v[lag-1] += fabs(y[i]-y[j]);
			}
	}
}
for (lagm1=0; lagm1< thismaxlag; lagm1++){
	/* Estimate of sigma is sqrt(pi)/2 times mean range. */
	if(count[lagm1] > 0){
		/* Estimate of sigma is sqrt(pi)/2 times mean range. */
		temp=v[lagm1]*.886226925452758/count[lagm1];
		v[lagm1]=temp*temp;
		}
	}
return;
}
/* ==================================================================== */
EXPORT robustvariog(double *x,double *v,INT *n,INT *maxlag)
{
/*	To compute robust estimate of variogram of x (length n) into v. */
INT lag,i;
double sum,temp;
if(*maxlag<=0)*maxlag=(*n-1);
for (lag=1; lag<=*maxlag; lag++){
	sum=0.;
	for (i=0; i+lag < *n; i++){sum=sum+fabs(x[i]-x[i+lag]);}
	/* Estimate of sigma is sqrt(pi)/2 times mean range. */
	temp=sum*.886226925452758/(*n-lag);
	v[lag-1]=temp*temp;
	}
return;
}
/* ==================================================================== */
EXPORT robustvariog2(double *x,double *v,INT *n,INT *maxlag)
{
/*	Compute robust estimate (type 2) of variogram of x (length n) into v. */
INT lag,i;
double sum,temp;
if(*maxlag<=0)*maxlag=(*n-1);
for (lag=1; lag<=*maxlag; lag++){
	sum=0.;
	for (i=0; i+lag < *n; i++){sum=sum+sqrt(fabs(x[i]-x[i+lag]));}
	/* If the data is i.i.d. with standard deviation sigma,
	expected root range is (2^.75)*gamma(.75)/sqrt(pi) * sqrt(sigma). */
	temp=sum/(*n-lag);
	v[lag-1]=1.0942198076132*pow(temp,4);
	}
return;
}
/* ==================================================================== */
EXPORT robustvariog4(double *x,double *v,INT *n,INT *maxlag)
{
/*	Compute robust estimate (type 4) of variogram of x (length n) into v. */
INT lag,i;
double sum,temp;
if(*maxlag<=0)*maxlag=(*n-1);
for (lag=1; lag<=*maxlag; lag++){
	sum=0.;
	for (i=0; i+lag < *n; i++){sum=sum+sqrt(sqrt(sqrt(fabs(x[i]-x[i+lag]))));}
	temp=sum/(*n-lag);
	v[lag-1]=1.54196463645237*pow(temp,16);
	}
return;
}
/*====================================================================== */
EXPORT SIMbyMAT(INT *nrowsS,INT *ncolsS,INT *nnzS,INT *irow,INT *icol,
	INT *ncolmatS,double *matS,double *resS)
{
/*	Multiply sparse indicator matrix by full-stored matrix matS, into resS.
	Sparse matrix must be nrowsS x ncolsS, with nnzS unit elements.
	Row and column stored in irow and icol, respectively.
	Elements may be in any order.
	Full-stored matrix must have ncolsS rows and ncolmatS columns.
	Result resS will be a full-stored matrix with nrowsS rows and
	ncolmatS columns. */
int i,j,nrows,ncols,nnz,ncolmat;
double *mat,*res;	/* with first element numbered 0, in C style */
mat=matS-1;
res=resS-1;

nrows=*nrowsS;
ncols=*ncolsS;
nnz=*nnzS;
ncolmat=*ncolmatS;

for (j=0; j<ncolmat; j++){
	for (i=1; i<=nrows; i++){res[i]=0.;}
	for (i=0; i<nnz; i++){res[irow[i]] += mat[icol[i]];}
	mat += ncols;
	res += nrows;
	}
return;
}
/*====================================================================== */
EXPORT SMbyMAT(INT *nrowsS,INT *ncolsS,INT *nnzS,INT *irow,INT *icol,
	INT *ncolmatS,double *ientry,double *matS,double *resS)
{
/*	Multiply sparse matrix by full-stored matrix matS, into resS.
	Sparse matrix must be nrowsS x ncolsS, with nnzS non-zero elements.
	Row, column and entry stored in irow, icol and ientry, respectively.
	Elements may be in any order.
	Full-stored matrix must have ncolsS rows and ncolmatS columns.
	Result resS will be a full-stored matrix with nrowsS rows and
	ncolmatS columns. */
int i,j,nrows,ncols,nnz,ncolmat;
double *mat,*res;	/* with first element numbered 0, in C style */
mat=matS-1;
res=resS-1;

nrows=*nrowsS;
ncols=*ncolsS;
nnz=*nnzS;
ncolmat=*ncolmatS;

for (j=0; j<ncolmat; j++){
	for (i=1; i<=nrows; i++){res[i]=0.;}
	for (i=0; i<nnz; i++){res[irow[i]] += ientry[i]*mat[icol[i]];}
	mat += ncols;
	res += nrows;
	}
return;
}
/*====================================================================== */
EXPORT TSMbyMAT(double *x,INT *n,INT *m,double *matrix,double *result,
	INT *errorcode)
{
/*	Multiply tridiagonal symmetric matrix by matrix
	Errorcode returned as 1 if dimensions don't cohere */
INT i,j,last,here,next;
double *below,*diag;
if(*n != x[0]){
	PRINT("TSMbyMAT: Matrix sizes incompatible\n");
	*errorcode=1;
        return;};
diag=&x[1];
below=diag+ *n-1;
*errorcode=0;

for (j=0; j<*m; j++){
	here=j* *n;
	next=here+1;
	result[here]=diag[0]*matrix[here]+below[1]*matrix[next];
	for (i=1; i<*n-1; i++){
		last=here;
		here=next;
		next += 1;
		result[here]=below[i]*matrix[last]+
			diag[i]*matrix[here]+below[i+1]*matrix[next];
		}
	last=here;
	here=next;
	result[here]=below[*n-1]*matrix[last]+diag[*n-1]*matrix[here];
	}
return;
}
/*====================================================================== */
EXPORT TSMinvMATldr(double *x,INT *n,INT *m,double *matrix,double *result,
	INT *errorcode)
{
/*	Tridiagonal symmetric matrix inverse times matrix using LDR
	Errorcode returned as 1 if dimensions don't cohere */
INT i,j,last,here,next;
double *below,*diag;
if(*n != x[0]){
	PRINT("TSMbyMAT: Matrix sizes incompatible\n");
	*errorcode=1;
        return;};
diag=&x[1];
below=diag+ *n-1;
*errorcode=0;

for (j=0; j<*m; j++){
	/* First compute L inverse times vector. */
	here=j* *n;
	result[here]=matrix[here];
	for (i=1; i<*n; i++){
		last=here;
		here += 1;
		result[here]=(matrix[here]-below[i]*result[last]);
		}
	/* Now compute R inverse times D inverse times vector. */
	result[here]=result[here]*diag[*n-1];
	for (i=*n-2; i>=0; i=i-1){
		next=here;
		here=here-1;
		result[here]=(result[here]*diag[i]-below[i+1]*result[next]);
		}
	}
return;
}
/*====================================================================== */
EXPORT variog(double *x,double *v,INT *n,INT *maxlag)
{
/*	To compute variogram of x (length n) into v. */
INT lag,i;
double sum,diff;
if(*maxlag<=0)*maxlag=(*n-1);
for (lag=1; lag<=*maxlag; lag++){
	sum=0.;
	for (i=0; i+lag < *n; i++){
		diff=x[i]-x[i+lag];
		sum=sum+diff*diff;
		}
	v[lag-1]=sum*.5/(*n-lag);
	}
return;
}
/*====================================================================== */
double millsratio(double x)
{
/* Computes Mill's ratio, which is the ratio of the right tail */
/*      probability to the density for a standard normal distribution. */
/* Note that this function should not be called for negative arguments. */

double val,t;
double rhodiff[26]={1.,3.94766755064487339,
        .262699373075432297,-4.79530994886257864,-.472387694632788044,
        -4.64178612251027508,.888029797897972282E-1,3.34506626993363627,
        -.253018663557510120,-3.22271919537673200,.367907414060725404,
        .515135649489030960,-.782205852853700406,-1.66658903030166903,
        -.317783299853388699,3.44266858016113165,.110601063267646032,
        -113.675951225046565,-.475674678528123486E-2,891.504346794291090,
        .155839071476740257E-3,-52392.5264048076612,-.145796580254323180E-4,
        5331522.33187226848,.169353718097630593E-6,-83341339.7482781260};
double xdata[25]={0.,1.,.2704,.5776,.0784,.04,.0576,.0256,.1024,
        .1296,.16,.1936,.2304,.0064,.3136,.36,.4096,.4624,
        .5184,.0144,.64,.7056,.7744,.8464,.9216};
if(x<0.){PRINT("(millsratio) Argument x should not be negative");
        return 0.;
}

/* Thiele interpolation by Abramowitz and Stegun 25.2.50*/
t=1./(x+1.);
val=t*(rhodiff[0]+(t-xdata[0])/(rhodiff[1]+(t-xdata[1])/(rhodiff[2]+
        (t-xdata[2])/(rhodiff[3]+(t-xdata[3])/(rhodiff[4]+
        (t-xdata[4])/(rhodiff[5]+(t-xdata[5])/(rhodiff[6]+
        (t-xdata[6])/(rhodiff[7]+(t-xdata[7])/(rhodiff[8]+
        (t-xdata[8])/(rhodiff[9]+(t-xdata[9])/(rhodiff[10]+
        (t-xdata[10])/(rhodiff[11]+(t-xdata[11])/(rhodiff[12]+
        (t-xdata[12])/(rhodiff[13]+(t-xdata[13])/(rhodiff[14]+
        (t-xdata[14])/(rhodiff[15]+(t-xdata[15])/(rhodiff[16]+
        (t-xdata[16])/(rhodiff[17]+(t-xdata[17])/(rhodiff[18]+
        (t-xdata[18])/(rhodiff[19]+(t-xdata[19])/(rhodiff[20]+
        (t-xdata[20])/(rhodiff[21]+(t-xdata[21])/(rhodiff[22]+
        (t-xdata[22])/(rhodiff[23]+(t-xdata[23])/(rhodiff[24]+
        (t-xdata[24])/(rhodiff[25]))))))))))))))))))))))))));
return val;
}
/*=====================================================================*/
void rightnormaltail(double z, double *derivative, double *tail)
{
/* Finds derivative and right tail probability for standard normal distribution */
/* Note that precision of the tail probability is poor if z is negative. */

double const dnorm0=.39894228040143267793994605993438187;
*derivative=-dnorm0*exp(-.5*z*z);
if(z<0.) *tail=1.+ *derivative * millsratio(-z);
else *tail=- *derivative * millsratio(z);
}
/*=====================================================================*/
void leftnormaltail(double z, double *derivative, double *tail)
{
/* Finds derivative and left tail probability for standard normal distribution */
/* Note that precision of the tail probability is poor if z is negative. */

double const dnorm0=.39894228040143267793994605993438187;
*derivative=dnorm0*exp(-.5*z*z);
if(z>0.) *tail=1.- *derivative * millsratio(z);
else *tail= *derivative * millsratio(-z);
return;
}
/*=====================================================================*/
EXPORT webmseplus(INT *nx,  INT *nL, INT *uL, double *x, double *y,
	double *nobs, double *B, double *V, double *W,
	double *minx, double *maxx, double *yL,
	double *webm, double *varebm, double *PL, double *precPL,
	INT *errorcode)
/*
Variables describing data:
	nx = number of distinct points
	x[i] = i'th value of x, which defines the ascending sequence of the data
	y[i] = i'th observed value

Variables describing precision of data:
	nobs[i] = number of observations averaged to give y[i].  May be zero.
	B = variance of Brownian motion for unit change in x.
	W = variance of white noise in true values of variable when nobs[i]=1.
	V = variance of measurement error when nobs[i]=1.

Variables describing computations required:
	minx = lower range of x for which probabilities to be calculated
	maxx = upper range of x for which probabilities to be calculated
	nL = number of limits
	yL = vector of limits
	uL = whether limits are lower or upper
	To compute probability of being lower: non-0 from C or T from Splus
	To compute probability of being higher: 0 from C or F from Splus

Output variables:
	webm = (weighted) estimate of Brownian motion
	varwebm = estimation variance of (weighted) estimate of Brownian motion
	PL = estimated probability that Y < or > yL[i] for x[j] < x < x[j+1]
	(These probabilities refer to the measured Y, not the true value)
	precPL = estimated precision of PL[i][j] as a standard deviation

Error code may take the following values:
	0: Apparently OK
	1: Number of data points, nx, not a positive integer.
	2: Locations x not in ascending order.
	3: Numbers of observations nobs not all positive (they may be zero).
	4: Brownian motion variance B must be positive.
	5: White noise variance W and measurement variance V must be positive.
		Their sum must be strictly positive.
	6: Number of limits, nL must be positive (may be zero).

*/
{
INT i,j,seq;
double *precdf,*precy,*below,*diag,*pleft,*pright,temp;
/*	Precisions of data in precy */
/*	Precisions of fertility differences in precdf */
double h,Bh,C,E,D,A,CC;
double lin,rQ,leftsqrtrQ,rightsqrtrQ,Q3,Q5,Q3l,Q5l,Q5l2,Q7l2,Q7l3,AC,AD,ACD;
double z1left0,z1left1,z2left0,z2left1,z3left0,z3left1,z4left0,z4left1;
double z1right0,z1right1,z2right0,z2right1,z3right0,z3right1,z4right0,z4right1;
double Amidx,Arightx,ACx,G,zleft,zright,z1,z2,z3,z4,z1sq,z1cubed,Fleft,Fright,F1,F2,F3,F4;
double F1left,F2left,F3left,F4left,F1right,F2right,F3right,F4right;
double halfh,z1mid0,z1mid1,z2mid0,z2mid1,z3mid0,z3mid1,z4mid0,z4mid1;
double zmid,midsqrtrQ,Fmid,F1mid,F2mid,F3mid,F4mid;
double error1,error2;

*errorcode=0;
if(*nx <= 0){*errorcode=1;return;}
if((*W<0) | (*V<0) | ((*W+*V) <= 0) ){*errorcode=5;return;}
if(*B <= 0){*errorcode=4;return;}
if(*nL < 0){*errorcode=6;return;}

precdf = MAKEARRAY(*nx,double);
precy  = MAKEARRAY(*nx,double);
below  = MAKEARRAY(*nx,double);
diag   = MAKEARRAY(*nx,double);
pleft = MAKEARRAY(*nx,double);
pright= MAKEARRAY(*nx,double);

if (precdf && precy && below && diag && pleft && pright) {
	/*	Compute precisions for WEBM calculations */
	temp=1./(*W+*V);
	for (i=0; i<*nx-1; i++){
		if(nobs[i]<0){*errorcode=3;return;}
		precy[i]=nobs[i]*temp;
		precdf[i]=1/(*B *(x[i+1]-x[i]));
		if(precdf[i]<0){*errorcode=2;return;}
	}
	i=*nx-1; precy[i]=nobs[i]*temp;

	/*	Factorize matrix as LD(L)transpose */
	diag[0]=precdf[0]+precy[0];
	for (i=0; i<*nx-1; i++){
		below[i]=-precdf[i]/diag[i];
		diag[i+1]=precdf[i]+precdf[i+1]+precy[i+1]+below[i]*precdf[i];
	}
	below[*nx-2]=-precdf[*nx-2]/diag[*nx-2];
	diag[*nx-1]=precdf[*nx-2]+precy[*nx-1]+below[*nx-2]*precdf[*nx-2];

	/*	Multiply y by precy then by L inverse */
	webm[0]=y[0]*precy[0];
	pleft[0]=precy[0];
	for (i=1; i<*nx; i++){
		webm[i]=y[i]*precy[i]-webm[i-1]*below[i-1];
		pleft[i]=precy[i]+(pleft[i-1]*precdf[i-1])/(pleft[i-1]+precdf[i-1]);
	}
	/*	Multiply webm by D inverse then by L transpose inverse */
	temp=1/diag[*nx-1];
	varebm[*nx-1]=temp;
	webm[*nx-1]=webm[*nx-1]*temp;
	pright[*nx-1]=precy[*nx-1];
	for (i=*nx-2; i>=0; i=i-1){
		webm[i]=webm[i]/diag[i]-webm[i+1]*below[i];
		temp=1/diag[i]+below[i]*below[i]*temp;
		varebm[i]=temp;
		pright[i]=precy[i]+(pright[i+1]*precdf[i])/(pright[i+1]+precdf[i]);
	}

#ifdef CHECK
for (i=0; i<*nx; i++){
	PRINT("webm[IFMT]=%g\n",i,webm[i]);
	PRINT("varebm[IFMT]=%g\n",i,varebm[i]);
	PRINT("pleft[IFMT]=%g\n",i,pleft[i]);
	PRINT("pright[IFMT]=%g\n",i,pright[i]);
}
#endif

	/* Deal with intervals within specified range */
	seq=0;
	for (i=1; i<*nx && x[i] <= *maxx; i++){
		if(x[i-1] < *minx)continue;
		/* Note that code could be made slightly more efficient by
			storing normal tail probabilities */
#ifdef CHECK
PRINT("\ni = IFMT\n",i);
#endif
		h= x[i]-x[i-1];
		Bh=*B * h;
		A=(webm[i]-webm[i-1])/h;

		temp=pleft[i-1]*pright[i];
		C=-*B * *B * temp/(pleft[i-1]+pright[i]+Bh*temp);
		E=*W+varebm[i-1];
		D=(varebm[i]-varebm[i-1])/h-C*h;
#ifdef CHECK
PRINT("A = %f\n",A);
PRINT("C = %f\n",C);
PRINT("D = %f\n",D);
PRINT("E = %f\n",E);
#endif

		/* Prepare to calculate derivatives of z as linear in G */
		/* At left hand side of interval */
		lin=D;
		rQ=1/(*W+varebm[i-1]);
		leftsqrtrQ=sqrt(rQ);
		Q3=rQ*leftsqrtrQ;
		Q5=rQ*Q3;
		Q3l=Q3*lin;
		Q5l=Q5*lin;
		Q5l2=Q5l*lin;
		Q7l2=Q5l2*rQ;
		Q7l3=Q7l2*lin;
		z1left0=leftsqrtrQ*A;
		z1left1=-.5*Q3l;
		CC=C*C;
		AC=A*C;
		AD=A*D;
		ACD=AC*D;
		z2left0=-AD*Q3;
		z2left1=-C*Q3+.75*Q5l2;
		z3left0=-3*AC*Q3+2.25*AD*Q5l;
		z3left1=4.5*C*Q5l-1.875*Q7l3;
		z4left0=18.*ACD*Q5-7.5*AD*Q7l2;
		z4left1=9*CC*Q5-22.5*C*Q7l2+6.5625*lin*Q7l3*rQ;

		/* At mid point of interval */
		halfh=.5*h;
		Amidx=A*halfh;
		lin=C*h+D;
		rQ=1/(E+halfh*(D+halfh*C));
		midsqrtrQ=sqrt(rQ);
		Q3=rQ*midsqrtrQ;
		Q5=rQ*Q3;
		Q3l=Q3*lin;
		Q5l=Q5*lin;
		Q5l2=Q5l*lin;
		Q7l2=Q5l2*rQ;
		Q7l3=Q7l2*lin;
		z1mid0=midsqrtrQ*A-.5*Amidx*Q3l;
		z1mid1=-.5*Q3l;
		ACx=AC*halfh;
		z2mid0=-(3*ACx+AD)*Q3+.75*Amidx*Q5l2;
		z2mid1=-C*Q3+.75*Q5l2;
		z3mid0=-3*AC*Q3+2.25*(4*ACx+AD)*Q5l-1.875*Amidx*Q7l3;
		z3mid1=4.5*C*Q5l-1.875*Q7l3;
		temp=6.5625*lin*Q7l3*rQ;
		z4mid0=(18.*ACD+45*ACx*C)*Q5+(-7.5*AD-37.5*ACx)*Q7l2+Amidx*temp;
		z4mid1=9*CC*Q5-22.5*C*Q7l2+temp;

		/* At right hand side of interval */
		Arightx=A*h;
		lin=2*C*h+D;
		rQ=1/(*W+varebm[i]);
		rightsqrtrQ=sqrt(rQ);
		Q3=rQ*rightsqrtrQ;
		Q5=rQ*Q3;
		Q3l=Q3*lin;
		Q5l=Q5*lin;
		Q5l2=Q5l*lin;
		Q7l2=Q5l2*rQ;
		Q7l3=Q7l2*lin;
		z1right0=rightsqrtrQ*A-.5*Arightx*Q3l;
		z1right1=-.5*Q3l;
		ACx=AC*h;
		z2right0=-(3*ACx+AD)*Q3+.75*Arightx*Q5l2;
		z2right1=-C*Q3+.75*Q5l2;
		z3right0=-3*AC*Q3+2.25*(4*ACx+AD)*Q5l-1.875*Arightx*Q7l3;
		z3right1=4.5*C*Q5l-1.875*Q7l3;
		temp=6.5625*lin*Q7l3*rQ;
		z4right0=(18.*ACD+45*ACx*C)*Q5+(-7.5*AD-37.5*ACx)*Q7l2+Arightx*temp;
		z4right1=9*CC*Q5-22.5*C*Q7l2+temp;

		/* loop over nL limits */
		for (j=0; j<*nL; j++){
			/* Find derivatives on left of region */
			G=webm[i-1]-yL[j];
			zleft=G*leftsqrtrQ;
			z1=z1left0+G*z1left1;
			z2=z2left0+G*z2left1;
			z3=z3left0+G*z3left1;
			z4=z4left0+G*z4left1;
			if (uL[j]){rightnormaltail(zleft,&F1,&Fleft);}
			else{leftnormaltail(zleft,&F1,&Fleft);}
			F2=-zleft*F1;
			F3=-F1-zleft*F2;
			F4=-2*F2-zleft*F3;
			F1left=F1*z1;
			z1sq=z1*z1;
			z1cubed=z1*z1sq;
			F2left=F2*z1sq+F1*z2;
			F3left=F3*z1cubed+3*F2*z2*z1+F1*z3;
			F4left=F4*z1*z1cubed+6*F3*z2*z1sq+F2*(3*z2*z2+4*z3*z1)+F1*z4;
#ifdef CHECK
PRINT("\nj = IFMT\n",j);
PRINT("G = %f\n",G);
PRINT("zleft and derivatives %f %f %f %f %f\n",zleft,z1,z2,z3,z4);
PRINT("Fleft and deriv. wrt z %f %f %f %f %f\n",Fleft,F1,F2,F3,F4);
PRINT("Derivatives wrt x %f %f %f %f\n",F1left,F2left,F3left,F4left);
#endif

			/* Find derivatives at mid point of region */
			zmid=(Amidx+G)*midsqrtrQ;
			z1=z1mid0+G*z1mid1;
			z2=z2mid0+G*z2mid1;
			z3=z3mid0+G*z3mid1;
			z4=z4mid0+G*z4mid1;
			if (uL[j]){rightnormaltail(zmid,&F1,&Fmid);}
			else{leftnormaltail(zmid,&F1,&Fmid);}
			F2=-zmid*F1;
			F3=-F1-zmid*F2;
			F4=-2*F2-zmid*F3;
			F1mid=F1*z1;
			z1sq=z1*z1;
			z1cubed=z1*z1sq;
			F2mid=F2*z1sq+F1*z2;
			F3mid=F3*z1cubed+3*F2*z2*z1+F1*z3;
			F4mid=F4*z1*z1cubed+6*F3*z2*z1sq+F2*(3*z2*z2+4*z3*z1)+F1*z4;
#ifdef CHECK
PRINT("\nj = IFMT\n",j);
PRINT("G = %f\n",G);
PRINT("zmid and derivatives %f %f %f %f %f\n",zmid,z1,z2,z3,z4);
PRINT("Fmid and deriv. wrt z %f %f %f %f %f\n",Fmid,F1,F2,F3,F4);
PRINT("Derivatives wrt x %f %f %f %f\n",F1mid,F2mid,F3mid,F4mid);
#endif

			/* Find derivatives on right of region */
			zright=(Arightx+G)*rightsqrtrQ;
			z1=z1right0+G*z1right1;
			z2=z2right0+G*z2right1;
			z3=z3right0+G*z3right1;
			z4=z4right0+G*z4right1;
			if (uL[j]){rightnormaltail(zright,&F1,&Fright);}
			else{leftnormaltail(zright,&F1,&Fright);}
			F2=-zright*F1;
			F3=-F1-zright*F2;
			F4=-2*F2-zright*F3;
			F1right=F1*z1;
			z1sq=z1*z1;
			z1cubed=z1*z1sq;
			F2right=F2*z1sq+F1*z2;
			F3right=F3*z1cubed+3*F2*z2*z1+F1*z3;
			F4right=F4*z1*z1cubed+6*F3*z2*z1sq+F2*(3*z2*z2+4*z3*z1)+F1*z4;
#ifdef CHECK
PRINT("\nj = IFMT\n",j);
PRINT("G = %f\n",G);
PRINT("zright and derivatives %f %f %f %f %f\n",zright,z1,z2,z3,z4);
PRINT("Fright and deriv. wrt z %f %f %f %f %f\n",Fright,F1,F2,F3,F4);
PRINT("Derivatives wrt x %f %f %f %f\n",F1right,F2right,F3right,F4right);
#endif

			/* Now compute integral and measure of precision */
PL[seq]=(54495/270270.)*(Fleft+Fright)+(161280/270270.)*Fmid+
	halfh*( (9450/270270.)*(F1left-F1right)+
	halfh*( (885/270270.)*(F2left+F2right)+(7680/270270.)*F2mid+
	halfh*( (45/270270.)*(F3left-F3right)+
	halfh*( (1/270270.)*(F4left+F4right)+(64/270270.)*F4mid))));
error1=fabs((5760/270270.)*(Fleft+Fright)+(-11520/270270.)*Fmid+
	halfh*( (2295/270270.)*(F1left-F1right)+
	halfh*( (375/270270.)*(F2left+F2right)+(-1920/270270.)*F2mid+
	halfh*( (30/270270.)*(F3left-F3right)+
	halfh*( (1/270270.)*(F4left+F4right)+(-32/270270.)*F4mid)))));
error2=.1*fabs((-21915/270270.)*Fleft+(23040/270270.)*Fmid+(-1125/270270.)*Fright+
	halfh*( (-9225/270270.)*F1left+(-11520/270270.)*F1mid+(-45/270270.)*F1right+
	halfh*( (-1605/270270.)*F2left+(3840/270270.)*F2mid+(105/270270.)*F2right+
	halfh*( (-138/270270.)*F3left+(-384/270270.)*F3mid+(-18/270270.)*F3right+
	halfh*( (-5/270270.)*F4left+(64/270270.)*F4mid+(1/270270.)*F4right)))));
	if(error1>error2){precPL[seq]=error1;}
	else{precPL[seq]=error2;}
		seq++;

		}

	}
}
if (precdf) free (precdf);
if (precy) free (precy);
if (below) free (below);
if (diag) free (diag);
if (pleft) free (pleft);
if (pright) free (pright);
return;
}
/*====================================================================== */

<html>
<body>
<p align="center"><font size="-2">
© Copyright 2012, CSIRO Australia<br>
Use of this web site and information available from<br>
it is subject to our </font>
<a href="http://www.csiro.au/legalnotices/disclaimer.html">
<font size="-2">Legal Notice and Disclaimer </font></a> and <a href="http://www.csiro.au/csiro/channel/pchb1,,.html"><font size="-2">Privacy
Statement</font></a>
</p>
<!-- Note: this is sourced from 'footer2.inc' -->
<!-- Clicky code for visitation statistices follows -->
<a title="Clicky Web Analytics" href="http://getclicky.com/16636"><!-- <img alt="Clicky Web Analytics" src="http://static.getclicky.com/media/links/badge.gif" border="0" /> --></a>
<script src="http://static.getclicky.com/26219.js" type="text/javascript"></script>
<!-- <noscript><p><img alt="Clicky" src="http://in.getclicky.com/26219-db2.gif" /></p></noscript> -->
</body>
</html>