stefan2904/simple_ecc_vis2.ipynb

## simple_ecc_vis2.ipynb
{
 "metadata": {
  "name": "",
  "signature": "sha256:63323a1853ec0ef49b5d7f7445c21469d7383ab2e504e8f7e6fdaa45fe2d8820"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Simple Elliptic Curve Visualization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Curve Equation:\n",
      " y^2 = x^3 + a*x + b"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Curve Parameter:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = 2\n",
      "b = 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Field Order:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# pick a nice number from https://en.wikipedia.org/wiki/List_of_prime_numbers\n",
      "p = 3571 # 941"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Setup:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from math import sqrt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Draw Elliptic Curve:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y, x = np.ogrid[-5:5:100j, -5:5:100j]\n",
      "plt.contour(x.ravel(), y.ravel(), pow(y, 2) - pow(x, 3) - x * a - b, [0])\n",
      "plt.grid()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD7CAYAAABKfn7LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHpFJREFUeJzt3Xd8VGW+BvDnl4QqTUGliDeuCFLUQCwgIEGpUap0JA4l\nAdELuHrXBWVxxYLreukoSDAgSxEIgSDEggQQpAkRlWK7KFVElN6See8fAQ8ihJCZOe87532+nw+f\nz04ymfnts2d/HJ45mRGlFIiIKHxE6B6AiIiuDBc3EVGY4eImIgozXNxERGGGi5uIKMxwcRMRhZmo\nUD+BiPB6QyKiAlBKycW+7soZt1JK+59hw4Zpn8GUP8yCWYQ6iwkT1qFu3cnIzs7R/t9LdxYF/ZOX\nkJ9xm2LHjh26RzAGs3AwC0ewsti16zD+8Y9MLF/uQ2RkeLaxph8X4ZkqERlJKYV+/RbhiSfuQo0a\n1+oex7OsOeP2+Xy6RzAGs3AwC0cwsnj77Szs3n0EqakNAx9II9OPC7lclxLwE4ioUD8HEen344+H\nEBs7CR9/nIDbbrte9zhhT0SgdL44aYLMzEzdIxiDWTiYhSOQLJRSSExMx6BB93hiaZt+XFizuIko\ndN5+OwsHDhzH3/5WX/coVmBVQkQB2bkztyL56KME3H57+J9tm4JVCRGFhN+v0LPnAgwaVJdL20XW\nLG7TOys3MQsHs3AUJIuxY9fixIlsz1Ukph8X1lwOSETBtX37AQwfvgJr1vRBVJQ154BGYMdNRFcs\nO9uPBg2mICHhDvTvf5fucTyJHTcRBdW//rUKV11VGP363al7FCtZs7hN76zcxCwczMKR3yzWr9+N\n0aPXIiWlDSIiLnpCGPZMPy6sWdxEFLijR0+je/dUjB8fj8qVS+sex1rsuIko3/r2TcepUzlISWmr\nexTPy6vj5lUlRJQvixZ9jQ8++B6ff95P9yjWC0pVIiKRIrJJRNKD8XihYHpn5SZm4WAWjryy2Lfv\nKJKS0pGS0galShVxbyhNTD8ugtVxDwSwBQA7ESKP8fsVevSYjz596qBRo2jd4xCC0HGLyA0AUgC8\nBOCvSqlWF3yfHTdRGPvXv1YhPf1rLFv2KH/RxkWh7rhHAvgfAKWC8FhEZJANG/bg3/9ejQ0bkri0\nDRLQ4haRhwDsV0ptEpG4S93P5/MhOjoaAFCmTBnExMQgLi737ue6pFDfPvc1t57P5NtZWVkYNGiQ\nMfPovD1q1Cgtx6OJty/8/8qxY6fRtu0IPPZYbdx4Y2nt87l5+9zX3M4/JSUFAH7fl5cSUFUiIi8D\n6AEgG0BR5J51z1NKJZx3HyOqkszMzN/Dsh2zcDALx/lZKKWQkJCGqKgIvP12G72DaWDCcZFXVRK0\n67hFpBGAp9lxE4W/yZM3YvTotVi7tg+KFy+kexwruXkdNzc0UZj78sv9GDx4KVas8HFpGyporzYo\npZYrpVoH6/GC7fzuynbMwsEsHJmZmTh+/Aw6dZqD115riurVr9U9kjamHxd8mZiIfjdgwBLUqVMB\nPl+M7lEoD3yvEiICALz11mcYNWot1qzpjZIlvf/bkabje5UQUZ42btyLIUM+xief9OTSDgPWVCWm\nd1ZuYhYOZgH89ttJdOw4B/37l0O1auV0j2ME048LnnETWUwphV69FiA+vgoaNy6uexzKJ3bcRBZ7\n+eWVWLBgO1as8KFIEZ7HmYQdNxH9yeLF32D8+PVYvz6RSzvMsOO2ELNw2JrFd98dRM+eC/Duux1Q\nsWJJAPZmcTGmZ2HN4iaiXMePn0H79u9i6ND7UL/+jbrHoQJgx01kEaUUunadh0KFIjFtWluIePNT\n2r2AHTcRAQBeeGE5fvjhEJYte5RLO4xZU5WY3lm5iVk4bMoiNXUrkpM3Yf78ziha9M/nbDZlcTmm\nZ8EzbiILfPHFT+jbdxGWLOmO8uVL6B6HAsSOm8jjfvrpKOrWTcaLLzZG9+636x6H8imvjtuaqoTI\nRsePn0GrVjORkHA7l7aHWLO4Te+s3MQsHF7Owu9XSEiYj2rVyuH55+Mue38vZ3GlTM+CHTeRRw0d\n+jH27TuKpUsTeAWJx7DjJvKg5OSNeOmllVi7tg+uvfYq3eNQAfA6biKLLFr0NZ57bhlWrPBxaXsU\nO24LMQuH17JYt243evVagAULuuCWW8pe0c96LYtAmJ6FNYubyOu+//5XtG07C1OmtMHdd1fSPQ6F\nEDtuIg/Ys+cIGjZ8G089VQ/9+9+lexwKAl7HTeRhBw+eQLNm76B379pc2pawZnGb3lm5iVk4wj2L\nY8dO46GHZqBFiyoYPLhBQI8V7lkEk+lZWLO4ibzm9OkcdOw4B1WrlsVrrzXltdoWYcdNFIbOnMld\n2iKCd9/tgEKFInWPREHG67iJPCQ7249u3VKhFLi0LWVNVWJ6Z+UmZuEItyz8foU+fRbi8OFTePfd\nDihcOHhLO9yyCCXTs+AZN1GYyMnxIzExHd9//yuWLOnOT2a3GDtuojCQk+OHz7cAe/YcwcKFXXDV\nVYV1j0Qhxo6bKIxlZ/uRkDAfP/98HOnpXVG8eCHdI5Fm7LgtxCwcpmdx+nQOHnkkFb/8cgILF3YJ\n6dI2PQs3mZ6FNYubKNwcO3YarVvPxIkT2UhL64xixXimTbnYcRMZ6NChk4iPn4GqVcvirbdaISqK\n51i24XuVEIWRn346igcemIY6dcojObk1lzb9ScBHhIhUFpFlIvKViHwpIgOCMViwmd5ZuYlZOEzL\nYtu2A6hXLxkPPVQVY8a0RESEe7/GbloWOpmeRTCuKjkD4EmlVJaIlADwmYh8qJTaGoTHJrLGp5/u\nRLt2szFiRBP4fDG6xyGDBb3jFpE0AGOVUkvP3mbHTXQZ8+dvRVLSIkyb1hYtW96iexwygGvXcYtI\nNIDaANYG83GJvEophVde+QQTJqzHkiXdceedFXWPRGEgaK96nK1J5gIYqJQ6GqzHDRbTOys3MQuH\nzixOnsxGjx7zkZq6FWvX9tG+tHlcOEzPIihn3CJSCMA8ANOVUmkXft/n8yE6OhoAUKZMGcTExCAu\nLg6AE1Cob5/j1vOZfDsrK8uoeXTezsrK0vL81arFomPHOYiM/AEvvtgQlSqVMiIP3s69fY6bz5+Z\nmYmUlBQA+H1fXkrAHbfkvnv7VAC/KKWevMj32XETnWflyh/Qpcs8JCXVwdChjVy9coTCR14ddzAW\ndwMAKwBsBnDuwQYrpTLOfp+Lmwi5ffYbb2zAP/+5HFOntkWLFlV0j0QGC+kv4CilPlFKRSilYpRS\ntc/+yQj0cYPtwn8C2YxZONzK4vDhU+jePRUTJqzH6tW9jFzaPC4cpmfBX8kiCrHPPtuD2rUnonTp\nIli3LhE333yN7pEozPG9SohC5Fw18vzzmRg/Ph4dO9bUPRKFEb4fN5HLdu8+jJ49F+DQoVP45JNe\nqFq1rO6RyEOsqUpM76zcxCwcocgiLW0b6tSZhIYNb8SqVeGztHlcOEzPgmfcREGyf/8xDByYgXXr\ndiMtrTPq1auseyTyKHbcREEwf/5WPPbYe+jR43b885+N+fFiFDB23EQhsmvXYfz1r+9j06Z9SE3t\njHvv5Vk2hR47bgsxC0dBs/D7FSZMWI+YmDdx663lsHlzv7Bf2jwuHKZnwTNuoiu0adNe/Pd/L4Hf\nr7ByZU9Ur36t7pHIMuy4ifLp0KGTGDJkKebO3Yrhwxujd+/aiIy05h+t5DJ+5iRRAJRS+M9/NqNm\nzQk4c8aPrVsfR1JSLJc2aWPNkWd6Z+UmZuG4XBYbN+5F/fpTMHLkGsye3QGTJrXCNdcUc2c4l/G4\ncJieBTtuoov49dcTeP75TMye/RVefvkB+HwxfPtVMgY7bqLz5OT4kZy8CUOHLkP79rdi+PD7Ua5c\ncd1jkYV4HTdRPmRl7cNjj72HiAhBRkZ31K5dQfdIRBfFjttCzMKRmZmJw4dPYdCgDDRr9g569ozB\nypU9rVzaPC4cpmdhzeImupg1a3ahZs0JOHLkFL76qj+SkmLZZZPx2HGTlQ4ePIEnn3wfK1f+gMmT\nW+P++2/SPRLRH/A6bqLzvP/+t7jttjdQunQRbN78GJc2hR1rFrfpnZWbbM0iO9uPwYM/Qu/eCzF9\nejuMGdMSGzas1j2WMWw9Li7G9Cx4VQlZYc+eI+jceS6KFy+EjRv74rrrrtI9ElGBseMmz/v8831o\n1WomEhPr4Nln7+OLjxQWeB03WeuDD77DI4+kYty4eHTqxA/rJW9gx20hW7JIT9+OHj3mIzW18yWX\nti1Z5AezcJieBc+4yZMWL/4GffqkY9Girrjrrkq6xyEKKnbc5DkbN+5F8+bTkZ7eFXXr3qB7HKIC\n4XXcZI2ffz6G9u1n4803H+TSJs+yZnGb3lm5yatZ+P0KXbrMQ/fut+Hhh2vk62e8mkVBMAuH6VlY\ns7jJ+8aPX4eTJ7PxwguNdY9CFFLsuMkTdu8+jDvueBOrVvVCtWrldI9DFLC8Om4ubvKEJ55YjGLF\novDaa810j0IUFHxxEuZ3Vm7yWhZ79x7BjBlf4Omn773in/VaFoFgFg7Ts7BmcZN3TZz4Gbp1uw3X\nX19C9yhErmBVQmFNKYWbbx6DOXM6Ija2ou5xiIKGVQl51rp1u1GkSBTq1LHvo8bIXtYsbtM7Kzd5\nKYvMzB1o3vxmiBTsHf+8lEWgmIXD9CwCXtwi0kJEtonINyLyTDCGIsqvlSt/RMOGN+oeg8hVAXXc\nIhIJYDuAJgB2A1gPoKtSaut592HHTSFTpcoYvPdeN167TZ4Tyo77bgDfKqV2KKXOAJgFoE2Aj0mU\nL36/wq5dh3HjjaV1j0LkqkAXdyUAO8+7vevs14xjemflJq9kceLEGURECIoVK1Tgx/BKFsHALBym\nZxHo+3HnqwPx+XyIjo4GAJQpUwYxMTGIi4sD4AQU6tvnuPV8Jt/Oysoyap5Abvv9/4fMzMwC/3xW\nVpZR/31424zb57j5/JmZmUhJSQGA3/flpQTacdcF8LxSqsXZ24MB+JVSr553H3bcFBKnT+egRImX\ncfLkc/wcSfKcUHbcGwDcIiLRIlIYQGcACwN8TKJ8KVw4EtdcUwz79h3VPQqRqwJa3EqpbABPAHgf\nwBYAs8+/osQkF/4TyGZeyuKmm67Gd98dLPDPeymLQDELh+lZBHwdt1JqiVKqmlKqilLqlWAMRZRf\n99xTCatX77z8HYk8hO9VQmEtNXUrJk36DBkZj+gehSio+F4l5FlNmvwFa9bswoEDx3WPQuQaaxa3\n6Z2Vm7yURalSRRAffwtmzfqyQD/vpSwCxSwcpmdhzeIm70pMrIOxY9chJ8evexQiV7DjprCnlEL9\n+lMwcOA96Ny5lu5xiIKCHTd5mohg2LBGGDp0GU6dytY9DlHIWbO4Te+s3OTFLJo3r4KaNa/D8OEr\nrujnvJhFQTELh+lZWLO4yfsmTIjHW29txKpVP+oehSik2HGTpyxe/A0SE9Oxfn0iKlYsqXscogJj\nx03WiI+/Bf3734n27Wfj+PEzuschCglrFrfpnZWbvJ7FkCENUbVqWbRvPxsnT+b9YqXXs7gSzMJh\nehbWLG6yh4hgypQ2KF26KNq1m40TJ3jmTd7Cjps8KzvbD58vDd9//ysWLuyKcuWK6x6JKN/YcZOV\noqIiMG1aOzRq9F+oVy8ZX3/9i+6RiILCmsVtemflJpuyiIgQvPJKEzzzTH00aDAFc+du+cP3bcri\ncpiFw/QsrFncZLc+fepg8eLueOaZj9C//3uXfdGSyGTsuMkqhw6dRGJiOrZtO4CpU9uidu0Kukci\nuqi8Om4ubrKOUgrTp2/G009/iKSkOnjuuftQpEiU7rGI/oAvTsL8zspNtmchIujR4w5kZfXFsmWZ\nuPPOt/Dpp/z4M9uPi/OZngVPM8haFSqUxPDhjfHTT9eiQ4c5iI+vghEjmqBsWV42SGZjVUKE3O77\nH/9YhlmzvsKrrzZBQsIdiIi46L9SiVzBjpsonzZu3Iu+fRchKioCY8e2xJ13VtQ9ElmKHTfM76zc\nxCwcF2ZRp04FrF3bB0lJddCq1UwkJaXj4METeoZzGY8Lh+lZWLO4ifIrIkLQs2dtbNv2OAoXjkSN\nGuMxceIGfqYlGYNVCdFlZGXtw8CBGTh06CQmTnwI99xzg+6RyALsuIkCpJTCzJlf4qmnPkDbttXw\n4ov38+oTCil23DC/s3ITs3DkNwsRQbdut2HLlv6IiopA9erjkZy8EV46KeFx4TA9C2sWN1EwXH11\nMYwdG48PPuiBN97YgGbNpuPbbw/qHossw6qEqICys/0YNWoNRoz4BAMG3IPBgxugUKFI3WORR7Dj\nJgqhnTsPISlpEfbvP4apU9uiVq3rdI9EHsCOG+Z3Vm5iFo5gZFG5cmksXtwNffvGonHjqRg2bBlO\nn84JfDiX8bhwmJ6FNYubKJREBElJscjK6osNG/aifv0p7L4pZFiVEAWZUgrjxq3DCy+swIsvNkZS\nUixE+L4ndGXYcRNpsGXLz0hImI+KFUti6tS2uPrqYrpHojDCjhvmd1ZuYhaOUGZRo8a1WL26N/7y\nl6sRGzsJa9fuCtlzBQOPC4fpWQS0uEXkNRHZKiKfi0iqiJQO1mBEXlC4cCRGjWqB119vhtatZ+HV\nVz/x1C/tkB4BVSUi0hTAUqWUX0RGAIBS6u8X3IdVCRFyLxvs0GEOoqPLIDm5NUqUKKx7JDJYyKoS\npdSHSqlzb5m2FgDffYfoEipXLo3ly3246qpCqFt3Mr7++hfdI1GYCmbH3QvA4iA+XlCZ3lm5iVk4\n3M6iaNEoJCe3xsCB96BBgylYvPgbV58/LzwuHKZncdnFLSIfisgXF/nT6rz7PAvgtFJqRkinJfIA\nEUFiYiwWLOiC3r0XYvToNey96Ypc9sOClVJN8/q+iPgAxAN44FL38fl8iI6OBgCUKVMGMTExiIuL\nA+D8zcbb7t4+x5R5dN0+9zUdz1+vXmWMHFkNgwfPxpYtP2PcuHisWrVSWx5xcXHa//ew+XZmZiZS\nUlIA4Pd9eSmBvjjZAsDrABoppQ5c4j58cZIoD0eOnEK3bqk4fvwMUlM7oXTporpHIgOE8jrusQBK\nAPhQRDaJyIQAHy9kLjzTtBmzcJiQRcmSRZCW1hm33loWjRqlYO/eI1rmMCELU5ieRaBXldyilPov\npVTts3/6B2swIptERkZg3Lh4dOpUE/feOwVffbVf90hkMP7KO5Fhpk/fjL/+9X2kpXXBvfdW1j0O\nacL3KiEKMxkZ36JHj/mYPbsD7r//Jt3jkAZ8rxKY31m5iVk4TM2iRYsqmDOnI7p0mYvU1K2uPKep\nWehgehaXvRyQiPSIi4tGRsYjePDBGTh69DQSEu7QPRIZglUJkeG2bv0ZTZq8g5dfvh+PPhqjexxy\nSV5VCc+4iQxXvfq1WLo0AU2aTENOjkKvXrV1j0SaseO2ELNwhEsWt95aDh9//CiGDcvEm29uCMlz\nhEsWbjA9C55xE4WJqlXLYvlyH+LiUlC0aBR8PtYmtmLHTRRmtm8/gMaNp2LkyObo3LmW7nEoRNhx\nE3lItWrlkJHxCJo1ewdRURF4+OEaukcil7HjthCzcIRrFrfffj0yMh7B448vxpIlwXlP73DNIhRM\nz8KaxU3kNTEx5TF/fmckJKTh00936h6HXMSOmyjMLVnyDXy+BVi6NAG1al2nexwKEv7KO5GHtWx5\nC0aNao7mzafj228P6h6HXGDN4ja9s3ITs3B4JYuuXW/DsGGN0KLFdPz887ECPYZXsggG07OwZnET\neV1SUiw6daqJ1q1n4fjxM7rHoRBix03kIUopPPpoGn777SRSUzsjKornZuGKHTeRJUQEkye3xqlT\nOejf/z1+erxHWbO4Te+s3MQsHF7MonDhSMyd2xHr1+/B669/mu+f82IWBWV6FtYsbiKblCxZBAsX\ndsHIkWuwcOF23eNQkLHjJvKwdet248EHZyAjoztiYyvqHoeuADtuIkvdfXclTJr0EFq3noUffzyk\nexwKEmsWt+mdlZuYhcOGLNq1q44nn6yLdu1m53mZoA1Z5JfpWVizuIls9tRT9VC9ejn07r2QV5p4\nADtuIkucOHEG992XgnbtbsWQIQ11j0OXwY6biFCsWCGkpXXGG29sQGrqVt3jUACsWdymd1ZuYhYO\n27KoVKkU0tI6o2/fRdi27cAfvmdbFnkxPQtrFjcR5YqNrYhXX22Cdu1m4/DhU7rHoQJgx01kqX79\nFmH//mOYO7cTIiIuWqWSRuy4iehPRo9ugX37juKll1boHoWukDWL2/TOyk3MwmFzFkWKRGHevE6Y\nOPEzZGR8a3UWFzI9C2sWNxH9WYUKJfGf/7SHz5eG/fuP6h6H8okdNxFhxIhPsHDhdmRm+lC4cKTu\ncQh5d9xc3EQEv1+hTZtZuOmmMhgzpqXucQh8cRKA+Z2Vm5iFg1nkiogQJCVdg8WLv8HMmV/oHkc7\n048LaxY3EeWtZMkimDOnIwYMyMB33/HT4k0WcFUiIk8BeA1AOaXUn/7XZlVCFF5Gj16Dd97ZjFWr\neqFIkSjd41grZFWJiFQG0BTAD4E8DhGZY8CAe3DDDaXwzDMf6R6FLiHQquR/AfwtGIOEmumdlZuY\nhYNZOM5lISJ4++02WLBgu7VvRmX6cVHgxS0ibQDsUkptDuI8RGSAq68uhnff7YB+/Rbxk3MMlGeB\nJSIfAih/kW89C2AwgGbn3/1Sj+Pz+RAdHQ0AKFOmDGJiYhAXFwfA+ZuNt929fY4p8+i6fe5rpsyj\n83ZcXNwfbt91VyW0bVsUrVq9go0bX0ZkZIRR83rtdmZmJlJSUgDg9315KQV6cVJEagFYCuD42S/d\nAGA3gLuVUvsvuC9fnCQKU36/QtOm76Bx42g899x9usexStBfnFRKfamUul4pdZNS6iYAuwDUuXBp\nm+TCM02bMQsHs3BcLIuICMG0aW0xbtw6rFr1o/tDaWL6cRGs67h5Sk3kUZUqlcLkya3RrVsqfv31\nhO5xCPyVdyLKpwEDluDAgeOYMeNh3aNYgb/yTkQBGzGiCTZt2odZs77UPYr1rFncpndWbmIWDmbh\nuFwWxYsXwvTp7TBgwBLs3OntSwRNPy6sWdxEFLjY2IoYNKguEhLSkJPj1z2OtdhxE9EVycnx44EH\npqFFiyr4+98b6B7Hs/h+3EQUVD/+eAixsZOwbNmjqFXrOt3jeBJfnIT5nZWbmIWDWTiuJIsbbyyN\nESMeQELCfJw5kxO6oTQx/biwZnETUXD16lUbFSqUxPDh/JR4t7EqIaIC27v3CGrXnoi0tC6oW/cG\n3eN4CqsSIgqJChVKYvz4ePToMR/Hj5/RPY41rFncpndWbmIWDmbhKGgWDz9cA3fdVRHPPfdxcAfS\nyPTjwprFTUShM2ZMS8yc+SVWr96pexQrsOMmoqCYN28Lhgz5GJs29UXx4oV0jxP2eB03Ebmia9d5\nKF/+Kowc2UL3KGGPL07C/M7KTczCwSwcwchi3LiWmD37KyxfviPgx9LJ9OPCmsWdlZWlewRjMAsH\ns3AEI4uyZYtjwoQHkZiYjpMns4MwlR6mHxfWLO7ffvtN9wjGYBYOZuEIVhZt296KO+4ojxdeWB6U\nx9PB9OPCmsVNRO4ZO7YlkpM3YePGvbpH8SRrFveOHTt0j2AMZuFgFo5gZlG+fAn8+99N8dZbnwXt\nMd1k+nHhylUlIX0CIiKP0nY5IBERBZc1VQkRkVdwcRMRhRnrFreIPCUifhG5RvcsuojIayKyVUQ+\nF5FUESmteya3iUgLEdkmIt+IyDO659FFRCqLyDIR+UpEvhSRAbpn0k1EIkVkk4ik657lUqxa3CJS\nGUBTAD/onkWzDwDUVErdAeBrAIM1z+MqEYkEMA5ACwA1AHQVkep6p9LmDIAnlVI1AdQF8LjFWZwz\nEMAWAMa+AGjV4gbwvwD+pnsI3ZRSHyqlzn1E91oAtr0D/t0AvlVK7VBKnQEwC0AbzTNpoZTap5TK\nOvufjwLYCqCi3qn0EZEbAMQDmAzgold0mMCaxS0ibQDsUkpt1j2LYXoBWKx7CJdVAnD++4/uOvs1\nq4lINIDayP3L3FYjAfwPAP/l7qhTlO4BgklEPgRQ/iLfeha5dUCz8+/uylCa5JHFEKVU+tn7PAvg\ntFJqhqvD6WfsP4F1EZESAOYCGHj2zNs6IvIQgP1KqU0iEqd7nrx4anErpZpe7OsiUgvATQA+FxEg\ntxr4TETuVkrtd3FE11wqi3NExIfcfxI+4MpAZtkNoPJ5tysj96zbSiJSCMA8ANOVUmm659HoXgCt\nRSQeQFEApURkmlIqQfNcf2LlL+CIyP8BiFVKHdQ9iw4i0gLA6wAaKaUO6J7HbSISBWA7cv/S2gNg\nHYCuSqmtWgfTQHLPZKYC+EUp9aTueUwhIo0APK2UaqV7louxpuO+gH1/W/3RWAAlAHx49rKnCboH\ncpNSKhvAEwDeR+7VA7NtXNpn1QfwCIDGZ4+FTWf/YieD94SVZ9xEROHM1jNuIqKwxcVNRBRmuLiJ\niMIMFzcRUZjh4iYiCjNc3EREYYaLm4gozHBxExGFmf8HY9oiPUvgyNwAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f661d5bd750>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Draw Elliptic Curve over Finite field of (prime) order n:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(0, p):\n",
      "    y2 = (pow(i, 3) + i * a + b) % p\n",
      "    y = sqrt(y2)\n",
      "    if y % 1 == 0:\n",
      "        y =  prime_mod_sqrt(y2, p)\n",
      "        for yy in y:\n",
      "            #print '(%d, %d)' % (i, yy)\n",
      "            plt.plot(i, yy, 'ro')\n",
      "#print '(Inf)'\n",
      "\n",
      "#plt.axis([-1, 5, -1, 5])\n",
      "plt.grid()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+QVOWd7/H3F5phENBBkjigJsNOcIPRWogx3CQbZ2J0\n0Jj1R+1e0KtEIkYIi7Cmdm+ACTWzUYh6K+uKoCarqUBc3Vi7NxaCP9AsDMld42RdRiDIBSYZbwYd\n4rIQRQEd+d4/+pzmTE/3/Og5PdPH+ryqujjnOb8+/Zynz9N9nu7B3B0RERGAYUMdQERESoc6BRER\nyVCnICIiGeoUREQkQ52CiIhkqFMQEZGMPnUKZjbczLaZ2ZPB/Olm9pyZ7TGzTWZWEVl3qZntNbPd\nZlYXKb/AzHYEy+6N/6mIiMhA9fWTwmJgFxD+qGEJ8Jy7nwP8LJjHzM4FZgHnApcB95uZBds8AMx1\n98nAZDO7LJ6nICIicem1UzCzs4AvAw8B4QX+SmBtML0WuDqYvgp4zN3fc/c2YB8w3cwmAGPdvTlY\nb11kGxERKRF9+aRwD/A3wIlI2RnufiCYPgCcEUxPBNoj67UDZ+Yo3x+Ui4hICemxUzCzrwC/d/dt\nnPyU0IWn/06G/laGiMgHQKqX5Z8DrjSzLwPlwKlm9mPggJlVuntHcGvo98H6+4GzI9ufRfoTwv5g\nOlq+P9cBzUwdjIhIP7l7zjfu/dXjJwV3X+buZ7v7JOBa4F/dfTawHrgxWO1G4Ilgej1wrZmVmdkk\nYDLQ7O4dwJtmNj0YeJ4d2SbXcUv60dDQMOQZlFM5lVMZw0ecevuk0O16Hfx7J/C4mc0F2oCZAO6+\ny8weJ/1NpU5ggZ9MvAD4ETAKeMrdnxlY9KHT1tY21BH6RDnjpZzxSkLOJGSMW587BXdvApqC6f8C\nLsmz3kpgZY7yl4DzC4spIiKDQb9oLsCcOXOGOkKfKGe8lDNeSciZhIxxs7jvRw2UmXmpZRIRKWVm\nhg/GQLPktmXLlqGO0CfKGS/ljFcSciYhY9zUKYiISIZuH4mIJJxuH4mISFGoUyhAUu4zKme8lDNe\nSciZhIxxU6cgIiIZGlMQEUk4jSmIiEhRqFMoQFLuMypnvJQzXknImYSMcVOnICIiGRpTEBFJOI0p\niIhIUahTKEBS7jMqZ7yUM15JyJmEjHFTpyAiIhkaUxARSbg4xxT6+99xDopZH/oQNQsXsqCxsajH\n2bpxI5tWrSJ1/DidI0dSt2gRF11xRWKPE4etGzeybvlyjrS1MdKd0ZMmcd6VV/LaCy/wxv79HO7o\nYMKECYyZOLHfz6Mv9VCMurq/sZGm1asZ1dnJ0VSqoLYVxz5yGcy20Z9j5WoH195+e871C30OcbeH\nJL3OSllJdgqTDx5ky+23A3DehRf2+0T3tbE9u3gxK1pbM2X1wXRvDf9X+/fzsbFj+cipp9I5ciQT\nP/tZXnvhhW7HC19YI155hQeOHev3cQbasLds2UJtbW3OZfku/i99//tUdnTwULjetm384/btXP/+\n+zwLfB/g4EHYubPH55HrePnqe9v27byxZQtv7N+P/eY3PHj0aLd1CqmDrRs3cu+iRZz+m9/wk0j5\n/BUruB/6fFG/v7GR7StW8I3OTmoL3Ee+fP1pg7m2z24rAP/Q2Ej16NFd2mZ/6nbrxo2svfnmTDvY\nCmzato0fXHMN/3TeeV06h/sbG9l+993d9rvzV7/iP9avz9upbN24kR/ccguPvPZazjy5XjtbgTU/\n/zn/u7q625uSgdZlPj29hj6w3D3vAygHXgRagF3Ad4PyRqAd2BY8Lo9ssxTYC+wG6iLlFwA7gmX3\n9nBMd/CbwGvA55eXuwdlTeCXl5X5n48Z41+tqPBvTJvmTRs2eNSahgafN2pUZhsHX1Zd3WW9NQ0N\n/pVUqss64WPutGk+d9o0nzVuXOYYaxoafFl1dSbD9ZH1m8DnZe1rWXV1Zpv6HMdw8G/PmOHZmjZs\n6HKcevDZ5eU5n2dfbN68OWd504YNflNlpS/LyvTnZWXd8tZn/duX55FLfV1d3vq+fuLEvMdoAp85\nfrw31NR4fV1dph6aNmzw+rq6buXZdTkzT+5Z48f3uo/QzPHj3cE359hHX7L0t076UqfRthJ9bXx9\nxIhMzmjbzK7bsH1dOXq0zxw/3hefd14me31dXWb9JujWTsLXU9OGDT4z67UWbvPfy8q6bXdD8Hpq\nqKnxq8aO7Vaf4XMPn1t91j7z5RhoXfZ0DvO9hkpN+lKe/1ren0fvK8Apwb8p4JfAnwINwDdzrHtu\n0IGMAKqAfZwct2gGPhNMPwVclud4mQaQ3ShuytEwbqus7HKhyNVIo41jTUODz0ulvCHfBWjEiJwX\ny+yLZGb9PBed8ELSkLV+fVA2a9y4bhePsGH39gIYqOiLPvpoyMobzd+Qlb8+mG+oqenTMRtqanJf\nnMeN63as3i4E0U46X/2EdXltnvPz5bKyXvcRuvG003Lu48bTTut2ce7PucpXJw01Nb6mocFnjh/v\nN552ms8cP97XNDR02fYb06Z1q6d8HXqudrish/pdfN55mfV7ejNQX1eX85zNzLFd9FhN4LPz7De8\nKGdn7u1NSU912Zu+nMNCO/7BEmen0Ou3j9z9nWCyDBgOHArmcw1qXAU85u7vuXtb0ClMN7MJwFh3\nbw7WWwdcne+Ym4AVdL23tQmYEJRH/V1HB8/dd196nVWrmBL5GBs1PPgI2rR6NQ92dtKZ57iT33uv\n2zHOe/fdzHSYaSvwLDAlz3MY1Zk+QmfW+neQ/pj1T4cO8ezixWzduPHkvo8fz+TIzrCitTXzPAcq\ndfx4zvuGnZG80TJIfyyM5r8jmG9/880+HbNz5Mic5dHS7GPnq4em1au73CYIy6P1E9blkTx5/N13\ne91H6Ggq913WY6kUm1at6vN+suWrkx2/+x3bV6zgJwcP8qM//IGfHDzI9hUruD+4VbV140aOvPJK\nZv1crxey5qN1G66fr35ff/31zPr57i8PP3YsfdsqUhZ9TWRvFz3WJuCjefb7fnl55txF991TDshf\nl++Xl+fZMpKtl3MY3pq6Y9MmGpuauGPTpm6v3Q+SXjsFMxtmZi3AAWCzu/86WHSrmb1sZg+bWUVQ\nNpH09SPUDpyZo3x/UJ5T2ACyG0VvDSO7kUaFjSO8WNcB9Vnr/L/y8rwXy+j0Fk428nzHCy8k4XH6\ncqEPG3Zvz7Ov8n3HunPkyJy564A9ZWVd6qUOmDd8OGXkyA+MtL594aFu0SLqq6u7lC2rrmZ0VRVh\nyuxzkq8ewnOYLVo/YV2eSvfzvAwYk+cikquOaxYuZH4qlckJMC+V4qKFCzMXsL7sJ1u+OvE33uDB\nrOf4YGcnW1evBtIXsY9G9p/9etmSNQ9d6zaV9W+2ispKXq+spJ787fv98vL0OEZkv9HXRPZ2qazp\nOuCGrHXml5dz6a23Zs5ddN+9va7z1eWlt96aZ8tInh7O4ZYtWwbU8SdRrwPN7n4CmGpmpwHPmlkt\n8ADwnWCV24HvAXPjChU2gLBR9HTxhZMNI9pIoxeweaNGcX3QOMKL9UXBsuWkP/68lEpx9pQpdG7b\n1m3/dcD8UaN48OhR6oAfAB/PkTG0rLqamhtuoP6RRzKN6eE82aMXj7pFi6hvbcWyGmD28xyoukWL\nWLt9O/UdHV1y/7Syktp58/iP9eu5rq2NMmDMpEn8yZ/9GS+uWgWHDnXb14fHju3TMcPBvuX33cfw\nY8d4v7ycy4Jz8oNbbqH2tdcy52TWqFFMqK7m9ddfTw9qZ8n3zj1aP2FdTmptpY6T5/l94DJgzZgx\nkONikKuOFzQ2cj/wnXvu4UdmHAs6hAWNjXx7xoxes+STr05eu/76nOuXBx1F6vhxLqb7ayNsi5dG\n5uenUjzY2dmlbt8fPhyOHMn7mvrIWWdx6V13sW75cl7fs4eb33mHh9wzy5dVV2fO3bOtrcxobWU5\nJ9/11QFr6fq6yH5jdRHpwcjoeemcMiVTJ/WtrZnXznJg+4gR3GLGDyKf2qM58tVlXwaZe/uUMZCO\nP4n6/O0jd/+DmW0EPu3uW8JyM3sIeDKY3Q+cHdnsLNJtZX8wHS3fn+9YPxk3jp+9+y4z3n6bPwCz\ngT+MGMHIEye44f33uRky3wL5i3Hj+EJNDZC+EPxg504ufO21TGNrKivjT2bOzDSOCV/5Clf++Mes\nP3GCi4ATwPeGDePy+nrOu/BC7pw9mxsOHeKRYP9bgNXjxnHxokUs/+Uv+V1HB4fefpvtBw7AW29x\nAvgIJxv31lNP5aq5c1mwdClbL7yQ2X/7twx/912OtbfDwYOZd3Fh/tZ33sl8w+GiK65g2/btbPrh\nD/lGezsPHDuWWX9T8AII3/2H34joab62tjb38tGjufGhh1i3fDkX793LCODj55zDrO98hxOjR3Nu\nsG10f6+98AJs2tRj/t7yXHTFFZwYPbrb8gsXLmR5UxPDjx2j9Z13+Nw117B46VK2btzIDbfcws2v\nvZY53vUTJzLh0kup/8UvWNHamr9+gvN911/9FdtffZX1772XOZ/fnziRmq9/nfpHHuHS4MJTS/oi\nc0ZNTc7ns6CxkQWNjd2e34dra7lh587Mt2i2AP8wcSLzgotVb+frxOjRfGnJki7LX41cgKP1fSyV\nYsuWLbS+/TaNQfls4DAwH3iQ9IX2YeCH5eWMPfdcys8/n9m7dlF9yim8X17O54LXSv3DDzOjtZUb\noMvr6fqJE/lMTQ0XXXEFF11xBVu2bOHlF17ocn4+fc01mdfTtu3b+Yef/pTqU07h6M6dbAk68RtJ\n3yO+GDgxbBiVVVV88513uLKjgw9zssOInr+v3n57pn3OuPdelt93H7/r6OD9sjL+uqEh/XyD19PZ\nlZVcduutnBg9usvrJ1f76m3+w7W1mU4omueyW2/lBND69ts5z8f75eX9ej3GOR9Ot7W1EbueBhyA\nDwEVwfQo0rcNvwRURta5DXg0mA4HmsuASUArJweaXwSmkx6L6HGg2T09sPPtGTO8oaYm842Epg0b\nfO60aX7tuHH+1XHjfMGnPpXzWyfZ22Vb09Dgs4JBvFlZg3h9OUa4Xvbg1NIeBhgLWb+35zGY+ps/\nzuPmqof+1E8c+ygkY6HCL0NE6/qWVCrTTnOdi69VVvqCT32qzxnCzLd88pM+K/j20UCz99ZGovU0\nd9q0fuUttp7O4VC1/f4gxoHmHn/RbGbnk/4kOCx4/Njd/5eZrQOmAg78Fpjn7geCbZYBN5H+lLjY\n3Z8Nyi8AfhR0Lk+5+6I8x/SeMpWC8J3J1o0beS7ycfXSXj6u9nf9uHLGpVj5k/Jd8MHMeX9jI1tX\nr6a8s7PL7apQT+diKOuzP20kCee90Nf6YIvzF836MxcFSEJjBuWMm3LGKwk5k5AR1CmIiEiE/j8F\nEREpCnUKBUjK31hXzngpZ7ySkDMJGeOmTkFERDI0piAiknAaUxARkaJQp1CApNxnVM54KWe8kpAz\nCRnjpk5BREQyNKYgIpJwGlMQEZGiUKdQgKTcZ1TOeClnvJKQMwkZ46ZOQUREMjSmICKScBpTEBGR\nolCnUICk3GdUzngpZ7ySkDMJGeOmTkFERDI0piAiknAaUxARkaJQp1CApNxnVM54KWe8kpAzCRnj\n1mOnYGblZvaimbWY2S4z+25QfrqZPWdme8xsk5lVRLZZamZ7zWy3mdVFyi8wsx3BsnuL95RERKRQ\nvY4pmNkp7v6OmaWAXwB/DVwJ/Ke7321m3wLGufsSMzsXeBS4EDgTeB6Y7O5uZs3AQndvNrOngFXu\n/kyO42lMQUSkHwZ1TMHd3wkmy4DhwCHSncLaoHwtcHUwfRXwmLu/5+5twD5guplNAMa6e3Ow3rrI\nNiIiUiJ67RTMbJiZtQAHgM3u/mvgDHc/EKxyADgjmJ4ItEc2byf9iSG7fH9QnkhJuc+onPFSzngl\nIWcSMsYt1dsK7n4CmGpmpwHPmtkXs5a7mcV6v2fOnDlUVVUBUFFRwdSpU6mtrQVOnqShnG9paSmp\nPEmfV32qPkt1vqWlpaTyhPPhdFtbG3Hr1+8UzGw5cBS4Gah1947g1tBmd/+EmS0BcPc7g/WfARqA\nV4N1pgTl1wE17j4/xzE0piAi0g+DNqZgZh8Kv1lkZqOAS4FtwHrgxmC1G4Engun1wLVmVmZmk4DJ\nQLO7dwBvmtl0MzNgdmQbEREpEb2NKUwA/jUYU3gReNLdfwbcCVxqZnuAi4N53H0X8DiwC3gaWBB5\n278AeAjYC+zL9c2jpIh+hCtlyhkv5YxXEnImIWPcehxTcPcdwKdylP8XcEmebVYCK3OUvwScX1hM\nEREZDPrbRyIiCae/fSQiIkWhTqEASbnPqJzxUs54JSFnEjLGTZ2CiIhkaExBRCThNKYgIiJFoU6h\nAEm5z6ic8VLOeCUhZxIyxk2dgoiIZGhMQUQk4TSmICIiRaFOoQBJuc+onPFSznglIWcSMsZNnYKI\niGRoTEFEJOE0piAiIkWhTqEASbnPqJzxUs54JSFnEjLGTZ2CiIhkaExBRCThNKYgIiJFoU6hAEm5\nz6ic8VLOeCUhZxIyxq3XTsHMzjazzWb2azPbaWaLgvJGM2s3s23B4/LINkvNbK+Z7Tazukj5BWa2\nI1h2b3GekoiIFKrXMQUzqwQq3b3FzMYALwFXAzOBt9z977LWPxd4FLgQOBN4Hpjs7m5mzcBCd282\ns6eAVe7+TNb2GlMQEemHQR1TcPcOd28Jpo8Ar5C+2APkCnEV8Ji7v+fubcA+YLqZTQDGuntzsN46\n0p2LiIiUiH6NKZhZFTAN+GVQdKuZvWxmD5tZRVA2EWiPbNZOuhPJLt/Pyc4lUZJyn1E546Wc8UpC\nziRkjFuqrysGt47+GVjs7kfM7AHgO8Hi24HvAXPjCDVnzhyqqqoAqKioYOrUqdTW1gInT9JQzre0\ntJRUnqTPqz5Vn6U639LSUlJ5wvlwuq2tjbj16XcKZjYC2AA87e5/n2N5FfCku59vZksA3P3OYNkz\nQAPwKrDZ3acE5dcBNe4+P2tfGlMQEemHQR1TMDMDHgZ2RTuEYIwgdA2wI5heD1xrZmVmNgmYDDS7\newfwpplND/Y5G3gijichIiLx6MuYwueBG4AvZn399C4z225mLwM1wG0A7r4LeBzYBTwNLIi89V8A\nPATsBfZlf/MoKaIf4UqZcsZLOeOVhJxJyBi3XscU3P0X5O48nu5hm5XAyhzlLwHn9yegiIgMHv3t\nIxGRhNPfPhIRkaJQp1CApNxnVM54KWe8kpAzCRnjpk5BREQyNKYgIpJwGlMQEZGiUKdQgKTcZ1TO\neClnvJKQMwkZ46ZOQUREMjSmICKScBpTEBGRolCnUICk3GdUzngpZ7ySkDMJGeOmTkFERDI0piAi\nknAaUxARkaJQp1CApNxnVM54KWe8kpAzCRnjpk5BREQyNKYgIpJwGlMQEZGiUKdQgKTcZ1TOeCln\nvJKQMwkZ49Zrp2BmZ5vZZjP7tZntNLNFQfnpZvacme0xs01mVhHZZqmZ7TWz3WZWFym/wMx2BMvu\nLc5TEhGRQvU6pmBmlUClu7eY2RjgJeBq4GvAf7r73Wb2LWCcuy8xs3OBR4ELgTOB54HJ7u5m1gws\ndPdmM3sKWOXuz2QdT2MKIiL9MKhjCu7e4e4twfQR4BXSF/srgbXBamtJdxQAVwGPuft77t4G7AOm\nm9kEYKy7NwfrrYtsIyIiJaBfYwpmVgVMA14EznD3A8GiA8AZwfREoD2yWTvpTiS7fH9QnjhJuc+o\nnPFSznglIWcSMsYt1dcVg1tH/wIsdve3zE5+UgluDcV2z2fOnDlUVVUBUFFRwdSpU6mtrQVOnqSh\nnG9paSmpPEmfV32qPkt1vqWlpaTyhPPhdFtbG3Hr0+8UzGwEsAF42t3/PijbDdS6e0dwa2izu3/C\nzJYAuPudwXrPAA3Aq8E6U4Ly64Aad5+fdSyNKYiI9MOgjilY+iPBw8CusEMIrAduDKZvBJ6IlF9r\nZmVmNgmYDDS7ewfwpplND/Y5O7KNiIiUgL6MKXweuAH4opltCx6XAXcCl5rZHuDiYB533wU8DuwC\nngYWRN76LwAeAvYC+7K/eZQU0Y9wpUw546Wc8UpCziRkjFuvYwru/gvydx6X5NlmJbAyR/lLwPn9\nCSgiIoNHf/tIRCTh9LePRESkKNQpFCAp9xmVM17KGa8k5ExCxripUxARkQyNKYiIJJzGFEREpCjU\nKRQgKfcZlTNeyhmvJORMQsa4qVMQEZEMjSmIiCScxhRERKQo1CkUICn3GZUzXsoZryTkTELGuKlT\nEBGRDI0piIgknMYURESkKNQpFCAp9xmVM17KGa8k5ExCxripUxARkQyNKYiIJJzGFEREpCjUKRQg\nKfcZlTNeyhmvJORMQsa49dopmNkPzeyAme2IlDWaWbuZbQsel0eWLTWzvWa228zqIuUXmNmOYNm9\n8T8VEREZqF7HFMzsC8ARYJ27nx+UNQBvufvfZa17LvAocCFwJvA8MNnd3cyagYXu3mxmTwGr3P2Z\nHMfTmIKISD8M6piCu/8cOJQrR46yq4DH3P09d28D9gHTzWwCMNbdm4P11gFXFxZZRESKZSBjCrea\n2ctm9rCZVQRlE4H2yDrtpD8xZJfvD8oTKSn3GZUzXsoZryTkTELGuKUK3O4B4DvB9O3A94C5sSQC\n5syZQ1VVFQAVFRVMnTqV2tpa4ORJGsr5lpaWksqT9HnVp+qzVOdbWlpKKk84H063tbURtz79TsHM\nqoAnwzGFfMvMbAmAu98ZLHsGaABeBTa7+5Sg/Dqgxt3n59ifxhRERPphyH+nEIwRhK4Bwm8mrQeu\nNbMyM5sETAaa3b0DeNPMppuZAbOBJwaQW0REiqAvX0l9DPg34I/N7HdmdhNwl5ltN7OXgRrgNgB3\n3wU8DuwCngYWRN72LwAeAvYC+3J98ygpoh/hSplyxks545WEnEnIGLdexxTc/bocxT/sYf2VwMoc\n5S8B3W4/iYhI6dDfPhIRSbghH1MQEZEPJnUKBUjKfUbljJdyxisJOZOQMW7qFEREJENjCiIiCacx\nBRERKQp1CgVIyn1G5YyXcsYrCTmTkDFu6hRERCRDYwoiIgmnMQURESkKdQoFSMp9RuWMl3LGKwk5\nk5AxbuoUREQkQ2MKIiIJpzEFEREpCnUKBUjKfUbljJdyxisJOZOQMW7qFEREJENjCiIiCacxBRER\nKQp1CgVIyn1G5YyXcsYrCTmTkDFuvXYKZvZDMztgZjsiZaeb2XNmtsfMNplZRWTZUjPba2a7zawu\nUn6Bme0Ilt0b/1MREZGB6nVMwcy+ABwB1rn7+UHZ3cB/uvvdZvYtYJy7LzGzc4FHgQuBM4Hngcnu\n7mbWDCx092YzewpY5e7P5DiexhRERPphUMcU3P3nwKGs4iuBtcH0WuDqYPoq4DF3f8/d24B9wHQz\nmwCMdffmYL11kW1ERKREFDqmcIa7HwimDwBnBNMTgfbIeu2kPzFkl+8PyhMpKfcZlTNeyhmvJORM\nQsa4pQa6g+DWUKz3e+bMmUNVVRUAFRUVTJ06ldraWuDkSRrK+ZaWlpLKk/R51afqs1TnW1paSipP\nOB9Ot7W1Ebc+/U7BzKqAJyNjCruBWnfvCG4NbXb3T5jZEgB3vzNY7xmgAXg1WGdKUH4dUOPu83Mc\nS2MKIiL9UAq/U1gP3BhM3wg8ESm/1szKzGwSMBlodvcO4E0zm25mBsyObCMiIiWiL19JfQz4N+CP\nzex3ZvY14E7gUjPbA1wczOPuu4DHgV3A08CCyNv+BcBDwF5gX65vHiVF9CNcKVPOeClnvJKQMwkZ\n49brmIK7X5dn0SV51l8JrMxR/hJwfr/SiYjIoNLfPhIRSbhSGFMQEZEPIHUKBUjKfUbljJdyxisJ\nOZOQMW7qFEREJENjCiIiCacxBRERKQp1CgVIyn1G5YyXcsYrCTmTkDFu6hRERCRDYwoiIgmnMQUR\nESkKdQoFSMp9RuWMl3LGKwk5k5AxbuoUREQkQ2MKIiIJpzEFEREpCnUKBUjKfUbljJdyxisJOZOQ\nMW7qFEREJENjCiIiCacxBRERKYoBdQpm1mZm281sm5k1B2Wnm9lzZrbHzDaZWUVk/aVmttfMdptZ\n3UDDD5Wk3GdUzngpZ7ySkDMJGeM20E8KDtS6+zR3/0xQtgR4zt3PAX4WzGNm5wKzgHOBy4D7zUyf\nVERESsiAxhTM7LfAp939YKRsN1Dj7gfMrBLY4u6fMLOlwAl3vytY7xmg0d1/mbVPjSmIiPRDKY0p\nOPC8mf27mX09KDvD3Q8E0weAM4LpiUB7ZNt24MwBHl9ERGI00E7h8+4+Dbgc+Esz+0J0YfCWv6e3\n/Yn8SJCU+4zKGS/ljFcSciYhY9xSA9nY3V8P/n3DzH4KfAY4YGaV7t5hZhOA3wer7wfOjmx+VlDW\nzZw5c6iqqgKgoqKCqVOnUltbC5w8SUM539LSUlJ5kj6v+lR9lup8S0tLSeUJ58PptrY24lbwmIKZ\nnQIMd/e3zGw0sAn4W+AS4KC732VmS4AKd18SDDQ/SrrjOBN4Hvh49gCCxhRERPonzjGFgXxSOAP4\nqZmF+/lHd99kZv8OPG5mc4E2YCaAu+8ys8eBXUAnsEBXfxGR0lLwmIK7/9bdpwaP89z9u0H5f7n7\nJe5+jrvXufvhyDYr3f3j7v4Jd382jicwFKIf4UqZcsZLOeOVhJxJyBg3/U5AREQy9LePREQSrpR+\npyAiIh8g6hQKkJT7jMoZL+WMVxJyJiFj3NQpiIhIhsYUREQSTmMKIiJSFOoUCpCU+4zKGS/ljFcS\nciYhY9zUKYiISIbGFEREEk5jCiIiUhTqFAqQlPuMyhkv5YxXEnImIWPc1CmIiEiGxhRERBJOYwoi\nIlIU6hQKkJT7jMoZL+WMVxJyJiFj3NQpiIhIhsYUREQSrlT+j+ZBt3XjRtYtX86RtjZGujN60iSu\nvf12Lrriii7rbFq1itTx43SOHEndokVdlgPc39hI0+rVjOrs5GgqRc3ChSxobOzXPvqzXiHr93ff\ng6GUMsVGOjpYAAAKcUlEQVRRl6X0fLL1t41O/Oxnee2FF/rdFt/Yv5/DHR1MmDCBMRMnDrgOeqrT\n6LL2N9+kDPjIqaeWRN331hZKua3Ezt0H9QFcBuwG9gLfyrHc3d2bNmzw+ro6b6ip8fq6Ol/T0OA3\nVVb6MnCPPG6rrPSmDRsy2yyrrvYm8HrwBvCZo0b5moYGD61paPB5qVSXfcxLpTLrrGlo8HmjRnVZ\nvqy62tc0NGTyfPnjH/eZf/RHPr+8vNt6YZZsYbbe1m/asMHnTpvWr33ns3nz5rzLsuu3t333NX8h\n/n7lyrxZcuXsT5Z8665paOj388lXn/2ty9701kaz2/kt4DcF6zWB3wA+u7zcvzFtWpd2m11/TdDt\n9dTfcxp97nOnTfPbKitz7i96HsLjbu7nceOu5+h+87WFzZs3F7XtxyW4bsZzjY5rR306GAwH9gFV\nwAigBZiStU6XxhU2/CvNvD6rAYePb8+Y4e6ebig5Gvq8UaMyJ3Dm+PGZ8mjn8ZXgRTczq0MI1ws7\niibwS4Ltcu1n5vjxXTqp+ro6v+WTn/SvZL3Is7OH6y+rru71efbVPffck7O8p84z7JRmjRvnX62o\n8G9Mm5Z5HgPNlO8Cf0nknOS7kESXfWPatD5ladqwIX0+Is+zPjhfM7OO2dvzadqwwS8555xuF6SB\nXjBy1Um+bLPGj3f37u08bC9h2T3Rdhtpd03Bef7zMWO6bJerDqLt4MpTTvFrxo71xeedl/O5h/U7\nK+s1MTco+4tUyq8ZOzazLDzuPTmOm10vt3zykz5z/Hj/Hx/7WM43a729geiLntr2PffcE0vbL7Y4\nO4XBvn30GWCfu7cBmNk/AVcBr0RXqty2jTuArcCzwAqgMd1h5DT82DEAUsePsylYP+rBo0dZft99\nXHTFFYzq7ISsfQPQ2cmsu+9mytGj3fa/KdhHOP35yLJu+zl4kPrFi9n5q1+x/5FHmNHayrPAhF6y\nA2xatYoVra009mHdvjh8+HDO8k2rVmVyZXIfPcrM736XF045hY8ePsxDYfm2bcy+4QaO5agXgLf2\n7+9Tlq0bN/Ls4sWsaG3NlNW3tnLo1FP5/MGDJ9cDrLWVh2fP5tiwYfwksgxgRWsrN44bl/MY0foJ\nj5c6eLDr8wTqAT9ypNd9RPe19uabOd7RAXv2ALB2+3Z46KHMOcvOGLa3nuSrkxF56ro8aLvZ7Tx8\nEYdljZH5B7Pa+0+OHs0sz/fi/317O2tvvpnKjg6+SqR979wJO3dS39rKzl/9iqbVq/nLSP2G+90K\nrAUqId2OOjtpfOutzP7D474MfDuY7wTeaG/vUi9hG/0+8O2DB7kjK2e0nvPVJdDreUgdP56zfPix\nYxw+fLjH5R9Eg/3tozOB30Xm24OyLrIbOaQbTWeenb5fXp5eZ+TIvA09PIFHU6lu+w5NOXo05zFS\nOabD9XLtZ0VrK02rV7OitTWzvLfscLJx9mXdgcjXeZ7z7rt87PDhLuVbgY8ePsywPC+Mjtdf79Mx\n8108325r63KsZ4E7gLWHDjElq0MI5U7StX7C471JjvMDHMnzfHLV8brly6ns6OBi0he+O4DKjg7W\nLV8+oAtGvjr5w7vv5lz/WNB2s9t52F6y2350Pvu1FP032+GODiZ0dLCC3O17Rmsr2+++mykHD+bc\n7ybSb4Ki23VmTW8lfcvgDk7Wqf3mN5l799HXTq7nFgrrOV9dPnfffXm2jOQZOTJnefS60tPyD5rB\n7hTyv92PyNXI64DXSb/Li7qtspJLb701vc6iRbwyalTOfYYnsGbhQuanUjkbWWdwnOxjvDxsWJd1\n2iLr5Wus4SeScHmu/S6rrs5kh5ONry/r9kVb5IIbla/zTNH9+YQvzIpcmYDTKiv7lCXfxfM46fqM\nHiuTM8++xlRVUV9d3TVLVv2ExxuTZx9WVtbrPkJH2tpYEclJkPPttrYBXTDy1Un5xInMT3U9E/NS\nKS5auBDo3s7D9hLWV5gzWn/Zr6V68rezCRMmZNbP1U7CT86defabqx1Fj1UHrAGqs9Z58OhRnrvv\nvky95Or4soX1PJDOuW7Rorxtoa2trcflH0SD+pVUM/tvQKO7XxbMLwVOuPtdkXUGL5CIyAeEx/SV\n1MHuFFLA/wW+BLwGNAPXufsrPW4oIiKDYlAHmt2908wWkr51PBx4WB2CiEjpKLlfNIuIyNApmb99\nZGaXmdluM9trZt8qgTxtZrbdzLaZWXNQdrqZPWdme8xsk5lVRNZfGmTfbWZ1Rcr0QzM7YGY7ImX9\nzmRmF5jZjmDZvYOUs9HM2oP63GZml5dAzrPNbLOZ/drMdprZoqC8pOq0h5wlVadmVm5mL5pZi5nt\nMrPvBuUlU589ZCypuowcY3iQ58lgvvh1GdcPHgbyoA8/ahuCTL8FTs8quxv4n8H0t4A7g+lzg8wj\nguewDxhWhExfAKYBOwrMFH4ybAY+E0w/BVw2CDkbgG/mWHcoc1YCU4PpMaTHu6aUWp32kLMU6/SU\n4N8U8EvgT0uwPnNlLLm6DPb7TeAfgfXBfNHrslQ+KWR+1Obu7wHhj9qGWvZo/pWkf5dD8O/VwfRV\nwGPu/p6nf5i3j/RzipW7/xw4NIBM081sAjDW3ZuD9dZFtilmTuhen0Ods8PdW4LpI6R/RHkmJVan\nPeSE0qvTd4LJMtJv9g5RevWZKyOUWF2a2VnAl0n/BjDMVvS6LJVOoU8/ahtkDjxvZv9uZl8Pys5w\n9wPB9AHgjGB6IunMocHM399M2eX7Gbyst5rZy2b2cORjb0nkNLMq0p9uXqSE6zSS85dBUUnVqZkN\nM7MW0vW22d1/TYnVZ56MUGJ1CdwD/A1wIlJW9LoslU6hFEe7P+/u04DLgb80sy9EF3r6s1hPuQf9\nOfUh01B6AJgETCX9O8TvDW2ck8xsDPAvwGJ3fyu6rJTqNMj5z6RzHqEE69TdT7j7VOAs4CIz+2LW\n8iGvzxwZaymxujSzrwC/d/dt5P4EU7S6LJVOYT9wdmT+bLr2boPO3V8P/n0D+Cnp20EHzKwSIPhY\n9vtg9ez8ZwVlg6E/mdqD8rOyyoue1d1/7wHSH4fD22tDmtPMRpDuEH7s7k8ExSVXp5Gcj4Q5S7VO\ng2x/ADYCF1CC9ZmV8dMlWJefA640s98CjwEXm9mPGYy6jHtgpJAH6QGfVtIDJGUM8UAzcArp+3AA\no4H/Q/rX+XcT/LlvYAndB3nKSL/baCUY5ClCtiq6DzT3KxPpWyTTSb8DKdYAWXbOCZHp24BHhzpn\nsN91wD1Z5SVVpz3kLKk6BT4EVATTo0j/iaMvlVJ99pCxspTqMitzDfDkYLXN2J/AAJ745aS/VbEP\nWDrEWSYFFdwC7AzzAKcDzwN7SP8JmIrINsuC7LuBGUXK9RjpX4K/S3oM5muFZCL97m1HsGzVIOS8\nifRFbTvpP475BOl7o0Od809J369tAbYFj8tKrU7z5Ly81OoUOB/4jyDnduBvCn3dFCtnDxlLqi6z\nMtdw8ttHRa9L/XhNREQySmVMQURESoA6BRERyVCnICIiGeoUREQkQ52CiIhkqFMQEZEMdQoiIpKh\nTkFERDL+P54ppAHEn60NAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f661d4dc850>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Helper function:\n",
      "SQRT mod p is nontrivial, so we ask the Internet:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# via http://codereview.stackexchange.com/questions/43210/tonelli-shanks-algorithm-implementation-of-prime-modular-square-root\n",
      "def legendre_symbol(a, p):\n",
      "    \"\"\"\n",
      "    Legendre symbol\n",
      "    Define if a is a quadratic residue modulo odd prime\n",
      "    http://en.wikipedia.org/wiki/Legendre_symbol\n",
      "    \"\"\"\n",
      "    ls = pow(a, (p - 1)/2, p)\n",
      "    if ls == p - 1:\n",
      "        return -1\n",
      "    return ls\n",
      "\n",
      "def prime_mod_sqrt(a, p):\n",
      "    \"\"\"\n",
      "    Square root modulo prime number\n",
      "    Solve the equation\n",
      "        x^2 = a mod p\n",
      "    and return list of x solution\n",
      "    http://en.wikipedia.org/wiki/Tonelli-Shanks_algorithm\n",
      "    \"\"\"\n",
      "    a %= p\n",
      "\n",
      "    # Simple case\n",
      "    if a == 0:\n",
      "        return [0]\n",
      "    if p == 2:\n",
      "        return [a]\n",
      "\n",
      "    # Check solution existence on odd prime\n",
      "    if legendre_symbol(a, p) != 1:\n",
      "        return []\n",
      "\n",
      "    # Simple case\n",
      "    if p % 4 == 3:\n",
      "        x = pow(a, (p + 1)/4, p)\n",
      "        return [x, p-x]\n",
      "\n",
      "    # Factor p-1 on the form q * 2^s (with Q odd)\n",
      "    q, s = p - 1, 0\n",
      "    while q % 2 == 0:\n",
      "        s += 1\n",
      "        q //= 2\n",
      "\n",
      "    # Select a z which is a quadratic non resudue modulo p\n",
      "    z = 1\n",
      "    while legendre_symbol(z, p) != -1:\n",
      "        z += 1\n",
      "    c = pow(z, q, p)\n",
      "\n",
      "    # Search for a solution\n",
      "    x = pow(a, (q + 1)/2, p)\n",
      "    t = pow(a, q, p)\n",
      "    m = s\n",
      "    while t != 1:\n",
      "        # Find the lowest i such that t^(2^i) = 1\n",
      "        i, e = 0, 2\n",
      "        for i in xrange(1, m):\n",
      "            if pow(t, e, p) == 1:\n",
      "                break\n",
      "            e *= 2\n",
      "\n",
      "        # Update next value to iterate\n",
      "        b = pow(c, 2**(m - i - 1), p)\n",
      "        x = (x * b) % p\n",
      "        t = (t * b * b) % p\n",
      "        c = (b * b) % p\n",
      "        m = i\n",
      "\n",
      "    return [x, p-x]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}