/*------------------------------------------------------------------------------*
 * File Name: OCN_g04.h															*
 * Creation: SDB 5/22/2001														*
 * Purpose: Origin C Header for NAG functions									*
 * Copyright (c) OriginLab Corp.	2001										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *																				*
 *------------------------------------------------------------------------------*/

#ifndef _O_NAG_G04_H
#define _O_NAG_G04_H

//#importdll "ONAG" // NAG DLL prepared by OriginLab
#pragma dll(ONAG)	
#include <NAG\nag_types.h>


/*	begin proto */

/**	g04bbc
		computes the analysis of variance and treatment means and standard
		errors for a randomized block or completely randomized design

Example1:
	Assume "Data1" Worksheet contains 7 columns and 36 rows and it is double type, "Data2" Worksheet 
	contains 2 columns and 36 rows and it is integer type. The first column contains 30 data in every
	worksheet. The result will write in rest columns.
	
	int n = 30, nt = 6, iblock = 10, tdc = 6, irdf = 0;
	double tol = 5e-6;
	double gmean;
	//Attach 7 Datasets to 7 columns of "Data1".
	Dataset dy("data1",0);
	Dataset dbmean("data1",1);
	Dataset dtmean("data1",2);
	Dataset dtable("data1",3);
	Dataset dc("data1",4);
	Dataset dr("data1",5);
	Dataset def("data1",6);
	//Attach 2 Datasets to 2 columns of "Data2".
	Dataset dit("data2",0);
	Dataset direp("data2",1);
	
	dy.SetSize(30);
	dbmean.SetSize(10);
	dtmean.SetSize(6);
	dtable.SetSize(20);
	dc.SetSize(36);
	dit.SetSize(30);
	dr.SetSize(30);
	def.SetSize(6);
	direp.SetSize(6);
	
	vector y = dy, bmean = dbmean, tmean = dtmean, table = dtable, c = dc, r = dr, ef = def;
	vector<int> it=dit, irep=direp;
	BOOL sucess = nag_anova_random(	n, y, Nag_SerialBlocks, iblock, nt, it, &gmean, bmean, tmean, 
									table,  c, tdc, irep, r, ef, tol, irdf );
	//put the result to the Worksheet "data1" and Worksheet "data2".
	dbmean = bmean;
	dtmean = tmean;
	dtable = table;
	dc = c;
	dr = r;
	def = ef;
	dit = it;
	direp = irep;
	
Example2:
	The data, given by John and Quenouille (1977), is for a balanced incomplete block design with
	10 blocks and 6 treatments and with 3 plots per block. The observations are the degree of pain
	experienced and the treatments are penicillin of different potency. 

void test_nag_anova_random()
{

    double tmean[6];
	double table[20];
	double ef[6], r[30];
	double gmean;
	double tol = 5e-6;
	double bmean[10];
	double c[36]; 
	int irep[6];
	int irdf = 0;
	int i, j;
	int tdc = 6;
	int sucess;
	int n = 30, nt = 6, iblock = 10;
	double y[30] = {1.00, 5.00, 4.00, 5.00, 10.00, 6.00, 2.00, 9.00, 3.00, 4.00, 
	8.00, 6.00, 2.00, 4.00, 7.00, 6.00, 7.00, 5.00, 5.00, 7.00, 2.00, 7.00, 
	2.00, 4.00, 8.00, 4.00, 2.00, 10.00, 8.00, 7.00};	
	int it[30] = {1, 2, 3, 1, 2, 4, 1, 3, 5, 1, 4, 6, 1, 5, 6, 2, 3, 6, 2, 4, 5, 
	2, 5, 6, 3, 4, 5, 3, 4, 6};

	printf("the input data is following\n");
	printf("n = 30, nt = 6, iblock = 10\n");
	printf("y:  \n");
	for( i = 0; i < 30; i++)
	{
		printf("%f ",y[i]);
		if((i + 1) % 10 == 0)
		printf("\n");
	}
	printf("\nit:\n");
	for(i = 0; i < 30; i++)
	{
		printf("%d ",it[i]);
		if((i + 1) % 10 == 0)
		printf("\n");
	}		  
	
	sucess = nag_anova_random(n, y, Nag_SerialBlocks, iblock, nt, it, &gmean, bmean, tmean, 
								table,  c, tdc, irep, r, ef, tol, irdf);
		
	if(sucess == 0)
	{
		printf("the function is called sucessfully\n\n\n");
		printf("the results are the following:\n\n");
		printf("ANOVA table\n");
		printf(" Source       df       SS  	    MS     F         Prob\n");
		printf("Blocks     ");
		for(i = 0; i < 5; i++)
		{
			printf("%10.4f",table[i]);
		}
		printf("\nTreatments ");
		for(i = 5; i < 10; i++)
		{
			printf("%10.4f",table[i]);
		}
		
		printf("\nResidual   ");
		for(i = 11; i < 14; i++)
		{
			printf("%10.4f",table[i]);
		}
		
		printf("\nTotal      ");
		for(i = 15; i < 17; i++)
		{
			printf("%10.4f",table[i]);
		}
		
		printf("\n\nEfficency Fators     \n");
		for(i = 0; i < 6; i++)
			printf("%10.5f",ef[i]);
		printf("\n\nGrand Mean%10.5f\n",gmean);
		printf("\nTreatment Means\n");
		for(i = 0; i < 6; i++)
			printf("%10.5f", tmean[i]);
		printf("\nthe standard erros of difference between means\n");
		for( i = 1; i < 6; i++)
		{
			for( j = 0; j < i; j++)
				printf("%10.5f",c[i*6+j]);
			printf("\n");		
		} 		
				
	}
	else
		printf("the function has not been called successfully");	 
}

	The output is following:

the input data is following
n = 30, nt = 6, iblock = 10
y:  
1.000000 5.000000 4.000000 5.000000 10.000000 6.000000 2.000000 9.000000 3.000000 4.000000 
8.000000 6.000000 2.000000 4.000000 7.000000 6.000000 7.000000 5.000000 5.000000 7.000000 
2.000000 7.000000 2.000000 4.000000 8.000000 4.000000 2.000000 10.000000 8.000000 7.000000 

it:
1 2 3 1 2 4 1 3 5 1 
4 6 1 5 6 2 3 6 2 4 
5 2 5 6 3 4 5 3 4 6 
the function is called sucessfully


the results are the following:

ANOVA table
 Source       df       SS  	    MS     F         Prob
Blocks         9.0000   60.0000    6.6667    4.7872    0.0039
Treatments     5.0000  101.7778   20.3556   14.6170    0.0000
Residual      20.8889    1.3926    0.0000
Total         29.0000  182.6667

Efficency Fators     
   0.00000   0.80000   0.80000   0.80000   0.80000   0.80000

Grand Mean   5.33333

Treatment Means
   2.50000   7.25000   8.08333   5.91667   2.91667   5.33333
the standard erros of difference between means
   0.83444
   0.83444   0.83444
   0.83444   0.83444   0.83444
   0.83444   0.83444   0.83444   0.83444
   0.83444   0.83444   0.83444   0.83444   0.83444


Return:
	The function returns NAG error code, 0 if no error.
	
	NE_BAD_PARAM(70): On entry, parameter blocks has illegal value.
	NE_INT_ARG_LT(11): On entry, n must not be less than 2: n = _value_.  On entry, nt must not be less than 1: nt = _value_.  On entry, irdf must not be less than 0: irdf = _value_.
	NE_REAL_ARG_LT(5): On entry, tol must not be less than 0.0: tol = _value_.
	NE_INT(90): On entry, iblock = _value_.  Constraint: iblock = 2when blocks = Nag NoBlocks.  On entry, nt = _value_.  Constraint: nt = 2when blocks = Nag NoBlocks.
	NE_2_INT_ARG_LT(17): On entry, tdc = _value_ while nt = _value_.  These parameters must satisfy tdc = nt.
	NE_INT_2(91): On entry, n = _value_, iblock = _value_.  Constraint: when iblock = 2, n must be a multiple of iblock.
	NE_INTARR(111): On entry, it[_value_] = _value_.  Constraint: 1 = it[i - 1] = nt, for i = 1, 2, . . . ,n.
	NE_IT_ARRAY(533): No value of it[j - 1] = j for some j = 1, 2, . . . ,nt.
	NE_ARRAY_CONSTANT(534): On entry, the elements of the array y are constant.
	NE_G04BB_STDERR(535): A computed standard error is zero due to rounding errors.  This is an unlikely error exit.
	NE_G04BB_DESIGN(537): The design is disconnected; the standard errors may not be valid. The design may be nested.
	NE_G04BB_TREAT(538): The treatments are totally confounded with blocks, so the treatment sum of squares and degrees of freedom are zero. The analysis of variance table is not computed, except for block and total sums of squares and degrees of freedom.  
	NE_G04BB_CONV(536): The eigenvalue computation has failed to converge.  This is an unlikely error exit.
	NE_G04BB_RES_DF(539): The residual degrees of freedom or the residual sum of squares are zero, columns 3, 4 and 5 of the analysis of variance table will not be computed and the matrix of standard errors and covariances, C, will not be scaled by s or the square of s.
	NE_ALLOC_FAIL(73): Memory allocation failed.
	
	successfully call of the nag_anova_random function.
*/

