#include <Stdio.h>
#include<stdlib.h>
#include <math.h>
#include <gsl/gsl_matrix.h>

# define N 4  /* array size */

int main()
{			  
		double tol=0.000000001;
		int sum;
		int iteration=0;
		double err[N];
		int max_iteration=100;
		double theta0[N]={0.6,0.3,0.05,0.05};
		double old_theta[N];
		double new_theta[N];
		double answer[N]={0,0,0,0};
		int i,j;
		int count1, count2,count3,count4,count5, count6,count7,count8,count9;

		/*define storage for input file */
		
		gsl_matrix *s =gsl_matrix_alloc(1842240,10);
		
		/* read the file in*/
  		FILE *fp;
  		fp =fopen("ant my.tab.dat","r");
  		gsl_matrix_fscanf (fp,s);
  		fclose(fp);

	 /* result file */
		FILE *f;
	  f =fopen("integer ant em.dat","w");
		

	 for(i=0;i<1842240;i++) {
		 if((i %100000)==0) printf("Pair finished %d\n",i);
	 
	 	count1=gsl_matrix_get(s,i,0);
		count2=gsl_matrix_get(s,i,1);
		count3=gsl_matrix_get(s,i,2);
		count4=gsl_matrix_get(s,i,3);
		count5=gsl_matrix_get(s,i,4);
		count6=gsl_matrix_get(s,i,5);
		count7=gsl_matrix_get(s,i,6);
		count8=gsl_matrix_get(s,i,7);
		count9=gsl_matrix_get(s,i,8);
		sum=gsl_matrix_get(s,i,9);
					
		iteration=0;					
		 err[0] = tol+1;
		 err[1] = tol+1;
		 err[2] = tol+1;
		 err[3] = tol+1;
	
		old_theta[0]=theta0[0]=0.6; 
		old_theta[1]=theta0[1]=0.3;
		old_theta[2]=theta0[2]=0.05;
		old_theta[3]=theta0[3]=0.05;
	 	
	 	new_theta[0]=theta0[0]=0.6;
 	   new_theta[1]=theta0[1]=0.3;
		new_theta[2]=theta0[2]=0.05;
		new_theta[3]=theta0[3]=0.05;
	
		
		while((err[0] > tol) && (err[1] > tol) && (err[2] > tol) && (err[3] > tol))

		{
				iteration = iteration +1;
				
				old_theta[0] = new_theta[0];
				old_theta[1] = new_theta[1];
				old_theta[2] = new_theta[2];
				old_theta[3] = new_theta[3];
		  
		
		new_theta[0] =((2*count1)+count2+count4+(count5*(old_theta[0]*old_theta[3])/((old_theta[1]*old_theta[2])+(old_theta[0]*old_theta[3]))))/(2*sum);
		
		new_theta[1] =(count2+(2*count3)+count6+(count5*(old_theta[1]*old_theta[2])/((old_theta[1]*old_theta[2])+(old_theta[0]*old_theta[3]))))/(2*sum);
		
		new_theta[2] =(count4+(2*count7)+count8+(count5*(old_theta[1]*old_theta[2])/((old_theta[1]*old_theta[2])+(old_theta[0]*old_theta[3]))))/(2*sum);
		
		new_theta[3] =((count5)*((old_theta[0]*old_theta[3])/((old_theta[1]*old_theta[2])+(old_theta[0]*old_theta[3])))+count8+(2*count9)+count6)/(2*sum);
			
				err[0] = fabs(old_theta[0]-new_theta[0]);
				err[1] = fabs(old_theta[1]-new_theta[1]);
				err[2] = fabs(old_theta[2]-new_theta[2]);
				err[3] = fabs(old_theta[3]-new_theta[3]);
		
		
				if (iteration >max_iteration) break;
		 
		}

 	 answer[0]= new_theta[0]*sum;
	 answer[1]= new_theta[1]*sum;
	 answer[2]= new_theta[2]*sum;
	 answer[3]= sum-answer[0]-answer[1]-answer[2];
  
	fprintf(f,"%.0f %.0f %.0f %.0f\n ", answer[0], answer[1], answer[2], answer[3]);
	
	}
	}