#include <math.h>
#include <iostream.h>
#include <fstream>
#include "contdist.h"

// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// global variable:
//   nsize = sample size
//   nq = rows of qv, that is, the number of q's to be computed
//   np = rows of pv, that is, the number of probability values 
//        to be computed for critical values
// other variables:
//   pv = vector of probability values, on which critical values
//        should be computed
//   qv = vector of q's. Pls note that nsize/q should be an integer
// Output:
//   res.txt = vector of values index be res_{ij} (i=1,2,...,nq)
//             (j=1,2,...,np), which can be reshaped to a nq by np
//             matrix.
// How To Use:
//   1. specify nsize, nq, and np;
//   2. choose probability values and write down for pv[]; and
//   3. choose values for qv[].
// Users are recommended to check the ooutput with the built-in function
// betap(x, a, b).
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#define nsize 30
#define nq 3
#define np 4


int main(void)
{
    ofstream output("res.txt"); 

    int iq, ip;
    int qv[nq];
    double pv[np]; 
    double q, n, mean, v, s, varia, xa, xb, p, zp;

	pv[0] = 0.025;
	pv[1] = 0.050;
	pv[2] = 0.950;
	pv[3] = 0.975;

	qv[0] = 2;
	qv[1] = 3;
	qv[2] = 5;
	//qv[3] = 5;
	//qv[4] = 6;
	//qv[5] = 10;

	for(iq=0; iq<nq; iq++)
	{
		q = qv[iq]; 
		n = (double) nsize/q;
		mean = (n-1)/(q*n-1);
		v = 2*(2*q-1.0)*(q-1.0)/(3*q*q*n);
        s = (n - 1.0)/(n*q - 1.0);
        varia = v*s*s;
		xa = (mean*mean - mean*mean*mean - mean*varia)/varia;
		xb = xa*(1-mean)/mean;

		for(ip=0; ip<np; ip++)
		{
			p = pv[ip];
			zp = betav(p,xa,xb);
			output << zp << "\n";
		}
	}
	return 0;
}