Buffonin neula.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simo Kivelä [pyysi laskinvalmistajia simuloimaan Buffonin neulaongelmaa](http://simokivela.blogspot.fi/2015/05/buffon-laskimet-ja-koulu.html) omilla tuotteillaan. En ole laskinvalmistaja, mutta olen jossain yhteydessä suositellut Pythonia, joten tässä pientä testailua IPython Notebookilla."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "rcParams['figure.figsize'] = (8,8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Piirretään muutama yhdensuuntainen viiva. Sovitaan, että neulan toinen pää osuu yksikköneliön sisään, kuvassa keltaisella. Tehdään kuvan piirtämisestä saman tien funktio, jotta sitä voidaan hyödyntää myöhemmin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAHaCAYAAAAkILGUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE9xJREFUeJzt3V+MpXd93/HP19r6IkGygASDbOxU/C2RIpdKjhGVPFLV\nYFtRlkpWgFYioVJkIVmJlEhxmiB5uQs3lUJI5DpyIlwFmSgSsIBp7AimiIu6Fs4WQ9awqMEYF7at\nwFH5c+E6317McZiuZ2YnnLNzvnvm9ZJGPs85vzm/n595dt7zPGf3THV3AIB5rlj3AgCAvYk0AAwl\n0gAwlEgDwFAiDQBDnVj3Ai5UVf66OQDHSnfXXvePi3Sy/2IPUlWnuvvUJVjOZcs+2Zv9sjf7ZW/2\nywvZJ3v7UffLQSenLncDwFAiDQBDbVKkt9e9gIG2172AobbXvYChtte9gKG2172AgbbXvYChtlf9\nhDXtbUGrqn+U16QB4HJ0UPc26UwaADaKSAPAUCINAEOJNAAMtXSkq+raqvp0VX2pqh6vql/dZ9z7\nq+pcVZ2pqhuWnRcANt0q3nHs/yb59e4+U1UvSvL5qnqou594fkBV3ZrkVd39mqr62ST3JLlpBXMD\nwMZa+ky6u7/V3WcWt7+b5GySay4YdjLJ/YsxjyS5qqquXnZuANhkK31Nuqp+KskNSR654KFrkjy1\na/vpvDDkAMAuK/sFG4tL3X+e5NcWZ9TLPNepXZvb3b19yM+5e5l5AWAJ7z3ML9ioqq0kW4d5wpW8\n41hVnUjyiSSf6u7f2+Pxe5J8prs/vNh+IsnN3X1+j7HecQyAY+Mo3nHsj5P89V6BXjid5J2LxdyU\n5Jm9Ag0A/NDSZ9JV9eYkn03yeJJefPx2kuuTdHffuxj3gSS3JPleknd192P7PJ8zaQCOjYO65xds\nAMAa+QUbAHAZEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFgKJEGgKFE\nGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkA\nGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAo\nkQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQa\nAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFgKJEGgKFEGgCGWkmkq+q+qjpfVV/Y5/Gbq+qZqnps\n8fGeVcwLAJvsxIqe50+S/H6S+w8Y89nu/oUVzQcAG28lZ9Ld/bkk37nIsFrFXABwXBzla9Jvqqoz\nVfXJqnrDEc4LAJelVV3uvpjPJ7muu79fVbcm+WiS1+43uKpO7drc7u7tS7s8ADgaVbWVZOtQY7t7\nVZNen+Tj3f0zhxj7N0n+WXd/e4/HurtdGgfgWDioe6u83F3Z53Xnqrp61+0bs/PDwQsCDQD80Eou\nd1fVh7Jz6v7Sqvp6kruTXJmku/veJLdX1buTPJvkB0netop5AWCTrexy96q43A3AcXJUl7sBgBUS\naQAYSqQBYCiRBoChRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQB\nYCiRBoChRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoCh\nRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJp\nABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFg\nKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIZaSaSr6r6qOl9VXzhgzPur6lxVnamqG1YxLwBsslWd\nSf9Jkrfs92BV3ZrkVd39miR3JLlnRfMCwMZaSaS7+3NJvnPAkJNJ7l+MfSTJVVV19SrmBoBNdVSv\nSV+T5Kld208v7gMA9nFi3QvYS1Wd2rW53d3bh/ycuy/RkgDgYt7b3acuNqiqtpJsHeYJq7uXW9IP\nJ70+yce7+2f2eOyeJJ/p7g8vtp9IcnN3n99jbHd3rWRRsIfbbqu3PPhg/ve618Hl4bbb8hMPPth/\nse51sLkO6t4qL3fX4mMvp5O8c7GYm5I8s1egAYAfWsnl7qr6UHZO3V9aVV/PzmXnK5N0d9/b3Q9W\n1W1V9dUk30vyrlXMCwCbbCWR7u5/fYgxd65iLgA4LrzjGAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAw\nlEgDwFAiDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAi\nDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQA\nDCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCU\nSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCUSAPAUCIN\nAEOtJNJVdUtVPVFVX6mqu/Z4/OaqeqaqHlt8vGcV8wLAJjux7BNU1RVJPpDkXyT5H0keraqPdfcT\nFwz9bHf/wrLzAcBxsYoz6RuTnOvuJ7v72SQPJDm5x7hawVwAcGysItLXJHlq1/Y3Fvdd6E1Vdaaq\nPllVb1jBvACw0Za+3H1In09yXXd/v6puTfLRJK89orkB4LK0ikg/neS6XdvXLu77e9393V23P1VV\nf1hVL+nub+/1hFV1atfmdndvr2CdALB2VbWVZOswY1cR6UeTvLqqrk/yzSRvT/KOCxZ0dXefX9y+\nMUntF+gk6e5TK1gXAIyzOPHcfn67qu7eb+zSke7u56rqziQPZec17vu6+2xV3bHzcN+b5PaqeneS\nZ5P8IMnblp0XADbdSl6T7u7/lOR1F9z3H3bd/oMkf7CKuQDguPCOYwAwlEgDwFAiDQBDiTQADCXS\nADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCUSAPA\nUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCUSAPAUCINAEOJ\nNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIA\nMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQ\nIg0AQ4k0AAwl0gAw1EoiXVW3VNUTVfWVqrprnzHvr6pzVXWmqm5YxbwAsMmWjnRVXZHkA0nekuSn\nk7yjql5/wZhbk7yqu1+T5I4k9yw7LwBsulWcSd+Y5Fx3P9ndzyZ5IMnJC8acTHJ/knT3I0muqqqr\nVzA3AGysVUT6miRP7dr+xuK+g8Y8vccYAGCXE+tewF6q6tSuze3u3j7k59x9iZbEhqla9wq4nJQD\nhsN5b3efutigqtpKsnWYJ1xFpJ9Oct2u7WsX91045pUXGfP3DvM/uc/n/IM/DwCO0uLEc/v57ara\n9wRzFZe7H03y6qq6vqquTPL2JKcvGHM6yTsXi7kpyTPdfX4FcwPAxlr6TLq7n6uqO5M8lJ3o39fd\nZ6vqjp2H+97ufrCqbquqryb5XpJ3LTsvAGy66u51r+H/U1Xd3V4AAuBYOKh73nEMAIYSaQAYSqQB\nYCiRBoChRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoCh\nRBoAhhJpABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJp\nABhKpAFgKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFg\nKJEGgKFEGgCGEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYCiRBoChRBoAhhJpABhKpAFgKJEGgKFE\nGgCGEmkAGEqkAWAokQaAoUQaAIY6scwnV9WLk3w4yfVJvpbkF7v7b/cY97Ukf5vk75I82903LjMv\nABwHy55J/1aSv+zu1yX5dJJ/t8+4v0uy1d3/VKAB4HCWjfTJJB9c3P5gkrfuM65WMBcAHCvLhvNl\n3X0+Sbr7W0lets+4TvJwVT1aVb+y5JwAcCxc9DXpqno4ydW778pOdN+zx/De52ne3N3frKqfzE6s\nz3b35w6Y89Suze3u3r7YOgHgclBVW0m2DjW2e7+uHmqis9l5rfl8Vb08yWe6+59c5HPuTvJ/uvvf\n7/N4d3f9yIsCgMvIQd1b9nL36SS/vLj9S0k+tsfkP1ZVL1rc/vEkP5fki0vOCwAbb9kz6Zck+bMk\nr0zyZHb+CdYzVfWKJH/U3T9fVf84yUeycyn8RJI/7e7fPeA5nUkDcGwc1L2lIn0piDQAx8mlvNwN\nAFwiIg0AQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0A\nQ4k0AAwl0gAwlEgDwFAiDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl\n0gAwlEgDwFAiDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgD\nwFAiDQBDiTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQIg0AQ4k0AAwl0gAwlEgDwFAiDQBD\niTQADCXSADCUSAPAUCINAEOJNAAMJdIAMJRIA8BQS0W6qm6vqi9W1XNV9cYDxt1SVU9U1Veq6q5l\n5gSA42LZM+nHk/yrJP95vwFVdUWSDyR5S5KfTvKOqnr9kvMCwMY7scwnd/eXk6Sq6oBhNyY5191P\nLsY+kORkkieWmRsANt1RvCZ9TZKndm1/Y3EfAHCAi55JV9XDSa7efVeSTvI73f3xS7Goqjq1a3O7\nu7cP+Tl3X4r1AMAhvLe7T11sUFVtJdk6zBNWdy+3pJ0JP5PkN7r7sT0euynJqe6+ZbH9W0m6u9+3\nz3N1dx90+RwANsZB3Vvl5e79wvpokldX1fVVdWWStyc5vcJ5AWAjLftPsN5aVU8luSnJJ6rqU4v7\nX1FVn0iS7n4uyZ1JHkrypSQPdPfZ5ZYNAJtvJZe7V8nlbgCOk6O63A0ArJBIA8BQIg0AQ4k0AAwl\n0gAw1MZEevEOLuxin+zNftmb/bI3++WF7JO9XYr9sjGRziHfYu2Y2Vr3AobaWvcChtpa9wKG2lr3\nAgbaWvcChtpa9RNuUqQBYKOINAAMNfIdx9a9BgA4Svu949i4SAMAO1zuBoChRBoAhhJpABjqsox0\nVd1eVV+squeq6o0HjLulqp6oqq9U1V1HucZ1qKoXV9VDVfXlqvqLqrpqn3Ffq6r/VlV/VVX/9ajX\neVQO8/WvqvdX1bmqOlNVNxz1Go/axfZJVd1cVc9U1WOLj/esY51Hraruq6rzVfWFA8Yct2PlwH1y\njI+Va6vq01X1pap6vKp+dZ9xqzleuvuy+0jyuiSvSfLpJG/cZ8wVSb6a5Pok/yjJmSSvX/faL/F+\neV+S31zcvivJ7+4z7r8nefG613uJ98VFv/5Jbk3yycXtn03yX9a97gH75OYkp9e91jXsm3+e5IYk\nX9jn8WN1rBxynxzXY+XlSW5Y3H5Rki9fyu8tl+WZdHd/ubvPJdnzr6wv3JjkXHc/2d3PJnkgyckj\nWeD6nEzywcXtDyZ56z7jKpfpVZR/gMN8/U8muT9JuvuRJFdV1dVHu8wjddg/Ewf9udpI3f25JN85\nYMhxO1YOs0+S43msfKu7zyxufzfJ2STXXDBsZcfLJn+jvibJU7u2v5EX7shN87LuPp/sHEhJXrbP\nuE7ycFU9WlW/cmSrO1qH+fpfOObpPcZsksP+mXjT4hLdJ6vqDUeztPGO27FyWMf6WKmqn8rO1YZH\nLnhoZcfLiR/lk45CVT2cZPdPHpWduPxOd398PatavwP2y16vB+33j+Df3N3frKqfzE6szy5+aobP\nJ7muu79fVbcm+WiS1655Tcx0rI+VqnpRkj9P8muLM+pLYmyku/tfLvkUTye5btf2tYv7LmsH7ZfF\nX/K4urvPV9XLk/zPfZ7jm4v//q+q+kh2LoNuWqQP8/V/OskrLzJmk1x0n+z+ZtPdn6qqP6yql3T3\nt49ojVMdt2Ploo7zsVJVJ7IT6P/Y3R/bY8jKjpdNuNy932sijyZ5dVVdX1VXJnl7ktNHt6y1OJ3k\nlxe3fynJCw6eqvqxxU+AqaofT/JzSb54VAs8Qof5+p9O8s4kqaqbkjzz/MsFG+qi+2T362ZVdWN2\n3pVw47/pLlT2/35y3I6V5+27T475sfLHSf66u39vn8dXdryMPZM+SFW9NcnvJ/mJJJ+oqjPdfWtV\nvSLJH3X3z3f3c1V1Z5KHsvPDyH3dfXaNyz4K70vyZ1X1b5M8meQXk2T3fsnOpfKPLN4j/USSP+3u\nh9a14Etlv69/Vd2x83Df290PVtVtVfXVJN9L8q51rvlSO8w+SXJ7Vb07ybNJfpDkbetb8dGpqg9l\n59cMvrSqvp7k7iRX5pgeK8nF90mO77Hy5iT/JsnjVfVX2XlZ8bez868mVn68eO9uABhqEy53A8BG\nEmkAGEqkAWAokQaAoUQaAIYSaQAYSqQBYKj/B3RNwXpNr3L6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109794d50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def pohja():\n",
    "    hlines([-1, 0, 1, 2], -1, 2)\n",
    "    bar(0, 1, 1, color='yellow', alpha=0.2)\n",
    "    xlim([-1.1, 2.1])\n",
    "    ylim([-1.1, 2.1])\n",
    "pohja()\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Piirretään kuvaan myös muutama mahdollinen neula. Jos alkupään koordinaatit ovat vaikka $(x, y)$ ja suuntakulma on $\\theta$, loppupään koordinaatit ovat $(x + \\cos\\theta, y + \\sin\\theta)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAHaCAYAAAAkILGUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucVXW9//HXBwa8kagYoihoopGZmZqXrOBoXuCoaMKx\nPOeYdjGLyvJSZpZYJ6tf3io100yzU1neEVGxEu1iZhreEpXyBiqad8XLAJ/fH2tzZsSZYXT2zFqz\n9+v5eMzDvS6z18c1m3nP97vW+n4jM5EkSdUzoOwCJElSxwxpSZIqypCWJKmiDGlJkirKkJYkqaJa\nyi5geRHh7eaSpKaSmdHR+sqFNHRebFciYlpmTuuFcvotz0nHPC8d87x0zPPyWp6Tjr3R89JV49Tu\nbkmSKsqQliSpohoppGeXXUAFzS67gIqaXXYBFTW77AIqanbZBVTQ7LILqKjZ9X7DqNqwoBGRb+Sa\ntCRJ/VFXuddILWlJkhqKIS1JUkUZ0pIkVZQhLUlSRfU4pCNi/Yj4XUTcGRG3R8TnOtnv+xFxb0TM\niYgte3pcSZIaXT1GHFsMHJaZcyJiCHBzRMzKzLnLdoiICcDGmblJRGwHnAFsX4djS5LUsHrcks7M\nRzNzTu3188BdwMjldpsEnFfb50ZgaESs09NjS5LUyOp6TToiNgS2BG5cbtNI4KF2ywt4bZBLkqR2\n6jbBRq2r+0Lg0FqLuifvNa3d4uzMnN3N7zm2J8eVJKkHjuvOBBsRMR4Y3503rMuIYxHRAswArszM\n73Ww/Qzg2sz8VW15LjAuMxd2sK8jjkmSmkZfjDj2E+DvHQV0zXTggFox2wNPdxTQkiSpTY9b0hGx\nI3A9cDuQta+jgdFAZuaZtf1OBXYHXgAOysxbOnk/W9KSpKbRVe45wYYkSSVygg1JkvohQ1qSpIoy\npCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQl\nSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmq\nKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihD\nWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qS\npIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSK\nMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaqouoR0RJwdEQsj4rZOto+LiKcj\n4pba1zH1OK4kSY2spU7vcw7wA+C8Lva5PjP3qtPxJElqeHVpSWfmH4CnVrBb1ONYkiQ1i768Jr1D\nRMyJiCsiYrM+PK4kSf1Svbq7V+RmYFRmLoqICcClwKad7RwR09otzs7M2b1bniRJfSMixgPju7Vv\nZtbroKOByzNzi27sex+wdWY+2cG2zEy7xiVJTaGr3Ktnd3fQyXXniFin3ettKf44eE1AS5KkNnXp\n7o6IX1A03YdFxIPAscBgIDPzTGByRHwKaAVeBParx3ElSWpkdevurhe7uyVJzaSvurslSVIdGdKS\nJFWUIS1JUkUZ0pIkVZQhLUlSRRnSkiRVlCEtSVJFGdKS1AMRrBHB0LLrUGMypCXpDYhgYASfAO4B\nvlF2PWpMfTULliQ1jAjeD3wP2LK26h0RDMxkSYllqQHZkpakbopgdAS/Bq6jCOiHgA8BOxnQ6g22\npCVpBSJYDfgScCSwMsVEQd8BvpvJojJrU2MzpCWpExEEsD9FII+srf4l8KVMHiqtMDUNQ1qSOhDB\nuymuO+9QW3ULcGgmfyivKjUbr0lLUjsRrBvBOcBfKAJ6IfAx4N0GtPqaLWlJAiJYGfg88BVgCNAK\nnAx8M5Nny6xNzcuQltTUatedJwEnAm+prZ4OHJHJvaUVJmFIS2piEWwOnALsXFv1d+DzmVxTXlVS\nG69JS2o6EQyL4FTgVoqAfgr4LPBOA1pVYktaUtOIYBBwCHAcsCawBDgNODaTJ8qsTeqIIS2pKUSw\nC0XX9ma1Vb+l6Nq+o7yqpK7Z3S2poUUwJoLLgFkUAf1PYG9gFwNaVWdIS2pIEawewf+juBlsL+B5\n4Chgs0wuyyRLLVDqBru7JTWUCAYABwLfAobXVp8DfCWTR8qqS3ojDGlJDSOC91IM5blVbdUNFEN5\n3lReVdIbZ3e3pH4vglER/BL4PUVALwD+E9jRgFZ/ZktaUr8VwarAF2tfqwAvAd8FvpPJC2XWJtWD\nIS2p36kN5bkf8P+ADWqrfw18MZMHSitMqjNDWlK/EsHWFNedd6ytmkNx3fn68qqSeofXpCX1CxGs\nE8GPgZsoAvpx4GBgGwNajcqWtKRKi2Al4HPAV4E3AYspWtLfyOSZMmuTepshLamSated9wBOAsbU\nVl8BHJ7J3aUVJvUhu7slVU4EmwFXUczrPAaYC0zIZI9GC+gIDolg59ofJdKrRGa1RsaLiMxMP6xS\nE4pgLWAa8GlgIPAMcCxweiatJZbWKyJYk+KZ7lWATTKZV3JJKkFXuWd3t6TSRdBCcRPYN4C1gKXA\nGcDXMnm8zNp62UcoAvoaA1odMaQllSqCnSmmkNy8tmo2xSNVt5VWVB+ojTH+6dri6WXWourymrSk\nUkSwcQSXAL+hCOj7gX2BnRo9oGt2BjYB5gMzSq5FFWVLWlKfiuBNwNHAYcBg4AXgm8DJmbxUZm19\nbFkr+keZLC61ElWWN45J6hO17t3/Br4NjKitPg/4ciYPl1ZYCSLYgKLnYCmwQSaPlluRyuSNY5JK\nFcEOFAOQvLu26kaK6843lldVqQ6muNx4gQGtrnhNWlKviWD9CP4X+BNFQD8CHAC8p1kDOoLBwCdq\ni94wpi7ZkpZUdxGsAhwBHAWsCrwMnAB8O5Pny6ytAvYB1gHuoJj/WuqUIS2pbmqjZk2mmNN5dG31\nRcCRmdxXWmHV8n+PXWVSrZuCVDneOCapLiLYkuK68/trq24DPp/JteVVVS0RbA7cDjwPrJfJcyWX\npAroKve8Ji2pRyIYHsGPgFsoAvoJ4BBgKwP6NT5V++95BrS6w5a0pDekdgPUZyjG1l6dYgrJU4Gv\nZ/JUmbVVUe358IeBIcA7Mrmj5JJUET6CJamuIpgInAxsWlt1FfCFTOaWV1Xl/RdFQF9vQKu7DGlJ\n3RbBWIpw3r226h6KcJ5ZXlXVV7uhznG69bp5TVrSCkWwRgQnU9z0tDvwLHA4RbetAb1i76UYn3wh\ncEnJtagfsSUtqVMRDAQ+DvwPsDaQwFnAMZk8VmZt/cyyVvRZmbxSaiXqV7xxTFKHIhhP8UjVFrVV\n11M8UvW30orqhyIYATwIDAQ2zOShkktSxfgIlqRui2CjCC4ErqUI6AeB/wDGG9BvyMeAQcB0A1qv\nl93dkgCIYAjFMJ5HACsBi4BvASdm8mKZtfVXEbQAn6wtesOYXjdDWmpytSkk9we+A6xXW/1z4KhM\n5pdWWGP4d2AD4F7gtyXXon7IkJaaWATbUVx33q626q8UU0j+qbyqGsqyG8Z+mMnSUitRv+SNY1IT\nimA9iq7sA2qrHqXo6v6ZYVIfEWxC8Rz5i8BIR2FTZxxxTBIAEawMHAYcDawGvAKcBBzvWNJ1d0jt\nv780oPVG2ZKWmkBtxKt9gBOBDWurLwWOyOQfZdXVqCJYFZgPrAlsk8nNJZekCrMlLTWxCLYATgH+\nrbbqTornnX9TXlUNbz+KgP6LAa2e8DlpqUFFsHYEPwT+RhHQTwJTgS0N6F7nON2qC7u7pQYTwSCK\nkJgGrAEsoQiLaZk8WWJpTSGCdwN/ofijaH2fMdeK9PqIYxFxdkQsjIjbutjn+xFxb0TMiYgt63Fc\nSa8Wwe7AbRTd22sA1wDvzORzBnSfWdaK/okBrZ6qV3f3OcBunW2MiAnAxpm5CcXoO2fU6biSgAg2\njWAGcCUwFpgH7AXslsmdpRbXRCIYBnyotujvOfVYXUI6M/8AXT5iMAk4r7bvjcDQiFinHseWmlkE\nQyM4AbiDYnSr54AvAptncnkm1bqe1fgOBFYGrvKuedVDX93dPRJeNbD8gtq6hX10fKmh1KaQ/Cjw\nTeDNFFNIng18JdN/V2WoDa/6qdqiN4ypLir5CFZETGu3ODszZ3fze47tpZKkCnkfxUie76ot/wE4\nNOCWjwEfC2+7LMluwFXAA8Bbpkc4cFsTOi4zp61op4gYD4zvzhvW7e7uiBgNXJ6ZW3Sw7Qzg2sz8\nVW15LjAuM1/zF793d6u3TZwYu82cyb/KruP1+uEPDxlx3HHHHrpw4YhdAAYPfnnh+PGzvzdjxh6z\nBg1aXHZ5DWviRNaeOTOvXtF+EUwH9gSOzuRbvV+ZGkVfDWYSta+OTKd4PvNXEbE98HRHAS3ptW69\ndYuVP/zhX35k7tyxB2QOWGnAgCUvjR0796e/+tV+P9t88ztfKrs+QQQbAntQDLN6drnVqJHUJaQj\n4hcUTfdhEfEgRbfzYCAz88zMnBkREyNiHvACcFA9jis1stbWFiZOnLn79de//3OvvLLScIARIx65\n+rjjjv3BwQef9WjZ9elVDqZopFyQyWNlF6PGUZeQzsz9u7HPZ+pxLKkZfP7zJ7/tnHMOOuLZZ4e+\nE2DIkOfmfuhD559w1lkHzym7Nr1aBCsBH68tesOY6qqSN45Jzer88/cbduSR3506f/4GewG0tLQ+\nucMON5x6+eV7zhg69FnvRKqmfSnusL8VuKHkWtRgDGmpAu6/f/Sgvfe+dP/bb3/Hx5YuHbhqxNLF\nY8bM+8W55x549nvec8MLZdenLv3fON0+l656M6SlErW2tjB58oXjZs3a9QsvvbTK+gBrr/349Uce\n+d2Tv/jF7z60ou9XuSJ4J7Aj8Czwi5LLUQMypKWSfPWrX9/4tNOmHv7UU2ttC7DKKov+OWnSZSf9\n8pf7/7ns2tRtywYv+Wkmz5daiRqSIS31sZkzJwydOvW0T95//4aTIQYMHLj4ua23vvmMyy6bdNGI\nEQt94LmfiGAo8F+1xR+WWYsalyEt9ZEnnlhr4MSJM/e9+eatD1mypGV1yKWjR99/wQ9+8Nkf7bnn\njKfLrk+v238DqwHXZnJX2cWoMRnSUh/Yf/+fb3fZZZMOW7RotY0B1ljjqZumTj3txP/5n6/OK7s2\nvX4RBO1uGCuzFjU2Q1rqRSeeeNj63/72UV/417/ePA5gpZVemr/rrrNOueiifWc7lGe/Ng54G/AI\ncFnJtaiBGdJSL5k69dTNf/jDT/04c0DLgAFLFm2++R0/ueiifX8xZsw/Xim7NvXYslb0mZm0llqJ\nGlpd5pOW9FrHH3/031ddddE/R46cf/nPfvbfH7z11i3PNaD7vwjWA/YBlgBnlVyOGpwtaamXDB36\n7NK//nWbg8aOvfvlsmtRXX2c4nfnRZksKLsYNTZb0lIvMqAbSwSDgE/WFr1hTL3OkJak7tsTWA+4\nG7i25FrUBAxpSeo+x+lWnzKkJakbIhgL7AwsAs4ruRw1CUNakrrnkNp/f56JI8SpTxjSkrQCEawG\nHFhb9IYx9RlDWpJW7MPAUOCGTOaUXYyahyEtSV1YunQgwNTaoq1o9SlDWpK68I9/fHUssCXwL+DC\nkstRkzGkJakLDz/8n3vWXp6dyUulFqOmY0hLUicuv3yPNRYtesv7gQR+VHY9aj6GtCR14mtf+/pe\nMGAQMDOT+8quR83HkJakDjzzzOoD7rrrbZNri6eVWoyaliEtSR346Ed/8p6XX155vQEDXngUuLrs\netScDGlJ6sB1142bDLDGGjfNyGRp2fWoORnSkrSc733vcyOfeGLYjhFLX9lkk2NmlV2PmpchLUnL\nOe20qR+EiHXXfWTWWmv98dmy61HzMqQlqZ158zYefN99G00C2GefSxy8RKUypCWpnYMPPvMDixcP\nWmPIkOfmnnzyF+4oux41N0Naktr561+3mQLwrnf97YJBgxaXXY6anCEtSTWHH37C2OeeW/0dAwcu\nfv700z/tY1cqnSEtSTUXXDBlMsDo0Q9M33zzOx2nW6UzpCUJmD173JAFC0buDvDxj//4orLrkcCQ\nliQAjjzyu3ssXTpw5TXXfPLGL3/52w+UXY8EhrQk0drawh13bD4FYMcd/3hB2fVIyxjSkpreQQed\n8+6XXlpl9ODBLz927rkH/r7seqRlDGlJTe+aa3aZArDppvdcPGzYk0vKrkdaxpCW1NR+8pODhj/2\n2PBxEUsXf/nL37q07Hqk9gxpSU3thBOO2Adi4PDhj127//6//FfZ9UjtGdKSmtajj67TMm/emH0A\nJky40hvGVDmGtKSmddBB54xvbR289iqrLPrnmWcefEvZ9UjLM6QlNa0//3n7yQBbbHGb43Srkgxp\nSU3puOO+ttHTT6+5zYABS1486aTDZpZdj9QRQ1pSUzr33AMnA2ywwUNXvOc9N7xQdj1SRwxpSU3n\n5pu3WuXBB0ftAbD//r+4sOx6pM4Y0pKazuc+9/0JS5cOXG311Z+Zc/zxX5lXdj1SZwxpSU2ltbWF\nOXO2nAKw3XY3+tiVKs2QltRUpk497Z2LFq22SUtL65Nnnnnw78quR+qKIS2pqcyYsccUgI03/sel\nG274QGvZ9UhdMaQlNY2LL95nzUcfHfEByKWHHXbSxWXXI62IIS2paUybNm3vzAEtb37z49cffPBZ\nj5Zdj7QihrSkpvDMM6sPuPvut+4LsPPOv/WxK/ULhrSkpvCRj/z0va+8stKIlVd+8aGzz/7YjWXX\nI3WHIS2pKfz+9++bArDZZn+/cNVVX8yy65G6w5CW1PBOPPGw9Z98ctgOEUtfPv74oy8vux6puwxp\nSQ3vjDMOmQwwcuSCq3fbbdazZdcjdZchLamhzZ371pXuu2+jvQD23fciRxhTv2JIS2pohxxyxi5L\nlrSsPmTIc3eecsoX7iq7Hun1MKQlNbSbb956CsDWW99sK1r9Tl1COiJ2j4i5EXFPRHypg+3jIuLp\niLil9nVMPY4rSV059NBTNnv++Te9feDAxc+eccYh15Rdj/R6tfT0DSJiAHAqsDPwMHBTRFyWmXOX\n2/X6zNyrp8eTpO66+OIPTgbYaKP7po8de/fLZdcjvV71aElvC9ybmQ9kZitwPjCpg/2iDseSpG65\n+updV1+wYORuAIcccoYjjKlfqkdIjwQearc8v7ZueTtExJyIuCIiNqvDcSWpU0cfffyemQNWWmut\nJ244/PCT5pddj/RG9Li7u5tuBkZl5qKImABcCmzaR8eWtAIR+deya+gtxSAmPfv/C/sBXyPT3tG+\nUI+QXgCMare8fm3d/8nM59u9vjIiTo+ItTLzyY7eMCKmtVucnZmz61CnJEmli4jxwPju7FuPkL4J\nGBMRo4FHgA8BH16uoHUyc2Ht9bZAdBbQAJk5rQ51SeqmzNim7BrqafjwhSc8/vjw8e94x22n3Xbb\nO8/pyXtNnMjaM2fm1fWqTao1PGcvW46IYzvbt8chnZlLIuIzwCyKa9xnZ+ZdEfHJYnOeCUyOiE8B\nrcCLwH49Pa4kdeTMMz8x4vHH3/z+iKWLv/rVb1xWdj1ST9TlmnRmXgW8dbl1P2r3+jTgtHocS5K6\ncvLJX9gHYsA66zw6a8qUCzvtsZP6A0cck9Qw5s8f2TJv3ph9APbYY4aPXanfM6QlNYyPfezsnRYv\nHrTWqqu+MO/00z89p+x6pJ4ypCU1jBtv3G4KwDvfeesFgwYtLrscqccMaUkN4ZhjvjHmmWfWeNeA\nAUteOOWUz19Zdj1SPRjSkhrC//7vf00GGDXqwSu23famRWXXI9WDIS2p3/vTn3ZY7aGHNpgIcMAB\n53nDmBqGIS2p3zv88BMnLF06cNWhQ5+++bjjpv2z7HqkejGkJfVrra0t3HbbFlMAdtjhhgvKrkeq\nJ0NaUr/2yU/+6F2LFq228aBBr/zrnHMOml12PVI9GdKS+rUrr5wwBWDMmHmXjBix0Oeu1FAMaUn9\n1vnn7zds4cJ1doJccsQRJ1xSdj1SvRnSkvqtb37zK/tkDmgZPvyx6z760XMeK7seqd4MaUn90hNP\nrDXwnns2/SDALrtc4w1jakiGtKR+6cADz33fK6+sNHyVVRbdf845B91Udj1SbzCkJfVLf/zjjlMA\n3v72Oy90nG41KkNaUr/zrW8dNfqpp9babsCAJS9997tHzii7Hqm3GNKS+p0f//jj+wKMHLngqvHj\nr3u+7Hqk3mJIS+pXbr11i5UfeGD0XgBTplzgON1qaIa0pH7ls5/9wW5LlrQMedObnr39xBOPmFt2\nPVJvMqQl9RutrS3ccstWUwC22eavPnalhmdIS+o3vvCFkzd/4YUhY1taWp8+88yDf1N2PVJvM6Ql\n9RuXXLLPZICNNrrvsjFj/vFK2fVIvc2QltQvXH75Hms88si6u0Lm1KmnXVx2PVJfMKQl9Qtf+9rX\n98wcMHjYsCf+eOih319Qdj1SXzCkJfULH//4j3+z6aZ3n7vbblf/ouxapL7SUnYBktQdU6ee/sjU\nqaefWnYdUl+yJS1JUkUZ0pIkVZQhLUlSRRnSkiRVlCEtSVJFGdKSJFWUIS1JUkUZ0pIkVZQhLUlS\nRRnSkiRVlCEtSVJFGdKSJFWUIS1JUkUZ0pIkVZQhLUlSRRnSkiRVlCEtSVJFGdKSJFWUIS1JUkUZ\n0pIkVZQhLUlSRRnSkiRVlCEtSVJFGdKSJFWUIS11obW1hWOO+caYt7/9joM33PC+o8quR1JzaSm7\nAKlqWltb+NKXvjP2iiv+/QMPPjhqp5deWmVUsSWXzJw54YcTJ175TLkVSmoWhrQELFq0Shx22Env\n+M1vPrDT/Pnr7/Tyyyuvt2xbS0vr0+uu+8jsHXa44bdbbXXLC2XWKam5GNJqWs88s/qAQw/93ruu\nu27cTgsWjNyptXXwm5dtGzTolX+tt97D144bd91vTzrpsL8NG/bkkjJrldScDGk1lQgGjRp16taj\nRu25zSOPrDt+8eJBay7bNnjwy49usMFDv91pp9/97pRTPn/bqqu+mGXWKkmGtBpeBCsDuwD7ApMe\nfHDqGsu2rbzyiw+NGvXgb3ff/arfnXDCEX8fNGhxaXVK0vIMaTWkCFYDJlAE8x7AkGXbBg589oFN\nN50/a++9L/3tcccdO89gllRVhrQaRgSrUwTyvhQBvUq7zX8DLgIu2nXXoaNnzuRfJZQoSa+LIa1+\nLYK1gEkUwbwLMLjd5j9TBPPFmfxz2cqJExndp0VK0htkSKvfiWAdYG+KYN4JGFjblMD1tAXz/HIq\nlKT6MKTVL0QwEvggMBl4L22j5S0BrqEI5kszWVhOhZJUf3UJ6YjYHTiF4hfn2Zn5nQ72+T7FdcIX\ngAMzc049jq3GFcFGFK3lfYHt2216hbZgnp7JEyWUJ0m9rschHREDgFOBnYGHgZsi4rLMnNtunwnA\nxpm5SURsB5zBq3/pSgBE8FaKUJ4MvKvdpheBqyiCeUYmDs0pqeHVoyW9LXBvZj4AEBHnU9zIM7fd\nPpOA8wAy88aIGBoR62SmXZNNLoIA3kFbi/nt7TY/D8ygCOYrM3FITklNpR4hPRJ4qN3yfIrg7mqf\nBbV1hnQTqgXz1rS1mMe02/w0MJ0imGdl8lLfVyhJ1VDJG8ciYlq7xdmZObub33NsL5WkHgtgB4pc\n/iCwYbttjwGXUuTytWtA6wHAAQARvVRNL72vGlP4gVH3HJeZ01a0U0SMB8Z35w0js2fDE0fE9sC0\nzNy9tnwUkO1vHouIM4BrM/NXteW5wLiOursjIjPTfxENIIIW4H20JfO67TY/AlxMkcy/z8RhvyQ1\npa5yrx4t6ZuAMRExmuIX74eADy+3z3RgKvCrWqg/7fXoxhTBYODfKLqx9wbWbrf5AWqjfgF/zmRp\n31coSf1Hj0M6M5dExGeAWbQ9gnVXRHyy2JxnZubMiJgYEfMoHsE6qKfHVXXUJrDYlaLFvBewRrvN\n99IWzDdn4sxSktRNPe7urje7u/uHdhNYTAb+nXYTWAB3UoTyhcAdBrMkda63u7vVJCIYStsEFrvz\n6gksbqHWYs7k7hLKk6SGY0irSxEMo+jCngx8gI4nsLgok/tKKE+SGpohrdeoTWCxD0WL+d9om8Bi\nKXAdRTBf4gQWktS7DGkBEMH6vHoCi2XXR5zAQpJKYkg3sQjeQttwnNu12/QKxd36yyaweLKE8iSp\n6RnSTSaCsbQF8/ITWFxJ2wQWz5ZQniSpHUO6wbWbwGIyRTBv1m7z88DlFMF8lRNYSFK1GNINqBbM\n29DWYl5+AovLKIL5GiewkKTqMqQbRAQDKGawmExxA9iodpsfp5jB4kLg2kxa+75CSdLrZUj3Y7UJ\nLN5P0Vreh1dPYPEwr57AYknfVyhJ6glDup+pTWCxE0WLeRIdT2BxIXCjE1hIUv9mSPcDEazCqyew\nGNpu8z20TWBxi+NkS1LjMKQrKoIhvHoCi9Xabb6DthbznQazJDUmQ7pCahNY7EnbBBYrt9t8M23j\nZN9TQnmSpD5mSJesNoHFJIpg3gUY1G7zDRTBfLETWEhS8zGkSxDBCNomsBjPayewuJBiAosFpRQo\nSaoEQ7qPRLABxfPL+/LqCSwW0zZO9qWZPFZOhZKkqjGke1EEG9M26te27TYtm8DiQuByJ7CQJHXE\nkK6zCN5GWzBv2W7Ti8BMihbzFU5gIUlaEUO6h2rjZG9B2wQWb2u3+TlgBkWL+apMFvV9hZKk/sqQ\nfgNqwfxu2lrMG7fb/BRtE1j8xgksJElvlCHdTbUJLN5DWzBv0G7zY7RNYDHbCSwkSfVgSHehNoHF\nONomsBjRbvPDtA3H+QcnsJAk1Zsh3YEIdgOmAHsDw9ptup+2YHYCC0lSr4rMag37HBGZmbHiPXuz\nBq6jmAISigksLqQI5r85TrYkqZ66yj1DusMa2A8YSxHMTmAhSeo1hrQkSRXVVe4N6OtiJElS9xjS\nkiRVlCEgnm3wAAAIk0lEQVQtSVJFGdKSJFWUIS1JUkUZ0pIkVZQhLUlSRRnSkiRVlCEtSVJFGdKS\nJFWUIS1JUkUZ0pIkVZQhLUlSRRnSkiRVlCEtSVJFGdKSJFWUIS1JUkUZ0pIkVZQhLUlSRRnSkiRV\nlCEtSVJFGdKSJFWUIS1JUkUZ0pIkVZQhLUlSRRnSkiRVlCEtSVJFGdKSJFWUIS1JUkUZ0pIkVZQh\nLUlSRRnSkiRVlCEtSVJFGdKSJFVUS0++OSLWBH4FjAbuB/4jM5/pYL/7gWeApUBrZm7bk+NKktQM\netqSPgr4TWa+Ffgd8OVO9lsKjM/MdxnQkiR1T09DehLw09rrnwJ7d7Jf1OFYkiQ1lZ4G5/DMXAiQ\nmY8CwzvZL4FrIuKmiPhED48pSVJTWOE16Yi4Blin/SqK0D2mg92zk7fZMTMfiYg3U4T1XZn5hy6O\nOa3d4uzMnL2iOiVJ6g8iYjwwvlv7ZnaWq9060F0U15oXRsQI4NrMfNsKvudY4LnMPKmT7ZmZ8YaL\nkiSpH+kq93ra3T0dOLD2+iPAZR0cfNWIGFJ7vRqwK3BHD48rSVLD62lLei3g18AGwAMUj2A9HRHr\nAmdl5h4RsRFwCUVXeAvw88z8dhfvaUtaktQ0usq9HoV0bzCkJUnNpDe7uyVJUi8xpCVJqihDWpKk\nijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoy\npCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQl\nSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmq\nKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihD\nWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qS\npIoypCVJqihDWpKkijKkJUmqKENakqSKMqQlSaooQ1qSpIrqUUhHxOSIuCMilkTEVl3st3tEzI2I\neyLiSz05piRJzaKnLenbgX2A6zrbISIGAKcCuwFvBz4cEWN7eFxJkhpeS0++OTPvBoiI6GK3bYF7\nM/OB2r7nA5OAuT05tiRJja4vrkmPBB5qtzy/tk6SJHVhhS3piLgGWKf9KiCBr2Tm5b1RVERMa7c4\nOzNnd/N7ju2NeiRJ6objMnPainaKiPHA+O68YWRmz0oqDngtcHhm3tLBtu2BaZm5e235KCAz8zud\nvFdmZlfd55IkNYyucq+e3d2dBetNwJiIGB0Rg4EPAdPreFxJkhpSTx/B2jsiHgK2B2ZExJW19etG\nxAyAzFwCfAaYBdwJnJ+Zd/WsbEmSGl9durvrye5uSVIz6avubkmSVEeGtCRJFWVIS5JUUYa0JEkV\nZUhLklRRDRPStRFc1I7npGOel455XjrmeXktz0nHeuO8NExI080h1prM+LILqKjxZRdQUePLLqCi\nxpddQAWNL7uAihpf7zdspJCWJKmhGNKSJFVUJUccK7sGSZL6UmcjjlUupCVJUsHubkmSKsqQliSp\nogxpSZIqql+GdERMjog7ImJJRGzVxX67R8TciLgnIr7UlzWWISLWjIhZEXF3RFwdEUM72e/+iLg1\nIv4WEX/p6zr7Snd+/hHx/Yi4NyLmRMSWfV1jX1vROYmIcRHxdETcUvs6pow6+1pEnB0RCyPiti72\nabbPSpfnpIk/K+tHxO8i4s6IuD0iPtfJfvX5vGRmv/sC3gpsAvwO2KqTfQYA84DRwCBgDjC27Np7\n+bx8B/hi7fWXgG93st8/gTXLrreXz8UKf/7ABOCK2uvtgD+XXXcFzsk4YHrZtZZwbt4LbAnc1sn2\npvqsdPOcNOtnZQSwZe31EODu3vzd0i9b0pl5d2beC3R4y3rNtsC9mflAZrYC5wOT+qTA8kwCflp7\n/VNg7072C/ppL8rr0J2f/yTgPIDMvBEYGhHr9G2Zfaq7/ya6+nfVkDLzD8BTXezSbJ+V7pwTaM7P\nyqOZOaf2+nngLmDkcrvV7fPSyL+oRwIPtVuez2tPZKMZnpkLofggAcM72S+BayLipoj4RJ9V17e6\n8/Nffp8FHezTSLr7b2KHWhfdFRGxWd+UVnnN9lnprqb+rETEhhS9DTcut6lun5eWN/JNfSEirgHa\n/+URFOHylcy8vJyqytfFeenoelBnD8HvmJmPRMSbKcL6rtpfzdLNwKjMXBQRE4BLgU1LrknV1NSf\nlYgYAlwIHFprUfeKyoZ0Zu7Sw7dYAIxqt7x+bV2/1tV5qd3ksU5mLoyIEcBjnbzHI7X/Ph4Rl1B0\ngzZaSHfn578A2GAF+zSSFZ6T9r9sMvPKiDg9ItbKzCf7qMaqarbPygo182clIlooAvpnmXlZB7vU\n7fPSCN3dnV0TuQkYExGjI2Iw8CFget+VVYrpwIG11x8BXvPhiYhVa38BEhGrAbsCd/RVgX2oOz//\n6cABABGxPfD0sssFDWqF56T9dbOI2JZiVMKG/6VbE3T++6TZPivLdHpOmvyz8hPg75n5vU621+3z\nUtmWdFciYm/gB8DawIyImJOZEyJiXeCszNwjM5dExGeAWRR/jJydmXeVWHZf+A7w64j4KPAA8B8A\n7c8LRVf5JbUx0luAn2fmrLIK7i2d/fwj4pPF5jwzM2dGxMSImAe8ABxUZs29rTvnBJgcEZ8CWoEX\ngf3Kq7jvRMQvKKYZHBYRDwLHAoNp0s8KrPic0LyflR2B/wRuj4i/UVxWPJriqYm6f14cu1uSpIpq\nhO5uSZIakiEtSVJFGdKSJFWUIS1JUkUZ0pIkVZQhLUlSRRnSkiRV1P8H/5N7kYnmQI0AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10af3b6d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def neula(x, y, theta):\n",
    "    plot([x, x + cos(theta)],\n",
    "         [y, y + sin(theta)],\n",
    "         color='blue', lw=2)\n",
    "\n",
    "pohja()\n",
    "neula(0.5, 0.5, 0)\n",
    "neula(0.1, 0.1, 1.1 * pi)\n",
    "neula(0.2, 0.9, pi/6)\n",
    "neula(0.8, 0.3, 3*pi/8)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ilmeisesti neula pysyy kahden viivan välissä, jos $0 < y + \\sin\\theta < 1$. Tehdään tästä funktio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def valissa(y, theta):\n",
    "    return 0 < y + sin(theta) < 1\n",
    "\n",
    "def ylittaa(y, theta):\n",
    "    return not valissa(y, theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sitten vain toistetaan monta kertaa. Tehdään tästäkin funktio, jotta voidaan kohta mitata aikaa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.63641\n",
      "0.637155\n",
      "0.636206\n"
     ]
    }
   ],
   "source": [
    "def simulaatio(n):\n",
    "    ylityksia = 0\n",
    "    n = 1000000\n",
    "    for i in range(n):\n",
    "        y = random.uniform(0, 1)\n",
    "        theta = random.uniform(0, 2*pi)\n",
    "        if ylittaa(y, theta):\n",
    "            ylityksia += 1\n",
    "    return 1.0 * ylityksia / n\n",
    "\n",
    "print simulaatio(1000)\n",
    "print simulaatio(10000)\n",
    "print simulaatio(1e6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aikaa menee minun koneellani pari sekuntia:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 1.84 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit simulaatio(1e6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nopeutta saa lisää vektoroimalla. Seuraava tuhlaa muistia mutta välttää silmukoiden ajamisen Python-tulkissa. Vielä vähän hienompi ratkaisu pilkkoisi tehtävän sopivan mittaisiin paloihin, joista kukin olisi erikseen vektoroitu, ja yhdistäisi tulokset Python-silmukalla."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.63663671\n",
      "1 loops, best of 3: 6.37 s per loop\n"
     ]
    }
   ],
   "source": [
    "def simulaatio2(n):\n",
    "    y = random.uniform(0, 1, size=n)\n",
    "    theta = random.uniform(0, 2*pi, size=n)\n",
    "    sin(theta, out=theta)\n",
    "    add(y, theta, out=y)\n",
    "    return logical_or(greater(y, 1), less(y, 0)).mean()\n",
    "print simulaatio2(1e8)\n",
    "%timeit simulaatio2(1e8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}