polynomial interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import operator\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import base64\n",
    "from io import BytesIO\n",
    "\n",
    "\n",
    "def generate_points(ymin=-3, ymax=4, n=6):\n",
    "    X = [random.randrange(ymin, ymax) for i in range(n)]\n",
    "    x = np.array(X)\n",
    "    \n",
    "    exps = np.arange(6, dtype=int) \n",
    "    \n",
    "    vandermonde = np.ones((n+1, n+1))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "generate_points() got an unexpected keyword argument 'amin'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-53-79516cd083e7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgenerate_points\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mamin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mamax\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: generate_points() got an unexpected keyword argument 'amin'"
     ]
    }
   ],
   "source": [
    "generate_points(amin=-3, amax=4, n=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_points(amin=-3, amax=4, n=5):\n",
    "\n",
    "    A = [random.randrange(amin, amax) for i in range(n-1)]\n",
    "    non_zero = [n for n in range(ymin, ymax) if n!=0]\n",
    "    A_top = random.choice(non_zero)\n",
    "    A.append(A_top)\n",
    "    a = np.array(A, dtype=int)\n",
    "    x  = np.arange(-n//2 + 1, n//2 + 1)\n",
    "\n",
    "    exps = np.arange(n, dtype=int) \n",
    "\n",
    "    vandermonde = ((np.ones((n, n))*x).T**exps).astype(int)\n",
    "    y = vandermonde.dot(a)\n",
    "    \n",
    "    points = list(zip(x, y))\n",
    "    \n",
    "    return points, list(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "points, a = generate_points()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, -1), (1, -2), (2, 0), (3, 0), (4, 1)]"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(enumerate(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAF2CAYAAAB6XrNlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8XNWB9vHfUbclW7IlWZbl3guuGOMYDBbVlNCThWwIZGEdQrIhCS+G5N0syWbT4E2ymwRISExZQiAOHUIA44YL7gV3W3KvsmSr2eo67x8zHo1syWqjuXPvPN/Pxx/fOzOSHl+Jh6sz555rrLWIiIh3xDgdQEREQkvFLiLiMSp2ERGPUbGLiHiMil1ExGNU7CIiHqNiFxHxGBW7iIjHqNhFRDwmLpxfLC0tzQ4dOjScX7LNTp06RXJystMxWqScoeOGjOCOnG7ICO7IuXbt2kJrbWZ7PjasxZ6VlcWaNWvC+SXbbNGiRcyYMcPpGC1SztBxQ0ZwR043ZAR35DTG7Gvvx2ooRkTEY1TsIiIeo2IXEfGYFovdGNPPGLPQGLPNGLPFGPOQ//Gexph5xphd/r97dH5cERFpSWvO2GuBh621o4CpwDeMMaOBx4D51tphwHz/voiIOKzFYrfWHrHWrvNvlwHbgBzgZuBF/8teBG7prJAiItJ6bRpjN8YMBCYCK4Esa+0R8JU/0CvU4UREpO1Ma2+NZ4xJARYDP7HWvmGMKbbWpgU9f9Jae844uzFmFjALIDMz88K5c+eGJnknKS8vJyUlxekYLVLO0HFDRnBHTjdkBHfkzM3NXWutndyuD7bWtvgHiAc+BL4b9NgOINu/nQ3saOnzDB8+3Ea6hQsXOh2hVZQzdNyQ0Vp35HRDRmvdkRNYY1vRz039ac2sGAPMAbZZa38V9NQ7wD3+7XuAt9v1fxYREQmp1iwpcAlwN7DJGLPB/9j3gZ8Dc40x9wH7gS90TkQREWmLFovdWrsUMM08fWVo44iIOKeu3hJjwDdQ4V668lRExG9pXiGTfjyP+19cw1vrDzkdp91U7CIifmv2nuDk6Ro+3naMDQeKnY7Tbip2ERG/VXtOBLYnD3TvKikqdhERoKq2rtFZ+pSBPR1M0zEqdhERYPOhEqpq6wEYkN6VXt2THE7Ufip2ERFg9d6Tge3JA9x7tg4qdhERAFYHja9PGeTe8XVQsYuIUF9vWbOv4Yz9IhePr4OKXUSEXQXllFTUAJCRksCgjGSHE3WMil1Eot6qvUHTHAf01JWnIiJuFzy+ftEgdw/DgIpdRKKctZbVQWfsF7n4wqQzVOwiEtUOFVdwpKQSgOSEWEZnd3c4Ucep2EUkqgUvIzBpQA/iYt1fi+7/F4iIdEBwsbt5GYFgKnYRiWrBxX7x4HQHk4SOil1EolZBWSW7C08BkBAXw7i+qQ4nCg0Vu4hEreCz9Qn90kiKj3UwTeio2EUkajUahvHA/PUzVOwiErUavXGqYhcRcbfi09VsP1oGQGyMYVJ/91+YdIaKXUSiUvD662NzUklOjHMwTWip2EUkKq3aUxTY9tL4OqjYRSRKrfTo+Dqo2EUkClXUWjYfKgHAGPffCu9sKnYRiTo7T9ZRb33bo7O7k9o13tlAIaZiF5Gos+NEfWD74kHeWEYgmIpdRKLO9hN1ge2pg701DAMqdhGJMuVVtewt9Z2xG+O9N05BxS4iUWbN3hOB8fWRvbuT1jXB2UCdQMUuIlFlxe6GaY5eHIYBFbuIRJmVQRcmTfXI+utn8841tCIi51NeTvXPf8GcX/2GtIoyTnbpRteYf4PvPwYpKU6nCykVu4h4X3k5TJ1KbF4ePauqAEivKIVf/xLeeQtWrPBUuWsoRkS878knIT+fWH+pB1RWQn6+73kPUbGLiPc9/bSvxJtSWQnPPBPePJ1MxS4i3ldU1LHnXUbFLiLel97C7JeWnncZFbuIeN+DD1KTkNj0c0lJ8PWvhzdPJ1Oxi4j3PfIIh3r2oTLurKtMk5JgyBB45BFncnUSFbuIeF5JTCI3fOkJfj/ldoq6pmKNgcxMmD3bc1MdQfPYRSQKrNxTxKn4Lvz39H9m4V0P8p0LapkxY4bTsTqNzthFxPM+3R20jMAQb71R2hQVu4h43qf5DcX+OY+uDxNMxS4inlZUXsX2o2UAxMUYLhrozRUdg6nYRcTTgpfpHdc3leRE77+1qGIXEU9bll8Y2L5kaIaDScJHxS4inhY8vj5tiIpdRMTVDhdXsKfwFACJcTFM7J/mcKLwULGLiGcty2sYhrloYE+S4mMdTBM+KnYR8azlwcMwQ70/zfEMFbuIeJK1luXBb5xGyfg6qNhFxKPyj5/iWKnvjkndk+K4ICfV4UTho2IXEU8KPlufOjid2BjjYJrwUrGLiCcFv3E6LQrWhwmmYhcRz6mrt42uOI2WC5POULGLiOdsPlRCSUUNAL26JTK0l7fWW2+Jil1EPGdp0DDMpUMzMCZ6xtdBxS4iHrR0V1CxD4uuYRhQsYuIx1RU17F238nAfrSNr0Mrit0Y85wxpsAYsznosR8aYw4ZYzb4/1zfuTFFRFpn1d4TVNfVAzA8K4Ws7kkOJwq/1pyxvwDMbOLxX1trJ/j/vB/aWCIi7bN01/HAdjSerUMrit1a+wlwoqXXiYhEgqV5DevDTI/C8XUAY61t+UXGDATes9Ze4N//IXAvUAqsAR621p5s5mNnAbMAMjMzL5w7d24IYnee8vJyUlIif2qUcoaOGzKCO3I6nbGkyvLQwtMAxBr43ZVd6RJ37owYp3O2Rm5u7lpr7eR2fbC1tsU/wEBgc9B+FhCL74z/J8Bzrfk8w4cPt5Fu4cKFTkdoFeUMHTdktNYdOZ3O+Nb6g3bAo+/ZAY++Z7/wzPJmX+d0ztYA1thW9GpTf9o1K8Zae8xaW2etrQf+CExp1/9VRERCKHiaY7SOr0M7pzsaY7KDdm8FNjf3WhGRcLDWsiTK56+f0eLtuo0xrwAzgAxjzEHgcWCGMWYCYIG9wNc6MaOISIvyCso5WloJ+JbpHd83epbpPVuLxW6tvauJh+d0QhYRkXZbvLPxNMe42Oi9/jJ6/+Ui4inBwzDTh2U6mMR5KnYRcb3KmjpW7tH89TNU7CLiemv2nqSyxreMwOCMZPr17OpwImep2EXE9ZYELSMQ7WfroGIXEQ/4ROPrjajYRcTVCsoq2XakFID4WMPnouz+pk1RsYuIqwVfbTqxfw+SE1ucxe15KnYRcbXg+euXD9cwDKjYRcTF6uotn6jYz6FiFxHX2nyohJOnawDI7JbImD7dHU4UGVTsIuJawcMwlw3LxJhz116PRip2EXGtRuPrIzQMc4aKXURcqeR0Dev3+27cZgxMj+L118+mYhcRV1qSd5x6/509x/dNo0dygrOBIoiKXURcafEOzYZpjopdRFzHWssnuzS+3hwVu4i4zrYjZRwrrQIgtUs84/umOZwosqjYRcR1Fu4oCGxfPjyT2BhNcwymYhcR11kUVOy5IzUMczYVu4i4SsnpGtbua5jmeJmW6T2Hil1EXOWTXY2nOaanJDobKAKp2EXEVYLH13NH9HIwSeRSsYuIa9TX20bz1zW+3jQVu4i4xqZDJRSdqgYgIyWBC/qkOpwoMqnYRcQ1Gk9z7EWMpjk2ScUuIq6xcHtDsc/Q1abNUrGLiCsUlFWy8WAJALExhsu0PkyzVOwi4gqLtje8aXrRwB6kdol3ME1kU7GLiCvM334ssH3lyCwHk0Q+FbuIRLyq2jqW7CoM7F8xSvPXz0fFLiIRb+XuE5yurgNgYHpXBmckO5wosqnYRSTiLQiaDXPFyCzdtLoFKnYRiWjW2sbj6xqGaZGKXUQiWl5BOQdOVACQkhjHRQN7Opwo8qnYRSSifbytYRhm+rAMEuJUWy3RERKRiPbxtoZhmKtGaZpja6jYRSRiFZZXsW6/76YaMQZyR2p8vTVU7CISsRZsL8D6b6oxeUBPeiYnOBvIJVTsIhKxPt4aNAwzWmfrraViF5GIVFnT+GpTja+3nopdRCLS8vxCKmp8V5sOzkxmcGaKw4ncQ8UuIhFp3taGaY5X62y9TVTsIhJx6ust84OnOY5WsbeFil1EIs7Gg8UUlFUB0KNrPJP693A4kbuo2EUk4ny0tfFFSbG6t2mbqNhFJOJ8tOVoYPuaMb0dTOJOKnYRiSh5BeXkHz8FQJf4WKYPy3A4kfuo2EUkoswLGoa5bHgGSfGxDqZxJxW7iESUj7YGDcOM1jBMe6jYRSRiFJRWsn5/MQCxMYYrtOhXu6jYRSRizAuauz5lYE96aNGvdlGxi0jE+HBLQ7FfM0YXJbWXil1EIkJJRQ3L8xoW/dI0x/ZTsYtIRJi/7Ri19b7F18f1TSUnrYvDidxLxS4iEeGDzQ2zYWZeoLP1jlCxi4jjTlXVsnjn8cD+TA3DdIiKXUQct3jncapq6wEYkdVNa693kIpdRBz3j6BhmGs1DNNhKnYRcVRlTR0LguavX6di77AWi90Y85wxpsAYsznosZ7GmHnGmF3+v7VYsoi0y7K8Qk5V+26BNyC9KyN7d3M4kfu15oz9BWDmWY89Bsy31g4D5vv3RUTa7P1NjWfDGKO11zuqxWK31n4CnDjr4ZuBF/3bLwK3hDiXiESB6tp65gUt+nX9BdkOpvGO9o6xZ1lrjwD4/9ZKPSLSZsvyCymtrAWgb48ujOub6nAibzDW2pZfZMxA4D1r7QX+/WJrbVrQ8yettU2OsxtjZgGzADIzMy+cO3duCGJ3nvLyclJSIn+qlXKGjhsygjtytjXjnE1VLDnkK/aZA+O5c2R4Fv1yw7HMzc1da62d3J6PjWvn1zxmjMm21h4xxmQDBc290Fr7LPAswIgRI+yMGTPa+SXDY9GiRUR6RlDOUHJDRnBHzrZkrKmr56HFHwf2v3b9RUwM002r3XAsO6K9QzHvAPf4t+8B3g5NHBGJFsvyCimpqAEgJ60LE/qltfAR0lqtme74CvApMMIYc9AYcx/wc+BqY8wu4Gr/vohIq/0jaDbMdZoNE1ItDsVYa+9q5qkrQ5xFRKJETV09HwbPhhmn2TChpCtPRSTsluUVUnzaNwyTnZrEhL4ahgklFbuIhN17nx0JbN8wNpuYGA3DhJKKXUTCqqq2jg+3NAzD3Di+j4NpvEnFLiJhtWRnIWX+i5L69ezCeF2UFHIqdhEJq/c+OxzYvmFsH82G6QQqdhEJm8qaOuZtbVii90bNhukUKnYRCZtFOwoCS/QOykhmTJ/uDifyJhW7iITNuxsbZsPcOC5bwzCdRMUuImFRXlXLx9uCh2E0G6azqNhFJCzmbT0auGH1yN7dGKE7JXUaFbuIhMXbGxpmw9w0QWfrnUnFLiKdrqi8iiW7CgP7n9cwTKdSsYtIp3t/0xHq6n039Zk8oAf9enZ1OJG3qdhFpNO9s7FhGOZmDcN0OhW7iHSqgydPs3rvSQBiYwzXj9VFSZ1NxS4inSp47vqlQzNIT0l0ME10ULGLSKex1vLm+oOBfQ3DhIeKXUQ6zdYjpew8Vg5Al/hYrh3T2+FE0UHFLiKd5s11hwLb147JIjmxxbtxSgio2EWkU9TVW94Omg1z66S+DqaJLip2EekUy/IKOV5WBUBGSiKXDEl3OFH0ULGLSKd4a33DMMxN4/sQF6u6CRcdaREJudPVtXwQdF/T2yblOJgm+qjYRSTkPth8lNP+G2oM7ZWiG2qEmYpdRELu9XUNc9dvm5SjG2qEmYpdRELqUHEFy/OLAIgxcNtEzYYJNxW7iITUm+sOYn0LOXLJ0Ax6pyY5GygKqdhFJGSstbwedFHSHRfqbN0JugxMREImv7iePYWVAHRLjOOa0VpCwAk6YxeRkFl6qDawfcO4bLokxDqYJnqp2EUkJCqq61h5tKHYb9cwjGNU7CISEh9uOUqFv9cHpndl8oAezgaKYip2EQmJv64+ENj+wuR+mrvuIBW7iHTY/qLTfLq7Ye66ZsM4S8UuIh322tqGs/UZI3qR1V1z152kYheRDqmrt7y2tmEJgS9O1tm601TsItIhS/MKOVzin7ueAFeMzHI4kajYRaRD/rp6f2B7Wp84EuJUK07Td0BE2q2wvIp5W48F9i/LiXcwjZyhYheRdnt97UFq6nwrfl04oAc53VQpkUDfBRFpF2strwbNXb9rSn8H00gwFbuItMuK3SfYU3gKgG5JcdwwNtvhRHKGil1E2uXVoDdNb5mQowW/IoiKXUTa7OSpav6xqeFm1RqGiSwqdhFps9fWHqS6rh6A8X1TGa2bVUcUFbuItEl9veUvqxqGYb50sc7WI42KXUTaZHl+UaM3TT8/vo/DieRsKnYRaZM/r9gX2L59Ul+6JugOm5FGxS4irXastJJ52xquNP1nDcNEJBW7iLTaq6sOUFfvu9J0yqCeDMvq5nAiaYqKXURapbaunleC3jTV2XrkUrGLSKvM23qMo6W+5XnTkxOYeUFvhxNJc1TsItIqL366N7D9pYv7kxinK00jlYpdRFq042gZK3afACA2xmjueoRTsYtIi15asTewfe2YLLJTuzgXRlqkYheR8yqtrOGNdYcC+1/53EDnwkirqNhF5LxeX3uQ09V1AIzI6sbFg3o6nEhaomIXkWbV11teWL43sH/35wZgjHEukLSKil1EmrVgewH7ik4D0D0pjtsm5TicSFpDxS4izXpu2Z7A9l0X99e6MC7Roe+SMWYvUAbUAbXW2smhCCUiztt+tJTl+UWAb4qj3jR1j1D87zfXWlsYgs8jIhHk+aV7A9szx/QmJ01THN1CQzEico6i8ire3NAwxfFfLh3oXBhps44WuwU+MsasNcbMCkUgEXHen1fsp7rWd+u7cX1TmdS/h8OJpC2Mtbb9H2xMH2vtYWNML2Ae8G/W2k/Oes0sYBZAZmbmhXPnzu1I3k5XXl5OSkqK0zFapJyh44aMEL6c1XWWhxefpqzat//AuESm9mndqK2OZejk5uaube/7lh0aY7fWHvb/XWCMeROYAnxy1mueBZ4FGDFihJ0xY0ZHvmSnW7RoEZGeEZQzlNyQEcKX85VV+ymr3gRAn9QkHv6nXOJjW/fLvY5lZGj3UIwxJtkY0+3MNnANsDlUwUQk/OrrLX9asjuw/9VLBrW61CVydOSMPQt4038VWhzwF2vtByFJJSKOWLSzgPzj/htVJ8Zx55R+DieS9mh3sVtrdwPjQ5hFRBz27CcNZ+t3TulHt6R4B9NIe+l3LBEBYMOB4kZrrt97ySCHE0l7qdhFBIDfL8oPbN80vo8uSHIxFbuIkH+8nA+3Hg3sf+3ywQ6mkY5SsYsIf/xkN2cuackdkcnI3t2dDSQdomIXiXLHSisb3SHpgcuHOJhGQkHFLhLl5izdQ3Wdb/mASf3TmKI7JLmeil0kip08Vc2fV+wL7D9w+RDdIckDVOwiUez5ZXsa3c/0qlFZDieSUFCxi0Sp0sqaRvcz/cYVQ4mJ0dm6F6jYRaLUS5/uo7SyFoBBGcncMDbb4UQSKip2kShUUV3Hc0sb7mf69RlDiNXZumeo2EWi0Msr91F0yrfgek5aF26dmONwIgklFbtIlKmoruP3ixuWD3hgxhAtzesx+m6KRJmXV+6jsNx3tt4nNYkvTu7rcCIJNRW7SBQ5+2z9wdyhJMbFOphIOoOKXSSKnH22/gWdrXuSil0kSpyurtXZepRQsYtEieeX7dXZepRQsYtEgZKKGv4QdLb+rSuH6Wzdw1TsIlFgzpLdgatMB6R35fYLdbbuZSp2EY8rKq9iTtBVpt+5arjmrXucvrsiHvfMonxO+VdwHJ6VwufH93E4kXQ2FbuIhx0qruB/P21Yb/27V4/QmjBRQMUu4mG/nrczcHek8f3SuHaM1luPBip2EY/afrSU19cdDOx/77qRujtSlFCxi3jUkx/swFrfdu6ITKYOTnc2kISNil3Eg1buLmL+9gIAjIHZM0c6nEjCScUu4jH19Zb/+vu2wP4tE3IYld3dwUQSbip2EY95e+MhNh0qASAxLob/c+0IhxNJuKnYRTyksqaOJz/YEdi/79JB5KR1cTCROEHFLuIhc5bu4XBJJQAZKQl8fcYQhxOJE1TsIh5RUFrJ0wvzAvvfvmo43ZLiHUwkTlGxi3jELz7YEVg6YFivFO68qJ/DicQpKnYRD1i//2Sji5H+4/OjidNCX1FL33kRl6uvt/zw3a2B/atHZzF9WKaDicRpKnYRl3tj/SE2HigGICE2hn+/YZTDicRpKnYRFyupqOHn/2i4GOn+6YMYkJ7sYCKJBCp2ERf71Uc7Avcx7d09iW/kDnU4kUQCFbuIS20+VMJLKxrWWv/BjaNJToxzMJFEChW7iAvV11t+8PZm6v2rN04flsH1Y3s7G0oihopdxIX+uuYA6/c3vGH6o5vGaK11CVCxi7hMQVklP3u/4Q3TWZcNZnBmioOJJNKo2EVc5j/f3UppZS0AA9K78s0r9IapNKZiF3GRhdsLeO+zI4H9n946lqT4WAcTSSRSsYu4RGWt5d/f2hzYv21iDpcMzXAwkUQqFbtIJCsvh8cfh8xMrr36St75z1v59pKX6RNby//VFabSDE16FYlU5eUwdSrk50NlJQZIryjlgVWv89Vja0n93pVAotMpJQLpjF0kUj35ZKDUgyXVVtP90H7f8yJNULGLRKqnnz6n1M8wlZXwzDNhDiRuoWIXiVRFRR17XqKWil0kQtn09PO/oKXnJWqp2EUi1Pzc26mMS2j6yaQk+PrXwxtIXEPFLhKBPt56jG/1vYp9ab3PLfekJBgyBB55xJlwEvFU7CIRpqC0ktmvf8bphC7cevcvmf/5e7CZmVhjIDMTZs+GFSsgRevDSNM0j10kgtTXWx7+20ZOnPLdPKNbRhrTfvw7TPKzLFq0iBkzZjgbUFxBZ+wiEWTO0j0s2VUIgDHw6y9OoEdyM+PsIs1QsYtEiM8OFvPEh9sD+1+7bAjTtBaMtIOKXSQClFTU8I2/rKOmzndLpHF9U/nu1cMdTiVupWIXcZi1ltmvbeTAiQoAuiXG8du7JpIQp/88pX30kyPisOeX7eXDLccC+0/cMY4B6ckOJhK361CxG2NmGmN2GGPyjDGPhSqUSLRYvfcEPw26zd290wZy3dhsBxOJF7S72I0xscBTwHXAaOAuY8zoUAUT8bqC0koefHkdtfUN4+rfu36kw6nECzpyxj4FyLPW7rbWVgOvAjeHJpaIt1XX1vPgy+s4XlYFQM/kBJ758oUkxuk2d9Jxxlrbvg805g5gprX2fv/+3cDF1tpvnvW6WcAsgMzMzAvnzp3bscSdrLy8nBQXXNGnnKHjRMYXt1Sx8IDvhtQGeOSiJEann7/UdSxDxw05c3Nz11prJ7fnYzty5alp4rFz/i9hrX0WeBZgxIgRNtKvnHPL1X3KGTrhzvjSp3tZeGBLYP/R60bywOVDWvw4HcvQcUvO9urIUMxBoF/Qfl/gcMfiiHjbsrxCfvju1sD+TeP78LXLBjuYSLyoI8W+GhhmjBlkjEkA7gTeCU0sEe/JP17Ogy+voy7ozdIn7hiHMU398ivSfu0eirHW1hpjvgl8CMQCz1lrt7TwYSJRqbC8iq8+v5qSihoAenVL5Nm7J5MUrzdLJfQ6tLqjtfZ94P0QZRHxpMqaOu5/cQ37T5wGoEt8LH/8ymR6pyY5nEy8SleeinSiunrLQ6+uZ8OBYsC3YuNv7prI+H5pDicTL4uqYj9VVet0BIki1lr+/a1NjZYLePzG0Vw9OsvBVBINoqbY6+otNz+1jPtfXMNnB4udjiNR4Jcf7eSVVQcC+/dfOoh7LxnkYCKJFlFzB6X3Nx0hr6CcvIJyPt52jNwRmfzblcOY1L+H09HEg/60ZDe/W5gX2L9tYg7fv36Ug4kkmkTNGfvWI6UEzypbuOM4tz29nLvnrGTdlv3w+OOQmcnlV1zhu6/k449DeblzgcW1Xvp0L//194aFvXJHZPKLO8YRE6NpjRIeUVPsj84cyYffvoybxvdpVPBrtxwgecZ0qn/2cygsxFgLhYXwxBMwdarKXdpk7uoD/ODthlm/Fw3swdP/fCHxsVHzn5pEgKj6aRue1Y3f3DWRed+5nFsn5hBjYNbKNxhQfJSEmurGL66shPx8ePJJZ8KK6/xtzQEefeOzwP6Efmk8d+9FdEnQXHUJr6gq9jOG9krh1/80gfkPz+D+Tf8gqba66RdWVmKfeSa84cSVXl21n9mvf8aZNfVGZ3fnxa9OoVtSvLPBJCpFZbGfMSgjmZTykvO+xhYWsmTXcdq7CqZ430sr9vHYG5salfqf77+Y1K4qdXFGVBc7AOnp5336ZFJ37p6zitufWc6iHQUqeAmw1vLUwjx+8NbmwGNjc1L5y79eTM/kBAeTSbRTsT/4ICQ1fWl3ZVwCL028HoB1+4u59/nV3PL0cuZvO6aCj3LWWn7y9208+eGOwGPj+6Xx5/svJq2rSl2cpWJ/5BEYMuTcck9KInboEEq/+RAJQTMaNh4o5r4X13Djb5fyweYj1Ner4KNNdW093527kT8t3RN4bNqQdF6+/2JSu2j4RZynYk9JgRUrYPZsyMzEGuObxz57NvGrV/Efd13M4tkzuHfaQBLiGg7XlsOlPPDndVz3P0t4Z+PhwFKs4m0lFTXc89wq3lx/KPDYtWOyeO7ei0hJjJrr/STCqdjBV+4/+hEUFLB4wQIoKPDt+2+dlZ3ahR/eNIYls3O579JBJMU3HLYdx8r41ivruepXi5m7+gDVtfVO/Sukk+0vOs0dzyzn091FgcfumtKPp740ScvvSkRRsbdBVvckfnDjaJY+egVfu3wwyUHzk/cUnmL2658x48mFvLBsDxXVdQ4mlVBbnl/ITU8tZVdBwwVrs2eO4Ke3jiVOFx9JhNFPZDtkpCTyvetGseyxK/jWlcPontTwK/jhkkp++O5WLv3FAn63YFfgxgriTtZaXli2h7vnrKL4tO97mRAXw//cOYEHZwzV3Y8kIqnYOyCtawLfvXo4yx67gtkzR5AeNMWt6FQ1/++jnUz72Xx++v42jpVWOphU2uNUVS3fenUDP3x3a+A9lF7dEvnrrKncPCHH4XQizdOvnHsRAAAM0UlEQVS7PSHQLSmeB2cM5avTBjF3zQGe/WQ3h4orADhVXcezn+zm+WV7uGVCDrMuG8ywrG4OJ5aW7DhaxoMvryX/+KnAY+P6pvLs3brzkUQ+FXsIdUmI5Z5pA/nSxf15e8Nh/rA4PzAmW1Nn+dvag/xt7UGuGNmLf50+mKmDe+pX+QhjreV/P93HT97f1uiN8Lum9Ofxz4/Wm6TiCir2ThAfG8MdF/bltok5zN9ewO8X57N238nA8wu2F7BgewEX5HTnvksHccPYPo2mUoozCsoq+d7rm5i/vSDwWFJ8DD+5ZSy3X9jXwWQibaNi70QxMYarR2dx9egs1uw9wR8+2c3H244F1hTZfKiU7/x1Iz99fztfmTqAuy7uT0ZKorOho5C1lrfWH+Lxd7Y0erN7ZO9u/PauiRo6E9dRsYfJ5IE9mTywJ7uPlzNn6R5eW3uQKv+v+sfLqvjlvJ38dkEeN47P5p7PDdTNjsPkUHEFv1lfxfqCDY0e/+olA3l05kgNvYgrqdjDbHBmCj+5dSwPXzOCV1bt58XleykoqwKguq6eN9Yd4o11hxjfN5UvTx3A58f3Ubl0gpq6el5Ytpdff7yT00HXHOSkdeGJO8ZxydAMB9OJdIyK3SE9kxP4Ru5Q/nX6YN7fdIQXlu9lw4GGm2xvPFjCxtc+48fvbeW2SX350sX9Ga4hgZBYvPM4P35vK3kFje+O9aWL+/P960dpaQBxPf0EOywhLoZbJuZwy8Qc1u8/yUsr9vHeZ0cCMzJKK2t5YfleXli+l4n90/inyf24cXwfh1O707YjpTzxwXYW7jje6PG+KYb//vJUJg/s6VAykdBSsUeQif17MLF/D/79htHMXXOAV1btZ1/R6cDz6/cXs35/MT96dysTMw3xfQuZOjidWN0k+bz2FJ7i1/N28u5nhwlebTklMY5vXjGUoXX7VeriKSr2CNQzOYEHLh/CrOmDWZZfyF9W7ufjbceoqfO1UkVNHcsPw/I/rSSreyI3je/DzRNyGNOnu+bFB9l+tJSnFubz988OE7z4pjHwxQv78fC1w+nVLYlFiw44F1KkE6jYI1hMjGH6sEymD8ukqLyKN9cf4m9rDrLjWFngNcdKq/jjkj38cckeBqZ35YZx2Vw/NpvR2dFZ8vX1lk92Hef5ZXtZvPP4Oc9fNSqLh68Zzqjs7g6kEwkPFbtLpKckcv/0wdx36SC2HC7lv99ewfqiGIpONdyIe2/RaZ5amM9TC/Pp17MLM8f05qpRWVw4oIfnVyAsKq/ijXWHeGX1fnYHLQNwxmXDM3noymFcOKCHA+lEwkvF7jLGGC7ISeXLoxN5ZvplLM0r5O31h5i39RingqbtHThRETiTT+saz+XDM7l8uO/sP7ObNy6CqqypY+H2At7ecJj52xuGqs4wBq4ZncU3cocyrq+uC5DooWJ3sfjYGHJH9CJ3RC8qa+pYtOM47286woLtBZRX1QZeV3y6hrc3HObtDYcBGJXdnWlD0pk2JJ3JA3qS2tU9t3Mrraxh0Y7jzNt6jIVn/TvPSEmM4wuT+3LvtIEMSE92IKWIs1TsHpEUH8vMC3oz84LeVNXWsTyviHnbjrFgWwFHz1oyeNuRUrYdKWXO0j0YAyOyujFpQA8m9EtjQr80hmSmdGymTXk5PPkkPP00lxcVQXq676bhjzwSuCtVa1VU17HxYDGf5hexNK+QDQeKm70N4cT+adx5UT9uHNeHZM1Flyimn34PSoyLJXdkL3JH9sLeYtlyuJTFO4+zeOdx1u07SW1QMVoL24+Wsf1oGX9ZuR/wLXw1snd3RmV3Z1ivFIZlpTA4M4Xs7knEtFT45eUwdSrk50NlJQagsBCeeAJef913f9lmyr2koobdx8vZdqSMLYdL2HyohC2HSxvlPdvA9K7cNL4PN03ow9BeuoBLBFTsnndmTP6CnFS+kTuU8qpaVu89wfK8QlbtOcHmw6XnnAFX1tSz4UBxoythARJiY+jbswt9UrvQOzWJrO6JpCcn0jM5gdQu8XRNiGXQb58kIy+fmKqzbixSWUldXh77Hn2c9f/yHU6cqqawvIrDJZUcLq5gX9FpCsurWvVvGpuTylWjfIurjcruFpWzf0TOR8UeZVIS4wLj8gCnq2tZv784UOQbDxQH1q45W3VdPbuPn2py1skZa//wzLml7hdbVUXq83N4uNsVbco8rFcKkwf2YNqQDKYNSSddK2CKnJeKPcp1TYjjkqEZjRa9KiyvYtuRUrYfKSOvoJxdBWXsKzrdaGplc3pUlLXwfGmzzyXExTAoPZmhWSmM6dOdMX1SGd83lbSuCc1+jIicS8Uu58hISQxcGBWsvKqWAydOc7SkkiMllRwrraT4dDVFp6opq6zldHUtZcndST1V0uznLu+Wyk3j+9AzOYGMlASyU7uQnZZE37Su5PToouURREJAxS6tlpIYx6js7ue/avPIQ743SiubGI5JSqL7d77Fb+6a2HkhRQRvX44o4ffIIzBkCCSddcPnpCTf44884kwukSiiYpfQSknxTWmcPRsyM7HGQGamb/88Ux1FJHRU7BJ6KSnwox9BQQGLFyyAggLfvkpdJCxU7CIiHqNiFxHxGBW7iIjHqNhFRDxGxS4i4jEqdhERj1Gxi4h4jIpdRMRjVOwiIh6jYhcR8RhjbfO3HQv5FzOmDNgRti/YPhlAodMhWkE5Q8cNGcEdOd2QEdyRc4S1tl33ewz3sr07rLWTw/w128QYsybSM4JyhpIbMoI7crohI7gjpzFmTXs/VkMxIiIeo2IXEfGYcBf7s2H+eu3hhoygnKHkhozgjpxuyAjuyNnujGF981RERDqfhmJERDymU4vdGPMFY8wWY0y9MabZd6CNMXuNMZuMMRs68k5wJ2ecaYzZYYzJM8Y8Fs6M/q/f0xgzzxizy/93j2ZeV+c/jhuMMe+EKdt5j40xJtEY81f/8yuNMQPDkauJHC3lvNcYczzo+N3vQMbnjDEFxpjNzTxvjDG/8f8bPjPGTIrAjDOMMSVBx/E/HMjYzxiz0Bizzf/f90NNvCYSjmVrcrb9eFprO+0PMAoYASwCJp/ndXuBjM7M0pGMQCyQDwwGEoCNwOgw53wCeMy//Rjwi2ZeVx7mXC0eG+BB4Pf+7TuBvzrwfW5NznuB34U721kZLgMmAZubef564B+AAaYCKyMw4wzgPYePYzYwyb/dDdjZxPc7Eo5la3K2+Xh26hm7tXabtTaiL0hqZcYpQJ61dre1thp4Fbi589M1cjPwon/7ReCWMH/95rTm2ARnfw240hhjwpgRIuN72CJr7SfAifO85Gbgf63PCiDNGJMdnnQ+rcjoOGvtEWvtOv92GbANyDnrZZFwLFuTs80iZYzdAh8ZY9YaY2Y5HaYJOcCBoP2DhODgt1GWtfYI+H4YgF7NvC7JGLPGGLPCGBOO8m/NsQm8xlpbC5QA6WHI1mQGv+a+h7f7fy1/zRjTLzzR2iQSfhZb43PGmI3GmH8YY8Y4GcQ/9DcRWHnWUxF1LM+TE9p4PDt85akx5mOgdxNP/V9r7dut/DSXWGsPG2N6AfOMMdv9ZwUhEYKMTZ1dhnw60flytuHT9Pcfy8HAAmPMJmttfmgSNqk1xyYsx68FrcnwLvCKtbbKGPMAvt8yruj0ZG0TCceyJeuAAdbacmPM9cBbwDAnghhjUoDXgW9ba0vPfrqJD3HkWLaQs83Hs8PFbq29KgSf47D/7wJjzJv4fm0OWbGHIONBIPjsrS9wuIOf8xzny2mMOWaMybbWHvH/uljQzOc4cyx3G2MW4TsD6Mxib82xOfOag8aYOCCV8P8q32JOa21R0O4fgV+EIVdbheVnsSOCi8la+74x5mljTIa1Nqxrsxhj4vGV5cvW2jeaeElEHMuWcrbneDo+FGOMSTbGdDuzDVwDNPluu4NWA8OMMYOMMQn43gAMy4yTIO8A9/i37wHO+U3DGNPDGJPo384ALgG2dnKu1hyb4Ox3AAus/12hMGox51njqzfhG++MNO8AX/HP6JgKlJwZoosUxpjeZ95DMcZMwdczRef/qJBnMMAcYJu19lfNvMzxY9manO06np38ju+t+P6vWAUcAz70P94HeN+/PRjfDIWNwBZ8wyPhfFe6xYy24R30nfjOfsOa0f/104H5wC7/3z39j08G/uTfngZs8h/LTcB9Ycp2zrEB/hO4yb+dBPwNyANWAYPDffxamfNn/p/BjcBCYKQDGV8BjgA1/p/L+4AHgAf8zxvgKf+/YRPnmW3mYMZvBh3HFcA0BzJeim9Y5TNgg//P9RF4LFuTs83HU1eeioh4jONDMSIiEloqdhERj1Gxi4h4jIpdRMRjVOwiIh6jYhcR8RgVu4iIx6jYRUQ85v8D4xXrG9odGwEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import operator\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import base64\n",
    "from io import BytesIO\n",
    "\n",
    "\n",
    "n = 4\n",
    "\n",
    "def generate_points(amin=-3, amax=4, n=4):\n",
    "\n",
    "    A = [random.randrange(amin, amax) for i in range(n-1)]\n",
    "    non_zero = [i for i in range(amin, amax) if i !=0]\n",
    "    A_top = random.choice(non_zero)\n",
    "    A.append(A_top)\n",
    "    a = np.array(A, dtype=int)\n",
    "    x  = np.arange(-n//2 + 1, n//2 + 1)\n",
    "\n",
    "    exps = np.arange(n, dtype=int) \n",
    "\n",
    "    vandermonde = ((np.ones((n, n))*x).T**exps).astype(int)\n",
    "    a = a.reshape(-1, 1)\n",
    "    y = vandermonde.dot(a)\n",
    "    \n",
    "    points = (x, y)\n",
    "    \n",
    "    return points, list(a)\n",
    "\n",
    "def poly(x, coeff):\n",
    "    y = 0\n",
    "    for (i, a) in enumerate(coeff):\n",
    "        y += a*(x**i)\n",
    "    return(y)\n",
    "\n",
    "points, a = generate_points(n=n)  \n",
    "  \n",
    "#print(points)  \n",
    "#print(a)\n",
    "\n",
    "#plotting   \n",
    "\n",
    "x = np.linspace(-n//2 + 0.9 , n//2 +0.1, 100)\n",
    "y = poly(x, a)\n",
    "fig = plt.figure(figsize=(5,5))  \n",
    "xmin, xmax =  np.min(points[0])-0.5, np.max(points[0])+0.5\n",
    "plt.xlim(xmin, xmax)\n",
    "plt.grid()\n",
    "plt.plot(x, y, lw=3)\n",
    "plt.plot(points[0], points[1], \"ro\", ms=8)\n",
    "plt.tight_layout(pad=0)    \n",
    "tmp = BytesIO()\n",
    "fig.savefig(tmp, format=\"png\", dpi=200)\n",
    "enc = base64.b64encode(tmp.getvalue())\n",
    "tmp.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1,  0,  1])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 3\n",
    "x  = np.arange(-n//2 + 1, n//2 + 1)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-3//2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 4\n",
    "np.arange(n, dtype=int) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0 4\n"
     ]
    }
   ],
   "source": [
    "N = 4\n",
    "\n",
    "def generate_points(amin=-3, amax=4, n=N):\n",
    "  \n",
    "    A = [random.randrange(amin, amax) for i in range(n-1)]\n",
    "    non_zero = [n for n in range(amin, amax) if n!=0]\n",
    "    A_top = random.choice(non_zero)\n",
    "    A.append(A_top)\n",
    "    a = np.array(A, dtype=int)\n",
    "    x  = np.arange(-(n//2), n/2 + 1)\n",
    "\n",
    "    exps = np.arange(n, dtype=int) \n",
    "  \n",
    "    return float(exps[-1]), n\n",
    "\n",
    "\n",
    "a, n = generate_points()\n",
    "\n",
    "print(a, n)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 4\n",
    "\n",
    "a = [n for n in range(10) if n != 7 ]\n",
    "n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "int(\"-2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2, -1,  1,  2])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 4\n",
    "x  = np.arange(-n//2, n//2 + 1)\n",
    "x  = x[x != 0]\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}1 & -2 & 4 & -8 & 16 & -32\\\\1 & -1 & 1 & -1 & 1 & -1\\\\1 & 1 & 1 & 1 & 1 & 1\\\\1 & 2 & 4 & 8 & 16 & 32\\\\1 & 3 & 9 & 27 & 81 & 243\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1  -2  4  -8  16  -32⎤\n",
       "⎢                     ⎥\n",
       "⎢1  -1  1  -1  1   -1 ⎥\n",
       "⎢                     ⎥\n",
       "⎢1  1   1  1   1    1 ⎥\n",
       "⎢                     ⎥\n",
       "⎢1  2   4  8   16  32 ⎥\n",
       "⎢                     ⎥\n",
       "⎣1  3   9  27  81  243⎦"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy\n",
    "sympy.init_printing(use_latex='mathjax')\n",
    "A = sympy.Matrix([[1, k, (k)**2, (k)**3, (k)**4, (k)**5] for k in [-2, -1, 1, 2, 3] ])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}2\\\\-3\\\\-1\\\\-6\\\\-23\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡ 2 ⎤\n",
       "⎢   ⎥\n",
       "⎢-3 ⎥\n",
       "⎢   ⎥\n",
       "⎢-1 ⎥\n",
       "⎢   ⎥\n",
       "⎢-6 ⎥\n",
       "⎢   ⎥\n",
       "⎣-23⎦"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = sympy.Matrix([2, -3, -1, -6, -23])\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = A.col_insert(6, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}1 & -2 & 4 & -8 & 16 & -32 & 2\\\\1 & -1 & 1 & -1 & 1 & -1 & -2\\\\1 & 1 & 1 & 1 & 1 & 1 & -1\\\\1 & 2 & 4 & 8 & 16 & 32 & -6\\\\1 & 3 & 9 & 27 & 81 & 243 & -23\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1  -2  4  -8  16  -32   2 ⎤\n",
       "⎢                          ⎥\n",
       "⎢1  -1  1  -1  1   -1   -2 ⎥\n",
       "⎢                          ⎥\n",
       "⎢1  1   1  1   1    1   -1 ⎥\n",
       "⎢                          ⎥\n",
       "⎢1  2   4  8   16  32   -6 ⎥\n",
       "⎢                          ⎥\n",
       "⎣1  3   9  27  81  243  -23⎦"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left ( \\left[\\begin{matrix}1 & 0 & 0 & 0 & 0 & 12 & -2\\\\0 & 1 & 0 & 0 & 0 & -4 & 2\\\\0 & 0 & 1 & 0 & 0 & -15 & 0\\\\0 & 0 & 0 & 1 & 0 & 5 & -1\\\\0 & 0 & 0 & 0 & 1 & 3 & 0\\end{matrix}\\right], \\quad \\left ( 0, \\quad 1, \\quad 2, \\quad 3, \\quad 4\\right )\\right )$$"
      ],
      "text/plain": [
       "⎛⎡1  0  0  0  0  12   -2⎤                 ⎞\n",
       "⎜⎢                      ⎥                 ⎟\n",
       "⎜⎢0  1  0  0  0  -4   2 ⎥                 ⎟\n",
       "⎜⎢                      ⎥                 ⎟\n",
       "⎜⎢0  0  1  0  0  -15  0 ⎥, (0, 1, 2, 3, 4)⎟\n",
       "⎜⎢                      ⎥                 ⎟\n",
       "⎜⎢0  0  0  1  0   5   -1⎥                 ⎟\n",
       "⎜⎢                      ⎥                 ⎟\n",
       "⎝⎣0  0  0  0  1   3   0 ⎦                 ⎠"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B.rref()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left\\{\\left ( - 12 \\tau_{0} - 2, \\quad 4 \\tau_{0} + 2, \\quad 15 \\tau_{0}, \\quad - 5 \\tau_{0} - 1, \\quad - 3 \\tau_{0}, \\quad \\tau_{0}\\right )\\right\\}$$"
      ],
      "text/plain": [
       "{(-12⋅τ₀ - 2, 4⋅τ₀ + 2, 15⋅τ₀, -5⋅τ₀ - 1, -3⋅τ₀, τ₀)}"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.linsolve((A, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = sympy.Matrix([-2, 2, 0, 0, 1, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}-2\\\\2\\\\0\\\\0\\\\1\\\\0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡-2⎤\n",
       "⎢  ⎥\n",
       "⎢2 ⎥\n",
       "⎢  ⎥\n",
       "⎢0 ⎥\n",
       "⎢  ⎥\n",
       "⎢0 ⎥\n",
       "⎢  ⎥\n",
       "⎢1 ⎥\n",
       "⎢  ⎥\n",
       "⎣0 ⎦"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}10\\\\-3\\\\1\\\\18\\\\85\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡10⎤\n",
       "⎢  ⎥\n",
       "⎢-3⎥\n",
       "⎢  ⎥\n",
       "⎢1 ⎥\n",
       "⎢  ⎥\n",
       "⎢18⎥\n",
       "⎢  ⎥\n",
       "⎣85⎦"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A*c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\operatorname{Poly}{\\left( x^{5} - 3 x^{4} - 5 x^{3} + 15 x^{2} + 4 x - 12, x, domain=\\mathbb{Z} \\right)}$$"
      ],
      "text/plain": [
       "Poly(x**5 - 3*x**4 - 5*x**3 + 15*x**2 + 4*x - 12, x, domain='ZZ')"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = sympy.Symbol(\"x\")\n",
    "\n",
    "m = [-2, -1, 1, 2, 3]\n",
    "\n",
    "y = 1\n",
    "for i in m:\n",
    "    y*= (x-i)\n",
    "    \n",
    "z = sympy.Poly(y)\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "dd = z.as_dict()\n",
    "u = np.array([dd[i] for i in dd], dtype=int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.int64"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(u[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5, 5) (5,)\n",
      "[-12   4  15  -5  -3   1]\n",
      "[-1  2  2  1  0  0]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAF2CAYAAAB6XrNlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt81fWd5/HX5+ROuASScI8ECIaL3DQgYEEbtbUOVltqa207TnXLzowz62xnZ6ZTZ9wZW7vdnd122d1uO8zY1k7dsR2xRWnVqlFABUEFRC6BQBQCCUmAhNwv53z3j4SYAwm5neR3zu+8n48Hj+T7O79zfp8vgTc/vr/f7/s15xwiIuIfAa8LEBGRyFKwi4j4jIJdRMRnFOwiIj6jYBcR8RkFu4iIzyjYRUR8RsEuIuIzCnYREZ9J9LqA7rKyslxubq7XZXQJBoMkJCR4XUZE+Kkv4K/+qC/RK9r6884771Q757L72i+qgj03N5e3337b6zK6VFdXk5WV5XUZEeGnvoC/+qO+RK9o64+Zfdif/TQUIyLiMwp2ERGfUbCLiPiMgl1ExGcU7CIiPqNgFxHxGQW7iIjPKNhFRHxGwS4i4jMKdhERn1Gwi4j4jIJdRMRnFOwiIj6jYBcR8RkFu4iIzyjYRUR8RsEuIuIzEQl2M8sws6fN7LCZHTKzlWY2wcxeMrOjnV/HR+JYIiJyZZE6Y98AvOCcmwssBg4B3wBecc7NAV7pbIuIyDAbcrCb2VhgDfA4gHOu1TlXA9wJPNG52xPAXUM9loiI9C0SZ+yzgCrgJ2a2x8z+2czSgUnOuXKAzq8TI3AsERHpQ2KEPuNa4E+dc2+Z2QYGMOxiZuuB9QA5OTlUV1dHoKTIqK2t9bqEiPFTX8Bf/VFfoles9icSwV4GlDnn3upsP01HsJ8xsynOuXIzmwJU9vRm59xGYCNAQUGBy8rKikBJkRNt9QyFn/oC/uqP+hK9YrE/Qx6Kcc5VACfNLL9z083AQeBZ4L7ObfcBm4d6LBER6VskztgB/hR40sySgePAV+n4R+OXZvYAcAK4O0LHEhGRK4hIsDvn9gIFPbx0cyQ+X0RE+k9PnoqI+IyCXURkEMrO1rO3tJrW9qDXpVwmUmPsIiJx5YU9J/m3HcdJTgwwP2c8n1sxi2V50fG4js7YRUQGYU9pxzM3re0h9paepaGl3eOKPqJgFxEZoAtNrRyruBC2bUlupkfVXE7BLiIyQPs+OIvr1p41aSwZ6Sme1XMpBbuIyADtLQ2f+mTJzOg5WwcFu4jIgO0pPRvWXpobXdMOKNhFRAagsraJU+cautoJAWPhjAkeVnQ5BbuIyADs/SB8GGbutAzSkqPrznEFu4jIAOy9dBhmZnQNw4CCXUSk35xzXfevX7REwS4iErs+qKzjXH1LVzstOYF50zI8rKhnCnYRkX5653j42fqiGZkkJkRfjEZfRSIiUerdS4Zhrp0VfcMwoGAXEemXYMhRUh6+Buq1UTi+DprdUUSkXxICxs8fKuRg2XnePV5N6ZkL5GSN9rqsHinYRUT6KTkxgSW5WSyJsidNL6WhGBERn1Gwi4j4jIJdRMRnFOwiIn1wzvW9UxTRxVMRkStwzvFHG7czPXM0BbOzuG52Ntlj07wu64oU7CIiV/BhVT2llXWUVtax/VA5KUkJbPqLT5AUhU+cXhS9lYmIRIG3j1WFtRfPmBDVoQ4KdhGRK3rneHiwXzc726NK+k/BLiLSi5a2IPs/PBe2rUDBLiISuw6X19EWDHW1J2WkMW1CuocV9Y+CXUSkF++dCJ/067pZ2ZiZR9X0n4JdRKQX752sCWvHwjAMKNhFRHpUWdvE6ZrmrnbAjCW5mR5W1H8KdhGRHlx6m+O86RmkpyZ5VM3A6AElEZFOzjmKT9fw3LbDbC8+gwUSSAy2Y6EQM7JH45yLiTF2BbuICNAeDLFh02727S9l7dZN3L+3iIyGGmrSM3h5cSFb3Dq+19DEQ+uWReU6p90p2EUk7jnn2LBpN2eLtrPxx39DaltL12uZdef4wutPc+dbz/HoA4+xAfj63cuj+sw9Iv/smNkHZrbfzPaa2dud2yaY2UtmdrTz6/hIHEtEJNKKT9ewb38pj1wS6t2ltrXwyOMPs29/KcWna3vcJ1pE8v8TH3fOLXHOFXS2vwG84pybA7zS2RYRiTpbth/mjq2beg31i1LbWli77Rm2bD80QpUNznAOFN0JPNH5/RPAXcN4LBGRQdtRUkXh3qJ+7Vu4p4idJVV97+ihSI2xO+B3ZuaAf3TObQQmOefKAZxz5WY2sac3mtl6YD1ATk4O1dXVESpp6Gpro/u/WwPhp76Av/qjvnivMRQgo6Gm7x2BjIYaGkKBqMqqS0Uq2G9wzp3uDO+XzOxwf9/Y+Y/ARoCCggKXlRVdq39HWz1D4ae+gL/6o754a1QgRE16Bpl15/rctyY9g/RAKKr7GZGhGOfc6c6vlcCvgOXAGTObAtD5tTISxxIRibSVedkULSns175FSwtZkRfdUwsMOdjNLN3Mxlz8HvgE8D7wLHBf5273AZuHeiwRkeGwdvVcnrtxHc1JKVfcryk5lS1r1rF29bwRqmxwInHGPgl43cz2AbuA3zjnXgC+C9xqZkeBWzvbIiJRJ39qBosXzuQ/f/XbvYZ7U3Iq37r/2yxemEv+1HEjXOHADHmM3Tl3HFjcw/azwM1D/XwRkeFmZjy0bhlfPHKG3/+Pj/PZHb/mlr2vdj15WrS0kC1r1rF4YS4PrVsW1Q8ngZ48FREBoKKmkQsuAUZn8OSN9/Dzj99LMCGJ9ECIFXnZPLxmHvlTM7wus18U7CIiwFtHP7q/ozU5lWtnZfHnn5wd1Xe/9Ca6Z7IRERkhu46G37i3PK/HR29igoJdROJeQ3Mb+0+E38N+/RwFu4hIzHr3eDXBkOtqT89MZ2oMLFrdGwW7iMS9ty4Zhonls3VQsItInAuGHLtKLhlfV7CLiMSu4tM11Da2drVHpSRyTc4EDysaOgW7iMS1ncVnwtrLZmdH/dJ3fYnt6kVEhmjHkfBgX3H1JI8qiRwFu4jErdPnGjhRXd/VDpixLIbvX79IwS4icWvnJXfDXHPVeMakJXlUTeQo2EUkrk0Y/dFsjit9MAwDmitGROLYZ6+fyV3Lcykpr2XHkTOsyp/sdUkRoWAXkbgWMOPqqRlcHSMzN/aHhmJERHxGwS4i4jMKdhERn1Gwi0jcOVpeS0NLm9dlDBtdPBWRuBJyjkee2k1dUxtLZmayKn8yNy+cRkpSgtelRYzO2EUkrhSfquFcfQttwRC7S6rY+NJBonxt6gFTsItIXHnz0km/8iaSnOifs3VQsItInHmzuCKsvSrfH0+bdqdgF5G4caK6nrKzDV3txIDF9KLVvVGwi0jcePNw+Nn64txM0lNjf9KvSynYRSRuvHFJsK/0ydwwl1Kwi0hcOFPTyJHy2q624c/xdVCwi0icuPRsfX7OeDLHpHpUzfBSsItIXHj9kmD/2Fx/DsOAgl1E4sC5+mYOnjwftm2Vgl1EJHa9WXwG1609Z8o4JmeM8qye4aZgFxHfe/1Q/AzDgIJdRHyutT1I2dn6sG03+DzYNbujiPhacmICP/sPhRwqO8/rhyo4ebaenKzRXpc1rCIW7GaWALwNnHLOrTWzmcBTwATgXeArzrnWSB1PRKS/AmYsyJnAgpwJXpcyIiI5FPMQcKhb+78C33fOzQHOAw9E8FgiItKLiAS7mU0Hfg/45862AYXA0527PAHcFYljiYjIlUXqjP1/An8JhDrbmUCNc669s10GTIvQsURE5AqGPMZuZmuBSufcO2Z208XNPezqetiGma0H1gPk5ORQXV091JIipra2tu+dYoSf+gL+6o/6Mjxa2oIEQ45RKYOPuWjqz0BE4uLpDcCnzex2IBUYS8cZfIaZJXaetU8HTvf0ZufcRmAjQEFBgcvKyopASZETbfUMhZ/6Av7qj/oSeS/uPcn//u37XDcri9Xzp7Dy6kmDmqI3WvozEEMeinHO/bVzbrpzLhe4Byhyzn0JeBX4XOdu9wGbh3osEZH+2nawnLZgiJ1HK/mHzft4eudxr0saMcP5gNJfAV83sxI6xtwfH8ZjiYh0udDUyp7S8GHdNfOmeFTNyIvoA0rOudeA1zq/Pw4sj+Tni4j0x5uHKwiGPrqsNz0zndyJYzysaGRpSgER8Z3XDpSHtW9aMJWOu7Djg4JdRHzlfH0L+z4IH4a5cX78DMOAgl1EfGb7oXK6jcIwa9JYrsqOn2EYULCLiM+8diD8zuqbFsTX2Too2EXERyprmzhwyUpJN86f6lE13lGwi4hvbDsYftF03rQMJo/370pJvVGwi4hvbL1kGObGBfF3tg4KdhHxibKz9Rwp/2huFwPWxNndMBcp2EXEF4r2h5+tL8rNJHNMqkfVeEvBLiK+cKDsXFi78Jr4HIYBrXkqIj7x3S9dz8Gy8xTtP8WOI2f4WBzNDXMpBbuI+IJ1W9f0wU9dQyCOphC4lIZiRMR34jnUQcEuIuI7CnYREZ9RsItIzGpuC1LT0OJ1GVFHF09FJGZtP1jO97e8x7LZ2dyyaDrXXz2R5MQEr8vynM7YRSRmvby/jGDIsfNoJd/e9C7/ur3E65KigoJdRGJSZW0T+0rPhm27KY4fSupOwS4iMenV90/RbT0Nrp4yjhlxtqBGbxTsIhJznHO8/N6psG23LJrmUTXRR8EuIjGnpOICJ6rru9oJAYvbKXp7omAXkZjzu30nw9rL8iaSkZ7iUTXRR8EuIjGltT3Iq++HT9F7y0INw3SnYBeRmLLzSCV1TW1d7TFpSVx/9UQPK4o+CnYRiSmXDsPcvHCaHkq6hIJdRGJG9YVm3jlWFbbtE4une1RN9FKwi0jMePm9MkLdbl6fPWkssyeP866gKKVgF5GY4Jzjd/vKwrZ9YonO1nuiYBeRmHC+oYXu62ckBoyPX6O7YXqi2R1FJCZMGJ3KP//RjRw4eZ7n95zAORg3KtnrsqKSgl1EYoaZcc1VE7jmqgk45/p+Q5zSUIyIxCSL83VNr0TBLiLiMwp2ERGfGXKwm1mqme0ys31mdsDM/r5z+0wze8vMjprZL8xMVzlEZMDeOV7F+XqtazoQkbh42gIUOufqzSwJeN3Mnge+DnzfOfeUmf0IeAD4YQSOJyJxorktyGNPv0tzW5AVcybyqWuv4tpZ2SQENL5+JUM+Y3cdLk6MnNT5ywGFwNOd258A7hrqsUQkvmw9cJqGlnaCIccbxWf478/uI6S7YfoUkdsdzSwBeAfIA34AHANqnHPtnbuUAT0+SWBm64H1ADk5OVRXV0eipIiora31uoSI8VNfwF/9UV969+u3joW1b5iTSe35cxE9xpXE6s8mIsHunAsCS8wsA/gVMK+n3Xp570ZgI0BBQYHLysqKREkRE231DIWf+gL+6o/6crmj5bUcr2wI2/a5G/LJmpAekc/vr1j82UT0rhjnXA3wGrACyDCzi/9wTAdO9/Y+EZFL/fbdE2Ht62ZlMXWEQz1WReKumOzOM3XMLA24BTgEvAp8rnO3+4DNQz2WiMSHhpY2ivaHL1a99roZHlUTeyIxFDMFeKJznD0A/NI5t8XMDgJPmdm3gT3A4xE4lojEgaL9p2luC3a1M8ekaJWkARhysDvn3gOW9rD9OLB8qJ8vIvHFOcdv3vkwbNunll5FQkDPU/aXfqdEJKq8f/I8pZV1Xe2AwW1LczysKPYo2EUkqjy7+4Ow9oqrJ5E9Ns2bYmKUgl1EosbZumbeOFwRtu3Ty3K9KSaGKdhFJGr89t0TBLstapqTmc6S3EwPK4pNCnYRiQptwdBl967fsSxX864PgoJdRKJCKOS4e+UspowfBUBacgK3LNKapoOhpfFEJCqkJCXw2RWzuOv6mbxzrIrK2ibSU5K8LismKdhFJKoEzFiWp4eRhkJDMSIiPqNgFxHxGQW7iIjPKNhFxDOt7UEe/Kft/OKNEi40tnpdjm/o4qmIeObV909TUnGBkooLPLntKJ9elsu/u6WndXpkIHTGLiKecM6xaefxrnZLe4ig1jONCAW7iHji7WNVfFhV39UOGNypeWEiQsEuIp54utvZOsDqeVOYnDHKo2r8RcEuIiPuWEUte0vPhm373MpZHlXjPwp2ERlxm3aWhrUXzZjA1VMzPKrGfxTsIjKiKmubeO3A6bBt61bobD2SFOwiMqJ+9VZp2Jzr0zPTWT5Hc8NEkoJdREbMhcbWy+ZcX7diFgHNuR5RCnYRGTHPvv0hzW3BrvaE0Smac30YKNhFZEQ0t7azeVf4RdPPrphJcmKCRxX5l4JdREbEC3tPcqGpras9OjWR26+9ysOK/EvBLiIjorGlneTEjyLnjoJcrZA0TDQJmIiMiHtXz+H2a69i864PeHHfSe5anut1Sb6lYBeREZORnsJ9H8/nyzfOISGgAYPhot9ZERlxCvXhpd9dERGfUbCLiPiMgl1EhkXIOf7br/fyZnEFTgtojChdPBWRYfHG4Qpe2X+KV/afYubEMXxp9RxWz5/idVlxQWfsIhJxIef4+dajXe3SyjpeP1zhYUXxRcEuIhG3/WA5H1TVdbUN+NLqPO8KijNDDnYzyzGzV83skJkdMLOHOrdPMLOXzOxo59fxQy9XRKJdMOT4+bajYds+fs1Ursoe41FF8ScSZ+ztwJ875+YBK4AHzWw+8A3gFefcHOCVzraI+Ny2g6c5UR2+SPWX1szxsKL4M+SLp865cqC88/s6MzsETAPuBG7q3O0J4DXgr4Z6PBGJHs45ik/X8Ny2w+w8VkVjKEBisJ3kUJDW5FQAChdOY3rmaI8rjS8RvSvGzHKBpcBbwKTO0Mc5V25mWiJFxEfagyE2bNrNvv2l3LF1E/fvLSKjoYaa9AxeXlzIplV30ZQyis+vmu11qXEnYsFuZqOBTcCfOecuWD9XRDGz9cB6gJycHKqrqyNV0pDV1tZ6XULE+Kkv4K/+xGJfnHP80wsHqX99Fxt/8rektrV0vZZZd44vvP40d771HH9736M8+cIEvnbbfPqbCdEkFn82EKFgN7MkOkL9SefcM52bz5jZlM6z9SlAZU/vdc5tBDYCFBQUuKysrEiUFDHRVs9Q+Kkv4K/+xFpfDp86z+Ej5ZeFenepbS1864lHWD/5p5wtXMzcaRkjXGVkxNrPBiJzV4wBjwOHnHPf6/bSs8B9nd/fB2we6rFEJDps2X6YO7Zu6jXUL0pta2HttmfYsv3QCFUmEJm7Ym4AvgIUmtnezl+3A98FbjWzo8CtnW0R8YEdJVUU7i3q176Fe4rYWVI1zBVJd5G4K+Z1Op4/6MnNQ/18EYk+jS5ARkNNv/bNaKihwelZyJGk320RGbBRFqImvX9j5jXpGaRbaJgrku4U7CIyYCvzsilaUtivfYuWFrIiL3uYK5LuFOwiMmBrV8/luRvX0ZyUcsX9mpJT2bJmHWtXzxuhygQU7CIyCPlTM1i8cCaPPvBYr+HelJzKt+7/NosX5pI/ddwIVxjfNB+7iAyYmfHQumVsANZn/ZS1256hcM9HT54WLS1ky5p1LF6Yy0PrlsXkw0mxTMEuIgNW29jKuFHJfP3u5RTfkM+WBTP4w5Iv0xAKkB4IsSIvm4fXzCN/amw+lBTrFOwiMiCNLe384T9uY9708XztlnnMnZbB3HtWAlBdXR2TT2r6jcbYRWRAfvFGCefqW3jjcAVf++FWfvbaEa9Lkkso2EWk38rPN7JpZ2lXuy0YIqSFqqOOgl1E+u1HvztIW/Cjh40mjE7hCzdoWt5oo2AXkX7ZdbSSnUfOhG27v3Auacm6VBdtFOwi0qfW9iA//N2BsG3zp4/n5kXTPKpIrkTBLiJ9emZnKafPNXa1DXjwtgUEdH96VFKwi8gVVdQ08v+2Hw3bdvt1V5E3RU+TRisFu4j0yjnHD144QEv7RxdMx6Ql8Qc35XtYlfRFwS4ivXrjcAW7joavavnvbp7L2FHJHlUk/aFgF5EeNba088MXD4ZtW5Aznk8syfGoIukvBbuI9OiJ14qprmvuaicEjP9w+0JdMI0BCnYRuUww5Dh9vjFs27oVs8idOMajimQgFOwicpmEgPHoFwr4yzsXMzYtiUkZaXxpzRyvy5J+0iNjItIjM+PmRdO5bnY21ReaSU1K8Lok6ScFu4hcUUZ6ChnpV14CT6KLhmJERHxGwS4iAJw61+B1CRIhCnYR4Wh5LV/74Vb+YfNe6pvbvC5HhkjBLhLnWtuD/I9n9xEMOV5+7xT//kfbePd4tddlyRAo2EXi3BOvHaG0sq6rXV3XTHNbu4cVyVAp2EXi2N7SajbtOB627eaF01iVP9mjiiQSFOwicepCUyv/sHkf3VcszR6byh/ftsCzmiQyFOwiccg5x//6zfthc8EY8Jd3LWF0apJ3hUlEKNhF4tDze06y/VB52La7V81m0YxMjyqSSFKwi8SZDyrr+OGL4euX5k0ey+/fdLVHFUmkKdhF4khzazuPbXqX1m4rIqUmJfCNzywlKUFx4Bf6SYrEkR+8cIAT1fVh2/709mvIyRrtUUUyHBTsInHi+T0n+N2+srBtty6azi2LpntUkQyXiAS7mf3YzCrN7P1u2yaY2UtmdrTz6/hIHEtEBq6lLcjPXjsStm16ZjoPfkq3NvpRpM7Yfwrcdsm2bwCvOOfmAK90tkXEAylJCXz/D1aRN3lsV/tvP3cdacmauduPIhLszrltwLlLNt8JPNH5/RPAXZE4logMzuTxo/j+V1dx29Ic/uz3FmqZOx8bzn+uJznnygGcc+VmNrGnncxsPbAeICcnh+rq6Jl8qLa21usSIsZPfQF/9Wek+/KVFVMBhuXvmp9+LhC7/fH8/2HOuY3ARoCCggKXlZXlcUXhoq2eofBTX8Bf/VFfolcs9mc474o5Y2ZTADq/Vg7jsUSkm6PltZypafS6DPHIcAb7s8B9nd/fB2wexmOJSKeKmkb+5l938dCP36T4dI3X5YgHInW7478CO4B8MyszsweA7wK3mtlR4NbOtogMo4bmNh55ajc1Da2cb2jhL57YcdmcMOJ/ERljd859sZeXbo7E54tI39qDIR7b9C4fVn30ZGlLe4ijp2tZPW+Kh5XJSNOTpyI+4Jzjf27ZzzuXLGm3et5k/qAw36OqxCsKdhEf+Mmrxbz0Xvh0AflTM/iLO5cQMPOoKvGKgl0kxm3eVcov3jgWtm1yRhp/94XrSElK8Kgq8ZKCXSSGvbSvjP/74sGwbeNGJfOde69nwuhUj6oSrynYRWLU64fK+d5z+8K2pSQl8K0vLmNaZrpHVUk0ULCLxKC3j1XxX57ZQ6jbStSJAeNvP3ct+VMzvCtMooKCXSTGvHO8ir//5du0d0v1gME3PrOUZXk9TskkcUbBLhJjntv9YdjSdgD/8Y5FrJ6ve9Wlg4JdJMb89WeXct3s7K72H9+2gE8szvGwIok2CnaRGJOSlMDfff46CmZn88efnM+dy3K9LkmijOfT9orIwCUndtz9ooePpCc6YxeJYmVn63t9TaEuvVGwi0Qh5xw/KTrM+h9tY3eJljKQgVGwi0SZtmCI7z33Hk+9cYxgyPGtp9/l8CnNqy79pzF2kRHmnKP4dA3PbTvMzmNVNIYCjAqEWJmXTeHyOTz1+lH2ffjR2vAtbUG+s+ldHn/wJpISdC4mfVOwi4yg9mCIDZt2s29/KXds3cT9e4vIaKihJj2DV5YU8l/3foaG5FRISul6z9i0JL65bqlCXfpNwS4yQpxzbNi0m7NF29n4478hta2l67XMunN8fvvTfHrnc3zz9x/lyNQ82pJSmJyRxmP3Lmd65mgPK5dYo1MAkRFSfLqGfftLeeSSUO8uta2F7/zsEUa1NDFz4hg23H+DQl0GTMEuMkK2bD/MHVs39RrqF6W2tfDZHZuZNTaRjPSUK+4r0hMFu8gI2VFSReHeon7te8veInYdq+57R5EeKNhFRkijC5DR0L/bFjMaamhw+uspg6M/OSIjZJSFqEnv31zpNekZpFuo7x1FeqBgFxkhK/OyKVpS2K99i5YWsiIvu+8dRXqgYBcZIWtXz+W5G9fRnHTlC6JNyalsWbOOtavnjVBl4jcKdpFh4Jzjxb0nqTjf2LUtf2oGixfO5NEHHus13JuSU/nW/d9m8cJc8qeOG6lyxWf0gJJIhH1YVcf/ef593vvwHMvzsnn0nmWYGWbGQ+uWsQFYn/VT1m57hsI9Hz15WrS0kC1r1rF4YS4Pret4j8hgKNhFIqS+uY2fbzvK5l0fEHId65HuKqni9UMVXcvWJSYE+Prdyym+IZ8tC2bwhyVfpiEUID0QYkVeNg+vmafFqGXIFOwiQ9QeDPH8nhP8y9aj1Da2Xvb640WHWTV3MgmBjjNwM2PutAzm3rMSgOrqarKyska0ZvE3BbvIIDnneP1wBT8pKubUuYYe9ymYnc2ffOqarlAXGQkKdpEBcs6x80gl/7L1CMfOXOhxn0nj0vjarfP42NzJGiuXEadgF+mnYMjxxuEKfvFGCSUVPQd6SlICn181m7tXziIlKWGEKxTpoGAX6af//dv9PL/nZI+vBcy4bWkOX14zh8wxqSNcmUg4BbtIP92yaHqPwb5m/hS+smYOV2WP8aAqkcsp2EW6aQ+GuNDUyoTRl591X3PVBBbkjOfAyfMAfGzuZL68Zg4zJ40d6TJFrmjYg93MbgM2AAnAPzvnvjvcxxQZCOccJRUXeGX/KYr2n2LOlHE8du/yHvf94sfy2HqwnM+vnKUzdIlawxrsZpYA/AC4FSgDdpvZs865g8N5XJG+OOcoraxj+8Fyth4sD7td8Z1jVVTUNDI5Y9Rl71uWN5FleRNHslSRARvuM/blQIlz7jiAmT0F3Ako2GXEtQVDHDhxjrdKKtlRfIbybvO4dOeA5989wVcL545sgSIRMtzBPg3ofrWpDLh+mI8pAnSclZ+srmfPB2fZc7yavR9U09Qa7PN96SmJuvdcYtpwB3tPfztc2A5m64H1ADk5OVRXR89yYLW1tV6XEDF+6gv0rz/vlJ7ney8c6dfnJQSMRTnjuOHqLK7LHU9yYmDE/iz66Wfjp75A7PYilnmOAAALwElEQVRnuIO9DMjp1p4OnO6+g3NuI7ARoKCgwEXbnBnRVs9Q+KkvbcEQte1JlFRcYMr4USyakXnZPtenjYErBHtSQoCls7L42NzJrMqfzJi0pOEs+Yr89LPxU18gNvsz3MG+G5hjZjOBU8A9wL3DfEzxkQtNrZw+18ips/WUnW3g5Nl6Pqyqp+xsPaHO//vdunh6j8GekZ7C9AnplHW7MDphdArXzc7m+ryJFORlk5asO37Ff4b1T7Vzrt3M/gR4kY7bHX/snDswnMeUDs45ik/X8Ny2w+w8VkVjKMCoQIiVedmsXT2P/KnjomocORgKsXn3h1TVNnGmtokzNY1U1DRS39ze53uP9fJ4P0BBXjaTz45iSW4mS2dmMWvyWAJR1G+R4TDspyvOud8Cvx3u48hH2oMhNmzazb79pdyxdRP37+22mMOSQr5zYB2LF87koXXLSEwY+iJazjla20PUN7d1/apr6vzV3EZdUyu1ja1caGxl1qSx3Lt6zmWfETDjp68W09LW98XNS31YVUdre5DkxMvnZvmjTy4YVJ9EYpn+H+ozzjk2bNrN2aLtbPzx35Da1tL1WmbdOe7e/jR37HyORx94jA3A1+9eTlswRF1TGw3NbTS0tFPf3EZDcztTJozqcdGHptZ2/tMTO2hoaaexpZ2G5jbaQ+6y/XpS39zOvasv325mZI9JDRs26cukcWnMnDiG2ZPH0RYM9RjsIvFIwe4zxadr2Le/9LJQ7y61rYVHHn+Y+8Y9ztqSs7QHew7lTy+b0WOwJycGep3dsC81DT3XBJA19vJgT0oIMGX8KKaOH0VO1uiuX6MDLVw1dfKgahDxOwV7DOq4oNhA+flGztR0jEnfung686ePZ8v2w9yxdVOvoX5RalsLn3nz1zx50xehl4WVG3oZ304IBEhJShjUsElNw+UrDF10y6LpLJmZxaRxaUwcl8akjDQyx6T2OCYeTbfFikQbBXuUcs5RXdfMB5V1lFbWcbK6vuPX2Qbqm9su2z934hjmTx/PjpIqvrq3qF/HuGXvq/z8473fpNTQw3EuSk9JDAv2pIQAo1ISGZOaxOi0JNJTkxiblsSYtCTGpCYzLj2ZcWnJZKQn9/qZty6e3q+6ReTKFOxRpKS8ltcOnOZoRS3HKi5Q19R7sF6qqrYJgEYXIKOhpl/vyWiooT0hkYSAMTo1idGpSaSnJpKeksTo1MQrLqr8nXuXk5TYEeajU5M0vi0SRRTsUaTsbAP/tuP4oN5bXdcMwCgLUZOeQWbduT7fU5OeweiA49+++akB3/qoqWpFopeCfYRUnG9k97Eq9pRWs27FTBbkTLhsn7nTej9D7i45seOC4pTx6UwZP4pJ49LImzIOgJV52RQtKeTu7U/3+TlFSwtZkZcdVfezi8jQKdiHSTDkOFR2njeLK9h1tJKTZz+622NG1ugeg31SRhoZ6cldFxjTkhPIzR5D7sQxzMge03FHSGY62ePSen3IZu3quXznwDru2PncFS+gNiWnsmXNOr65et4Qeyoi0UbBHkHBkOP9E+d47cBpdhSf4Xwvt/btKa3mvo/nX7bdzLi/cC4piQnMmTKOKRNGDfgpyfypGSxeOJNHH3iMRx5/uMdwb0pO5Vv3f5vFC3PJnzpuQJ8vItFPwR4Bx89c4KV9Zbx24DTn6q98myF03Gve0NxGeurlk059cklOD+/oPzPjoXXL2ACsz/opa7c9Q+Gebk+eLi1ky5p1LF6Yy0PrlmkYRsSHFOyD1BYM8cKek7yw50S/HtZJCBgLcsZTMDuba2dlk5YyfL/1iQkBvn73copvyGfLghn8YcmXaQgFSA+EWJGXzcNr5l3xjhcRiW0K9m4GMnFWQsB4escxKmqaev28USmJLM+byKr8SRTMzu7xDH24mBlzp2Uw956VQMcDPbE4/aiIDJyCvdNAJ84KmHH7tTP4cdHhsM9JS05gVf5kblwwhaUzs3R/t4iMOAU7g5s4y8z45JLp/Oy1YoIhx9JZWXxi8XRW5k8mNUlhLiLeUbAzsImz1mf9lOIb8pk7LYOM9BT+6jNLmTstg4nj0ka4ahGRng19Mm4fGMjEWWu3PcOW7Ye6tq2ZP0WhLiJRRcEO7CiporCfE2cV7iliZ0nVMFckIjJ4cR/sFTWNNIQGNnFWg4v73zYRiWJxO8be0hbk3948xi/ePEZisG1AE2elW2gEKhQRGZy4PPV853gV//4ft/Ev247S2h7CQiFeXlzYr/denDhLRCRaxVWwN7e2899+vZdvPrmL8vONXdtbk1PZtOoumntZSeiiixNnrdXEWSISxeIq2JOTEsIC/aJxo5KYMWsKf//AY72GuybOEpFYEVdj7AEzHvq9hTz4T9tpDzkMuGPZDO67KZ/UpAQ2bErXxFkiEvPiKtihY23Qz6+azY4jZ/iztQuZO21812uaOEtE/CDugh3gi6vz+NKaOSQmhI9EaeIsEfGDuAx2TcwlIn4WVxdPRUTigYJdRMRnFOwiIj6jYBcR8RkFu4iIzyjYRUR8RsEuIuIzCnYREZ9RsIuI+MyQgt3M7jazA2YWMrOCS177azMrMbNiM/vk0MoUEZH+GuqUAu8DnwX+sftGM5sP3AMsAKYCL5vZ1c654BCPJyIifRjSGbtz7pBzrriHl+4EnnLOtTjnSoESYPlQjiUiIv0zXJOATQN2dmuXdW67jJmtB9YD5OTkUF1dPUwlDVxtba3XJUSMn/oC/uqP+hK9YrU/fQa7mb0MTO7hpYedc5t7e1sP21xPOzrnNgIbAQoKCly0TZMbbfUMhZ/6Av7qj/oSvWKxP30Gu3PulkF8bhmQ0609HTg9iM8REZEBGq7bHZ8F7jGzFDObCcwBdg3TsUREpJuh3u74GTMrA1YCvzGzFwGccweAXwIHgReAB3VHjIjIyBjSxVPn3K+AX/Xy2mPAY0P5fBERGTg9eSoi4jMKdhERn1Gwi4j4jIJdRMRnFOwiIj6jYBcR8RkFu4iIz5hzPU7h4gkzqwI+9LqObrKA6JmVbGj81BfwV3/Ul+gVbf2Z4ZzL7munqAr2aGNmbzvnCvreM/r5qS/gr/6oL9ErVvujoRgREZ9RsIuI+IyC/co2el1ABPmpL+Cv/qgv0Ssm+6MxdhERn9EZu4iIzyjY+2Bm/2Bmh83sPTP7lZlleF3TYJnZ3WZ2wMxCZhZzV/oBzOw2Mys2sxIz+4bX9QyFmf3YzCrN7H2vaxkqM8sxs1fN7FDnn7GHvK5psMws1cx2mdm+zr78vdc1DZSCvW8vAdc45xYBR4C/9rieoXgf+CywzetCBsPMEoAfAJ8C5gNfNLP53lY1JD8FbvO6iAhpB/7cOTcPWAE8GMM/mxag0Dm3GFgC3GZmKzyuaUAU7H1wzv3OOdfe2dxJx/qtMck5d8g5V+x1HUOwHChxzh13zrUCTwF3elzToDnntgHnvK4jEpxz5c65dzu/rwMOAdO8rWpwXIf6zmZS56+YuhipYB+Y+4HnvS4ijk0DTnZrlxGj4eFnZpYLLAXe8raSwTOzBDPbC1QCLznnYqovQ1oazy/M7GVgcg8vPeyc29y5z8N0/HfzyZGsbaD605cYZj1si6kzKb8zs9HAJuDPnHMXvK5nsDrXaF7SeU3tV2Z2jXMuZq6FKNgB59wtV3rdzO4D1gI3uyi/P7SvvsS4MiCnW3s6cNqjWuQSZpZER6g/6Zx7xut6IsE5V2Nmr9FxLSRmgl1DMX0ws9uAvwI+7Zxr9LqeOLcbmGNmM80sGbgHeNbjmgQwMwMeBw45577ndT1DYWbZF+9+M7M04BbgsLdVDYyCvW//BxgDvGRme83sR14XNFhm9hkzKwNWAr8xsxe9rmkgOi9i/wnwIh0X537pnDvgbVWDZ2b/CuwA8s2szMwe8LqmIbgB+ApQ2Pn3ZK+Z3e51UYM0BXjVzN6j42TiJefcFo9rGhA9eSoi4jM6YxcR8RkFu4iIzyjYRUR8RsEuIuIzCnYREZ9RsIuI+IyCXUTEZxTsIiI+8/8BUlxfwNhznH4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import operator\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import base64\n",
    "from io import BytesIO\n",
    "import sympy\n",
    "\n",
    "n = 5\n",
    "m = 4\n",
    "\n",
    "def generate_points(amin=-2, amax=3, n=4):\n",
    "  \n",
    "    non_zero = [i for i in range(amin, amax) if i!=0]\n",
    "    A = [random.randrange(amin, amax) for i in range(m-1)]\n",
    "    top_coeff = random.choice([i for i in range(amin, amax) if i!=0])\n",
    "    A.extend([top_coeff] + [0]*(n-m))\n",
    "    a = np.array(A, dtype=int)\n",
    "    x  = np.arange(-n//2+1, n//2 + 2)\n",
    "    x = x[x != 0]\n",
    "    exps = np.arange(n, dtype=int) \n",
    "    vandermonde = ((np.ones((n, n))*x).T**exps).astype(int)\n",
    "    print(vandermonde.shape, a.shape)\n",
    "    y = vandermonde.dot(a)\n",
    "    points = (x, y)\n",
    "  \n",
    "    t = sympy.Symbol(\"t\")\n",
    "    b_poly = 1\n",
    "    for k in x:\n",
    "       b_poly *= (t-k)\n",
    "    b_dict = sympy.Poly(b_poly).as_dict()\n",
    "    b_coeff = np.array([b_dict[i] for i in b_dict], dtype=int)\n",
    "    print(b_coeff)\n",
    "    az = np.array(A + [0])\n",
    "    print(az)\n",
    "    b =  az + b_coeff\n",
    "    \n",
    "    return points, list(a), list(b)\n",
    "\n",
    "def poly(x, coeff):\n",
    "    y = np.zeros_like(x)\n",
    "    for (i, a) in enumerate(coeff):\n",
    "        y += a*(x**i)\n",
    "    return(y)\n",
    "\n",
    "points, a, b = generate_points(n=n)  \n",
    "  \n",
    "#print(points)  \n",
    "#print(a)\n",
    "#print(b)\n",
    "\n",
    "#plotting   \n",
    "x = np.linspace(-n//2 + 0.8 , n//2 +1.2, 100)\n",
    "y = poly(x, a)\n",
    "fig = plt.figure(figsize=(5,5))  \n",
    "xmin, xmax =  np.min(x) - 0.5, np.max(x) + 0.5\n",
    "ymin, ymax =  np.min(y) , np.max(y)\n",
    "ymin, ymax =  ymin - 0.1*(ymax-ymin), ymax + 0.1*(ymax-ymin)\n",
    "plt.xlim(xmin, xmax)\n",
    "plt.ylim(ymin, ymax)\n",
    "plt.grid(ls=\"-\", lw=0.5, alpha=0.5)\n",
    "plt.plot(x, y, '--', lw=4, color = \"steelblue\", alpha = 1, zorder = 10)\n",
    "\n",
    "plt.plot(points[0], points[1], \"o\", fc='r', ec=\"steelblue\",  ms=12, zorder = 20)\n",
    "plt.tight_layout(pad=0)    \n",
    "tmp = BytesIO()\n",
    "fig.savefig(tmp, format=\"png\", dpi=200)\n",
    "enc = base64.b64encode(tmp.getvalue())\n",
    "tmp.close()\n",
    "\n",
    "\n",
    "xx = [int(t) for t in points[0]]\n",
    "yy = [int(t) for t in points[1]] \n",
    "\n",
    "s = \"\"\n",
    "for i in range(len(xx)):\n",
    "  s += \"({}, {}), \\hskip 2mm\\n\".format(xx[i], yy[i])\n",
    "s = \"\\[ \" + s + \"\\]\" \n",
    "\n",
    "def format_num(n):\n",
    "  if n == 0:\n",
    "      return str(0)\n",
    "  elif n == 1:\n",
    "      return \"+\"\n",
    "  elif n == -1:\n",
    "      return \"-\"\n",
    "  elif n > 0:\n",
    "     return \"+\" + str(n)\n",
    "  else:\n",
    "      return str(n)\n",
    " \n",
    "ansb = \"{} \".format(a[0]) if a[0] != 0 else \"\"\n",
    "for i in range(1, 4):\n",
    "  fa = format_num(a[i])\n",
    "  if fa == \"0\":\n",
    "     continue\n",
    "  ansb += fa + (\"*\" if fa not in \"+-\" else \"\") + \"x^{} \".format(i)\n",
    "  \n",
    "ansa = \"{} \".format(b[0]) if b[0] != 0 else \"\"\n",
    "for i in range(1, 5):\n",
    "  fb = format_num(b[i])\n",
    "  if fb == \"0\":\n",
    "     continue\n",
    "  ansa += fb + (\"*\" if fb not in \"+-\" else \"\") + \"x^{} \".format(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}