#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<mpfr.h>


#define THRESHOLD 0.0000001

#define Np0   1
#define Dp0   262144

#define Np_r  1
#define Dp_r  2048

#define Np_b  1
#define Dp_b  1

#define Np_rb 3
#define Dp_rb 131072


#define DEG 3

#define GMP_PRECISION 657400



int binom[DEG+1][DEG+1];
int k;

mpfr_t p;
mpfr_t p_ERR;

mpfr_t r,b,w;
mpfr_t r_ERR,b_ERR,w_ERR;

mpfr_t qr,qb,qw;
mpfr_t qr_ERR,qb_ERR,qw_ERR;

mpfr_t P_r[DEG+1][DEG+1];
mpfr_t P_r_ERR[DEG+1][DEG+1];
mpfr_t P_b[DEG+1][DEG+1];
mpfr_t P_b_ERR[DEG+1][DEG+1];
mpfr_t P_w[DEG+1][DEG+1];
mpfr_t P_w_ERR[DEG+1][DEG+1];

mpfr_t w_rb[DEG+1][DEG+1];
mpfr_t w_rb_ERR[DEG+1][DEG+1];



void add(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR, mpfr_t op2, mpfr_t op2_ERR)
{
  mpfr_add(r_ERR, op1,op2, MPFR_RNDU);
  mpfr_add(r_ERR, r_ERR, op1_ERR, MPFR_RNDU);
  mpfr_add(r_ERR, r_ERR, op2_ERR, MPFR_RNDU);

  mpfr_add(r,op1,op2,MPFR_RNDN);
}

void add_(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR)
{
  mpfr_t x_ERR;

  mpfr_init(x_ERR);

  mpfr_add(x_ERR,r,op1,MPFR_RNDU);
  mpfr_add(x_ERR,x_ERR,r_ERR,MPFR_RNDU);
  mpfr_add(x_ERR,x_ERR,op1_ERR,MPFR_RNDU);

  mpfr_set(r_ERR, x_ERR, MPFR_RNDU);
  mpfr_clear(x_ERR);

  mpfr_add(r,r,op1,MPFR_RNDN);
}

void add_ui_(mpfr_t r, mpfr_t r_ERR, unsigned long int op1)
{
  mpfr_t x_ERR;

  mpfr_init(x_ERR);

  mpfr_add_ui(x_ERR,r,op1,MPFR_RNDU);
  mpfr_add(x_ERR,x_ERR,r_ERR,MPFR_RNDU);

  mpfr_set(r_ERR, x_ERR, MPFR_RNDU);
  mpfr_clear(x_ERR);

  mpfr_add_ui(r,r,op1,MPFR_RNDN);
}

void sub(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR, mpfr_t op2, mpfr_t op2_ERR)
{
  mpfr_sub(r_ERR, op1,op2, MPFR_RNDU);
  mpfr_add(r_ERR, r_ERR, op1_ERR, MPFR_RNDU);
  mpfr_add(r_ERR, r_ERR, op2_ERR, MPFR_RNDU);

  mpfr_sub(r,op1,op2,MPFR_RNDN);
}

void sub_(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR)
{
  mpfr_t x_ERR;

  mpfr_init(x_ERR);

  mpfr_sub(x_ERR,r,op1,MPFR_RNDU);
  mpfr_add(x_ERR,x_ERR,r_ERR,MPFR_RNDU);
  mpfr_add(x_ERR,x_ERR,op1_ERR,MPFR_RNDU);

  mpfr_set(r_ERR, x_ERR, MPFR_RNDU);
  mpfr_clear(x_ERR);

  mpfr_sub(r,r,op1,MPFR_RNDN);
}

void mul(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR, mpfr_t op2, mpfr_t op2_ERR)
{
  mpfr_t x;

  mpfr_init(x);

  mpfr_mul(r_ERR, op1,op2, MPFR_RNDU);

  mpfr_mul(x,op1,op2_ERR,MPFR_RNDU);
  mpfr_add(r_ERR, r_ERR, x, MPFR_RNDU);
  mpfr_mul(x,op1_ERR,op2,MPFR_RNDU);
  mpfr_add(r_ERR, r_ERR, x, MPFR_RNDU);
  mpfr_clear(x);

  mpfr_mul(r,op1,op2,MPFR_RNDN);
}

