{
 "metadata": {
  "name": "9932_05_05"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Chapter 5, example 5\n",
      "====================\n",
      "\n",
      "Here we show how to use Cython with NumPy to accelerate Python algorithms operating on arrays. The code consists in simulating a stochastic process."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 10000"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This function simulates a brownian motion with random, independent steps (up or down). We use a Python loop here."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def step():\n",
      "    return sign(rand(1) - .5)\n",
      "\n",
      "def sim1(n):\n",
      "    x = zeros(n)\n",
      "    dx = 1./n\n",
      "    for i in xrange(n - 1):\n",
      "        x[i+1] = x[i] + dx * step()\n",
      "    return x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(sim1(10000))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "[<matplotlib.lines.Line2D at 0x5a99b90>]"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAD9CAYAAAB5lZr/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX+P/DXmFimGZoJyVAoFwEviGG0dqMMNEuytMQt\nZc3KLNdcrTW3NGxToXLLckv9VirWqq3l5VdI6rp0U2RVrBRLLEhAoBQRb4TA5/fHx8OZMzcuZ4Zh\nZl7Px4PH53POfM45n5ls3nM+t2MQQggQERHp0M7VFSAiIvfHYEJERLoxmBARkW4MJkREpBuDCRER\n6cZgQkREuukOJpmZmQgPD0doaCjS0tKslpk2bRpCQ0MRFRWF3NxcAEBRURFuv/129O3bF/369cOb\nb77ZUD4lJQVGoxHR0dGIjo5GZmam3moSEZETtddzcF1dHaZOnYrt27cjICAAgwcPRmJiIiIiIhrK\nZGRk4MiRI8jPz8fu3bsxZcoUZGdnw8fHB6+//joGDhyIM2fO4Prrr0dCQgLCw8NhMBgwY8YMzJgx\nQ/cbJCIi59N1Z5KTk4OQkBAEBQXBx8cHSUlJ2LRpk6bM5s2bkZycDACIjY1FZWUlysvL4e/vj4ED\nBwIAOnfujIiICJSUlDQcx7mURETuQ1cwKSkpQWBgYMO20WjUBARbZYqLizVlCgsLkZubi9jY2IZ9\nb731FqKiojBp0iRUVlbqqSYRETmZrmYug8HQpHLmdxmmx505cwZjxozB4sWL0blzZwDAlClTMHfu\nXADAnDlzMHPmTLz33nstujYREWk5o+VH151JQEAAioqKGraLiopgNBrtlikuLkZAQAAA4MKFCxg9\nejQefvhhjBo1qqFMjx49YDAYYDAY8OijjyInJ8fq9YUQ/BMCL774osvr0Fb++Fnws+BnYf/PWXQF\nk5iYGOTn56OwsBA1NTVYt24dEhMTNWUSExORnp4OAMjOzoavry/8/PwghMCkSZMQGRmJ6dOna44p\nLS1tyG/YsAH9+/fXU00iInIyXc1c7du3x5IlSzBs2DDU1dVh0qRJiIiIwLJlywAAkydPxogRI5CR\nkYGQkBB06tQJK1asAAB88803+OCDDzBgwABER0cDABYuXIjhw4dj1qxZ2L9/PwwGA3r16tVwPiIi\napsMwpn3PU5kMBicesvmTrKyshAXF+fqarQJ/CxU/CxU/CxUzvruZDAhIvIizvru5HIqRESkG4MJ\nERHpxmBCRES6MZgQEZFuDCZERKQbgwkREenGYEJERLoxmBARkW4MJkREpBuDCRER6cZgQkREujGY\nEBGRbgwmRESkG4MJERHpxmBCRES6MZgQkde6/34gK8vVtfAMfDgWEXmlBQuA55+XeUd8lXTsCFRV\nAT4++s/lTHzSohkGEyLSw2BQ83q/SmprZRCpqAC6dtV3LmfjkxaJiJysqAi49lrgwoXmHXf2rEyr\nqhxfJ3fBYEJEXsc0WHTsKO8sAGDxYhlQXn65eefbulWmQUEOqV6D33937PmcicGEiLzOqVPAZZcB\nH30EnD8P3HKL3H/6tExfeql55/v4YzVfX++YOgoh61hc7JjzORuDCRF5nf/9D6iuBu64Q25nZ8t0\n7161THO6FYYMUfNKQGquZ55R75g++wx48kmZP3WqZedrbeyAJyKvM2gQkJsrA4bSEf/007KZS3H+\nvLwzaArTzvxXX5V9J3//uwwsnTs37xwHDgD9+qn7R4wAMjLkHY/pderqgJoa2UzXHBzNZYbBhIia\n6l//Am69FTAa5baPj+wnMQ0mijfeAKZPB8rKAD+/xs99/jxw+eW2Xz93rvEvfGU0mD1nz2qv89xz\nQFpa80eicTQXEVELPfQQEBiobisd7gDw+OPaskrA+e9/m3bu9HTb5wKA3r3tH79lS9PmplRWarfT\n0ho/pjUxmBCRR1u/Xs2//768QzGVmqrdPndOps88o+7LyQF69JD5U6e0dwOXXCLToiKgSxfL65eV\n2a9fYqL91xUBATI1GLR3U+fOtY1+FQYTIvJoDzyg5idNAr76SuYzMmRqPsmwUyeZ/uEP6r7sbOC3\n32Te11d7V3D+PBAaKu9odu9W97/4okwjIuzXz/QuydTp03LOiylrrVOdOsk6uZruYJKZmYnw8HCE\nhoYizcZ917Rp0xAaGoqoqCjk5uYCAIqKinD77bejb9++6NevH958882G8hUVFYiPj0dYWBgSEhJQ\naX5/R0Sk0113qXllKPDcuUBCghwqrAQdAHj3XZk+/bRMZ89Wg8C0aUB+vsybHvPUU/K8hw61rH6d\nOwP33afdd/So7fJLl7bsOg4jdKitrRXBwcGioKBA1NTUiKioKJGXl6cp89lnn4m77rpLCCFEdna2\niI2NFUIIUVpaKnJzc4UQQpw+fVqEhYWJQ4cOCSGEePbZZ0VaWpoQQojU1FQxa9Ysi2vrrDoReQn5\ne97yz1R9vRBr11oec/y4ED//bP348nJtWSGEmD9fe+7u3S2vZa5HD8tznzwpX1u5Urs/Otr2+2nq\nV6Kzvjt13Znk5OQgJCQEQUFB8PHxQVJSEjZt2qQps3nzZiQnJwMAYmNjUVlZifLycvj7+2PgwIEA\ngM6dOyMiIgIlJSUWxyQnJ2Pjxo16qklEXujXX4FVq2R+82b7ZQ0GYOxYdfuf/5Rp9+62O9CVpVMe\neABYu1bmZ8/WTlpUhvia98so6upkPQ8fBv70J7mvvl5ttrrpJpl+9plMLzbstEm6gklJSQkCTYZI\nGI3GhoBgr0yx2ZTOwsJC5ObmIjY2FgBQXl4Ov4tj8vz8/FBeXq6nmkTkhYYPV7+g77wTmDhRfa1b\nN/vH3n134+ffsUOmp06pX/7mneMbNsh09mzL45OTgQ4dZN7XVy7hsnGj9viQEHnPMWKE9tiSEtkc\n1xxffAEsWtS8Y5pDVzAxmA/QtkGY9RqZHnfmzBmMGTMGixcvRmcrs3sMBoPN66SkpDT8ZfGhBER0\nUVCQ9ld8x45yJNcPP8jtTz6xf7y1YKMsmfLOOzKdPFmmW7cCpaXWz2PaMd6nj6yDIj1dvYu58ko5\nWuvee+3XS9GzJ/D550Bhodx+8UWgf385q99cVlYWUlJSEBeXgmeeSWnaBVpCTxvZrl27xLBhwxq2\nFyxYIFJTUzVlJk+eLNasWdOw3adPH1FWViaEEKKmpkYkJCSI119/XXNMnz59RGlpqRBCiGPHjok+\nffpYXFtn1YnIg7WkL8FUfb32+DlzhKiuFqJ3b+35lfzXXzetLuPGyX21tc2vX2ysLHvkiOVre/fK\n18aObUo92mCfSUxMDPLz81FYWIiamhqsW7cOiWaDphMTE5F+cVZPdnY2fH194efnByEEJk2ahMjI\nSEyfPt3imFUXGztXrVqFUaNG6akmEXkR06VIAGD+/Oafw2CQD8/68EO5/dJLwKWXAj/9JLe7d5fp\nL7/IWekXu38bpTTSfPll8+u0fbu8c+nVy/I1ZTjzunXA3/4GvP1288+vm95olJGRIcLCwkRwcLBY\nsGCBEEKIpUuXiqVLlzaUeeqpp0RwcLAYMGCA2Lt3rxBCiK+++koYDAYRFRUlBg4cKAYOHCi2bNki\nhBDixIkTYujQoSI0NFTEx8eLk8rQBhMOqDoReSDTX/xPPaX/fD/9ZLkvOVmef+RIIS65RN7JNKU+\nytfWjh3q9oYN+utYXGx/dFdRkfPvTLg2FxF5jKoq2f+geOYZufCio5l2uvv6AidP2i9v2u1bXw+s\nWSOXeAGAgweByEh99TlxQr1bMlVfL8/fv7+mNlybi4jIHtNVfwFgzhznXMc0YFlbQsXcvn3A1Kky\nf+6cGkjGjJEd83p16wbMnGm5v7RUG0islXEUBhMi8himS8bv29e0L/qWioqSqb1Z6YroaEBZ5CMl\nRd3/73+ra3vpYTAAr72mXYcMAObNU/PdummDoKMxmBCRx7jySuCxx2TPQHS0c6/V3AmESlPXa6/J\n+SUDBji+TqNHawcgLF+u5isqnLuGF/tMiMhjdOsmv7RPnGid6ykBoqlfRUr5MWPkyKyHH3Z8nc6e\ntf1Arn37gEGD2GdCRGSVwQBcdZXsCK+oaN1rf/tt849Zv955TU7WHtRVVSX/nHm31t55pyYicj7l\n+SNKEFmzpvWu3dwf+H5+gLI6VFMf59tcBoNcDubaa9XZ+ldc4ZxrmeKdCRG5tbNntds1Na6pR1N8\n842aV56P4gyffqqu+7Vli/OuY4rBhIjcmnkwUdaraouCg9W8swcI3HabXO9r+HDnXkfBYEJEbm31\nau32hAmuqUdTzZ0rR4KFhjr3OvfdB5gt0O5UHM1FRG5NGSG1Zo18QqLyrHSyzlnfnQwmROTWmjs8\n19s567uTo7mIyK1dcYW8IyHXYp8JEbm106eBceNcXQtiMCEit3X8uEyNRtfWg9hnQkRubPBgYM8e\noK4OaMefxk3CDngzDCZE3qW6Wk7EU4JGXR3Q/mKvL78Kmo4d8ETk1Tp2BBYuBOLi5PbKla6sDZnj\nnQkRuQXTpxUCctJffr7M86ug6djMZYbBhMi7mAcTRUkJ0LNn69bFnTGYmGEwIfIutoIJvwaax1nf\nnRz/QERtnvlijgr2m7QdDCZE1OYdPmx9f2lp69aDbGMwIaI2z2AAwsPV7TvukKkyaZFcj0ODiahN\n27oVSE8Hrr4a+OQT4KefgHvuAfr3B/76V1fXjhTsgCeiNsu8053/y+vHDngiImqzGEyIyC0kJLi6\nBmQPm7mIqM0ybeZq106ux0X6sJmLiLyK+fddv36uqQc1je5gkpmZifDwcISGhiItLc1qmWnTpiE0\nNBRRUVHIzc1t2P/II4/Az88P/fv315RPSUmB0WhEdHQ0oqOjkZmZqbeaRORmbrpJzb/8MvDNN66r\nCzVOVzCpq6vD1KlTkZmZiby8PKxZswaHDh3SlMnIyMCRI0eQn5+P5cuXY8qUKQ2vTZw40WqgMBgM\nmDFjBnJzc5Gbm4vhw4frqSYRuaFdu2T688/A888DnTu7tj5kn65gkpOTg5CQEAQFBcHHxwdJSUnY\ntGmTpszmzZuRnJwMAIiNjUVlZSXKysoAALfccgu6du1q9dzsDyFyvfp6YP16116vV6/Wuz61nK5J\niyUlJQgMDGzYNhqN2L17d6NlSkpK4O/vb/fcb731FtLT0xETE4NFixbB19fXokxKSkpDPi4uDnHK\ngw6IyCEKCoAHHmi9+R0//CCv9/vvrXM9b5CVlYWsrCynX0dXMDHYWsbTjPldRmPHTZkyBXPnzgUA\nzJkzBzNnzsR7771nUc40mBCR4yn/69bWqk81dKaZM2W6b59Mu3Rx/jU9nfkP7Xnz5jnlOrqauQIC\nAlBUVNSwXVRUBKPRaLdMcXExAgIC7J63R48eMBgMMBgMePTRR5GTk6OnmkTUQsrDp6z8lnMKpQv1\nD3+Q6bJlrXNd0k9XMImJiUF+fj4KCwtRU1ODdevWITExUVMmMTER6enpAIDs7Gz4+vrCz8/P7nlL\nTZYC3bBhg8VoLyKyLTdXzs+oqtJ3nl9+AUaMkPknntB3rpgY4MSJ5h+XlKTvutR6dAWT9u3bY8mS\nJRg2bBgiIyMxduxYREREYNmyZVh28SfFiBEj0Lt3b4SEhGDy5Ml4++23G44fN24chgwZgsOHDyMw\nMBArVqwAAMyaNQsDBgxAVFQUvvjiC7z++ut6qknkVZQxMEePWr6Wn9+0/o+ffgK+/NJxddq7V73L\nIc/EGfBEHmbePCAlRd6hDByofc1gkMHGrAHBgrVuTT3/uxkMssnq8cftl+vZU/uMEv4v7nicAU9E\nzXLmjPX9lZXNO88ll1juy8sDbrtNu89gAD780LLsCy/IdPJk6+ffvRsYPBi47DIZSPbvB77+GmBX\nqXthMCHyMEeOyLSqSq5lZTAAhYXq6+Xl9o9XRlQpnn8e6NhRu2/HDm0zWE2NTM0HCgkBzJ8v8717\nq/vz8tS7jhtvBPbsUYcDDxggZ78PHmy/ntS2MJgQeZgPPpBpVZX67PReveRMcqDxB0r94x9q/uxZ\n2WT2++9yeLDCvBlMCRjm80MWLlTzyvUBoG9feUdiTRNnHFAbw2BC5EStvcqt6V3HuHHA99+r26tX\nq3nTZ6rX18s/helgy44d5Zd7p07a0WHK++reXaYFBTKdMEFbn+eft6yjckfy73/bfy/kXhhMiJxk\n0aKmTfQTAli3zjGzvpUvdcXNN6t50zm+Sgc9IPtEJk6U+bo6IDhYLafcJZw+Dcyape5/+mmZKsN9\nlUUulCB18KCsi/n7r66WnwugDW7k/hhMiJzkmWeaVu6jj+R8Ckcs6KBM9ouNtV9uzRpg0CC1Sezi\nVDB88gmwcycwfLj6muLdd22fT5mr/NFHMu3XT/aRKE1jSj/MqVPqMZ06ydR0FaSKCvv1praLwYTI\nyRobhZmRIdPdu+1PNOzUCfjtN9uvm66xOnKk9TLmM8rHj1fzv/8OrF0r81ddBTz0kOXx589b7jt9\nGnjySdv1AoDoaOC66wB/f/VuprBQzoUxXTbKxrqv5AYYTIic7NVXLfft2aPm77lHpv/9L3DlldaD\nT20tcO6cnEwIyLWrzPtjHn1Uph06ALNnA4sXW57n8ceBjz+23N+5sxyae/Kk3L64NF4Dpbnq1Cm1\nf2XvXllf8/Wz8vK028uWAX/8o5xRDwCvvaa+ZnoH9f77lvUi98FgQmRm6VJtB3VLXXONTM07ms+e\nlcNef/wRuOEG4LvvtK+HhMjAYUrp32h38f/Y66+XI7TuvFMtc/y4TKurZTkluJjP1+jQwbKuypwU\n5Tnr3bppX6+ulunp0+qdRO/e2mYrRd++2u3vv9eO0DIdFXbxaRSoqVH7bcg9tcI6oETu4ehRoH9/\n2dR0xx3Af/4j948dK+dC/OUv1o8LC5NLhfz5z8DVVwNz5sj9p0/L1LzpRhmuGx4u0//9T/v6zz/L\nfTfeCFx6qdz38svqsUozUVGR/PvwQ+0jbZUv7ssvl+nVV2vPb29pPGUQgPndxiWXyGtUVKiTHq+8\nUlumtFQNoKbMZ+Fb4+PTeBlq44SbcuOqUxu1bZsQspFJiMREdb+yT1FfL8TOnZavA0JcdZXcV1sr\nRLt26v4DB9Tyhw5pj7H3V1Ehj5k9W25v3CjEt99alnvnHZnOnq19T+fOyXTxYiGSkrT7IyPlMTfe\nqJ4nPl6ImTOtfz7m1zTfZ749f76sf3295Wvp6UI88IDlZ0vO56zvTjZzEUE2vUyfrm4rzVzWlh45\neBAYMsR6E49y13DmjOwwj4mR2xMnqrPETdY6baDcRZhTOraV4boLFlgvpzRTKU1bCmXm+rRpcgSX\n6f6DB4ElS+SwZMW2bcA//2n9GrYMG6YeExYm0969gb/9Td6VKXdKpsOEO3XSvxIxtTFOCVGtwI2r\nTm2QtbsCIYTYutXy1/Mrr2j3xcRoj/v9dyFWrpT5F15Q97/+uizfv7/ltYKDhdi3T4jycsvXfvtN\nu52RYVnmqqtkWlur//1/8on1Mq++avleACG6dlXL1NTIOydrKiqEiI6Wx9TVCbF9u8xnZbWsztQy\nzvru5J0Jeb1Dh6zvP3bMeme1+XIk5hPzNmyQfwAwerS6f+tWmSp9JoDax9Clixw+26OHZZ/Gr79q\nt5VnjAByPgig3hFZW5SxuWwNzzXtlzG9izOd9e7jI9fWsnVeZaJiu3bqhErzBSPJPTGYkNdTRkGZ\nCwiQTVl33SW/pM2H7P7pTzI1DyZJSeqcj4EDgchImd+yRab33Qe88oo8X06OHJmldPYDlnM27r1X\npmPGaPevXAm8+aa6PWmSrXfYONOmrs6drZexNmRZCOCNN5p+nbAwIDVV5i+9lEvMexIGE/J6pkNV\nb74ZuHBB3S4tle37l14KlJRoj1u5Ut4RfP219fPed59MTeddGAxyOZGdO+W20SjnnFi7G3jkEZkq\nqwAr62ApkpOB0FD1l72eR+s++KCav/5662X8/Vt+foWPj3ZZFvIcDCbk9UxXr/3yS+2dxhNPAOvX\ny3kfgYHqch8zZsjU3ox0panrpZcsX9uxw/ZxSjAzn5S4eLGcjAgAt96q7jcfxquXrVV7o6O1dSAy\nxSctktdTvjy/+kptx//0U9tLkgAyuFx+uWwKM79jUYSFyYmJx48DU6dqm5Ly8oCICOvH/fqrHFV1\nyy1yCRKF8s/988/lhEVl5NTPP8vVeV9/Xd/dw+bN6ppa5Lmc9d3JYEJeTwkmpv+cMjNlXwkg+wRM\nO5wBuaRIO5P7euWOJjRU9o1MmCBX1lX6E6qq1El+7dtrm9LsOXtW7cPgP3dyBGd9d3IGPBHkvBFT\npjOye/WyLG/eFHTLLWp+/HjtAooAcMUVat7Xt+n1UlbW5QOjqK1jMCGvN2SIHF1lyjQ4KF/oCmVC\nY7duTV8y3WBQm8RsjR6zZe9eBhNq+9gBT14vP98yYJjOLzHvE1GanbZvb951Dhxoft0A+dyR6OiW\nHUvUWnhnQl5NCDkiy9oTEfPy5ByRwYPl0iDffgt89pm69InyBa/Mm2iMry/wf/+nLvNO5EnYAU9e\n7dw5eVdy5Ij2cbXWVFfLvpCPPmKzE7kvjuYyw2BCjvDrr3L5Ev5TIm/hrO9O9pmQVzNdSZeIWo7B\nhLyaMvGPiPRhMCGvdu4ccP/9rq4FkfvTHUwyMzMRHh6O0NBQpKWlWS0zbdo0hIaGIioqCrnKw6wB\nPPLII/Dz80P//v015SsqKhAfH4+wsDAkJCSg0toTiogc4IknHPO8dyJvpyuY1NXVYerUqcjMzERe\nXh7WrFmDQ2YPh8jIyMCRI0eQn5+P5cuXY8qUKQ2vTZw4EZmZmRbnTU1NRXx8PA4fPoyhQ4citalj\nL4ma6fjxls//ICKVrmCSk5ODkJAQBAUFwcfHB0lJSdikPMjhos2bNyM5ORkAEBsbi8rKSpSVlQEA\nbrnlFnS1sva26THJycnYuHGjnmoS2dWOjb1EuumatFhSUoLAwMCGbaPRiN2m63nbKFNSUgJ/O8ub\nlpeXw+/i4+b8/PxQXl5utVxKSkpDPi4uDnFxcS14F+St6utlqjwvhMgTZWVlISsry+nX0RVMDE2c\nuWU+prmpxyllbZU3DSZETaX8c/zuO5laW8iRyFOY/9CeN2+eU66j6wY/ICAARUVFDdtFRUUwGo12\nyxQXFyMgIMDuef38/BqawkpLS9GjRw891SRqsGWLbNZ64QXgm29cXRsiz6ErmMTExCA/Px+FhYWo\nqanBunXrkJiYqCmTmJiI9PR0AEB2djZ8fX0bmrBsSUxMxKpVqwAAq1atwqhRo/RUk6jBiBEyXbBA\nPrCKiBxDVzBp3749lixZgmHDhiEyMhJjx45FREQEli1bhmXLlgEARowYgd69eyMkJASTJ0/G22+/\n3XD8uHHjMGTIEBw+fBiBgYFYsWIFAOC5557Dtm3bEBYWhh07duC5557TU03ycD/+CPTt6+paEHk3\nrs1Fbm/RIuCZZxpfX2v9euCBB7T7BgyQqwETeQuuzUVkwzPPNK2ctaChdMITkT4MJuQ1Xn5Zpm+9\npe779FPX1IXI0zCYkNu77TaZrlzZtKXk77tPzd99t1OqROR1GEzIrdXXA198IfMTJwJWVufRmD1b\nPoudiByLwYTapBMngNraxstVVWm37S3j1qcP8PDD6vbNN7esbkRkicGE2qTu3YFrr228nHkw+fJL\n+2WvvFLm77xT3qUQkWMwmFCbkZ8PvPqqut2vX+PHnDpluU9Zc8vUt98CpaVqMNm2TZ3ASET6MZhQ\nm7FwIfDXvwIXLsjt666zX/7UKWDsWJk/cAB4802ZT0gAzpzRlp05U6adOjmuvkSkYjChNqG4GLi4\nAIJF05Utvr6A8vicvn2BP/9Z5v/zH+CKK7RlL7lEps1YY5SImoHBhNoEk6cU4Phxmf7+u+3yTZ3A\n+/e/y+ayrVtbXjciahyDCbU54eEyXb1a9nOcPw88/TTw6KNqme3b1fzKldbPIwQwdy5w8KDcfvFF\np1SXiKDzeSZEetXVATt22H69Z0+ZXnKJLPvuu3I7IUEtc/GhnABk57vy5MQnntCe66qr9NeXiKzj\nnQm51M6d2sBgS12d9f3mExANBuCNN2R++XLta+b9KETkOAwm5FKmfR9NebKoafnly4H9+y3LTJtm\n/djx45tVNSJqBgYTchkhgFmz1O3YWJkGB8vX3n/f8pilS4Fdu2R+zBg5udGcrRFbyoguInI8Ps+E\nXOaXX4CgIHVbCCA7G+jaVS59cv48cPnlto+vrbUdIJSAEhMD7Nmjnp/I2/F5JuRxzp+33HfjjTKQ\nAEDHjsDPP9s+3t6dRkyMTJOTgaNHgZqalteTiBrHYEIuY7r8+/Tp1sv06iXvKJr7Q0pZd2vYMDmH\nxcenZXUkoqZhMxe5jNIUZTQCRUWNl//wQ2DCBHXtLf7nJ2o+NnORRzFtdjpwoGnHhIYCgwbJ/MmT\njq8TEbUcgwm5hLL+1jvvqCv5NqZzZ6CyUvaVNPUYImodnAFPLqEsHW86e70xnToBR47IPBdsJGpb\neGdCLhESItOOHZt+THPKElHrYjChVjd6dMuO69rVsfUgIsfhaC5qdaZNVM39T6gcy//0RC3D0Vzk\nEUz/DX/wQcvO0ZTH+RJR6+KdCbWq3Fx1eG9NDScTErU23pmQRzh9Ws0zkBB5Dt3BJDMzE+Hh4QgN\nDUVaWprVMtOmTUNoaCiioqKQm5vb6LEpKSkwGo2Ijo5GdHQ0MjMz9VaT2ggO6SXyTLrmmdTV1WHq\n1KnYvn07AgICMHjwYCQmJiIiIqKhTEZGBo4cOYL8/Hzs3r0bU6ZMQXZ2tt1jDQYDZsyYgRkzZuh+\ng9S2nDoF3HqrupgjEXkGXXcmOTk5CAkJQVBQEHx8fJCUlIRNmzZpymzevBnJF2emxcbGorKyEmVl\nZY0ey/4Qz3TyJODvb/kURCJyb7ruTEpKShAYGNiwbTQasXv37kbLlJSU4NixY3aPfeutt5Ceno6Y\nmBgsWrQIvr6+FtdPSUlpyMfFxSEuLk7P26FWMGGCTNetc209iLxFVlYWspryGFOddAUTQxMbwJt7\nlzFlyhQ7LeNRAAAVcUlEQVTMnTsXADBnzhzMnDkT7733nkU502BCRESWzH9oz5s3zynX0RVMAgIC\nUGSydnhRURGMRqPdMsXFxTAajbhw4YLNY3v06NGw/9FHH8XIkSP1VJPaCOU3xcKFrq0HETmerj6T\nmJgY5Ofno7CwEDU1NVi3bh0SExM1ZRITE5Geng4AyM7Ohq+vL/z8/OweW1pa2nD8hg0b0L9/fz3V\npDaiokKmwcGurQcROZ6uO5P27dtjyZIlGDZsGOrq6jBp0iRERERg2bJlAIDJkydjxIgRyMjIQEhI\nCDp16oQVK1bYPRYAZs2ahf3798NgMKBXr14N5yP3tnq1TO09bpeI3BNnwFOr2b4diI8Hdu8GbrjB\n1bUh8k6cAU9ur6oKGDGCgYTIEzGYUKupqgKuvtrVtSAiZ2AzlwN16CAXMuzb19U1aZu4fDyR67GZ\nyw1cuAB8/72ra9G2hYa6ugZE5AwMJg5WV+fqGrRtt9zi6hoQkTMwmDjIiy/K9OGHXVuPturECZk+\n/7xr60FEzsE+EwfR8yhab6B8PlVVwBVXuLYuRN6MfSZu5PhxV9eg7erc2dU1ICJnYDBxAKUJR9EK\nC3S6neRk4Oab+XAsIk+lazkVkrp3l+nOncCQIeoaVKRatcrVNSAiZ+KdiU7nzqn5mBiZ/uUvrqlL\nW+ekla+JqA1gMNGprEzN+/jIpwieOwecPu26OrVV7Hgn8lwMJjqZ3pmYbnfpAmzZAlh5ppfXUQaO\n8LE0RJ7LY4NJdTUQG+v861RVabejotT8iBHAo486vw5t3fnzwGWXASEhrq4JETmLxwaTjh2BnBzL\nOwdHO3JEu/3xx807XgigoMBx9WlrPvkE6NQJaM+hHkQezSODydGjan7UKOdeKzlZpkpTjjKyy9Su\nXbaP37kT6N3b8fVqK7Ztk+mZM66tBxE5l0cGE9OmJeXLrLUYDMDFh0k2GDLEdvnKSufWx9W4VhmR\nd/DIYGIaQObPb/3r/+lPlvv277de9p57ZOppS7CcOCH/lOHSROTZPK4lOyNDzRuNwNmzzr1e9+7A\nwYONl7PWd1Nfr+bPnvWspUZiYuQgCH9/uf3jj66tDxE5l8fdmdx9t0wnTQKefVY72mr2bODddx13\nLSFkM9WVV1p/PSZGDg+OjdUGDsWvv6r5U6ccV6+2oLBQzsHZv1+uFBwW5uoaEZEzeUww+fBD7bpP\nd90lh6MuWaJ+kaemAgsWOO6aV18N1NYCl15q+Vr37rKJbfhwoGtXyyHEgHam/AsvOK5ebc0nn7i6\nBkTkbB4TTMyfI2IwyC96ACgqUvcXFADffqtur1sHlJQ07RoFBcCNN6rb5gs8mvrtNyAhQea7dLEe\nTNauVfMrVzatDu5iwAA1f+iQ6+pBRK3DY4LJE0+o+ZtukrOtlZFE48drO7gHDgRuuEHuS0qSfSsG\nQ+PDVx98ENi9W+br64GICDmHojFdugDjxsnAZWrMGODxxxs/3h2ZjuJiExeR5/OYYBIQoOY//1yu\nkzV4sNwuKgJuv11b/n//s1zJNi/P/jX27FHzf/yj/MW9enXjdbtwQaZvvaXdf8UVMqj17Cm3t29v\n/Fzu4Lvv1EEJnTsDBw64tj5E5HweE0yWL1fzl18uU6VJqrAQ+OILy2MmTtRu21t+xfTOprZWfWaJ\nrc53U0rQ+uYb7f4VK4B27YCTJ+X24cONn8sd/Pe/av7tt2VgJyLP5hHBZMYMbb+IaUd8ly76z//A\nA/JLX1FVBZSXy/z5840ff+utlvv27ZPpnj3A++/LvLs+VKumRru9dKlMv/4aGDu29etDRK3PI4LJ\n66/L9P33LeeVKH0cAPCvf1nvFzFvAjO3fr1223Tk1TXXNF6/0FA1rzSlXX+9TO+4A7jqKpn/978b\nP5er/fST9uFfR4/K0WymQXXyZOCRR2TfVYcOrV9HImp9bh1MamuB555TtxMT1SYuRXi4mr/mGtlh\nLoQaIN55B9ixQx1ZVV2tPf733y2v+847aj46uvF6vvKKmu/bV/taYqJ2hFhrLD/y179an/fSFCEh\nwIQJMr9hg5xHAwBz5qhlTp9uWpAlIs+hO5hkZmYiPDwcoaGhSEtLs1pm2rRpCA0NRVRUFHJzcxs9\ntqKiAvHx8QgLC0NCQgIqbSxgdewYYHrJxvovbrtNzY8eLYOKMgpMaYoy77ew14/xzTdNe6Z5t25q\ns5g5Hx/ZEf+Pf8jtQYMaP5/C2gi0c+fs12nFCuDVV+WKyi2lDHO+/37181u0SM7tOX0a+OEHxzQv\nEpH70BVM6urqMHXqVGRmZiIvLw9r1qzBIbNJBRkZGThy5Ajy8/OxfPlyTJkypdFjU1NTER8fj8OH\nD2Po0KFITU21en3zuwhby5x36yZTe1+yShllbor6HtV8Xp52zoS9BRzN9ehh/ZwK5cv9u++adj5l\nZeQfftDuV443HXlmSunPMO0kb65LLrG+PzNTBpF//cv+SslE5Hl0BZOcnByEhIQgKCgIPj4+SEpK\nwqZNmzRlNm/ejOSL67THxsaisrISZWVldo81PSY5ORkbN260en3TJUhmzrRdz+Lixp8/fuml8u7E\nfFkT5Uu3slLOK3HEL25lSRfTu6rmriF23XUyHTxYG+D+8Ad1vzXKMOW//Q0w+0/VZE0ZKFBa2rJz\nE5F70hVMSkpKEBgY2LBtNBpRYjad3FaZY8eO2Ty2vLwcfn5+AAA/Pz+U22gjUtruAdl0Y0vHjsDc\nuY2/n0svBY4f1+6bMUOmShOaMidED6Vp6IYb1H2ms+GVO5evvmras07eeEPeUTXlQWAmrYzNftaL\nacd7Y6yNYCMiz6Vr1WBDUzoMAIgmrK8uhLB6PoPBYPM6P/yQ0pD/4os4xMXFNak+tmzbJv9KS9XV\nbh2pb18ZrJTYaFrdyy8Hhg2TEy5zc+WkylOn5BIutbXaJjzTh38BMvg8+6wMKs1VX68d9mzPb7+p\neWVoMyBHpr3zjhocg4K0d11E5DpZWVnIaoV5B7qCSUBAAIpMJngUFRXBaDTaLVNcXAyj0YgLFy5Y\n7A+4OI3dz88PZWVl8Pf3R2lpKXqYdjhopDTkdMYRAHKZlf375UikzEz55W6Nr2/LHmp14ADw0EOy\nT8GazEzZr6M0USnrjUVGagcCKE1cil9/Bd57T92eOVN2iAvR+ACBt94Cnn66afU/fVrNK0ObAdm0\nZtqs5smPISZyN3Fx2h/a8xpr828hXc1cMTExyM/PR2FhIWpqarBu3TokJiZqyiQmJiI9PR0AkJ2d\nDV9fX/j5+dk9NjExEasuThtftWoVRjn72bsXjR+v5pUvch8fy5FP+flNXxzSnNJcpjwUy54PPlCv\nZ41yJ/L//p92v7Iacbt26hDgmho58iskRC7b8umncv9HH9mvgzI7H7C9srHyn/ymm+yfi4g8mNAp\nIyNDhIWFieDgYLFgwQIhhBBLly4VS5cubSjz1FNPieDgYDFgwACxd+9eu8cKIcSJEyfE0KFDRWho\nqIiPjxcnT560uC4AAQjRq5cQs2bpfRfS+PFCyN/z8k8IIbp2FeL4ccecXwghxo2T5y4qsv666fVN\n/157TYiaGm2ZH3+0XraqSs0vWybEgQMy7+Mj09JSeZ6nnxbilVds1/WHH2T5Cxfk9i23WF6rrs6y\n7kTUdjnga98qw8WTux3ZjyKQmgrMmuWYcz7xBLBsmbpdXy/vTM6dc9xMbqXZqapKzi8xd/fd2qdF\nmtq8Wa6GrJxDCDkiy/zGra7O9vBdQL6fjh2Bl16S/TEvvWS/rsePy1n6s2bJZ7O8+KK6hIrpv56D\nB+WEyM8+s31tInItg8HQpH7s5nLrGfAAsHev4861aJF2++RJ+cXsyCVBvv5apo09otd05r7i3Dl1\naO/338vUdPb8+vXAxo2Nd6hfdplMu3SRgw0aG5ZcUSED6yuvyHMrqx/v3Kkt17cvAwmRt3L7YHLt\ntY47V6dO2mVGlDWzHCk6Wq7VZatj/Ikn5DNOlMAWEaG+Nm6cOg+mXz+Zdu2qvj56NHDvvTJv64dH\nt27qtX185JwXa4HNdJ7Iyy+rdzrp6bJ+110nH0tMRAS4eTD5+GM5+c6RDIaWjdRqqssvt79Ey8iR\nsqltxAgZEEyfsSKE9k4EkHdNH3xguRqAuaQkmZrOFXnzTTV/5ox878qMetPBCKYBduRImRYWcml5\nIlK5dTC5/351GRRHamvrSpneff30k+XrDz1k/Tn0pubPt9w3bpyaVyZ9btgg0//8R6aPPKKOKgMa\nX0mAiLyTWwcTZzFtgjL99e4qv/zSsuNMBxNYm0mfkqLmlU548+YxX181v3gxl5QnIusYTBrRlPkg\nraGpEwtNPf64fM6Istz+qlWW80qGD9d22D//vHbZ/f371fy0ac2vAxF5B7ceGuzMqr/8snxGx6+/\nAldf7bTLNNnWreqM/LffBi4uvqzbb79pVzQ2degQMGCAOoLMPf+lEJEpDg1uZQkJMu3UybX1UCQk\nqE8zvP9+x53XtBnL3DXXqAtkMpAQkT28M7GhulpO7Kuvb9oDsNyZrfdXWysD2Ndfy+YwInJ/zvru\nZDCxo6QEuLj2pEdTgsmNNwLZ2ep+9/yXQUT2sJnLBbwhkADAY4/J+SW7dslmvenTGUiIqHl4Z0JE\n5EV4Z0JERG0WgwkREenGYEJERLoxmBARkW4MJkREpBuDCRER6cZgQkREujGYEBGRbgwmRESkG4MJ\nERHpxmBCRES6MZgQEZFuDCZERKQbgwkREenGYEJERLoxmBARkW4tDiYVFRWIj49HWFgYEhISUFlZ\nabVcZmYmwsPDERoairS0tEaPLywsRMeOHREdHY3o6Gg8+eSTLa2i18jKynJ1FdoMfhYqfhYqfhbO\n1+Jgkpqaivj4eBw+fBhDhw5FamqqRZm6ujpMnToVmZmZyMvLw5o1a3Do0KFGjw8JCUFubi5yc3Px\n9ttvt7SKXoP/o6j4Waj4Waj4WThfi4PJ5s2bkZycDABITk7Gxo0bLcrk5OQgJCQEQUFB8PHxQVJS\nEjZt2tTk44mIyD20OJiUl5fDz88PAODn54fy8nKLMiUlJQgMDGzYNhqNKCkpafT4goICREdHIy4u\nDl9//XVLq0hERK2kvb0X4+PjUVZWZrF//vz5mm2DwQCDwWBRznyfEMJmOWV/z549UVRUhK5du2Lf\nvn0YNWoUDh48iCuuuKLR83uzefPmuboKbQY/CxU/CxU/C+eyG0y2bdtm8zU/Pz+UlZXB398fpaWl\n6NGjh0WZgIAAFBUVNWwXFxcjICDA7vEdOnRAhw4dAACDBg1CcHAw8vPzMWjQIM25hRBNfItERORs\nLW7mSkxMxKpVqwAAq1atwqhRoyzKxMTEID8/H4WFhaipqcG6deuQmJho9/jjx4+jrq4OAPDzzz8j\nPz8fvXv3bmk1iYioFRhEC3/iV1RU4MEHH8TRo0cRFBSEjz76CL6+vjh27Bgee+wxfPbZZwCALVu2\nYPr06airq8OkSZMwe/Zsu8d/8sknmDt3Lnx8fNCuXTu89NJLuPvuux33jomIyPGEG9qyZYvo06eP\nCAkJEampqa6ujsMdPXpUxMXFicjISNG3b1+xePFiIYQQJ06cEHfeeacIDQ0V8fHx4uTJkw3HLFiw\nQISEhIg+ffqIzz//vGH/nj17RL9+/URISIiYNm1aq78XR6mtrRUDBw4U99xzjxDCez+LkydPitGj\nR4vw8HAREREhsrOzvfazWLBggYiMjBT9+vUT48aNE9XV1V7zWUycOFH06NFD9OvXr2GfI997dXW1\nePDBB0VISIiIjY0VhYWFjdbJ7YJJbW2tCA4OFgUFBaKmpkZERUWJvLw8V1fLoUpLS0Vubq4QQojT\np0+LsLAwkZeXJ5599lmRlpYmhBAiNTVVzJo1SwghxMGDB0VUVJSoqakRBQUFIjg4WNTX1wshhBg8\neLDYvXu3EEKIu+66S2zZssUF70i/RYsWiT/+8Y9i5MiRQgjhtZ/FhAkTxHvvvSeEEOLChQuisrLS\nKz+LgoIC0atXL1FdXS2EEOLBBx8UK1eu9JrP4ssvvxT79u3TBBNHvvd//vOfYsqUKUIIIdauXSvG\njh3baJ3cLpjs3LlTDBs2rGF74cKFYuHChS6skfPde++9Ytu2baJPnz6irKxMCCEDTp8+fYQQ8leH\n6R3asGHDxK5du8SxY8dEeHh4w/41a9aIyZMnt27lHaCoqEgMHTpU7Nixo+HOxBs/i8rKStGrVy+L\n/d74WZw4cUKEhYWJiooKceHCBXHPPfeIrVu3etVnUVBQoAkmjnzvw4YNE9nZ2UII+aOle/fujdbH\n7dbmsjd3xRMVFhYiNzcXsbGxNufmHDt2DEajseEY5TMx3x8QEOCWn9Vf/vIXvPrqq2jXTv3n6o2f\nRUFBAa6++mpMnDgRgwYNwmOPPYazZ8965WfRrVs3zJw5E9deey169uwJX19fxMfHe+VnoXDkezf9\nnm3fvj2uvPJKVFRU2L2+2wUTb5pbcubMGYwePRqLFy+2mGdja26Pp/n000/Ro0cPREdH2xwO7i2f\nRW1tLfbt24cnn3wS+/btQ6dOnSyWMfKWz+Knn37CG2+8gcLCQhw7dgxnzpzBBx98oCnjLZ+FNa54\n724XTMznrhQVFWmiq6e4cOECRo8ejfHjxzcMm1bm5gDQzM2xNp/HaDQiICAAxcXFmv3KPB93sXPn\nTmzevBm9evXCuHHjsGPHDowfP94rPwuj0Qij0YjBgwcDAMaMGYN9+/bB39/f6z6LPXv2YMiQIbjq\nqqvQvn173H///di1a5dXfhYKR/w/oXyXBgQE4OjRowDkj5hTp06hW7dudq/vdsHE3twVTyGEwKRJ\nkxAZGYnp06c37Lc1NycxMRFr165FTU0NCgoKkJ+fjxtuuAH+/v7o0qULdu/eDSEEVq9ebXU+UFu2\nYMECFBUVoaCgAGvXrsUdd9yB1atXe+Vn4e/vj8DAQBw+fBgAsH37dvTt2xcjR470us8iPDwc2dnZ\nOH/+PIQQ2L59OyIjI73ys1A44v+Je++91+Jc69evx9ChQxuvQMu7f1wnIyNDhIWFieDgYLFgwQJX\nV8fhvvrqK2EwGERUVJQYOHCgGDhwoNiyZYs4ceKEGDp0qNWhf/PnzxfBwcGiT58+IjMzs2G/MvQv\nODhY/PnPf3bF23GYrKyshtFc3vpZ7N+/X8TExIgBAwaI++67T1RWVnrtZ5GWltYwNHjChAmipqbG\naz6LpKQkcc011wgfHx9hNBrF+++/79D3Xl1dLR544IGGocEFBQWN1qnFkxaJiIgUbtfMRUREbQ+D\nCRER6cZgQkREujGYEBGRbgwmRESkG4MJERHp9v8BSm6eBkoP+XYAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's evaluate the performance of this function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit sim1(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 216 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we convert it in Cython. We specify the type of each variable, and we rewrite the step function in C so that it only uses pure C functions."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "# We need to import NumPy with cimport\n",
      "cimport numpy as np\n",
      "DTYPE = np.double\n",
      "# We define the type of the array for Cython\n",
      "ctypedef np.double_t DTYPE_t\n",
      "\n",
      "# We rewrite the step function in C using the standard library\n",
      "from libc.stdlib cimport rand, RAND_MAX\n",
      "from libc.math cimport round\n",
      "\n",
      "cdef double step():\n",
      "    return 2 * round(float(rand()) / RAND_MAX) - 1\n",
      "\n",
      "def sim2(int n):\n",
      "    # here we specify the type of each variable with cdef, including the type of the array\n",
      "    # with a special syntax so that Cython knows at compile-time the number of dimensions of this array.\n",
      "    cdef int i\n",
      "    cdef double dx = 1. / n\n",
      "    cdef np.ndarray[DTYPE_t, ndim=1] x = np.zeros(n, dtype=DTYPE)\n",
      "    for i in range(n - 1):\n",
      "        x[i+1] = x[i] + dx * step()\n",
      "    return x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit sim2(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 900 us per loop\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We achieve a 240x speed improvement here!"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}