int nag_anova_random(
	int n, 		//the number of observations n
	const double y[], //the observations in the order as described by blocks and nt
	Nag_Blocks blocks, 
	int iblock, //indicates the number of blocks
	int nt, 	//indicates the number of treatments
	int it[], 	//indicates the treatments used
	double *gmean, //the grand mean
	double bmean[], //if blocks = Nag SerialBlocks or Nag ParallelBlocks, bmean contains the mean for the jth block
	double tmean[], //contains the mean of the treatments
	double table[], // the analysis of variance table
	double c[], 	// the variance/covariance matrix of the treatment effects
	int tdc, 	//the second dimension of the array c
	int irep[], //the treatment replications 
	double r[], //the residuals
	double ef[], //the canonical efficiency factors
	double tol, //the tolerance value used to check for zero eigenvalues of the matrix.
	int irdf 	//an adjustment to the degrees of freedom for the residual and total
);	// General block design or completely randomized design.


/**	g04bcc
		computes the analysis of variance for a general row and column design
		together with the treatment means and standard error
		
Example:
	The data for a 5 x 5 Latin square is inputted and the ANOVA and treatment means are computed and printed.

void test_nag_anova_row_col()
{
	double c[25], c_b20 = 1e-5, cmean[5], ef[5], r[25], gmean;
	double	rmean[5], rpmean[1], tmean[5], table[30];
	int irep[5], c_0 = 0, i, j;
	int nrep = 1, nt = 5, nrow = 5, ncol = 5;
	int sucess;
	double y[25] = {6.67, 7.15, 8.29, 8.95, 9.62, 5.40, 4.77, 5.40, 7.54, 6.93,
	7.32, 8.53, 8.50, 9.99, 9.68, 4.92, 5.00, 7.29, 7.85, 7.08, 4.88, 6.16,
	7.83, 5.38, 8.51};
	int it[25] = {5, 4, 1, 3, 2, 2, 5, 4, 1, 3, 3, 2, 5, 4, 1, 1, 3, 2, 5, 4, 4,
	1, 3, 2, 5};

	sucess = nag_anova_row_col(nrep, nrow, ncol, y, nt, it, &gmean, tmean, table, c,
							 nt, irep, rpmean, rmean, cmean, r, ef, c_b20, c_0);

	if(sucess == 0)
	{
		printf("\nthe input data are following\n");
		printf("nrep = 1, nt = 5, nrow = 5 , ncol = 5\n");
		printf("y:\n");
		for(i = 0; i < 25; i++)
		{
			printf("%3.2f   ",y[i]);
			if((i + 1) % 5 == 0)
				printf("\n");
		}

		printf("it:\n");
		for(i = 0; i < 25; i++)
		{
			printf("%d   ",it[i]);
			if((i + 1) % 5 == 0)
				printf("\n");
		}
        
		printf("\n\n\n");
		printf("ANOVA TABLE\n");
		

		printf("Rows       ");
		for(i = 5; i < 10; i++)
		{
			if(i == 5)
				printf("%3.0f   ", table[i]);
			else
				printf("%6.4f   ", table[i]);
		}

		printf("\nCloumns    ");
		for(i = 10; i < 15; i++)
		{
			if(i == 10)
				printf("%3.0f  ", table[i]);
			else
				printf("%6.4f  ", table[i]);
		}

		printf("\nTreatments ");
		for(i = 15; i <20; i++)
		{
			if(i == 15)
				printf("%3.0f  ", table[i]);
			else
				printf("%6.4f  ", table[i]);
		}

		printf("\nResidual   ");
		for(i = 20; i <23; i++)
		{
			if(i == 20)
				printf("%3.0f  ", table[i]);
			else
				printf("%6.4f  ", table[i]);
		}

		printf("\nTotal      ");
		for(i = 25; i < 27; i++)
		{
			if(i == 25)
				printf("%3.0f  ", table[i]);
			else
				printf("%6.4f  ", table[i]);
		}

		printf("\nTreatment means\n");
		for(i = 0; i < 5; i++)
			printf("%10.4f", tmean[i]);

		printf("\n\nS.E. of difference  (orthogonal design) = %10.4f\n", c[5]);
	}
	else
		printf("the function is not sucessfully called");
}
	
	The output is following:
		
the input data are following
nrep = 1, nt = 5, nrow = 5 , ncol = 5
y:
6.67   7.15   8.29   8.95   9.62   
5.40   4.77   5.40   7.54   6.93   
7.32   8.53   8.50   9.99   9.68   
4.92   5.00   7.29   7.85   7.08   
4.88   6.16   7.83   5.38   8.51   
it:
5   4   1   3   2   
2   5   4   1   3   
3   2   5   4   1   
1   3   2   5   4   
4   1   3   2   5   



ANOVA TABLE
Rows         4   29.4231   7.3558   9.0266   0.0013   
Cloumns      4  22.9950  5.7487  7.0545  0.0037  
Treatments   4  0.5423  0.1356  0.1664  0.9514  
Residual    12  9.7788  0.8149  
Total       24  62.7392  
Treatment means
    7.3180    7.2440    7.2060    6.9000    7.2600

S.E. of difference  (orthogonal design) =     0.5709

	
Return:

	The function returns NAG error code, 0 if no error.
	
	NE_INT_ARG_LT(11):  On entry, nrep must not be less than 1: nrep = <value>.  On entry, nrow must not be less than 2: nrow = <value>.  On entry, ncol must not be less than 2: ncol = <value>.  On entry, nt must not be less than 1: nt = <value>.  On entry, irdf must not be less than 0: irdf = <value>.
	NE_2_INT_ARG_LT(17):  On entry, tdc = <value> while nt = <value>.  These parameters must satisfy tdc >= nt.
	NE_REAL_ARG_LT(5):   On entry, tol must not be less than 0: tol = <value>.
	NE_ARRAY_CONS(58):  The contents of array it are not valid.  Constraint: if nt >= 2, 1 <= it[i] <= nt, i = 0, 1, 2, ..., nrep*nrow*ncol.  Constraint: some value of it = j for all j = 1, 2, ..., nt.
	NE_ALLOC_FAIL(73):  Memory allocation failed.
	NE_ARRAY_CONSTANT(534): On entry, the elements of the array y are constant.
	NE_G04BC_ST_ERR(543): A computed standard error is zero due to rounding errors, or the eigenvalue computation failed to converge.  Both are unlikely errors.
	NE_G04BC_REPS(544): The treatments are totally confounded with replicates, rows and columns, so the treatment sum of squares and degrees of freedom are zero.  The analysis of variance table is not computed, except for replicate, row, column, total sum of squares and degrees of freedom.
	NE_G04BC_RESD(545): The residual degrees of freedom or the residual sum of squares are zero, columns 3, 4, and 5 of the analysis of variance table will not be computed and the matrix of standard errors and covariances, c, will not be scaled.
	NE_G04BC_DISCON(546): The design is disconnected, the standard errors may not be valid.  The design may have a nested structure.
	NE_INTERNAL_ERROR(74):  An internal error has occurred in this function.  Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
		
	successfully calll of the nag_anova_row_col function.	
*/

	int nag_anova_row_col(
		int nrep, //the number of replicates.
		int nrow, //the number of rows per replicate.
		int ncol, //the number of columns per replicate.
		const double y[], // the observation of the k column of the jth row of the ith replicate.
		int nt,			//the number of treatments
		const int it[], // indicate which of the nt treatments unit i received, i = 1, 2, 3, ...., n.
		double * gmean, //the grand mean.
		double tmean[], //contains the (adjusted) mean for the lth treatment, l = 1, 2, 3, ..., nt.
		double table[], //the analysis of variance table
		double c[], 	// the variance/covariance matrix of the treatment effects.
		int tdc, 		//the second dimension of c.
		int irep[], 	//the treatment replications.
		double rpmean[], //contains the mean for the ith replicate, i = 1, 2, ..., b.
		double rmean[], //contains the mean for the jth row, j= 1, 2, ..., r.
		double cmean[], //contains the mean for the kth column, k = 1, 2, ..., c.
		double r[], 	//contains the residuals
		double ef[], 	//the canonical efficiancy factor
		double tol, 	//the tolerance value used to check for zero eigenvalue of the matrix Omega.
		int irdf 		//an adjustment to  the degrees of freedom for the residual and total.
	);	// Analysis of variance, general row and column design, treatment means and standard errors.


