#include<curses.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<gmp.h>

int n, i, N;
double x_t[20],x_e[20],val, sum;
char input[10], output[10];
double chi_2(double val, int N);

int main(){
    WINDOW *pop_ptr;

   initscr();
echo();
resize_term(35, 100);
    if(!has_colors()){
	endwin();
	exit(0);
    }
    if(start_color() != OK){
	endwin();
	exit(0);
    }
    init_pair(1, COLOR_RED, COLOR_BLACK);
    init_pair(2, COLOR_YELLOW, COLOR_BLACK);
    init_pair(3, COLOR_BLUE, COLOR_BLACK);
    init_pair(0, COLOR_WHITE, COLOR_BLACK);
    init_pair(6, COLOR_BLUE, COLOR_WHITE);
    init_pair(5, COLOR_CYAN, COLOR_BLACK);

 

    move(1, 10);
    attrset(COLOR_PAIR(6));
    printw("201 Chi^2 analysis");
    attrset(COLOR_PAIR(0));
    move(3,1);
    printw("Number of data points N=");
    attron(A_UNDERLINE);
    attrset(COLOR_PAIR(1));
    getstr(input);
    attrset(COLOR_PAIR(0));
    attroff(A_UNDERLINE);
    N=atoi(input);
    for(n=1;n<=N;n++){
f10:	move(5+n, 1);
         printw("x_{theor}[%d]= ",n);
	 attron(A_UNDERLINE);
	 attrset(COLOR_PAIR(1));
	 getstr(input);
	 attrset(COLOR_PAIR(0));
	 attroff(A_UNDERLINE);
	 x_t[n]=atof(input);
if(x_t[n]<0.0){
	     pop_ptr=newwin(4,15,10, 30);
	     box(pop_ptr,'|','-');
	     mvwprintw(pop_ptr,2,2,"x_{theor}[%d] must >0.0",n);
	     wrefresh(pop_ptr);  /* draw */
	     napms(1500);
	     wclear(pop_ptr);    /* clear */
	     wrefresh(pop_ptr);  /* redraw */
	     delwin(pop_ptr);    /* delete */
	     refresh();
	move(5+n,1);
	for(i=0;i<strlen(input);i++)
	delch();
	goto f10;
    }


	move(5+n, 30);
         printw("x_{exp}[%d]= ",n);
	 attron(A_UNDERLINE);
	 attrset(COLOR_PAIR(1));
	 getstr(input);
	 attrset(COLOR_PAIR(0));
	 attroff(A_UNDERLINE);
	 x_e[n]=atof(input);

	 


    }
    sum=0.0;
    for(n=1;n<=N;n++){
	sum=sum+(x_e[n]-x_t[n])*(x_e[n]-x_t[n])/(x_t[n]);
    }
	move(7+N, 1);
       
	printw("The value of ");
	attrset(COLOR_PAIR(5));
	printw("Chi^2 = %f", sum);
	move(9+N, 1);
	printw("P(Chi^2>=%f) = %f", sum, chi_2(sum, N));
	refresh();


move(29, 1);
attrset(COLOR_PAIR(0));
printw("Type anything <RET> to quit ");
attrset(COLOR_PAIR(2));
getstr(input);
   endwin();
   exit(0);

}



