#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <math.h>
using namespace std;

void assay(double [6][8][4],double);
double normal(double,double);
int poisson(double);
double myrand(double,double);
void space(double [6][8][4],double,double);
double grow(double,int,double,double,int);
double getp95(double [6][8][4]);

main(int argc, char *argv[]) {

  int myage,growgens,adultgens,numlambda,nreps,i,j,ii,jj,kk,succ,itime,seed;
  double p,q,lambda,lowlambda,inclambda,tarfreq,errfreq,tarp95,errp95,extra,qq,freq,p95,mtot,dtemp,muts,x,tot,stot;
  double testis[6][8][4];

  int age, thresh,half;
  double diegn,freq2,p952;
  cin >> seed;
  srand48(seed);
  growgens = 30;

  while (cin >> age >> lambda >> p >> diegn) {

    for (half=0;half<2;half++) {

    adultgens = 23*(age-13); if (adultgens > 1081) { adultgens = 1081; }
    mtot = pow(2.,((double)(growgens)));
    qq = diegn;

      // FORMAT
      extra = 0.;
      for (ii=0;ii<6;ii++) {
	for (jj=0;jj<8;jj++) {
	  for (kk=0;kk<4;kk++) {
	    testis[ii][jj][kk] = 0.;
	  }
	}
      }
      
      // GROWTH PHASE
      for (ii=0;ii<growgens;ii++) {
	dtemp = pow(2,((double)(ii+1)));
	muts = myrand(dtemp,lambda);
	if (muts > 0.) {
	  dtemp = pow(2.,((double)(growgens-ii-1)));
	  for (jj=0;jj<((int)(floor(muts)));jj++) {
	    itime = adultgens;
	    x = grow(dtemp,itime,p,qq,adultgens);
	    extra += x; 
	    space(testis,x,mtot);
	  }
	}
      }

      // ADULT PHASE
      for (ii=0;ii<adultgens;ii++) {
	muts = myrand(mtot,lambda);
	if (muts > 0.) {
	  itime = adultgens - ii;
	  for (jj=0;jj<((int)(floor(muts)));jj++) {
	    x = grow(1.,itime,p,qq,adultgens);
	    extra += x; 
	    space(testis,x,mtot);
	  }
	}
      }

      // OUTPUT
      if (extra > (.05*mtot)) { cerr << "BIG EXTRA" << endl; }

      assay(testis,mtot);
    
	freq = 0;
	for (ii=0;ii<6;ii++) {
	  for (jj=0;jj<8;jj++) {
	    for (kk=0;kk<4;kk++) {
	      freq += testis[ii][jj][kk];
	    }
	  }
	}
	freq = freq/192.; 
	p95 = getp95(testis);

	if (half == 0) {
	  freq2 = freq;
	  p952 = p95;
	}
    }
    
    freq = 0.5*(freq + freq2);
    cout << age << " " << freq << " " << p95 << " " << p952 << endl;
  }

  return 1;
}

void assay(double testis[6][8][4],double mtot) {
  int i,j,k,flag,flag2,tot,pos,tube;
  double dil,f;

  for (i=0;i<6;i++) {
    for (j=0;j<8;j++) {
      for (k=0;k<4;k++) {
	tot = 10;
	dil = 1.;
	flag = 1;
	while (flag) {
	  pos = 0;
	  for (tube=0;tube<tot;tube++) {
	    if (myrand(25000.,(192.*testis[i][j][k]/(dil*mtot))) >= 1 ) { pos += 1; }
	  }
	  flag2 = 0;
	  if (tot == 10) { if (pos <= 5) { flag2 = 1; } }
	  if (tot == 40) { 
	    if (pos <= 25) { flag2 = 1; }
	    if (pos == 0) {
	      flag2 = 0;
	      dil = dil*0.05;
	    }
	  }
	  if (flag2) {
	    flag = 0;
	    f = -log(1. - (((double)(pos))/((double)(tot))))*40*dil;
	    if (f <= 0.) { f = 0.; } 
	    testis[i][j][k] = f;
	  }
	  else {
	    tot = 40;
	    dil = dil*10.;
	  }
	}
      }
    }
  }
}



double getp95(double testis[6][8][4]) {
  // OUTPUT THE PERCENTAGE OF TESTIS PIECES REQUIRED TO COMPRISE 95% OF THE MUTANTS
  // INPUT 3-DIMENSIONAL TESTIS ARRAY

  int i,j,k,cnt,ibig;
  double big,tot,lin[192];

  tot = 0; // LATE ADDITION
  cnt = -1;
  for (i=0;i<6;i++) {
    for (j=0;j<8;j++) {
      for (k=0;k<4;k++) {
	cnt += 1;
	lin[cnt] = testis[i][j][k];
	tot += testis[i][j][k];
      }
    }
  }
  for (i=0;i<191;i++) {
    big = lin[i];
    for (j=(i+1);j<192;j++) {
      if (lin[j] > big) {
	big = lin[j];
	ibig = j;
      }
    }
    if (big > lin[i]) {
      lin[ibig] = lin[i];
      lin[i] = big;
    }
  }

  // THIS IS A CHANGE FROM THE PUBLISHED CODE
  cnt = -1;
  big = 0.;
  while ((big/tot) < 0.95) {
    cnt += 1; if (cnt > 191) { cnt = 191; break; }
    big += lin[cnt];
  }

  return ((double)(cnt+1))/192;
}