/**	g04cac
		computes an analysis of variance table and treatment means for a
		complete factorial design.
	Example:   
		These are the program data:
		The data, given by John and Quenouille (1977), is for the yield of turnips (54 observation) for a 
		factorial experiment with two factors, the amount of phosphate with 6 levels and the amount of 
		liming with 3 levels.  The design was replicated in 3 blocks. 

		void test_nag_anova_factorial()
		{
			double r[54];
			double *bmean;
			double *e;
			double *semean;
			double *table;
			double *tmean;
			int c__27 = 27;
			int lfac[2] = {6, 3};
			int *imean;
			int mterm = 6;
			int nfac, irdf;
			int i, j, k, l, n;
			int num;
			int inter, nblock;
			int itotal, ntreat;
			n = 54;
			nblock = 3;
			nfac = 2;
			inter = 2;
			double y[54] = {274, 361, 253, 325, 317, 339, 326, 402, 336, 379, 345, 361,
			352, 334, 318, 339, 393, 358, 350, 340, 203, 397, 356, 298, 382, 376, 355,
			418, 387, 379, 432, 339, 293, 322, 417, 342, 82, 297, 133, 306, 352, 361,
			220, 333, 270, 388, 379, 274, 336, 307, 266, 389, 333, 353};

			irdf = 0;
			nag_anova_factorial(n, y, nfac, lfac, nblock, inter, irdf, &mterm, &table, &tmean,
				&c__27, &e, &imean, &semean, &bmean, r);

			itotal = mterm;
			printf("\n ANOVA table\n\n");
			printf(" Source     df      SS      MS      F     	Prob\n\n");
			k = 0;
			if (nblock > 1)
			{
				++k;
				printf("%s"," Blocks ");
				for (j = 1; j <= 5; ++j)
				printf("%10.2f  ",table[((1)-1) * ( 5) + ((j)-1)]);
				printf("\n");
			}
			ntreat = mterm - 2 - k;
			for (i = 1; i <= ntreat; ++i)
			{
				printf("%s%2i  "," Effect ",i);
				for (j = 1; j <= 5; ++j)
				printf("%10.2f ",table[((k+i)-1) * ( 5) + ((j)-1)]);
				printf("\n");
			}
			printf("%s "," Residual ");
			for (j = 1; j <= 3; ++j)
				printf("%10.2f  ",table[((mterm -1)-1) * ( 5) + ((j)-1)]);
			printf("\n");
			printf("%s "," Total ");
			for (j = 1; j <= 2; ++j)
				printf("%10.2f  ",table[((mterm)-1) * ( 5) + ((j)-1)]);

			printf("\n");
			printf("\n");
			printf(" Treatment Means and Standard Errors ");
			printf("\n");
			printf("\n");
			k = 0;
			for (i = 0; i < ntreat; ++i)
			{
				l = imean[i];
				printf("%s%2i"," Effect ",i+1);
				printf("\n");
				printf("\n");
				num=1;
				for (j = k; j < l; ++j)
				{
					printf("%10.2f ",tmean[j]);
					if(num % 8 == 0)
						printf("\n");
					num++;
				}
				printf("\n");
				printf("\n%s%10.2f\n\n"," SE of difference in means = ",semean[i]);
				k = l;
			}

			nag_anova_factorial_free(&table, &tmean, &e, &imean, &semean, &bmean);
		}


		The output is following:

		ANOVA TABLE
		Source             df           SS         MS         F        Prob


		Blocks           2.0000    30118.7778     15059.3889  7.6848   0.0018   
		Effect0          5.00      73008.17       14601.63    7.45     0.00   
		Effect1          2.00      21596.33       10798.17    5.51     0.01   
		Effect2          10.00     31191.67       3119.17     1.59     0.15   
		Residual         34.0000   66627.8889     1959.6438  
		total            53.0000   222542.8333  

		Treatment Means and Standard Errors


		Effect 1
			254.78,   339.00,   333.33,   367.78,   330.78,   360.67

		SE of difference inm means =      20.87

		Effect 2
			334.28,	  353.78,   305.11

		SE of difference inm means =      14.76

		Effect 3
			235.33,   332.67,   196.33,   342.67,   341.67,   332.67,   309.33,   370.33,
			320.33,   395.00,   370.33,   338.00,   373.33,   326.67,   292.33,   350.00,
			381.00,   351.00

		SE of difference in means =      36.14


	Return:
		The function returns NAG error code, 0 if no error.
	
		NE_INT_ARG_LT(11):  On entry, n must not be less than 4: n = _value_.  On entry, nfac must not be less than 1: nfac = _value_.  On entry, nblock must not be less than 0: nblock = _value_.  On entry, inter must not be less than 0: inter = _value_.  On entry, irdf must not be less than 0: irdf = _value_.  
		NE_2_INT_ARG_GT(19):  On entry, inter = _value_ while nfac = _value_.  These parameters must satisfy inter = nfac.
		NE_INTARR(111): On entry, lfac[_value_] = _value_.  Constraint: lfac[i - 1] = 2 for i = 1, 2, . . . ,nfac.
		NE_INT_2(91):  On entry, nblock = _value_, n = _value_.  Constraint: n must be a multiple of nblock, when nblock > 1.
		NE_PLOT_TREAT(541): The number of plots per block is not a multiple of the number of treatment combinations.
		NE_ARRAY_CONSTANT(534): On entry, the elements of the array y are constant.
		NE_G04CA_RES_DF(542): There are no degrees of freedom for the residual or the residual sum of squares is zero. In either case the standard errors and F-statistics cannot be computed.
		NE_ALLOC_FAIL(73):  Memory allocation failed.

		successfully call of the nag_anova_factorial function.	
*/
int nag_anova_factorial(
	int n,	// the number of observations
	const double y[], // the observations in standard order
	int nfac,	// the number of factors
	int lfac[],	// the number of levels for the ith factor
	int nblock,	// the number of blocks
	int inter, // the maximum number of factors in an interaction term, inter=0 or 1 if no interaction terms are computed.
	int irdf,	// the adjustment to the residual and total degrees of freedom
	int *mterm,	// the number of terms in the analysis of variance table
	double **table,	// the table of analysis of variance.
	double **tmean,	// contains the treatment means. 
	int *maxt,	// the number of treatment means that have been computed.
	double **e,	// contains the estimated effects in the same order as for the means in tmean.
	int **imean,	// indicates the position of the e.ect means in tmean. 
 	double **semean, // the standard error of the di.erence between means corresponding to the ith treatment effect in the ANOVA table.
	double **bmean,	// bmean[0] contains the grand mean, bmean[1] up to bmean[nblock] contain the block means.
	double r[]	// the residuals
);	// Complete factorial design

