data step;
   array apow(k) finp1-finp3;
   array epow(k) finp4-finp6;
   array an(k) finn1-finn3;
   array en(k) finn4-finn6;
   za=1.96;
   zb=0.842;
   desired=0.80;
   nsimul=100000;
do mean=0.50,1.0,3.0;
  do rho=0.10,0.50,0.80;
    beta=rho/(1-rho);
    alpha=mean/beta;
    do eta=log(0.50),log(1.2),log(1.5);
     rmt=exp(eta);
     do k=1 to 3;
     if k eq 1 then gamma=log(0.9);
     if k eq 2 then gamma=log(1.0);
     if k eq 3 then gamma=log(1.1);
     rmp=exp(gamma);
     a=1+exp(eta+gamma);
     b=exp(eta)+exp(gamma);
     p1=exp(eta+gamma)/(1+exp(eta+gamma));
     p2=exp(gamma)/(exp(gamma)+exp(eta));
     pp=(exp(eta+gamma)+exp(gamma))/(1+exp(eta+gamma)+exp(eta)+exp(gamma));
      v0=(1/(mean*a)+1/(mean*b))/(4*pp*(1-pp));
     va=(1/(mean*a*p1*(1-p1))+1/(mean*b*p2*(1-p2)))/4;
     n=ceil((za*sqrt(v0)+zb*sqrt(va))**2/eta**2);
    sumrej=0;
    do sim=1 to nsimul;
     mu=beta*rangam(1432171,n*alpha);
     ydot1=ranpoi(1457151,mu*(1+exp(eta+gamma)));
     mu=beta*rangam(1432151,n*alpha);
     ydot2=ranpoi(1457151,mu*(exp(eta)+exp(gamma)));
     if ydot1 gt 0 and ydot2 gt 0 then do;
     lowtail=0;
      uptail=0;
     y12=ranbin(1432151,ydot1,p1);
     y22=ranbin(1432141,ydot2,p2);
      lrange=max(0,(y12+y22)-ydot2);
      urange=min(ydot1,(y12+y22));
      do y=lrange to y12;
       lowtail+exp(lgamma(ydot1+1)+lgamma(ydot2+1)+lgamma(y12+y22+1)+
               lgamma(ydot1-y12+ydot2-y22+1)-lgamma(y+1)-lgamma(ydot1-y+1)-
               lgamma(y12+y22-y+1)-lgamma(ydot2-(y12+y22)+y+1)-lgamma(ydot1+ydot2+1));
      end;  ** end of do lrange to y12;
      do y=y12 to urange;
       uptail+exp(lgamma(ydot1+1)+lgamma(ydot2+1)+lgamma(y12+y22+1)+
               lgamma(ydot1-y12+ydot2-y22+1)-lgamma(y+1)-lgamma(ydot1-y+1)-
               lgamma(y12+y22-y+1)-lgamma(ydot2-(y12+y22)+y+1)-lgamma(ydot1+ydot2+1));
      end;  ** end of do lrange to y12;
      if min(lowtail,uptail) le 0.025 then sumrej+1;
      end;  ** end of ydot1 gt 0 and ydot2 gt 0;
      end; ** end of do sim=1 to nsimul;
     power=sumrej/nsimul;
     apow=power; an=n;
     if power lt desired then do;
     do while (power lt desired);
      n+10;
    sumrej=0;
    do sim=1 to nsimul;
     mu=beta*rangam(1432151,n*alpha);
     ydot1=ranpoi(1457151,mu*(1+exp(eta+gamma)));
     mu=beta*rangam(1432151,n*alpha);
     ydot2=ranpoi(1457151,mu*(exp(eta)+exp(gamma)));
     if ydot1 gt 0 and ydot2 gt 0 then do;
     lowtail=0;
      uptail=0;
     y12=ranbin(1432151,ydot1,p1);
     y22=ranbin(1432141,ydot2,p2);
      lrange=max(0,(y12+y22)-ydot2);
      urange=min(ydot1,(y12+y22));
      do y=lrange to y12;
       lowtail+exp(lgamma(ydot1+1)+lgamma(ydot2+1)+lgamma(y12+y22+1)+
               lgamma(ydot1-y12+ydot2-y22+1)-lgamma(y+1)-lgamma(ydot1-y+1)-
               lgamma(y12+y22-y+1)-lgamma(ydot2-(y12+y22)+y+1)-lgamma(ydot1+ydot2+1));
      end;  ** end of do lrange to y12;
      do y=y12 to urange;
       uptail+exp(lgamma(ydot1+1)+lgamma(ydot2+1)+lgamma(y12+y22+1)+
               lgamma(ydot1-y12+ydot2-y22+1)-lgamma(y+1)-lgamma(ydot1-y+1)-
               lgamma(y12+y22-y+1)-lgamma(ydot2-(y12+y22)+y+1)-lgamma(ydot1+ydot2+1));
      end;  ** end of do lrange to y12;
      if min(lowtail,uptail) le 0.025 then sumrej+1;
      end;  ** end of ydot1 gt 0 and ydot2 gt 0;
      end;  ** end of do sim=1 to nsimul;
     power=sumrej/nsimul;
     end;  ** end of do while power lt desired;
     nlow=n-10;
     nup=n+10; ** do this is to assure that the final value falls between these two limits;
               ** otherwise the final sample estimate may be the previous result;
     do nsam=nlow to nup;
      n=nsam;
    sumrej=0;
    do sim=1 to nsimul;
     mu=beta*rangam(1432151,n*alpha);
     ydot1=ranpoi(1457151,mu*(1+exp(eta+gamma)));
     mu=beta*rangam(1432151,n*alpha);
     ydot2=ranpoi(1457151,mu*(exp(eta)+exp(gamma)));
     if ydot1 gt 0 and ydot2 gt 0 then do;
     lowtail=0;
      uptail=0;
     y12=ranbin(1432151,ydot1,p1);
     y22=ranbin(1432141,ydot2,p2);
      lrange=max(0,(y12+y22)-ydot2);
      urange=min(ydot1,(y12+y22));
      do y=lrange to y12;
       lowtail+exp(lgamma(ydot1+1)+lgamma(ydot2+1)+lgamma(y12+y22+1)+
               lgamma(ydot1-y12+ydot2-y22+1)-lgamma(y+1)-lgamma(ydot1-y+1)-
               lgamma(y12+y22-y+1)-lgamma(ydot2-(y12+y22)+y+1)-lgamma(ydot1+ydot2+1));
      end;  ** end of do lrange to y12;;
      do y=y12 to urange;
       uptail+exp(lgamma(ydot1+1)+lgamma(ydot2+1)+lgamma(y12+y22+1)+
               lgamma(ydot1-y12+ydot2-y22+1)-lgamma(y+1)-lgamma(ydot1-y+1)-
               lgamma(y12+y22-y+1)-lgamma(ydot2-(y12+y22)+y+1)-lgamma(ydot1+ydot2+1));
      end;  ** end of do lrange to y12;
      if min(lowtail,uptail) le 0.025 then sumrej+1;
      end;  ** end of ydot1 gt 0 and ydot2 gt 0;
      end;  ** end of do sim=1 to nsimul;
     power=sumrej/nsimul;
     if power ge desired then do;  fpower=power;
                                   nfinal=n;
                                   nsam+20;
                              end;
     end;  ** end of do nsam=nlow to nup;
     end;  ** end of do if power lt desired;
   epow=fpower; en=nfinal;
  if k eq 3 then do;
put mean 1-3 1 rho 4-8 2 rmt 9-13 1 finn1 16-20 finn4 26-30 finn2 36-40 finn5 46-50 
                                    finn3 56-60 finn6 66-70/
                              @15 "(" finp1 16-20 3 @21 ")" @25 "(" finp4 26-30 3 @31 ")"
                              @35 "(" finp2 36-40 3 @41 ")" @45 "(" finp5 46-50 3 @51 ")"
                              @55 "(" finp3 56-60 3 @61 ")" @65 "(" finp6 66-70 3 @71 ")";

  end;  ** end of if k eq 3;
  end;  ** end of do gamma;
  end;  ** end of do eta;
  end;  ** end of do rho;
  end;  ** end of do mean;