void mul_(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR)
{
  mpfr_t x;
  mpfr_t x_ERR;

  mpfr_init(x);
  mpfr_init(x_ERR);

  mpfr_mul(x_ERR, r,op1, MPFR_RNDU);

  mpfr_mul(x,r,op1_ERR,MPFR_RNDU);
  mpfr_add(x_ERR, x_ERR, x, MPFR_RNDU);
  mpfr_mul(x,r_ERR,op1,MPFR_RNDU);
  mpfr_add(x_ERR, x_ERR, x, MPFR_RNDU);
  mpfr_clear(x);

  mpfr_set(r_ERR, x_ERR, MPFR_RNDU);
  mpfr_clear(x_ERR);

  mpfr_mul(r,r,op1,MPFR_RNDN);
}

void div_(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR)
{
  mpfr_t x,y;
  mpfr_t x_ERR;

  mpfr_init(x);
  mpfr_init(y);
  mpfr_init(x_ERR);

  mpfr_div(x_ERR, r,op1, MPFR_RNDU);

  mpfr_mul(x,r,op1_ERR,MPFR_RNDU);
  mpfr_mul(y,r_ERR,op1,MPFR_RNDU);
  mpfr_add(y,y,x,MPFR_RNDU);
  mpfr_mul(x,op1,op1,MPFR_RNDD);
  mpfr_div(y,y,x,MPFR_RNDU);

  mpfr_add(x_ERR,x_ERR,y,MPFR_RNDU);
  mpfr_clear(x);
  mpfr_clear(y);

  mpfr_set(r_ERR, x_ERR, MPFR_RNDU);
  mpfr_clear(x_ERR);

  mpfr_div(r,r,op1,MPFR_RNDN);
}

void div_ui(mpfr_t r, mpfr_t r_ERR, mpfr_t op1, mpfr_t op1_ERR, unsigned long int op2)
{
  mpfr_t x;

  mpfr_init(x);

  mpfr_div_ui(r_ERR, op1, op2, MPFR_RNDU);

  mpfr_mul_ui(x,op1_ERR,op2,MPFR_RNDU);
  mpfr_div_ui(x,x,op2*op2,MPFR_RNDU);

  mpfr_add(r_ERR,r_ERR,x,MPFR_RNDU);
  mpfr_clear(x);

  mpfr_div_ui(r,op1,op2,MPFR_RNDN);
}

void div_ui_(mpfr_t r, mpfr_t r_ERR, unsigned long int op1)
{
  mpfr_t x;
  mpfr_t x_ERR;

  mpfr_init(x);
  mpfr_init(x_ERR);

  mpfr_div_ui(x_ERR, r, op1, MPFR_RNDU);

  mpfr_mul_ui(x,r_ERR,op1,MPFR_RNDU);
  mpfr_div_ui(x,x,op1*op1,MPFR_RNDU);

  mpfr_add(x_ERR,x_ERR,x,MPFR_RNDU);
  mpfr_clear(x);

  mpfr_set(r_ERR, x_ERR, MPFR_RNDU);
  mpfr_clear(x_ERR);

  mpfr_div_ui(r,r,op1,MPFR_RNDN);
}

void status() {
  int i,j;

  printf("%d:",k);

  mpfr_printf(" p=%.7RNf", p);
  mpfr_printf(" r=%.7RNf", r);
  mpfr_printf(" b=%.7RNf", b);
  mpfr_printf(" w=%.7RNf", w);

  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) 
    mpfr_printf(" %d%d=%.7RNf", i,j,w_rb[i][j]);

  printf("\t");
  printf(" p_E=");
  mpfr_out_str(stdout, 10, 6, p_ERR, GMP_RNDU);
  printf(" r_E=");
  mpfr_out_str(stdout, 10, 6, r_ERR, GMP_RNDU);
  printf(" b_E=");
  mpfr_out_str(stdout, 10, 6, b_ERR, GMP_RNDU);
  printf(" w_E=");
  mpfr_out_str(stdout, 10, 6, w_ERR, GMP_RNDU);
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) {
    printf(" %d%d_E=",i,j);
    mpfr_out_str(stdout, 10,6, w_rb_ERR[i][j], GMP_RNDU);
  }
  printf("\n");
}