/** g04czc
    frees Nag allocated memory to some parameters in nag anova factorial
	
	Example:   
		These are the program data:
		The data, given by John and Quenouille (1977), is for the yield of turnips (54 observation) for a 
		factorial experiment with two factors, the amount of phosphate with 6 levels and the amount of 
		liming with 3 levels.  The design was replicated in 3 blocks. 

		void test_nag_anova_factorial()
		{
			double r[54];
			double *bmean;
			double *e;
			double *semean;
			double *table;
			double *tmean;
			int c__27 = 27;
			int lfac[2] = {6, 3};
			int *imean;
			int mterm = 6;
			int nfac, irdf;
			int i, j, k, l, n;
			int num;
			int inter, nblock;
			int itotal, ntreat;
			n = 54;
			nblock = 3;
			nfac = 2;
			inter = 2;
			double y[54] = {274, 361, 253, 325, 317, 339, 326, 402, 336, 379, 345, 361,
			352, 334, 318, 339, 393, 358, 350, 340, 203, 397, 356, 298, 382, 376, 355,
			418, 387, 379, 432, 339, 293, 322, 417, 342, 82, 297, 133, 306, 352, 361,
			220, 333, 270, 388, 379, 274, 336, 307, 266, 389, 333, 353};

			irdf = 0;
			nag_anova_factorial(n, y, nfac, lfac, nblock, inter, irdf, &mterm, &table, &tmean,
				&c__27, &e, &imean, &semean, &bmean, r);

			itotal = mterm;
			printf("\n ANOVA table\n\n");
			printf(" Source     df      SS      MS      F     	Prob\n\n");
			k = 0;
			if (nblock > 1)
			{
				++k;
				printf("%s"," Blocks ");
				for (j = 1; j <= 5; ++j)
				printf("%10.2f  ",table[((1)-1) * ( 5) + ((j)-1)]);
				printf("\n");
			}
			ntreat = mterm - 2 - k;
			for (i = 1; i <= ntreat; ++i)
			{
				printf("%s%2i  "," Effect ",i);
				for (j = 1; j <= 5; ++j)
				printf("%10.2f ",table[((k+i)-1) * ( 5) + ((j)-1)]);
				printf("\n");
			}
			printf("%s "," Residual ");
			for (j = 1; j <= 3; ++j)
				printf("%10.2f  ",table[((mterm -1)-1) * ( 5) + ((j)-1)]);
			printf("\n");
			printf("%s "," Total ");
			for (j = 1; j <= 2; ++j)
				printf("%10.2f  ",table[((mterm)-1) * ( 5) + ((j)-1)]);

			printf("\n");
			printf("\n");
			printf(" Treatment Means and Standard Errors ");
			printf("\n");
			printf("\n");
			k = 0;
			for (i = 0; i < ntreat; ++i)
			{
				l = imean[i];
				printf("%s%2i"," Effect ",i+1);
				printf("\n");
				printf("\n");
				num=1;
				for (j = k; j < l; ++j)
				{
					printf("%10.2f ",tmean[j]);
					if(num % 8 == 0)
						printf("\n");
					num++;
				}
				printf("\n");
				printf("\n%s%10.2f\n\n"," SE of difference in means = ",semean[i]);
				k = l;
			}

			nag_anova_factorial_free(&table, &tmean, &e, &imean, &semean, &bmean);
		}


		The output is following:

		ANOVA TABLE
		Source             df           SS         MS         F        Prob


		Blocks           2.0000    30118.7778     15059.3889  7.6848   0.0018   
		Effect0          5.00      73008.17       14601.63    7.45     0.00   
		Effect1          2.00      21596.33       10798.17    5.51     0.01   
		Effect2          10.00     31191.67       3119.17     1.59     0.15   
		Residual         34.0000   66627.8889     1959.6438  
		total            53.0000   222542.8333  

		Treatment Means and Standard Errors


		Effect 1
			254.78,   339.00,   333.33,   367.78,   330.78,   360.67

		SE of difference inm means =      20.87

		Effect 2
			334.28,	  353.78,   305.11

		SE of difference inm means =      14.76

		Effect 3
			235.33,   332.67,   196.33,   342.67,   341.67,   332.67,   309.33,   370.33,
			320.33,   395.00,   370.33,   338.00,   373.33,   326.67,   292.33,   350.00,
			381.00,   351.00

		SE of difference in means =      36.14
		
	Parameter:
		Input:pointers to which memory has been allocated internally in nag anova factorial(g04cac).
        Output:The memory allocated to each of the pointers is freed and the pointers are set to NULL.
	Return:
		None
*/ 
void nag_anova_factorial_free(
	 double **table, 
	 double **tmean, 
	 double **e,
	 int **imean, 
	 double **semean, 
	 double **bmean
);


