/*

Bounds:
-10 <= X <= 10   Positive, Firm 2 sells to Firm 1.  Negative, the converse
  0 <= P <= 39
  0 <= Q1 <= B
  0 <= Q2 <= B


Profits:

PI1 = 39*Q1 - P*X - 3*(Q1 - 10 - X)^2   
PI2 = 39*Q2 + P*X -   (Q2 - 10 + X)^2 


Permit Demand and Supply:

Demand FOC: P = 6*Q1 - 60 - 6*X      
Supply FOC: P = 2*Q2 - 20 + 2*X

Equlibrium: P(Q1,Q2) = (3*Q1 + 3*Q2 - 60)/2      0 <= P <= 39
            X(Q1,Q2) = (3*Q1 - Q2 - 20)/4      -10 <= X <= 10

Derivatives: DP1 = 3/2
             DP2 = 3/2
             DX1 = 3/4
             DX2 = -1/4


Reaction Functions:
i.e. lhs of R1 FOC is the derivatve of PI1 wrt Q1; R2 is same for PI2 wrt Q2

R1 FOC  39 - DP1*X - P*DX1 - 6*(Q1 - 10 -X)(1 - DX1) = 0
R2 FOC  39 + DP2*X + P*DX2 - 2*(Q2 - 10 -X)(1 + DX2) = 0

Q1 to solve R1FOC given Q2 computed using a root finder over 0 <= Q1 <= B
Q2 to solve R2FOC given Q1 computed using a root finder over 0 <= Q2 <= B

*/

#include "libscl.h"

using namespace std;
using namespace scl;

namespace {

  REAL const B = 20;                   // Bound on Q1 and Q2

  REAL const DP1 = 3.0/2.0;
  REAL const DP2 = 3.0/2.0;
  REAL const DX1 = 3.0/4.0;
  REAL const DX2 = - 1.0/4.0;

  REAL P(REAL Q1, REAL Q2) 
  {
    REAL bound = 39.0;
    REAL price = (3*Q1 + 3*Q2 - 60.0)/2.0;
    if (price > bound) price = bound;
    else if (price < 0) price = 0.0;
    return price;
  } 

  REAL X(REAL Q1, REAL Q2) 
  {
    if (P(Q1,Q2) == 0.0) return 0.0;
    REAL bound = 10.0;
    REAL quantity = (3*Q1 - Q2 - 20.0)/4.0;
    if (quantity > bound) quantity = bound;
    else if (quantity < -bound) quantity = -bound;
    return quantity;
  } 

  REAL PI1(REAL Q1, REAL Q2)
  {
    return 39*Q1 - P(Q1,Q2)*X(Q1,Q2) - 3*pow(Q1 - 10 - X(Q1,Q2), 2);
  }

  REAL PI2(REAL Q1, REAL Q2)
  {
    return 39*Q2 + P(Q1,Q2)*X(Q1,Q2) - pow(Q2 - 10 + X(Q1,Q2), 2);
  }

  class R1FOC : public nleqn_base {
  private:
    REAL Q2;
  public:
    R1FOC(REAL q2) : Q2(q2) {}
    REAL operator() (REAL Q1) 
    {
      return 39 - DP1*X(Q1,Q2) - P(Q1,Q2)*DX1 - 6*(Q1-10-X(Q1,Q2))*(1 - DX1); 
    }
  };
    
  class R2FOC : public nleqn_base {
  private:
    REAL Q1;
  public:
    R2FOC(REAL q1) : Q1(q1) {}
    REAL operator() (REAL Q2) 
    {
      return 39 + DP2*X(Q1,Q2) + P(Q1,Q2)*DX2 - 2*(Q2-10+X(Q1,Q2))*(1 + DX2);
    }
  };

  REAL R1(REAL Q2)
  {
    R1FOC f(Q2);
    return nlroot(0.0,B,f,1.0e-5);
  }

  REAL R2(REAL Q1)
  {
    R2FOC f(Q1);
    return nlroot(0.0,B,f,1.0e-5);
  }
}

int main(int argc, char** argp, char** envp)
{

  stringstream header;
  header << '\t' << "            Best                             \n"
     << '\t' << "Quantity    Response     Permits          Profit \n"
     << '\t' << "--------  ----------  ----------   ------------- \n"
     << '\t' << "  Q1  Q2    Q1    Q2     X     P     PI1     PI2 \n";


  cout << starbox("/Equilibrium Solutions//") << header.str();
  for (INTEGER q1 = 1; q1<=B; ++q1) {
    for (INTEGER q2 = 1; q2<=B; ++q2) {
      REAL Q1 = q1;
      REAL Q2 = q2;
      if (fabs(Q1-R1(Q2)) <= 1.0 && fabs(Q2-R2(Q1)) <= 1.0 && P(Q1,Q2) > 0.0) {
        cout << '\t' << fmt('d',4,q1) << fmt('d',4,q2)
             << fmt('f',6,1,R1(Q2)) << fmt('f',6,1,R2(Q1))
             << fmt('f',6,1,X(Q1,Q2)) << fmt('f',6,1,P(Q1,Q2))
             << fmt('f',8,1,PI1(Q1,Q2)) << fmt('f',8,1,PI2(Q1,Q2))
             << '\n';
       }
     }
  }

  return 0;
}