eqlbrm35.f

      function eqlbrm_f_g_h(b, astar, bigb)
      external bigb
      eqlbrm_f_g_h = bigb(astar, b)
      end

      function eqlbrm_f_g(a, astar, biga, bigb, bstar, n)
      external biga
      bstar = argmax_2(astar, bigb, bstar, n)
      eqlbrm_f_g = biga(a, bstar)
      end

      function eqlbrm_f(astar, biga, bigb, bstar, n)
      eqlbrm_f = argmax_1(astar, biga, bigb, bstar, astar, n) - astar
      end

C     astar & bstar: guesses in, optimized values out
      subroutine eqlbrm(biga, bigb, astar, bstar, n)
      astar = root_2(biga, bigb, bstar, astar, n)
      end

      function root_1_1(astar, biga, bigb, bstar, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
 1669 x = x - y / yprime
      root_1_1 = x
      end

      function root_1_2(astar, bigb, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2_1_2(astar, bigb, x, y, yprime)
 1669 x = x - y / yprime
      root_1_2 = x
      end

      function root_2(biga, bigb, bstar, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2_2(biga, bigb, bstar, n, x, y, yprime)
 1669 x = x - y / yprime
      root_2 = x
      end

      subroutine deriv2_1_1_adf3(x, astar, biga, bigb, bstar, n, y)
      y = argmax_fprime_1(x, astar, biga, bigb, bstar, n)
      end

      subroutine deriv2_1_1(astar, biga, bigb, bstar, n, x, y, yprime)
      adf(x)
      call deriv2_1_1_adf3(x, astar, biga, bigb, bstar, n, y)
      end adf(yprime = tangent(y))
      end

      subroutine deriv2_1_2_adf4(x, astar, bigb, y)
      y = argmax_fprime_2(x, astar, bigb)
      end

      subroutine deriv2_1_2(astar, bigb, x, y, yprime)
      adf(x)
      call deriv2_1_2_adf4(x, astar, bigb, y)
      end adf(yprime = tangent(y))
      end

      subroutine deriv2_2(biga, bigb, bstar, n, x, y, yprime)
      adf(x)
      y = eqlbrm_f(x, biga, bigb, bstar, n)
      end adf(yprime = tangent(y))
      end

      function argmax_fprime_1(x, astar, biga, bigb, bstar, n)
      argmax_fprime_1 = deriv1_1(astar, biga, bigb, bstar, n, x)
      end

      function argmax_fprime_2(x, astar, bigb)
      argmax_fprime_2 = deriv1_2(astar, bigb, x)
      end

      function argmax_1(astar, biga, bigb, bstar, x0, n)
      argmax_1 = root_1_1(astar, biga, bigb, bstar, x0, n)
      end

      function argmax_2(astar, bigb, x0, n)
      argmax_2 = root_1_2(astar, bigb, x0, n)
      end

      subroutine deriv1_1_adf1(x, astar, biga, bigb, bstar, n, y)
      y = eqlbrm_f_g(x, astar, biga, bigb, bstar, n)
      end

      function deriv1_1(astar, biga, bigb, bstar, n, x)
      adf(x)
      call deriv1_1_adf1(x, astar, biga, bigb, bstar, n, y)
      end adf(deriv1_1 = tangent(y))
      end

      subroutine deriv1_2_adf2(x, astar, bigb, y)
      y = eqlbrm_f_g_h(x, astar, bigb)
      end

      function deriv1_2(astar, bigb, x)
      adf(x)
      call deriv1_2_adf2(x, astar, bigb, y)
      end adf(deriv1_2 = tangent(y))
      end

      function gmbiga(a, b)
      price = 20 - 0.1*a - 0.1*b
      costs = a*(10 - 0.05*a)
      gmbiga = a*price - costs
      end

      function gmbigb(a, b)
      price = 20 - 0.1*b - 0.0999*a
      costs = b*(10.005 - 0.05*b)
      gmbigb = b*price - costs
      end

      program main
      read *, astar
      read *, bstar
      read *, n
      call eqlbrm(gmbiga, gmbigb, astar, bstar, n)
      print *, astar, bstar
      end