/**	g04dbc
		Computes confidence intervals for differences between means computed
		by nag_anova_random (g04bbc) or nag_anova_row_col (g04bcc).
	
Example:
	In the example taken from Winer (1970) a completely randomised design with unequal treatment replication 
	is analysed using nag_anova_random (g04bbc) and the confidence intervals are computed by 
	nag_anova_confid_interval using the Tukey-Kramer method.  There are 26 observations and 4 treatments.

void test_nag_anova_confid_interval()
{
	double bmean[1], tol = 5e-6;
	double irdf = 0;
	int nbolck = 1;
	double c[16];
	double cil[6];
	double ciu[6];
	double clevel = 0.95;
	double ef[4];
	double gmean;
	double r[26];
	double rdf;
	matrix table;
	table.SetSize(4,5);
	double tmean[4];
	int i, ij;
	int irep[4];
	int isig[6];
	int success;
	int j, n, nblock, nt;
	Nag_IntervalType type_enum;
	n = 26;
	nt = 4;
	type_enum = Nag_TukeyInterval;	
	double y[26] = {3, 2, 4, 3, 1, 5, 7, 8, 4, 10, 6, 3, 2, 1, 2, 4, 2, 3,
	1, 10, 12, 8, 5, 12, 10, 9};
    int it[26] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3 ,3, 3, 3, 3, 3, 3, 4, 
    4, 4, 4, 4, 4, 4};
    
    success =nag_anova_random(n, y, Nag_NoBlocks, nblock, nt, it, &gmean, bmean, tmean, 
								table,  c, nt, irep, r, ef, tol, irdf);   
    
    printf("the input data and results as folowing\n\n");
    printf("input data:\n");
    printf("n = %d,  nt = %d  clevel = %4.2f type_enum = Nag_TurkeyInterval\n",n, nt, clevel);
    printf("\ny:\n");
    for(i = 0; i < 26; i ++)
  	{
  		printf("%2.0f,",y[i]);
  		if( (i + 1) % 13 == 0)
  		{
  			printf("\n");
  		}
  	}
  	
    printf("\nit:\n");
    for(i = 0; i < 26; i ++)
  	{
  		printf("%d, ",it[i]);
  		if( (i + 1) % 13 == 0)
  		{
  			printf("\n");
  		}
  	}
  	
    if(success ==0)
    {
    	printf("the results:\n");
    	printf("\nANOVA TABLE");
    	printf("\n\n  Source       df    SS         MS       F      Prob\n");
    	printf("Treatments   ");
    	for(j = 0; j < 5; j++)
    	{
    		if( j == 0)
    		printf("%3.0f", table[1][j]);
    		else if(j == 1 || j == 2)
    		printf("%10.1f", table[1][j]);
    		else if(j == 3)
    		printf("%10.3f", table[1][j]);
    		else
    		printf("%9.4f", table[1][j]);   		
    	} 
    	
    	printf("\nResidual    ");
    	for(j = 0; j < 3; j++)
    	{  	if( j == 0)
    		printf(" %3.0f", table[2][j]);
    		if(j == 1 || j == 2)
    		printf("%10.1f", table[2][j]);    		   	
    	}		
		
		printf("\nTotal   ");		
		for(j = 0; j < 2; j++)
    	{  	if( j == 0)
    		printf("     %3.0f", table[3][j]);
    		if(j == 1 || j == 2)
    		printf("%10.1f", table[3][j]);    		   	
    	}

    	printf("\nTreat means  ");		
		for(i = 0; i < nt; i ++)
  		{
  			printf("%6.3f  ",tmean[i]);
  			if( (i + 1) % 8 == 0)
  			{
  				printf("\n");
  			}
  		}
		
		printf("\nSimiltanous Confidence Intervals\n");
		rdf  = table[2][0];
		success = nag_anova_confid_interval(type_enum, nt, tmean, rdf, c, nt, clevel, cil, ciu, isig);
		if(success ==0)
		{
			ij = 0;
			for(i = 1; i <= nt; i++)
			{
				for(j = 1; j <= i-1; j++)
				{
					++ij;
					printf("%2ld%2ld %10.3f  %10.3f", i, j, cil[ij - 1], ciu[ij - 1]);
					if(isig[ij - 1] == 1)
					printf("    *");
					if(isig[ij - 1] == 2)
					printf("\0");
					if(isig[ij - 1] == 0)
					printf("  ");
					printf("\n");
				}								
			}			
		}
		else
		printf("there is some problem with function nag_anova_random");	    
    }
    else
    printf("there is some problem with the function of nag_anova_random"); 
}
		
	The output is following:
the input data and results as folowing

input data:
n = 26,  nt = 4  clevel = 0.95 type_enum = Nag_TurkeyInterval

y:
 3, 2, 4, 3, 1, 5, 7, 8, 4,10, 6, 3, 2,
 1, 2, 4, 2, 3, 1,10,12, 8, 5,12,10, 9,

it:
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 
3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 
the results:

ANOVA TABLE

  Source       df    SS         MS       F      Prob
Treatments     3     239.9      80.0    24.029   0.0000
Residual      22      73.2       3.3
Total         25     313.1
Treat means   3.000   7.000   2.250   9.429  
Similtanous Confidence Intervals
 2 1      0.933       7.067    *
 3 1     -3.486       1.986  
 3 2     -7.638      -1.862    *
 4 1      3.610       9.247    *
 4 2     -0.538       5.395  
 4 3      4.557       9.800    *


Return:
	The function returns NAG error code, 0 if no error.
	
	NE_INT_ARG_LT(11): On entry, nt must not be less than 2: nt = <value>.
	NE_2_INT_ARG_LT(17): On entry, tdc = <value> while nt = <value>.
	NE_REAL_ARG_LT(5):  On entry, rdf must not be less than 1: rdf = <value>.
	NE_REAL(96): On entry, clevel <value>.  Constrant 0 < clevel < 1.
	NE_BAD_PARAM(70): On entry, parameter type has an illegal value.
	NE_ALLOC_FAIL(73): Memory allocation failed.
	NE_INTERNAL_ERROR(74): An internal error has occurred in this function.  Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
	NE_STUDENTIZED_STAT(549): There has been a failure in the computation of the studentized range statistic.  Try using a smaller value of clevel.
			
	successfully call of the nag_anova_confid_interval function.	
*/

	int nag_anova_confid_interval(
		Nag_IntervalType type,
		int nt,				// the number of treatment means
		const double tmean[], // the treatment means
		double rdf,			// the residual degrees of freedom
		const double c[],	// the standard errors of the differences between the means
		int ldc,			// the second dimension of the array c
		double clevel,		// the confidence level for the computed intervals
		double cil[],		// lower limit to the confidence interval
		double ciu[],		// upper limit to the confidence interval
		int isig[]			// difference between the ith and jth means in tmean is significant.
	);	// Computes confidence intervals for differences between means computed 