void init() {
int i,j;
int x;

  mpfr_set_default_prec(GMP_PRECISION);

  for (i=0; i<=DEG; i++) {
    binom[i][0]=1;
    x=i;
    for (j=1; j<=i; j++) {
      binom[i][j]=(x*binom[i][j-1])/j;
      x--;
    }
  }
  for (i=0; i<=DEG; i++) for (j=i+1; j<=DEG; j++) binom[i][j]=0;

  mpfr_init(p);
  mpfr_init(p_ERR);
  mpfr_init(r);
  mpfr_init(r_ERR);
  mpfr_init(b);
  mpfr_init(b_ERR);
  mpfr_init(w);
  mpfr_init(w_ERR);
  mpfr_init(qr);
  mpfr_init(qr_ERR);
  mpfr_init(qb);
  mpfr_init(qb_ERR);
  mpfr_init(qw);
  mpfr_init(qw_ERR);

  for (i=0; i<=DEG; i++) for (j=0; j<=DEG; j++) {
    mpfr_init(w_rb[i][j]);
    mpfr_init(w_rb_ERR[i][j]);
    mpfr_init(P_w[i][j]);
    mpfr_init(P_w_ERR[i][j]);
    mpfr_init(P_r[i][j]);
    mpfr_init(P_r_ERR[i][j]);
    mpfr_init(P_b[i][j]);
    mpfr_init(P_b_ERR[i][j]);
  }
}


