{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "다음 조건을 만족시키는 두 자연수 $a$, $b$의 모든 순서쌍 $(a, b)$의 개수를 구하여라.\n",
    "\n",
    "- $1 \\le a \\le 10$, $1 \\le b \\le 100$\n",
    "- 곡선 $y=2^x$이 원 $(x-a)^2 + (y-b)^2 = 1$과 만나지 않는다\n",
    "- 곡선 $y=2^x$이 원 $(x-a)^2 + (y-b)^2 = 4$와 적어도 한 점에서 만난다"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "먼저 임의의 점 $(p, q)$와 $y=2^x$와의 거리를 측정해보도록 하자. 이는 $y=2^x$라는 constraint가 존재할 때 $(x-p)^2 + (y-q)^2$를 최소화시키는 문제와 같다.\n",
    "\n",
    "따라서 라그랑지안 $\\mathcal{L}$은 다음과 같이 정의된다\n",
    "\n",
    "$$\n",
    "\\mathcal{L}(x, y, \\lambda) = (x-p)^2 + (y-q)^2 + \\lambda(xlog2 - logy)\n",
    "$$\n",
    "\n",
    "$\\nabla\\mathcal{L} = 0$이어야 하므로 다음 식이 성립한다\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\nabla_{x}\\mathcal{L} &= 2x - 2p + \\lambda log2 = 0 \\\\\n",
    "\\nabla_{y}\\mathcal{L} &= 2y - 2q - \\frac{\\lambda}{y} = 0 \\\\\n",
    "\\nabla_{\\lambda}\\mathcal{L} &= xlog2 - logy = 0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "1, 2번식을 정리하면 $x$와 $y$를 $\\lambda$에 대해서 표현할 수 있다.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "x &= p - \\frac{\\lambda log2}{2} \\\\\n",
    "y &= \\frac{q + \\sqrt{q^2 + 2\\lambda}}{2}(\\because y \\gt 0)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "이 식을 3번에 대입하면 $\\lambda$에 대한 식으로 바뀌고, 이는 newton raphson method로 해를 구할 수 있다. 이렇게 구한 좌표 $(x, y)$와 $(p, q)$간의 거리가 $y=2^x$와 $(p, q)$ 사이의 거리가 된다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from scipy import optimize\n",
    "from functools import partial\n",
    "log2 = math.log(2)\n",
    "\n",
    "def x(p, l):\n",
    "    return p - l*log2/2\n",
    "    \n",
    "def y(q, l):\n",
    "    return (q + math.sqrt(q**2 + 2*l))/2\n",
    "\n",
    "def f(p, q, l):\n",
    "    return x(p, l)*log2 - math.log(y(q, l))\n",
    "\n",
    "def dist(p, q):\n",
    "    \"\"\"Get Euclidean distance between y=2^x and (p, q)\"\"\"\n",
    "    l = optimize.newton(lambda l: f(p, q, l), 1.0)\n",
    "    x_ = x(p, l)\n",
    "    y_ = 2**x_\n",
    "    return math.sqrt((x_-p)**2 + (y_-q)**2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7012892057857038e-15\n"
     ]
    }
   ],
   "source": [
    "# test\n",
    "print(dist(3, 8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 5)\n",
      "(1, 6)\n",
      "(1, 7)\n",
      "(1, 8)\n",
      "(2, 1)\n",
      "(2, 9)\n",
      "(2, 10)\n",
      "(2, 11)\n",
      "(2, 12)\n",
      "(2, 13)\n",
      "(2, 14)\n",
      "(2, 15)\n",
      "(2, 16)\n",
      "(3, 2)\n",
      "(3, 3)\n",
      "(3, 17)\n",
      "(3, 18)\n",
      "(3, 19)\n",
      "(3, 20)\n",
      "(3, 21)\n",
      "(3, 22)\n",
      "(3, 23)\n",
      "(3, 24)\n",
      "(3, 25)\n",
      "(3, 26)\n",
      "(3, 27)\n",
      "(3, 28)\n",
      "(3, 29)\n",
      "(3, 30)\n",
      "(3, 31)\n",
      "(3, 32)\n",
      "(4, 4)\n",
      "(4, 5)\n",
      "(4, 6)\n",
      "(4, 7)\n",
      "(4, 33)\n",
      "(4, 34)\n",
      "(4, 35)\n",
      "(4, 36)\n",
      "(4, 37)\n",
      "(4, 38)\n",
      "(4, 39)\n",
      "(4, 40)\n",
      "(4, 41)\n",
      "(4, 42)\n",
      "(4, 43)\n",
      "(4, 44)\n",
      "(4, 45)\n",
      "(4, 46)\n",
      "(4, 47)\n",
      "(4, 48)\n",
      "(4, 49)\n",
      "(4, 50)\n",
      "(4, 51)\n",
      "(4, 52)\n",
      "(4, 53)\n",
      "(4, 54)\n",
      "(4, 55)\n",
      "(4, 56)\n",
      "(4, 57)\n",
      "(4, 58)\n",
      "(4, 59)\n",
      "(4, 60)\n",
      "(4, 61)\n",
      "(4, 62)\n",
      "(4, 63)\n",
      "(4, 64)\n",
      "(5, 8)\n",
      "(5, 9)\n",
      "(5, 10)\n",
      "(5, 11)\n",
      "(5, 12)\n",
      "(5, 13)\n",
      "(5, 14)\n",
      "(5, 15)\n",
      "(5, 65)\n",
      "(5, 66)\n",
      "(5, 67)\n",
      "(5, 68)\n",
      "(5, 69)\n",
      "(5, 70)\n",
      "(5, 71)\n",
      "(5, 72)\n",
      "(5, 73)\n",
      "(5, 74)\n",
      "(5, 75)\n",
      "(5, 76)\n",
      "(5, 77)\n",
      "(5, 78)\n",
      "(5, 79)\n",
      "(5, 80)\n",
      "(5, 81)\n",
      "(5, 82)\n",
      "(5, 83)\n",
      "(5, 84)\n",
      "(5, 85)\n",
      "(5, 86)\n",
      "(5, 87)\n",
      "(5, 88)\n",
      "(5, 89)\n",
      "(5, 90)\n",
      "(5, 91)\n",
      "(5, 92)\n",
      "(5, 93)\n",
      "(5, 94)\n",
      "(5, 95)\n",
      "(5, 96)\n",
      "(5, 97)\n",
      "(5, 98)\n",
      "(5, 99)\n",
      "(5, 100)\n",
      "(6, 16)\n",
      "(6, 17)\n",
      "(6, 18)\n",
      "(6, 19)\n",
      "(6, 20)\n",
      "(6, 21)\n",
      "(6, 22)\n",
      "(6, 23)\n",
      "(6, 24)\n",
      "(6, 25)\n",
      "(6, 26)\n",
      "(6, 27)\n",
      "(6, 28)\n",
      "(6, 29)\n",
      "(6, 30)\n",
      "(6, 31)\n",
      "(7, 32)\n",
      "(7, 33)\n",
      "(7, 34)\n",
      "(7, 35)\n",
      "(7, 36)\n",
      "(7, 37)\n",
      "(7, 38)\n",
      "(7, 39)\n",
      "(7, 40)\n",
      "(7, 41)\n",
      "(7, 42)\n",
      "(7, 43)\n",
      "(7, 44)\n",
      "(7, 45)\n",
      "(7, 46)\n",
      "(7, 47)\n",
      "(7, 48)\n",
      "(7, 49)\n",
      "(7, 50)\n",
      "(7, 51)\n",
      "(7, 52)\n",
      "(7, 53)\n",
      "(7, 54)\n",
      "(7, 55)\n",
      "(7, 56)\n",
      "(7, 57)\n",
      "(7, 58)\n",
      "(7, 59)\n",
      "(7, 60)\n",
      "(7, 61)\n",
      "(7, 62)\n",
      "(7, 63)\n",
      "(8, 64)\n",
      "(8, 65)\n",
      "(8, 66)\n",
      "(8, 67)\n",
      "(8, 68)\n",
      "(8, 69)\n",
      "(8, 70)\n",
      "(8, 71)\n",
      "(8, 72)\n",
      "(8, 73)\n",
      "(8, 74)\n",
      "(8, 75)\n",
      "(8, 76)\n",
      "(8, 77)\n",
      "(8, 78)\n",
      "(8, 79)\n",
      "(8, 80)\n",
      "(8, 81)\n",
      "(8, 82)\n",
      "(8, 83)\n",
      "(8, 84)\n",
      "(8, 85)\n",
      "(8, 86)\n",
      "(8, 87)\n",
      "(8, 88)\n",
      "(8, 89)\n",
      "(8, 90)\n",
      "(8, 91)\n",
      "(8, 92)\n",
      "(8, 93)\n",
      "(8, 94)\n",
      "(8, 95)\n",
      "(8, 96)\n",
      "(8, 97)\n",
      "(8, 98)\n",
      "(8, 99)\n",
      "(8, 100)\n",
      "196\n"
     ]
    }
   ],
   "source": [
    "cnt = 0\n",
    "for a in range(1, 11):\n",
    "    for b in range(1, 101):\n",
    "        if 1 < dist(a, b) <= 2:\n",
    "            print((a, b))\n",
    "            cnt += 1\n",
    "            \n",
    "print(cnt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}