/* EQUIL.G: procedure to compute equilibria of the prisoner's dilemman game

	    Fixed point (equilibrium condition) is given by

	    p_a(x_a,x_b)=br2_a(x_a,x_b,p_a(x_a,x_b))  where

            br2_a(x_a,x_b,p_a)=br_a(br_b(p_a,x_b),x_a)

            is the "second order" best response function, and
	    br_a and br_b are the "first order" best response
	    functions for prisoners a and b.

	    John Rust, University of Maryland, November 2005 */

    proc (2)=equil(xa,xb);

       local pt,pt1,df,i;

       pt1=1;
       pt=0;

       i=1;

       if (sa_steps > 0);

          if debg_e;

	     "Equil.g: Initial successive approximation steps ";

          endif;

          do until (abs(pt-pt1) < sa_tol);
             
	     pt=pt1;
	     {pt1,df}=br2_a(xa,xb,pt1); 

	     if debg_e;
                "successive approximation "$+ftos(i,"*.*lf",3,0);; 
		pt;; pt1;; abs(pt1-pt);
	     endif;

	     i=i+1;

             if (i > sa_steps); 
               "warning: terminating successive approximation iterations at ";;
	       ""$+ftos(i,"*.*lf",3,0)$+" iterations.";
               break; 
             endif;

          endo;
	  
       endif;

       if debg_e;
         "starting equil with p=";; pt1;
       endif;

       i=0;

       do until (abs(pt-pt1) < ctol);

       i=i+1;

       if (i > maxnewts); 
         "warning: terminating Newton iterations at ";;
	 ""$+ftos(maxnewts,"*.*lf",3,0)$+" iterations.";
	 "xa="$+ftos(xa,"*.*lf",6,3)$+" xb="$+ftos(xb,"*.*lf",6,3);
         break; 
       endif;

       pt=pt1;

       {pt1,df}=br2_a(xa,xb,pt);


       if debg_e;
         "newton "$+ftos(i,"*.*lf",1,0);; 
	 pt;; (pt-pt1);; (1-df);; pt-(pt-pt1)/(1-df); 
       endif;

       pt1=pt-(pt-pt1)/(1-df);

       if (pt1 < 0);
          pt1=pt/2;
       endif;

       endo;

       {pt,df}=br_b(pt1,xb);

       retp(pt1,pt);

    endp;