void step() {
int i,j,i_,j_,l;
mpfr_t x,y,z,d;
mpfr_t x_ERR,y_ERR,z_ERR,d_ERR;
mpfr_t w_new[DEG+1][DEG+1];
mpfr_t w_new_ERR[DEG+1][DEG+1];

  k++;

  mpfr_init(x);
  mpfr_init(x_ERR);
  mpfr_init(y);
  mpfr_init(y_ERR);
  mpfr_init(z);
  mpfr_init(z_ERR);
  mpfr_init(d);
  mpfr_init(d_ERR);
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) {
    mpfr_init(w_new[i][j]);
    mpfr_init(w_new_ERR[i][j]);
  }


  //d=0
  mpfr_set_ui(d,0,MPFR_RNDN); mpfr_set_ui(d_ERR,0,MPFR_RNDU);

  for (i=0; i<DEG; i++) for (j=0; j<DEG-i; j++) {
    //x=DEG-i-j
    mpfr_set_ui(x, DEG-i-j,MPFR_RNDN); mpfr_set_ui(x_ERR, 0,MPFR_RNDU);

    //x*=w_rb[i][j]
    mul_(x, x_ERR, w_rb[i][j],w_rb_ERR[i][j]);

    //d+=x
    add_(d, d_ERR,x,x_ERR);
  }
  //d/=DEG
  div_ui_(d,d_ERR,DEG);

  mpfr_set_ui(qr,0,MPFR_RNDN); mpfr_set_ui(qr_ERR,0,MPFR_RNDU);
  mpfr_set_ui(qb,0,MPFR_RNDN); mpfr_set_ui(qb_ERR,0,MPFR_RNDU);
  mpfr_set_ui(qw,0,MPFR_RNDN); mpfr_set_ui(qw_ERR,0,MPFR_RNDU);

  for (i=0; i<DEG; i++) for (j=0; j<DEG-i; j++) {
    //x=DEG-i-j
    mpfr_set_ui(x, DEG-i-j,MPFR_RNDN); mpfr_set_ui(x_ERR, 0,MPFR_RNDU);
    
    //x*=w_rb[i][j]/DEG
    div_ui_(x,x_ERR,DEG);
    mul_(x, x_ERR, w_rb[i][j],w_rb_ERR[i][j]);

    //y=x*P_w[i][j]
    mul(y, y_ERR, x, x_ERR, P_w[i][j],P_w_ERR[i][j]);
    //qw+=y
    add_(qw,qw_ERR,y,y_ERR);

    //y=x*P_r[i][j]
    mul(y, y_ERR, x, x_ERR, P_r[i][j],P_r_ERR[i][j]);
    //qr+=y
    add_(qr,qr_ERR,y,y_ERR);

    //y=x*P_b[i][j]
    mul(y, y_ERR, x, x_ERR, P_b[i][j], P_b_ERR[i][j]);
    //qb+=y
    add_(qb,qb_ERR,y,y_ERR);
  }

  //qr/=d qb/=d qw/=d
  div_(qw, qw_ERR, d,d_ERR);
  div_(qr, qr_ERR, d,d_ERR);
  div_(qb, qb_ERR, d,d_ERR);


  //compute probability p_k
  mpfr_set_ui(x,0,MPFR_RNDN); mpfr_set_ui(x_ERR, 0,MPFR_RNDU);
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) {
    //y = DEG-i-j
    mpfr_set_ui(y, (DEG-i-j),MPFR_RNDN); mpfr_set_ui(y_ERR, 0,MPFR_RNDU);
    //y*=qb
    mul_(y,y_ERR,qb,qb_ERR);
    mpfr_div_ui(y,y,2,MPFR_RNDN);
    //y+=j
    add_ui_(y,y_ERR,j);
    //y*=P_r[i][j]
    mul_(y,y_ERR,P_r[i][j],P_r_ERR[i][j]);
    //y*=w_rb[i][j]
    mul_(y,y_ERR,w_rb[i][j],w_rb_ERR[i][j]);
    //x+=y
    add_(x,x_ERR,y,y_ERR);

    mpfr_set_ui(y, (DEG-i-j),MPFR_RNDN); mpfr_set_ui(y_ERR, 0,MPFR_RNDU);
    mul_(y,y_ERR,qr,qr_ERR);
    mpfr_div_ui(y,y,2,MPFR_RNDN);
    //y+=i
    add_ui_(y,y_ERR,i);
    //y*=P_b[i][j]
    mul_(y,y_ERR,P_b[i][j],P_b_ERR[i][j]);
    //y*=w_rb[i][j]
    mul_(y,y_ERR,w_rb[i][j],w_rb_ERR[i][j]);
    //x+=y
    add_(x,x_ERR,y,y_ERR);
  }
  //y=2*w*x/DEG
  mpfr_set_ui(y,2,MPFR_RNDN); mpfr_set_ui(y_ERR, 0,MPFR_RNDU);
  mul_(y,y_ERR,w,w_ERR);
  mul_(y,y_ERR,x,x_ERR);
  div_ui_(y,y_ERR,DEG);
  
  //p+=y
  add_(p,p_ERR,y,y_ERR);

  
  //compute probability r_k
  mpfr_set_ui(x,0,MPFR_RNDN); mpfr_set_ui(x_ERR, 0, MPFR_RNDU);
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) {
    //y=w_rb[i][j]*P_r[i][j]
    mul(y,y_ERR,w_rb[i][j],w_rb_ERR[i][j],P_r[i][j],P_r_ERR[i][j]);
    //x+=y
    add_(x,x_ERR,y,y_ERR);
  }
  //y=x*w
  mul(y,y_ERR,x,x_ERR,w,w_ERR);
  //r+=y
  add_(r,r_ERR,y,y_ERR);



  //compute probability b_k
  mpfr_set_ui(x,0,MPFR_RNDN); mpfr_set_ui(x_ERR, 0, MPFR_RNDU);
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) {
    //z=w_rb[i][j]*P_b[i][j]
    mul(z,z_ERR,w_rb[i][j],w_rb_ERR[i][j],P_b[i][j],P_b_ERR[i][j]);
    //x+=z
    add_(x,x_ERR,z,z_ERR);
  }
  //z=x*w
  mul(z,z_ERR,x,x_ERR,w,w_ERR);
  //r+=z
  add_(b,b_ERR,z,z_ERR);


  //compute probability p_k
  //w-=y+z
  sub_(w,w_ERR,y,y_ERR);
  sub_(w,w_ERR,z,z_ERR);
  

  mpfr_set_ui(d,0,MPFR_RNDN);  mpfr_set_ui(d_ERR, 0, MPFR_RNDU);
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i;j++) {
    //y=w_rb[i][j]*P_w[i][j]
    mul(y,y_ERR,w_rb[i][j],w_rb_ERR[i][j],P_w[i][j],P_w_ERR[i][j]);
    //d+=y
    add_(d,d_ERR,y,y_ERR);
  }
  
  for (i=0;i<=DEG;i++) for (j=0;j<=DEG-i;j++) {
    mpfr_set_ui(w_new[i][j],0,MPFR_RNDN); mpfr_set_ui(w_new_ERR[i][j],0,MPFR_RNDU);
    for (i_=0; i_<=i; i_++) for (j_=0; j_<=j; j_++) {
      //x=(DEG-i_-j_ \choose i-i_)*(DEG-i-j_ \choose j-j_)
      mpfr_set_ui(x,binom[DEG-i_-j_][i-i_]*binom[DEG-i-j_][j-j_],MPFR_RNDN); mpfr_set_ui(x_ERR,0,MPFR_RNDU);
      //x*=w_rb[i_][j_]*P_w[i_][j_]
      mul_(x,x_ERR,w_rb[i_][j_],w_rb_ERR[i_][j_]);
      mul_(x,x_ERR,P_w[i_][j_],P_w_ERR[i_][j_]);

      for (l=0; l<i-i_; l++) mul_(x,x_ERR,qr,qr_ERR);  //x*=qr
      for (l=0; l<j-j_; l++) mul_(x,x_ERR,qb,qb_ERR);  //x*=qb
      for (l=0; l<DEG-i-j; l++) mul_(x,x_ERR,qw,qw_ERR); //x*=qw
  
      //w_new[i][j]+=x
      add_(w_new[i][j],w_new_ERR[i][j],x,x_ERR);
    }
    //w_new[i][j]/=d
    div_(w_new[i][j],w_new_ERR[i][j],d,d_ERR);
  }

  for (i=0;i<=DEG;i++) for (j=0;j<=DEG-i;j++) {
    mpfr_set(w_rb[i][j],w_new[i][j],MPFR_RNDN); mpfr_set(w_rb_ERR[i][j], w_new_ERR[i][j],MPFR_RNDU);
  }
  
  mpfr_clear(d);
  mpfr_clear(d_ERR);
  mpfr_clear(x);
  mpfr_clear(x_ERR);
  mpfr_clear(y);
  mpfr_clear(y_ERR);
  mpfr_clear(z);
  mpfr_clear(z_ERR);
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) {
    mpfr_clear(w_new[i][j]);
    mpfr_clear(w_new_ERR[i][j]);
  }
}