double grow(double sta,int itime,double p,double qq,int ag) { 
// OUTPUT NUMBER OF MUTANTS IN THE CLUSTER
// INPUT INITIAL NUMBER MUTANTS (1 FOR ADULT PHASE, MORE FOR GROWTH PHASE),
// TIME REMAINING IN GENERATIONS, p AND qq PARAMETERS DISCUSSED ELSEWHERE

  int i,temp,loage,hiage;
  double x;
  loage = 43; hiage = 47;

  x = sta;
  for (i=0;i<itime;i++) {
    temp = i + (ag - itime);
    if ( (23*(loage-13)<=temp) && (temp<23*(hiage+1-13)) ) {
      x -= myrand(x,qq);
      if (x <= 0.) { 
	x = 0.;
	break;
      }
    }
    x += myrand(x,p);
  }

  return x;
}

void space(double testis[6][8][4],double mclump,double mtot) {
  // OUTPUT DISTRIBUTES MUTATION CLUSTER INTO 3-DIMENSIONAL TESTIS ARRAY
  // INPUT 3-DIMENSIONAL TESTIS ARRAY,  
  // NUMBER OF MUTANTS IN CLUMP, TOTAL NUMBER OF CELLS IN TESTIS

  double x,y,z,x1,x2,y1,y2,z1,z2,leng,dtemp;
  int i,j,k;

  x = 6.*drand48(); 
  y = 8.*drand48(); 
  z = 4.*drand48();
  leng = pow( 192.*mclump/mtot, 0.333);

  if (leng > 4) {
    dtemp = mclump/192.;
    for (i=0;i<6;i++) {
      for (j=0;j<8;j++) {
	for (k=0;k<4;k++) {
	  testis[i][j][k] += dtemp;
	}
      }
    }
  }
  else {
    if (x-(0.5*leng) < 0) { x += (0.5*leng)-x; }
    if (x+(0.5*leng) > 6) { x -= x+(0.5*leng)-6; }
    if (y-(0.5*leng) < 0) { y += (0.5*leng)-y; }
    if (y+(0.5*leng) > 8) { y -= y+(0.5*leng)-8; }
    if (z-(0.5*leng) < 0) { z += (0.5*leng)-z; }
    if (z+(0.5*leng) > 4) { z -= z+(0.5*leng)-4; }
    for (i=0;i<6;i++) {
      for (j=0;j<8;j++) {
	for (k=0;k<4;k++) {
	  x1 = max( ((double)(i)),  (x-(0.5*leng)) );
	  x2 = min( ((double)(i+1)),(x+(0.5*leng)) );
	  y1 = max( ((double)(j)),  (y-(0.5*leng)) );
	  y2 = min( ((double)(j+1)),(y+(0.5*leng)) );
	  z1 = max( ((double)(k)),  (z-(0.5*leng)) );
	  z2 = min( ((double)(k+1)),(z+(0.5*leng)) );
	  if ((x2 > x1) && (y2 > y1) && (z2 > z1)) {
	    dtemp = (x2-x1)*(y2-y1)*(z2-z1);
	    testis[i][j][k] += mclump*dtemp/pow(leng,3.);
	  }
	}
      }
    }
  }
}

double normal(double mean,double sd) {
  double y1,y2,x;
  int flag;
  flag = 1;
  while (flag) {
    y1 = -log(drand48());
    y2 = -log(drand48());
    if (y2 >= (.5*(y1-1)*(y1-1))) {
      if (drand48() < .5) { x = y1; }
      else { x = -y1; }
      flag = 0;
    }
  }
  x = mean + x*sd;
  return x;
}

int poisson(double x) {
  double temp;
  int flag, cnt;
  temp = 0;
  flag = 1;
  cnt = 0;
  while (flag) {
    cnt += 1;
    temp = temp + log(drand48());
    if (temp < (-x)) {
      flag = 0;
      cnt -= 1;
    }
  }
  return cnt;
}

double myrand(double num,double q) { 
// POISSON. IF NUM*Q HIGH THEN NORMAL APPROXIMATION TO SPEED UP THE ALGORITHM
  double dnum,temp;
  int i;

  if (num > 0.) {
    if (num < 100.) {
      temp = 0.;
      for (i=0;i<num;i++) {
	if (drand48() < q) { temp += 1.; }
      }
    }
    else {
      dnum = ((double)(num))*q;
      if (dnum < 50.) { temp = ((double)(poisson(dnum))); }
      else { 
	temp = (round(normal(dnum,(sqrt(dnum))))); 
	if (temp < 0) { temp = 0.; }
      }
    }
  }
  else { temp = 0.; }
  return temp;
}