double chi_2(double val, int N){

int m,p, num;
mpf_t w[21];
mpf_t y[21];
mpf_t D[21];
mpf_t factors[21];
mpf_t powers[21];
mpf_t chi, sum, tmp, term, two;
    
   
    for(m=0;m<=20;m++){
	mpf_init2(y[m], 256);
	mpf_init2(w[m],256);
	mpf_init2(D[m],256);
	mpf_init2(powers[m],256);
	mpf_init2(factors[m],256);
    }
    mpf_init_set_ui(sum,0);
    mpf_init_set_ui(two,2);
    mpf_init2(chi, 256);
    mpf_init2(term, 256);
    mpf_init2(sum, 256);
    mpf_init2(tmp, 256); 
        
 
  /* the abscissa points and weights */
 mpf_set_str(y[1],"0.705398896919887533667e-1",10);
 mpf_set_str(w[1],"0.168746801851113862149223899677e0",10);
 mpf_set_str(y[2],"0.372126818001611443794e0",10);
 mpf_set_str(w[2],"0.291254362006068281716795323806e0",10);
 mpf_set_str(y[3],"0.916582102483273564668e0",10);
 mpf_set_str(w[3],"0.266686102867001288549520868998e0",10);
 mpf_set_str(y[4],"0.170730653102834388068e1",10);
 mpf_set_str(w[4],"0.166002453269506840032140094428e0",10);
 mpf_set_str(y[5],"0.27491992553094321296e1",10);
 mpf_set_str(w[5],"0.748260646687923705422201683691e-1",10);
 mpf_set_str(y[6],"0.404892531385088692234e1",10);
 mpf_set_str(w[6],"0.249644173092832210734541468302e-1",10);
 mpf_set_str(y[7],"0.561517497086161651408e1",10);
 mpf_set_str(w[7],"0.620255084457223684757314111118e-2",10);
 mpf_set_str(y[8],"0.745901745367106330976e1",10);
 mpf_set_str(w[8],"0.114496238647690824204551915748e-2",10);
 mpf_set_str(y[9],"0.959439286958109677244e1",10);
 mpf_set_str(w[9],"0.155741773027811974784485571297e-3",10);
 mpf_set_str(y[10],"0.120388025469643163096e2",10);
 mpf_set_str(w[10],"0.15401440865224915690105126352e-4",10);
 mpf_set_str(y[11],"0.148142934426307399785e2",10);
 mpf_set_str(w[11],"0.108648636651798235149891019091e-5",10);
 mpf_set_str(y[12],"0.179488955205193760174e2",10);
 mpf_set_str(w[12],"0.533012090955671475092792098967e-7",10);
 mpf_set_str(y[13],"0.214787882402850109757e2",10);
 mpf_set_str(w[13],"0.17579811790505820036686240782e-8",10);
 mpf_set_str(y[14],"0.254517027931869055035e2",10);
 mpf_set_str(w[14],"0.372550240251232087264898344503e-10",10);
 mpf_set_str(y[15],"0.299325546317006120067e2",10);
 mpf_set_str(w[15],"0.476752925157819052455714099041e-12",10);
 mpf_set_str(y[16],"0.350134342404790000062e2",10);
 mpf_set_str(w[16],"0.337284424336243841253294925521e-14",10);
 mpf_set_str(y[17],"0.40833057056728571062e2",10);
 mpf_set_str(w[17],"0.115501433950039883096855598093e-16",10);
 mpf_set_str(y[18],"0.476199940473465021399e2",10);
 mpf_set_str(w[18],"0.153952214058234355350101868104e-19",10);
 mpf_set_str(y[19],"0.558107957500638988907e2",10);
 mpf_set_str(w[19],"0.528644272556915782906500664212e-23",10);
 mpf_set_str(y[20],"0.665244165256157538186e2",10);
 mpf_set_str(w[20],"0.165645661249902329596426420137e-27",10);
 /* denominators as function of N */
 mpf_set_str(D[2],"0.1e1", 10);
 mpf_set_str(D[3],"0.125331413731550025120788264241e1",10);
 mpf_set_str(D[4],"0.2e1",10);
 mpf_set_str(D[5],"0.375994241194650075362364792722e1",10); 
 mpf_set_str(D[6],"0.8e1",10);
 mpf_set_str(D[7],"0.187997120597325037681182396361e2",10); 
 mpf_set_str(D[8],"0.48e2",10);
 mpf_set_str(D[9],"0.131597984418127526376827677453e3",10);
 mpf_set_str(D[10],"0.384e3",10);
 mpf_set_str(D[11],"0.118438185976314773739144909707e4",10);

 
 mpf_set_d(chi, val);

 for(m=1;m<=20;m++){
     mpf_mul(factors[m], y[m], two);
     mpf_add(factors[m], factors[m],chi);
     mpf_sqrt(factors[m], factors[m]);
     mpf_pow_ui(powers[m], factors[m], N-2);
    
 }


 /* denominator=pow(2.0, N/2.0-1.0)*exp(gamma(N/2.0)); */
 /* compute exp(-chi/2) */
 /* need to use at least 1000000 to get 6 decimal place accuracy */
 mpf_div_ui(tmp, chi, 1000000);
 mpf_sub(tmp, two, tmp);
 mpf_div(tmp, tmp, two);
 mpf_pow_ui(tmp, tmp, 1000000);
  /* a=chi^2 */

for(m=1;m<=20;m++){
    mpf_mul(term, w[m], powers[m]);
    mpf_add(sum,sum,term);
}
 mpf_div(sum, sum, D[N]); 
 /* mpf_set_d(tmp, exp(-val/2.0)); */
 
 
 
 mpf_mul(sum, sum, tmp);
 /* mpf_out_str(stdout,10, 30, sum); */
 /* gmp_printf("\nProb(chi^2 >= %f)=%.*Ff\n", val, 5, sum); */
 

 for(m=0;m<21;m++){
	mpf_clear(y[m]);
	mpf_clear(w[m]);
	mpf_clear(D[m]);
	mpf_clear(powers[m]);
	mpf_clear(factors[m]);
	mpf_clear(two);
	mpf_clear(sum);
	mpf_clear(tmp);
	mpf_clear(term);
	mpf_clear(chi);

 }
 return(mpf_get_d(sum));
}