int main() {
  int i,j,l;
  mpfr_t x;

  init();


//the first step of the procedure
  k = 1;

  //r=p_0/2
  mpfr_set_ui(r,Np0,MPFR_RNDN);
  mpfr_div_ui(r,r,2*Dp0,MPFR_RNDN);
  mpfr_set_ui(r_ERR,0,MPFR_RNDU);

  //b=r
  mpfr_set(b,r,MPFR_RNDN);
  mpfr_set_ui(b_ERR,0,MPFR_RNDU);

  //w=1-r-b
  mpfr_set_ui(w,1,MPFR_RNDN);
  mpfr_sub(w,w,b,MPFR_RNDN);
  mpfr_sub(w,w,r,MPFR_RNDN);
  mpfr_set_ui(w_ERR,0,MPFR_RNDU);

  //p=2*r*b
  mpfr_set_ui(p,2,MPFR_RNDN); mpfr_set_ui(p_ERR,0,MPFR_RNDU);
  mul_(p,p_ERR,r,r_ERR);
  mul_(p,p_ERR,b,b_ERR);

  for (i=0; i<=DEG; i++) for (j=0; j<=DEG; j++) {
    //w_rb[i][j] = (DEG \choose i)*(DEG-i \choose j)
    mpfr_set_ui(w_rb[i][j], binom[DEG][i]*binom[DEG-i][j],MPFR_RNDN);
    mpfr_set_ui(w_rb_ERR[i][j],0,MPFR_RNDU);

    for (l=0;l<i;l++)  mul_(w_rb[i][j],w_rb_ERR[i][j],r,r_ERR); //w_rb[i][j]*=r
    for (l=0;l<j;l++)  mul_(w_rb[i][j],w_rb_ERR[i][j],b,b_ERR); //w_rb[i][j]*=b
    for (l=0;l<DEG-i-j;l++) mul_(w_rb[i][j],w_rb_ERR[i][j],w,w_ERR); //w_rb[i][j]*=w
  }
    

//the first phase -- set up the probabilities
  for (i=0; i<=DEG; i++) for (j=0; j<=DEG-i; j++) {
    mpfr_set_ui(P_r[i][j], 0,MPFR_RNDN);
    mpfr_set_ui(P_b[i][j], 0,MPFR_RNDN);
    mpfr_set_ui(P_w[i][j], 0,MPFR_RNDN);
    if      (i+j==0) mpfr_set_ui(P_w[i][j], 1,MPFR_RNDN);
    else if (i>=2)   mpfr_set_ui(P_b[i][j], 1,MPFR_RNDN); 
    else if (j>=2)   mpfr_set_ui(P_r[i][j], 1,MPFR_RNDN);
    else if (i+j==2) mpfr_set_ui(P_w[i][j], 1,MPFR_RNDN);
  
    else if (i) {
      mpfr_set_ui(P_b[i][j], Np_b,MPFR_RNDN);
      mpfr_div_ui(P_b[i][j], P_b[i][j], Dp_b,MPFR_RNDN);
      mpfr_set_ui(P_w[i][j], 1,MPFR_RNDN);
      mpfr_sub(P_w[i][j],P_w[i][j],P_b[i][j],MPFR_RNDN);
    }
    else if (j) { 
      mpfr_set_ui(P_r[i][j], Np_r,MPFR_RNDN);
      mpfr_div_ui(P_r[i][j], P_r[i][j], Dp_r,MPFR_RNDN);
      mpfr_set_ui(P_w[i][j], 1,MPFR_RNDN);
      mpfr_sub(P_w[i][j],P_w[i][j],P_r[i][j],MPFR_RNDN); 
    }

    mpfr_set_ui(P_w_ERR[i][j],0,MPFR_RNDU);
    mpfr_set_ui(P_r_ERR[i][j],0,MPFR_RNDU);
    mpfr_set_ui(P_b_ERR[i][j],0,MPFR_RNDU);
  }


//run the first phase until the probability that v is white and has exactly one non-white neighbor is < 10^-7
  mpfr_init(x);
  do {
    status();
    step();
    mpfr_set(x,w_rb[1][0],MPFR_RNDN);
    mpfr_add(x,x,w_rb[0][1],MPFR_RNDN); 
    mpfr_mul(x,x,w,MPFR_RNDN);
 } while (mpfr_cmp_d(x,THRESHOLD) > 0);

  mpfr_clear(x);

  printf("\n THE END OF THE FIRST PHASE (round %d) \n\n",k);


//the second phase -- change the probabilities P_r[1][1] P_b[1][1] and P_w[1][1]
  mpfr_set_ui(P_r[1][1],Np_rb,MPFR_RNDN);
  mpfr_div_ui(P_r[1][1],P_r[1][1],Dp_rb*2,MPFR_RNDN);

  mpfr_set(P_b[1][1],P_r[1][1],MPFR_RNDN);

  mpfr_set_ui(P_w[1][1],1,MPFR_RNDN);
  mpfr_sub(P_w[1][1],P_w[1][1],P_r[1][1],MPFR_RNDN);
  mpfr_sub(P_w[1][1],P_w[1][1],P_b[1][1],MPFR_RNDN);


//run the second phase until the probability that v is white is < 10^-7
  while (mpfr_cmp_d(w,THRESHOLD) > 0) {
    status();
    step();
  }

  status();
  printf("\nProbabilities were: p0=%d/%d p_r=%d/%d p_b=%d/%d p_rb=%d/%d\n", Np0, Dp0, Np_r,Dp_r,Np_b,Dp_b,Np_rb,Dp_rb);

  return 0;
}