/**	g04eac
		Computes orthogonal polynomials or dummy variables for factor/classification variable.
		
Example:
	Data are read in from an experiment with four treatments and three observations per 
	treatment with the treatment coded as a factor, nag_dummy_vars is used to compute the 
	required dummy variables and the model is then fitted by nag_regsn_mult_linear (g02dac).

void test_nag_dummy_vars()
{
	double b[5], cov[15], p[35], h[12], q[72], rep[4], res[12]; 
	double df, rss, tol;
	double se[5], v[1], com_ar[45], wt[12], x[48];
	double *wtptr;
	int i, j, ip, irank, levels, m, n, tdq, tdx;
	int isx[4];
	Boolean  svd;
	Nag_DummyType dum_type;
	Nag_IncludeMean  mean_enum;
	int success, success1;
	int ifact[12] = {1, 4, 2, 3, 4, 2, 4, 1, 3, 1, 3, 2};	
	double y[12] = {33.63, 39.62, 38.18, 41.46, 38.02, 35.83, 35.99, 36.58, 42.92,
	37.80, 40.43, 37.89};
	
	n = 12;
	levels = 4;
	dum_type = Nag_AllLevels;
	mean_enum = Nag_MeanInclude;
	tdx = levels;
	wtptr = NULL;
		
	printf("the input data are as following\n");
	printf("n = 12,	levels = 4,	dum_type = Nag_AllLevels,\nmean_enum = Nag_MeanInclude");
	printf(",tdx = levels, wtptr = NULL, tol = 1e-5, ip = 5, m = 4, tdx =4 \n");
	printf("ifact:\n");
	for(i = 0; i < n; i ++)
	{
		printf("%5.2f,  ",ifact[i]);
		if((i + 1) % 6 == 0)
			printf("\n");
	}
	
	printf("y:\n");
	for(i = 0; i < n; i ++)
	{
		printf("%5.2f,  ",y[i]);
		if((i + 1) % 6 == 0)
			printf("\n");
	}
	
	success = nag_dummy_vars(dum_type, n, levels, ifact, x, tdx, v, rep);
	
	if(success == 0)
	{
		m = tdx;
		ip = m;
		if(mean_enum == Nag_MeanInclude)
			ip++;
		for(j = 0; j < m ; j++)
			isx[j] = 1;
		tol = 1e-5;
		tdq = 6;
		success1 = nag_regsn_mult_linear(mean_enum, n, x, tdx, m, isx, ip, y, wtptr, &rss, &df, b,
										se, cov, res, h, q, tdq, &svd, &irank, p, tol, com_ar);
	
		if(success1 ==0)
		{
			if(svd)
				printf("Model not full of rank, rank = %4ld\n\n ", irank);
			printf("Residual sum of squares = %12.4e\n", rss);
			printf("Degrees of freedom = %3.1f\n\n", df);
			printf("Variable Parameter estimate Standard error\n\n");
			for (j=0; j<ip; j++)
				printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]);
			printf("\n");			
		}
		else
			printf("what the hell  is the pronblem with nag_regsn_mult_linear");
	}	
	else
		printf("what the hell  is the pronblem with the fuction");

}
	
	The output is following:
	
the input data are as following
n = 12,	levels = 4,	dum_type = Nag_AllLevels,
mean_enum = Nag_MeanInclude,tdx = levels, wtptr = NULL, tol = 1e-5, ip = 5, m = 4, tdx =4 
ifact:
 0.00,   0.00,   0.00,   0.00,   0.00,   0.00,  
 0.00,   0.00,   0.00,   0.00,   0.00,   0.00,  
y:
33.63,  39.62,  38.18,  41.46,  38.02,  35.83,  
35.99,  36.58,  42.92,  37.80,  40.43,  37.89,  
Model not full of rank, rank =    4

 Residual sum of squares =  2.2227e+001
Degrees of freedom = 8.0

Variable Parameter estimate Standard error

     1         3.0557e+001         3.8494e-001
     2         5.4467e+000         8.3896e-001
     3         6.7433e+000         8.3896e-001
     4         1.1047e+001         8.3896e-001
     5         7.3200e+000         8.3896e-001

Return:
	The function returns NAG error code, 0 if no error.
	
	NE_INT_ARG_LT(11): On entry, levels must not be les than 2: levels = <value>.
	NE_2_INT_ARG_LT(17): On entry, n = <value> while levels = <value>.  These parameters must satisfy n >= levels.  On entry, tdx = <value> while levels = <value>.  These parameters must satisfy tdx >= levels.  On entry, tdx = <value> while levels - 1 = <value>.  These parameters must satisfy tdx >= levels - 1.
	NE_BAD_PARAM(70): On entry, parameter type has an illegal value.
	NE_ALLOC_FAIL(73): Memory allocation failed.
	NE_ARRAY_CONS(58): The contents of array v are not valid.  Constraint: all values of v must be distinct.
	NE_INT_ARRAY_CONS(103): On entry, factor{0] = <value>.  Constraint: 1 <= factor[0] <= levels.
	NE_G04EA_LEVELS(547): All levels are not represented in array factor.
	NE_G04EA_ORTHO_POLY(548): An orthogonal polynomial has all values zero.  This will be due to some values of v being close together.  This can only occur if type = Nag_Poly.
	NE_INTERNAL_ERROR(74): An internal error has occurred in this function.  Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
		
	successfully call of the nag_dummy_vars function.	
*/

int nag_dummy_vars(
	Nag_DummyType type,
	int n,	// the number of observations for which the dummy variables are to be computed.
	int levels,	// the number of levels of the factor.
	const int ifact[], // the n values of the factor.
	double x[],	// the n by k matrix of dummy variables
	int tdx,	// the second dimension of the array x
	const double v[], // If type equal Nag_Poly, the k distinct values of the underlying variable for which 
					// the orthogonal polynomial is to be computed.  If type not equal Nag_Poly, v is not reference.
	double rep[]	// contains the number of replications for each level
);	// Computes orthogonal polynomials or dummy variables for factor/classification variable.

/* end proto */

#endif //!_O_NAG_G04_H