prngs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pseudo-random number generators and random sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Anyone who attempts to generate random numbers by deterministic means is, of course, living in a state of sin. --John von Neumann\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Properties of PRNGs\n",
    "\n",
    "+ dimension of output\n",
    "    - commonly 32 bits, but some have more\n",
    "    \n",
    "+ number of states\n",
    "    - dimension of state space in bits\n",
    "    - sometimes state = output, but better generators generally have output = f(state)\n",
    "\n",
    "+ period\n",
    "    - maximum over initial states of the number of states visited before repeating\n",
    "    - period &le; number of states\n",
    "    - if state has $s$ bits, period $\\le 2^s$\n",
    "    - for some PRNGs, period is much less than number of states\n",
    "    - for some seeds for some PRNGs, number of states visited is much less than period\n",
    "    \n",
    "+ $k$-distribution\n",
    "    - suppose $\\{X_i\\}$ is sequence of $P$ $w$-bit integers\n",
    "    - define $t_v(X_i)$ to be the first $v$ bits of $X_i$\n",
    "    - $\\{X_i\\}$ is $k$-distributed to $v$-bit accuracy if each of the $2^{kv}-1$ possible nonzero $kv$-bit vectors occurs equally often among the $P$ $kv$-bit vectors\n",
    "\\begin{equation*}  (t_v(X_i),\\,t_v(X_{i+1}), \\ldots ,t_v(X_{i+k-1}))\\quad (0\\le i<P),\\end{equation*}\n",
    "and the zero vector occurs once less often.\n",
    "    - amounts to a form of uniformity in $k$-dimensional space, over an entire cycle\n",
    "    - does not measure dependence or other \"serial\" properties\n",
    "\n",
    "+ sensitivity to initial state; burn-in\n",
    "    - many PRNGs don't do well if the seed has too many zeros\n",
    "    - some require many iterations before output behaves well\n",
    "    - for some seeds, some PRNGs repeat very quickly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some PRNGs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Middle Square\n",
    "\n",
    "Dates to Franciscan friar ca. 1240 (per [Wikipedia](https://en.wikipedia.org/wiki/Middle-square_method)); reinvented by von Neumann ca. 1949.\n",
    "\n",
    "Take $n$-digit number, square it, use middle $n$ digits as the \"random\" and the new seed.\n",
    "\n",
    "E.g., for $n=4$, take $X_0 = 1234$.\n",
    "\n",
    "$1234^2 = 1522756$, so $X_1 = 2275$.\n",
    "\n",
    "$2275^2 = 5175625$, so $X_2 = 7562$.\n",
    "\n",
    "+ $10^n$ possible states, but not all attainable from a given seed\n",
    "+ period at most $8^n$, but can be very short. E.g., for $n=4$,\n",
    "    - 0000, 0100, 2500, 3792, & 7600 repeat forever\n",
    "    - 0540 &rightarrow; 2916 &rightarrow; 5030 &rightarrow; 3009 &rightarrow; 0540"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Congruential Generators (LCGs) \n",
    "\n",
    "\\begin{equation*}  X_{n+1} = (aX_n +c)\\mod m.\\end{equation*}\n",
    "\n",
    "LCG period is at most $m$.\n",
    "\n",
    "#### Hull-Dobell Theorem: the period of an LCG is $m$ for all seeds $X_0$ iff\n",
    "\n",
    "+ $m$ and $c$ are relatively prime\n",
    "+ $a-1$  is divisible by all prime factors of $m$\n",
    "+ $a-1$  is divisible by 4 if $m$ is divisible by 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Marsaglia (PNAS, 1968): Random Numbers Fall Mainly in the Planes\n",
    "Multiplicative congruential generators ($c=0$), Lehmer (1949).\n",
    "\n",
    "![Marsaglia68](./RandPics/marsaglia68.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RANDU\n",
    "\n",
    "RANDU is a particularly bad [linear congruential generator](https://en.wikipedia.org/wiki/Linear_congruential_generator) promulgated in the 1960s and widely copied.\n",
    "\n",
    "RANDU is given by the recursion\n",
    "\n",
    "\\begin{equation*}  X_{j+1} = 65539 X_j\\mod 2^{31}.\\end{equation*}\n",
    "\n",
    "Period is ($2^{29}$); all outputs are odd integers. \n",
    "\n",
    "Triples of values from RANDU fall on 15 planes in 3-dimensional space, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.special import comb, factorial\n",
    "from scipy.optimize import brentq\n",
    "from scipy.stats import chisquare, norm\n",
    "import scipy.integrate\n",
    "from random import Random\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# LCG; defaults to RANDU, a particularly bad choice\n",
    "class lcgRandom: # defaults to RANDU: BEWARE!\n",
    "    def __init__(self, seed=1234567890, A=0, B=65539, M = 2**31): \n",
    "        self.state = seed\n",
    "        self.A = A\n",
    "        self.B = B\n",
    "        self.M = M\n",
    "        \n",
    "    def getState(self):\n",
    "        return self.state, self.A, self.B, self.M\n",
    "    \n",
    "    def setState(self,seed=1234567890, A=0, B=65539, M = 2**31):\n",
    "        self.state = seed\n",
    "        self.A = A\n",
    "        self.B = B\n",
    "        self.M = M\n",
    "\n",
    "    def nextRandom(self):\n",
    "        self.state = (self.A + self.B * self.state) % self.M\n",
    "        return self.state/self.M\n",
    "\n",
    "    def random(self, size=None):  # vector of rands\n",
    "        if size==None:\n",
    "            return self.nextRandom()\n",
    "        else: \n",
    "            return np.reshape(np.array([self.nextRandom() for i in np.arange(np.prod(size))]), size)\n",
    "    \n",
    "    def randint(self, low=0, high=None, size=None):  # integer between low (inclusive) and high (exclusive)\n",
    "        if high==None:  # numpy.random.randint()-like behavior\n",
    "            high, low = low, 0\n",
    "        if size==None:\n",
    "            return low + np.floor(self.nextRandom()*(high-low)) # NOT AN ACCURATE ALGORITHM! See below.\n",
    "        else:\n",
    "            return low + np.floor(self.random(size=size)*(high-low))\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# generate triples using RANDU\n",
    "reps = int(10**5)\n",
    "randu = lcgRandom(12345)\n",
    "xs = np.transpose(randu.random(size=(reps,3)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO8AAADrCAYAAACB4cGvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB+v0lEQVR4nO2dd3hTZfvHP+neu006GIVC6YQuOkCgTkQFRVG2ijhQFNfrQqsvCojjp6974EJlOXEwVGjZtLTQMgqllNW998w4vz9CQkOSNkkLUsjnury8yDl5ck6a+zzPc4/vLRIEATNmzPQ9LP7tCzBjxoxpmI3XjJk+itl4zZjpo5iN14yZPorZeM2Y6aOYjdeMmT6KVQ/fb44zmTFzYRHpO2Ceec2Y6aOYjdeMmT6K2XjNmOmjmI3XzCWLIAjI5XLa29sxp/Fq01OHlRkzvY5cLkcqldLR0YFCoQDA0tISKyvzz7Uz5m/DzCWBQqFQG6xMJkMkEtHR0UFFRQXt7e2Eh4djaWmJSKTX+XrFYTZeM/8agiDQ3t6uMcPK5XKqqqooLy9HEAR8fHxoaWmhvr4eLy8vs/F2wmy8Zi4qgiAgk8no6OhAKpVSXFyMSCTC1taW8vJy2tra8Pb2ZtiwYdjb26NQKLCzs6OkpAR3d3csLMxuGhVm4zVzwVE5nlQGq5plGxoaKC0tRS6XIxaLCQwMxMnJSev93t7enDhxAqlUipWVlXn2PYvZeM1cEARBUO9jVd5iQRBoaWmhvLycmpoa3NzcsLGxYciQITqNFkAkEmFpaYmXlxeVlZUEBARgaWl5ke/m0sRsvGZ6FYVCQUdHh3ofKwgCHR0dlJeXU1VVhb29PRKJhMGDB2NhYUFWVhZ2dnYaYwiCQH19vfo9UVFR+Pv7c+zYMSQSidl4z2I2XjM9prOnWC6XAyCVSqmsrKSiogIrKyvEYjFRUVFa4R65XK5+rampifLycqqrq3FyckIikeDl5UVRURHh4eFIpVJaW1uxsbExL50xG68ZE1E5ntrb25HJZOplcnV1NeXl5ep9bEREBDY2NnrHEASB06dPU1lZiY2NDRKJhIEDB2JpaalOzDhx4gRyuRw/Pz8OHz5MeHg4rq6uF/N2L0nMxmvGYFQGq5plVa/V1tZSXl5Oa2srXl5eBAcHY29vr3ccqVRKRUUFZWVltLe3Y21tzfDhw7G2ttY618LCArFYTHl5Ob6+vuTl5ak/90qffc3Ga6ZLzvcUq2bLxsZGysvLqa+vx8PDgwEDBuDk5KTXoORyOdXV1ZSVldHR0YGPjw/9+/enuroaPz8/nZ9bX1+Pra0t/v7+HDp0CD8/P0QiEc3NzTg7O5uN99++ADOXHp09xa2trZw8eZLAwEBaW1vVnmJnZ2ckEgnBwcF6jUihUFBXV0dZWRlNTU14enoyePBgHB0dAaioqNByVjU3N2s4txQKBXFxcQA0NjZia2tLSUkJ3t7eV3zM12y8ZtSc73hShXYqKyupq6vDzs4OiUTCoEGD9BqOalYuKyujtrYWNzc3/P39cXFx0TLytrY2bG1taW9vp7y8nIqKCqytrZFIJAwYMACRSEROTg6tra34+/tz+vRpXF1daWhooK2t7YqP+Yp6WK1hLvXo4+jzFFdVVVFRUQGATCYjOjq6y8KAlpYWysrKqKqqwtHREbFYjIeHh14jl8lkHDx4EJlMhqWlJWKxGG9vbw3nliAIVFZW0tbWRmBgINu3b6dfv34IgoCdnZ3asXWZo/fpZJ55r0B0eYoFQVB7iqVSKWKxmLCwMORyOSdOnNBpuO3t7VRUVFBeXo61tTVisZiYmBi9BqXyRpeVldHW1oZCoWDIkCF4eHjovVaxWExGRgaDBw/GxsYGQRDw9/fn4MGDV3zChtl4rxA6O55UnmKFQqFOhmhubsbLy4ugoCD1nhSUsdfOs6dMJqOyspLy8nIUCgU+Pj5ERkZqzJg5RfUs2Xic6hYpMf1c+M8YMZUVFTQ2NuLh4aFOg8zMzNQK+ai816p9cmJiIq6urtTU1GBtbU19fT1BQUGIRCKampqwtra+Yve+ZuO9jFEZrGpZrJphVckQdXV1uLu7069fP73eW7lcjoWFhdpgW1tb8fb21hsOyq9o4p5vc5Ap05dZn1vF+twq3rypHzfEDdP4DIVCoY7nqrzXqrRJX19fqqqqqK6uxt/fn8LCQjo6OhCJREilUvz9/SkpKcHFxcVsvGYuDzp7ilX7WJFIpPYUd85eGjJkSJeOp7q6Os6cOUNDQwNWVlbdhoNaW1v5JDVPbbid+c+fhWw73cqSSaEA6tTJkydPUllZqd4nq9ImFQoFtra2nDhxgqioKI4cOQKAr68vpaWl+Pv7c+LECYKCgq5Yx5XZeC8Tzi9mB+jo6KCyslJn9pIuOs/KNTU1uLi44OrqipOTE4MHD9b5HlXBfHl5uXJckf5Z8PdDVfx+aBsf3CTBQVpPR0cH9vb2ep1hzs7O6mW+p6cnVVVV+Pn5sX//fvr374+7uztVVVX4+/ubjddM30IQBHUOsYODAyKRCJlMpi5mB6XDR1/2kgrVrFxZWYm9vT1isVgdDqqoqKClpUXjfLlcrl5Gy2QyfHx8CA8Px9bWFiufBv7Kz+7yuuf/WcY1gQ48HCtGIpFoHVfFev39/fHz86OkpARnZ2dKSkqwtbXF1taWxsZGAgICOHnyJGKxWG8K5uWM2Xj7GOcXswMcP34ciURCVVUVHR0deHt7ExoaqpUA0ZnzZ0x9hQPbCup4bUsZcs5gASyItSfcXcDb25shQ4bg4OCgcX64nwuRvg4cKNU0+PPZfLKFLSdbyHlhCCKRSCPWa2Njg4ODA4WFhQwZMoS9e/fi6+uLjY0NTU1N+Pv7U1RUxLBhw2hpaaGtrQ1ra+srbvY1x3n7APpSFFWe4srKSiQSCf7+/nrrYkH3jOnj44Otra3O809UNTPp0yyt10PF9qyZG6fzOuvr6ykqKWXaL5UG39/LibYMclOGmnx8fLC2tkahUJCZmUliYiKHDh1CLpfj6uqKTCZjyJAh7Nq1i6SkJM6cOQPA4MGDL1eBOnPHhL6GymBbW1tpaGigqamJ9vZ2GhsbKSgoICsrS73f8/T0JCAgQKfhKhQKqqqqOHz4MFlZWbS2tjJkyBBiYmLo16+fTsNtbm6moKCAxb9oGy5Abnkrccu20dah3Fs3NTVRUFDA3r17KS0txU8i5uYQw6t+/ru7na21Lvj7+2Ntba0uOvD29qaiogJ/f3/q6+vx9/enqqoKUKprVFZW4ufnp1bjuNLkYS/LR1VfRpfsaVtbmzrfV+WVDQoKUnuKVT9eFYIg0NDQQFlZGXV1dXh4eHQZDlJ9RkVFhXrZKpFI8PHyhgrdM2ibDOLe3MVwT3h2lDtisVjDGbZ0sgd/LN5m8H1/uaeEr/aUsOfpJBTSNuRyOQEBARw5coSYmBikUinW1tZ4eHioH1p5eXmIxWJsbW2vSIE687L5EkCf7KmqmF2VveTt7a3TU5yfn4+Xlxc2NjaUlZVphIPc3d0RiUQIgsAfhyrYfrwGRxtL5iT1Q+JkpV5GC4KgsWwFKK5rZfyHe7u9fhcbEVufTMKq07UJgkBxRQ03Lj9s9PcxJ8yKJF9LRo0axd69ewkLCyMzM5OhQ4fi7OzM8ePHiYqKIj09neHDh9PQ0EBVVRXh4eFdOub6KHqfRmbj/Zc4P0URdMue+vj4dOlJbWtr48iRI7S1teHo6IhEIsHT01PLyFekF/HlrkJqW6Uozv7VEiTw4rUDEIvFep1b87/eztZiw/7M38wKJ8TbXsNzvWRXAwcrZQa9vzPhXhZ8eGc4MmkH9fX1SKVSpFIpcXFx7Nmzh+joaHUq58CBA9m1axeJiYnY29tfbrOv2XgvBXR5igVBoKamRkP2VCwWG1TMXl5ejkgkwsrKCk9PT511sarPuOnDdArrO7SOiYDtTyTg6nDuAaFQKNQF9tXV1Ty8zfA/83UDLPjPNYPx9vbGysoKhULB8KU7DH7/+TyZPIBgitW1v+Hh4VRWVqJQKAgICCAjI4PExESOHTuGs7Mz/fv3v9zync3G+2/RlexpeXk5DQ0NeHp6IhaLjSpmVxm5nZ0dRUVFWFhYaBjv+SmHC3dLKW/Wkfp0lvlX9WPaCE+NUj6JREJRURHevv5c8/EBg+9ZBKx7KIZAT2WO9BfbjvPu9hKD338+TpawatoAHB0caG1tZcCAAWRmZpKQkMDBgwfp168f1tbWHD16lNjYWL3e8z6K2XgvJufLnhYWFmJvb4+dnZ3amFxdXRGLxbi5uek12POT9D09PZFIJBqFA6B0WMlkMvr166eWVu2ccujh4cHyXYW8v/V0l9ctAv64dygBEh+1M+zQoUMMGjSIAyUt3L8m16jv4YF4b8YPsKKmpoZFu9o503Xot1synhnFoex9JCUlkZOTQ2BgIDKZjNLSUsLDw0lPTyciIqLL77QPYjbei4E+2dPjx4/T2NioNlhPT0+jitklEolWMfuWY1V8uauQhjYZowJsGd9PuSRXKTV6eXlpxD0FQSByyXaD7uOX+6MI8nEGICcnBz8/P2pra5n7Qxm1Rn4nImDnE3FYWoiIfzvDyHdrc8cwe568fhiCIFBRUUFISAi7du0iPj5eXWo4dOjQy8lxZTbeC0V3sqeWlpY4OTkhl8sJDg7WO44xxezb8qt55tcjNHdoLoM3PjgCfy8XnePLZDLSj57moXXFBt3XPbHe3DbElsLCQjw9PfH19cXV1ZWoZbsMev/5LL5pCCerGlieXm7S+zvzwAgnnrw1nl27dpGQkMDJkyfVaZ3p6ekkJiZiZ2d3ucy+ZuPtTQyRPfXx8VHn3La2tpKfn09kZKTGOLqK2fWFg+BcMfszvx0ns0yq85z7Evx5/JrBGud3LuX7Ma+FbzKrDLpPZ2sRb461JSH2XIH9r9nFvPRngaFflQZzE/1Yvtv0va8KRytIf24sp06dwsnJCTc3Nw4ePMjIkSM5ePAggiDoTPXso5iVNHpKd7KnLS0teHt7M3ToUK18Xzs7O9rb2wHDitnP/9y6ujoNpUY7e3tAt/F+saeYH/YV8+kET5rP7pM7l/I9PRC+yTQseaJRKvDQP208IStkzqiBKBQKrupnh68dlLYZ9r11ZvnuEuytRbRKe/bMF0RQVlaGv78/hw8fxtfXF5FIREtLi7rq6ErItjLPvF1gqOypWCzuMntJoVCwZ88enJ2dDQ4HNTU1UVZWpqHUqHLE7DtTz93f5nR7/V/PCCNmoKfW63K5nBGv7zT8iwBG+MD8SFvc3Nzw8vJi3CfGJ1/0FnMS/Bnn0UB8fDwZGRlERERQV1dHU1MTvr6+7NmzhzFjxnTpve9DmJfNhnK+p1ihUKif6p1lT8VisTp7Sd84qhmzoaEBqVRKSEiI1nsEQWD5zjN8uuMMUoWAn5MlL8db4ejgoHZu6VpGxyzdRof+yI+aSD8Hvr83Vn09ZWVlNDQ0INg5c/8fhiyfBTr/fl64biBTov2oa2kn+X3duc8XkuGesObRa9i3bx/BwcHU19fT2tpKYGAge/bsYeDAgZSUlKhrly+DpbPZeLtDl+ypVCpVh13s7Oy6NCbQXcyuCgfl5eWpnT6d+WlfEa9sOKE11kd3hXBVkLfW66oEjdLSUu7dZFjsxQL4+GprPNzdNcJTC389yG+HjfUfw4RhHrx153Be/eMQq/YZXj1kCqHetowc5IG3rIoZ143k5MmTeHt7I5fLqa+vZ/DgwaSnp5OUlMSRI0dob29HLBZTWFhIfHz85RDzNRuvLlQGW1dXR3t7Oy4uLuo9aUVFBSKRSO1EMraY/fxw0OnTp9UPgM5pkP/ZUkdJs+5x/Zyt2PhooroyqKysDKlUqnaG1bcruOb97nOPVbw+KZibwsUar0UYUTxwPjueSmL026Z5n/UhsYd7RwcyaogXg7yUVVKCIFBcXIxcLsfT05OTJ08SGRmpLgvMzc3F19cXKysrsrKyGD16NAcOHCA4OLhL+dk+gtl4VejyFNfX13P69GksLS01jKOrp7auYnZVSqAuVKmGoNzPenl5IRaLuWvFYU7XdO39WTZKxLD+/ojFYq0EjcxTNdz7/SGD73/0IFc+vCtCna3V3NzM/Zu10yYN5cZhnmw4Wm3y+89n+5Oj8HTSdN6pVkH79u0jMTGR9PR0YmJiKCgowNPTE2trawoLCwkPD+eff/5h3Lhx1NTUUF9fT0hISF+P+V7ZxqtL9rTznrSpqQmpVEpUVJSWcXTGkGL2lg45B0sasLa0INLPieZObStFIhEhISEaCRd/5Vbw1C9Hu72HaSM8eOGmcJ3X82NGPp8eVNDF31mL/13nRvTQ/jg4ONDaIWP0/zINfu/52FpCu7z78wxl0U1B3BHTT+M1hUJBbm4u/fv3p76+HpFIhIeHB8ePH2fEiBHs3r2bsLAwdbpkv379LpdihSvPeLuSPa2oqKC2thb3s3tAFxcX9u3bx/Dhw7VmToVCoS4cUGkbSyQSrXAQwOmaFu74PIs22bmv5ZFoG6aMHIKzszOHDx8mKipK630TP9rDyVrDZr+dj8chbWtRz5qqGfym97KoNsCB1ZkQsT2r5yhjuBtzy3n65+4fIhcLD3tLtj99ldroOhdwBAcHk5WVRWJiorrCqKioiJaWFqytrampqSExMZHc3Fw8PT3x9/fvy8UKV4bxni97qioCaG1tpaKigqqqKpycnNSe4s57ofz8fDw8PPD09NQqZnd3d0cikXQZDmpvb+emj/fqTP53s7Ng25PK2tSRI0eqz+/cn+epf+opNzB2+kKiPRNihmrM4Hv37uXpLa3UGF99x95nRuFoZ8NtH24jr7oXp9Be4K3bgpkQoSy4kMvl7N27l/j4eA4ePEhQUBA1NTUIgoBEImH37t2Eh4dTVFSkVrs8fvw4MTExfVmg7vI2Xl2yp51TFG1sbNT5vvqewKrZ2MrKSqOY3c3Nrct+O6p9ryAIzNnYhLyLb+SdqywIGRxIRUWFznrd6MXb9KReaONgBenPjlH/OyMjA19fXyasOGngCJrMu6ofjyYHEb4oFSMn8AuOu50lO/5zlVrn2dHREWtra6qrqxk8eDBZWVkkJCSwZcsWoqOjkUqlVFdXExISwu7du4mKitLZ6KyP0DPjFYlE44H/AZbAckEQXj97SADljLdgwQLWr1+Pg4MDX3/9NdHR0b1w3fpRKBTq2liVwapkT1XGcb4yhC5UEjPl5eW0t7cTHBzcZTjo/H47nRMukt7aQWN71z/92WG2zL8hUmeChjHFAyrWTB2ItKmWuro6+vfvj4OTC9d+ZloChY+TFesfTSC2B/W3F5JtTyZhL5KpS/9U3uYDBw4wcOBAsrOz8fX1ZciQIezevZvExEQKCwtRKBRqcfY+iOnGKxKJLIFjwHVAEbAXmCYIQi5njXf9+vW8//77rF+/nvT0dBYsWEB6enqvXb0KXcXsnfek7e3tak9xV7KnnYvZAbWR79+/n7i4OEQiEXnlTTz9cy4VjR14OdnwxoR+WLTWqetvJRKJluDb/sI6Zq8wrO51/7NJ6h9T5xlcKpVy39/G5R5OjfTgpn4yBg8ejL29PalHynjyN+3YsaEsnxbC3FVHTH7/heSaYHfuCZITHh5OYWEhrq6uWFlZqUNJLS0tJCUlkZ+fj6urKx4eHura3z5arNAj400EXhEE4Yaz/34eQBCEpZw13gcffJBx48Yxbdo0AIKDg0lLS8PX17fHV65yPJ3f0U4le9rY2Kh22nQne6qvmF3FoUOHGDhwIO1YM/6jvbRJNWfRET62rJg7UucPQOUMW/pnNr8XGLab+PzWABwVzVozeIdMTswy49IXnaxh9dRABgUOBODHjFOkbDRtCQ3gZgd1JuQvXwwcrGDVlAACAgLIy8sjOjqarVu34u/vj0wmw9PTEwcHB/Ly8oiJiSE7O5uBAwf21YbcPSpM8AcKO/27CIjvfEJxcTH9+p1z7QcEBFBcXGyy8aocT6rQTlNTE2fOnKFfv35UVFSom1H5+fnh6upqVDH7oEGD9Bq5q6srVVVV/HW8UctwAbIr2olcsp09TyXiaKdcip+foPHENcNwcipnVU5Nt/d5/69FeNlZsOXJURr3YGNlSfZzo4zKP26Sws3fnuTnB9wZJnHljpED+f1gCXuL2w0eozN1bWCNvvKHf5cWGUxaVcQHd7rifDaN1dbWVi2NoxKok8lktLW14e/vT3Fx8eWQsKGBITPvFOAGQRDmnt37fgXYA6+fnX256aabeP755zl9+jTLli3jxIkTDB48mBUrVjB8+HCDL6Yr2dMzZ86owzRd/RF0dZwTi8VdGrkq4aKkpISOjg5OCt68sa3rutN7hzsy1hcsLCyQSCRaCRonKxuZ+Nl+g+9908Ox+Lk7aJQJyoG5BqZAdiY6wJnv5sQCELoo1ej39yX8nS35/Lb+6trpkSNHsnv3bmJiYqioqKCjo4PAwEB27txJfHw8jo6OfW3p3PNlMzAB5d73J6ABmHL48OHI0NBQ9bJ5wIABhISEkJCQQEpKCu+99163e9+Ojg6ys7MJCQlRF7N3bpClUoaoq6vD29sbT0/tKhlAr/yLysgVgoBULmBrpfx3554+qoQLLy8vDh06RPiIaBLf3Nml5xjAwRrSnxmj81hrays/7snnjV11XQ/SicmDYOJQR3x8fPDw8MDGxobsU1XMWZtv8Bid2f5UEp6Otr1iwC420GB6ItYFZ0mCiH6ezkRERKgdlp0F6vLz83F0dGTAgAF9LebbI+O1Qmm0TwMPAWJgOjBxyZIlS55//nn+/PNPPvjgA7XD6rHHHmPTpk2Eh4dTXNy1ckNbWxuxsbH89ttvVFdX09raikKhUDuRVGEUVd3ssGHD1O9VzVIVFRV65V8A/jhUzuKNx2lql2MhgjuG2nJDf5F6n9k54SIrK4vIyEgqmqSM/8iwrKOdT8Tj4mCrnsE750VvO9PC4s2GF6CLgIynRmJtbY2lpSUikYjU3GIe/fm4wWN05pZwT5ZMCutRDnNf4dmxYsYOcKBfv37q8JEq48rGxobc3Fzi4uL6WrFCj0NFE4DlgBOwTBCExSKR6KuxY8fek5aWhiAIzJ8/n40bN+Lg4MBXX31FWloaR48eZfny5XrH7ejoYMOGDTz11FMIgsDLL7/M+PHj9YZRMjIyiI6OVs+Ycrlcy8jP51BJA3O+zaFVpn2fB164SmsJlZ+fj62tLa2trRw9Vc6LBmYNjvSGedEOOvOidxwtYd5Pxhnfnw+OIFDsrv73l9uO81ZaYRfv0I+jNbwx1pFH/tFTAXEZ8cnVNowbc5VaoE4ul1NcXExERATp6emEh4d3Wcp5CdLzJA3V3hf4EWXM1y0xMdFn1y7tqpLU1FTuvfdeCgsLWbNmDXfccYfOMR9//HGsra0ZMmQIO3bs4MMPP9R5niq2mp+fj0gkwtfXt9ti9ubmZsrKyvg5p5wVR/SnHf12fySBPm7qtpJlZWVYWloSFBSEu7s79a2tjHnX8LrVnOdHa+3Hm5qa2J5TwDNp9QaPAxDgas1fC0ar//3E6v1sOlZn1Bid+XhKCI/9cOSSdEL1FrMjHZmXHIxCodASqFMpngQHB/elYoVeMV7V3jcIZcx3ukQieXXz5s2Ehoaqzztw4AC33nqret82Z84cvcarQi6XM3z4cLZs2aI2yPPlXzw9PbGzs6O5uVmvkJvKudW5mXRONTzZTc6u2Bb+73o3tSJGXl6eRg6yscoTH9wxhPgB7ur0R1tbW3x8fNic38Dr28oMHkdF5jNJONgpZ/IHvklnx+nunFiaBfSdSRrghJsVrC9oMvo6+gJ3x/lwc3+IiIhQC9SdOnUKW1tbJBJJXxSo65UugXuBcKAYZbhoyl133cW6devUJ5w5c4bJkydzyy23MGvWLHx8fAwa2NLSkhtvvJENGzZw8uRJjh8/TkZGBmVlZYjFYkaOHElQUBB+fn7U1dVR39JBYW0rUrkyLbK4uJh9+/aRm5uLlZUVw4cPJzIyEh8fH5KDvbHq5m9U3g4zfq/Dzd0dBwcHteQNKGfw06dP88V1+pM+zmf+j/ncvTwDKysrwsLCCAkJwdPTk6lJg3l+rPHhs9g3dlHZoDTYz+6OZ1yg/ni2Ev03vOt0E+sLmlgxXb+SZV/mlhH9aWpqUosAqroMlpSUYGVlhaurqzofuq9jsPEKgiADvgaGA0eAtTExMWzatIlPPvkEgEWLFlFZWckXX3zBZ599xh9//GHQ2CdPnqS9vZ0XXniBRYsW4erqSlxcnJZsjIWFBTsrrLju/XQmfLSX6Nd3cNuHu5HJZISGhhIdHa1uE6nCykLE1gUjDbqO2Df38M2OAuzs7CgoKCAzM5Pjx4/j4OBATEwMBxfq9izrIrcexn9dQKtUjrW1NTY2NtjY2DBr7DAWGB49UzP23XRilqYB8O5dwxmu2+luMLNX5rH6nuGEeFwOcU+lIU4ZasFgD+WKSyVQV1xcjJ2dHRYWFjQ3N6tfU6XU9mWM/ctlAz8IgjAYyHrmmWfIzs6mrq4OgOXLl3P99dfzzz//kJOTQ2xsLE8//TRhYWGMHTtW54AHDhzgwQcfJCoqCrFYzFtvvaUzE0YQBHYdLWbloWYN59PpZrj5u9PUNLVpnV9TU8ORI0c4evgAv9ylu4/P+by1tZgZv9fR1tZGZGQkw4cPRyKRqMMLBxeOwdOIFNlxH+zn7b/zNO5nQnQQL17lZvggZ2mVCoQuSmVXxl6W3TKYR5MMW9noY+rXOUgtbPj8rr41C7vZQIi3HcM8LBnuAQ8Nt2X93FAeTh5GcXExfn5+FBcX4+DggEKhoLW1VW20bm5K/0ZxcXGfn32Nqio6P+b79NNPD3JxceHHH39k1apVhIaGEhgYqM6QKioqwt3dnS+++IKkpKRul9HLli3DxsaGuXPnAtoJFwcabHkvo0Hv+30cRfx63wi16qKuliLGhEx+mxtJoNhN57X8kN/Opq67h2hx+KVxiEQitSpEnWMAT/x8zLhBzjIhxJO3pkRyoqKOmz8xPBlEHytnhTL9W+PamejDCujNee2aIW68Py1KLR9UWlpKR0cHYrEYsVisTKw5eZLhw4ezZ88ekpKSNATq2traGDhwoPrY8ePHKSkp4eqrr+4LxQq9UxJ4fsw3MjLyupUrV/Lbb78B8Pzzz6vP/eijj/jkk09ISUnp1mGlorCwkGnTpvHmm29ibW1NfX29RsJFdnEjd6/oXvJ0ywMheluKCILA5Pe3c7zRoEtisJsV70zwpbKyEgcHZSjI1dUVCwsLlv2ezZojrYYNdJa0x0bi4+ZITk4OgwYNYvfpeh7/2bQkDCsg56VxAIx/M43CHuYiXzfMgyGe9ny007CuChcTaxG8n2yDn0SsVRSiUCjIzMwkOjqa48ePIxaLkclkOgXqvL29aW9v5/jx44wbNw4bG5tL3XHVKw4r1b53PvABkHDnnXcSFhbGsWPH2LRpk8a5x44do6Ojg5dffpmYmBhWrFjR5dhlZWX8/PPPHDt2jAULFtDa2kp0dDRhYWF4eXlhYWFBVIALbrbdX/LVnx3hrwNnOl83jY2NakfYkuu8efMGcRcjnKOgTsbElYWEhYUxdOhQ3N3d1XvYl6ck8ISRDqhx72Vw9ds78PX1pbS0lOvDA3hrUpBRY6iQAWGvptEqlbPpmWSeGuNv0jgq/j5aw0c7i/ltrgmbcj301rwmFeChLR3sq7fTyk0XiURIJBJKSkrUy2MfHx+qqqqwsLBQO6lUOfeNjY3Y2toqJXD78NLZaG+FIAjrgQWcdVgFBwezfv16tTaUCplMhpOTEwMGDKC9vZ3777+fpUuX6h136dKlWFlZ8dxzz3H99dczYsQIrSWNXC7nh5mqH3rXX/p//jxDxOJtnDp1iszMTE6dOqXhCBsfG8zfD48w+L4T3s3iQGE1NjY2WFtbq2f1+8cO47M7h3Xzbk3KmqWM++QwZeXlnDp1Cq/2UpaMNl3pIfb17dz4/i7uGzeU9P8kmjyOionLc0gZH4jEcAe7XnrbLfTKn/lc9eY2LaMTi8WUlpbi4uJCc3MzcrkcDw8Pqqur1Qbt4uJCa2sr1dXVan1nVUpuX8RUV2MR0O+RRx5hw4YNzJ8/n5MnT5Kbe27PFBAQgKurKxERERw6dIjbb7+dxYsXaxm5iv/973888sgjzJkzhz/++ENdmKCSPT18+DD79u2jtbWV9CdiERkotnbL92cIGjKEiIgILUeYj6sT/8w13Flz7+p8xrytvWcePcyXfx6NM3gcFQ+lyvhi92lGjBjBrVeP4tBL40yeqU7XthO6KBWpQkRuSjIhHiYOdJZFG0/SohDxzYzI7k++yFS3yol8LY2mNmW6iUgkwtraGnt7exobG9WrGn9/f4qKinBzc6OpqUmdw67KzKuqqtIIC/Y1TDXevUC4v78/AQEB/Pjjj5wf8500aRKnTp2ivr6e5uZmMjMz8fDw6NZB4ObmRlBQEN988w0//PADe/fupba2ln79+hEXF0dgYCAODg5svDeIZC/DLnbMhwd4dq0yz1GV/JGXl8fevXupr69n+yOG/0BrWpUe37pGzSQHP3cnMp8bredd+vkpX0Hs2dpdC5GIAynJPJBg2JJeF6Pf3sU7fx3lp/nJbJgXa/I4AA0dAnd/f4CvZoZ3f7KBOPZSTYBcgJFv7ODbPUqvoUgkws/Pj6KiIvz8/CgpKcHFxYW2tjakUikSiYTS0lJ16aCFhQWenp5UVlaqJ4q+hskaViKRaLGLi8sLqj60KqObOnUqDz30EACvvvoqS5cupa2tDUEQGDhwID/99JNeiZxTp06xfPlyvv/+e7y9vXnppZcYM2aMTsdTe3s7+/fvp6W1gwVG6H5/cZ1yz+Tj44OzszMWFhaIRCIsLCyY9N4OThmR/ps4wIEv7tYobUahUBD+2lbDB+nEP4+NxM9NKT1b19RC0v+ZrkZiDWS/NI6WlhYS387o8fJ1hJ8jCkHGgVLT6oMvJNH9nPj67hhEZ/PfExMTyc7OZujQoRoCdTk5OTg7O9PY2EhISAgikYj8/Hyio6Mv5WKF3nFYnUf2lClT+PPPP9m2TbmUjIiIUBsuQEhICNdffz033HADx44do62tjQceeEDvgFVVVURGRpKdnU1zczOjRo3SMFyVgsaxY8fIzs5GLpcTER7CgReuMvii7/u7DZm9G66urmrHk62tLdbW1qx/Kpn7ot27H+Qsu0+3aJXbtbS0sG5agMFjdOba9zJYmaGcSdycHMhNSWZgd8lUepCidGaVN8v5coIzz4w1LM6tj+ySZg6UtvPrfSN6NM6FYF9hE5GvbWV7fhXe3t5UVFRoCEKoZlyVuODAgQMpKirC2dmZjo4OWltb++TSuSfGW1RYWMiYMWPw8PBAKpXi56f5A/nqq69QKBTMnj2bIUOGEBISQmVlJaWlpToHjI2N5c4778TV1ZUxY8aQlpYGKA3ixIkT7N27l6KiIjw9PYmLi2PAgAG0tbUhEok4uHAMuttKazNj5XGS/2+3huNJxVM3j+DPh2OM+iJCF6WyKeMQ6enpFBQU4OLiwsGFY4j2Nj4E8drGE4QuSkV+dim3/slkVswM7eZd+rn5kyx+OaFgbH8bNhuYadYVt36RTcr4QTj3eKTe5+G1udy/rpTCwkK8vb2pqqrCysoKe3t7Ghoa1K1mfHx8qK2tRaFQ4OfnR2lpaZ90XPXEePfm5+dz8uRJOjo6qK+vZ+LEiRon9O/fn6NHj9KvXz/Ky8vJy8tjwIAB3db4AowfP57333+fqVOnkp+fr05RDAsLU8dwvb29qaysRC6XU1FRwce3uPFchGEXXy/TrzIR6OViVCokwBMbK3llTzvDhw9HLBZjaWnJd/PGsekR4x1ZABGvbeXPA0UAxA4Ss+dp073Ivxxt5qZvT+Fmb0NuSjLhPj3beC7aeAK5lYj/3jioR+NcCMqbZcz6s4HTlXV4enpSVVWl9jaLRCJU2zzVDO3r66vuAdXXZt8e6TavX79eePzxx2lra1OLiKvynB966CFKSkoYPnw4zs7OODg48Nxzz/HVV1/xxhtvEBOje3Zrampi8uTJ6kL7VatWqQW0O6PSpzp8+DDW1tZ4eXnh7e2tLqyPfmO3wffx4R1BJIf203nMFBWKgwvHaKg1CIJA2KtpRo8D4OdsyT9PnHuQ3P3lbvYWmZ6N8VRyf+67ajDHyuq41QiZHn18dUd/7v3xTPcn/gsM87Hn5UR7oqKi2LVrF66urkilUvz9/XF0dOTIkSPExsaSk5ND//798fHxuRQ1ri6s6PqpU6e4+eabOXRIu+GVKcqSRUVFBAQEsHDhQgYOHMjUqVOVH3ZWobFz+iOAvb09fn5+aseTpaUlFhYWRhmeDZCdkqx9g4JA9KtpGOummTDMjbfu1Gxt8upvB1iVbVpTrq1PJuDtpCyXzC2s4o6vDpo0DoCLrQW7nxlDu1RG/Os7elzfe2+cmJpWGesO9V7Dsd5k51MJ1JSVUFpayogRIzhx4gRRUVGkp6czfPhw9W8qMjLyUqzzvSAOK4OYOHEiK1asQBAE9uzZg6ura7eqkgEBSofPzJkzWbt2LXl5eZw4cUKdbOHi4sKIESMYPHgw/v7+VFdXazieVE/P3JRkbh1g2HV2oJxlpWflZRsaGsjLy2PPnj18PsGZm/obd9/rj9YRdt7D46WJkfz+oGnZS2P/bw/vb1amUYb28+LgwjH0069F0CUN7QrCXk2jorGdnJRkbhnWsx3sV3vL2X2qnlUzh/ZonM7Y9+Ivc9TbezjdZIlMJsPV1ZW2tjY6OjrUy2lPT091m9e+tHTu8cw7bdo00tLSqKqqQiwW89///lctiP7QQw/plMiJje0+/lhRUcGaNWtYtGgRAwcO5MMPP2TAgAHqJ6OlpaV6hs3MzCQ8PFyv0Hpraysxb+4x+KZGecCjY5TC6l5eXurc2Yjhw4kzYhwVfz0SS4CnpoHEL0rFwPRqLTovy7/ffYLFfxtZIdGJYG87fpmXSEldK9e+Z/y9nc9T0fDBPoxeqVwMbugHb8wcrW5P01mgTlX6eQkK1PW9XkXPPfccfn5+FBQUUFtby7vvvqteEquWxyqKioqQyWQMHDiwyzGN3b/mdlpGHzhwgAEDBuDq6so1y1IxNtx5fZAD707XjAn/knmchetN06X69YFohkqU24aWDimxr/esRUn2C2OwsbLkyZWZbDS0akMPiQE2jAvxZWkPHioXivdv8mPciCC1QN2hQ4fw9/fHzs6OQ4cOMXLkyEst5tv3jFfF0aNHueuuu9i5c6deZ4KqxG7kSN3dDDpz1aJUjNmZfTllEAkhA6isrKS6ulqtXvnPoSIeM6EaKPe8fXWHTMYII/sTqfB0sGT70+ecWXd9vI2DlaaHPG4dDLcOcaTRwolH13etW20I6x6KZtIn+3o8DiibZPVGMGfqUCtemT5WLVCnKl2NiIhg7969agGIS8hx1XeNF+DGG29k8eLFGrKvnWlokzL6jR0aWUTnG4mqp1FZWRmZRyt5w4jSVQvgwItjSU9PJz4+XuMPa4o3+s2Jg7hphOZm/Mlvd7LxpGnCyDkvXIX12bTTgsp6bvnYdINxs7dg13/G0iKV9krDsfsS/amqb2JdrnHie72PUtdrXowLs0cHaQjU7d69m5EjR1JRUUFTUxPDhg27lBxXfdt4v/vuOw4ePMjLL7+s8/jopanU6HGZ7l4QTWlpqbqZtq+vr7rdo7GG99MUiTok1ZlZH6eSVWnUUDhYQuZCzQdMQUUDt3xiuFJlZ55M7s/cq5QhNblCQYSJKZoqtj+ZgKeTPZ+mHuN/23tW32trAW9PCWP+GtO6F/YWfs7WrJgZSl1lGeHh4VoCdb6+vuzevZukpKRLSaCubxtvc3MziYmJ7Nq1S6c8Tncx1G0PR+htkXLLa6kUGJGX/nQ0zLlZO6TU0NxCwtvG5yLvfSYJR7tze6wzZ84w/usCo8cBsLaAnBfPXdvL63L4wYCeSfqYGi0m5eZQCsrrueXTni9/X74xkK93nOZ0Y+8UArxxlTXXjErAUiSQW9rE8YpGKhql7C6oorJFhreDJcX1HQhyOZMivYnzaOeqhDi1FOzJkydxcnLC3d2dnJwc4uPjOXz4MD4+Pvj5+V0qS+e+bbygDBvde++9JCZqZxoZMoMujoDbbtM2OoCGlhYS3jLO8M5flhtzLedzV6Q7L986AjjX/uVgmytvpxneaaEzu58ehb0VlJeXc6qwiAf/MT2pwxLYt3AMVhYW3PVhGodMfxYAEORlw0OjBvD0OtPUQ87nrYmDmDBCdzxQ1WFy69atXHPNNWRlZREREUFJSQm2trZ4eHiQm5tLbGwsmZmZhISEIJVKOXXqFFFRUXqF/C8y/16ct7eYMWMGq1ev1nhNJpNRUmLYD3zhQf2G5eLgoNcY9RG6KJW6Jm3t49yUZGYZUdwAsOZArfraVMX+U2MCyHjG+BJDgMS3drL8950IgsDImGhyU5KJkphWKSwHhi/exu6CatbOT+ajO0JMGkfF8aoOnl6XT/ozSfRCrT9P/3aCUW/o3iKo9rU2NjbU1NSo47qq/3cWqAsICKCoqAhXV1daWlrUlXCXMn1m5pXJZIwYMYJ//vmHuro66urqaGlpwcfHB4lEQqwRM2dXhrp2x0Fe2VJl8FgeNrDjOe3x2js6iDJCqF1FxjNJNNfX0djYSFCQUjXEMKF13RxYOAars3HLzJMVzP7W9H1npK8Dq++Pp665ldFv76Gni99nrgnkYGE1G47pFxU0hoxnRuFoa61u66oS63d1daW8vFzd8iQpKYmsrCxCQkKoq1MqhQYGBrJr1y6SkpI4ffo0FhYWDBo06FIQqOvby2aFQsGuXbt45JFHqK+v56WXXuKmm27C2dlZw6lgzJI12ga+02F0YFousr4HwtSP0zhQadzX5G4Lb42xIyEhQX1/B85UMfVr01Iiv50dScxApdBzc1sHcW8Y/1DpzO5nRuNqZ80Tq/ey6VjPOi84WsF7d4Zx38recWZdHwAPJChF6lSa33K5nMzMTGJjYzl27Bi+vr50dHTQ2NjIoEGD1AZ99OhRvLy8cHFxYf/+/cTHx2Nra/tvO676tvF++eWX7Nixg5iYGDZv3sx3332n99yffkvlpWzDx+5qFjZ2//q/WwdxXaT2/iuvvJbbPjXios6y/ZFIjZamLS0tRq0wOqMKAamYtyKdradMm80BHhwdwIKrh5CeX869q3ouGfvZtHCeXH2Ipl74RbnaWbL1yVHYWClXHAqFgsLCQiwtLXFxcaGwsJDw8HC1Z/nw4cP4+flhZWVFQUEBUVFRZGZmIpVK9YpBXET6tvGqUCgUxMTEsH79enVRgi4yMjK4Z6PhkhhdGfC2Y4U8tNq4Dn+96cxadMNARvpaUl5ejqWlJb6+vryddpr1+cZJzqpIXRCP2FVZebVrXy5z/zA9GcPaArIXjkOuUHD1G9uo6mGFwzVD3BnkYcvn6cb3c9LFkonB3DrCD0EQaG9v58CBAyQkJKjjuseOHVM3RS8sLCQyMlLdmPvEiRPU1tYyatSofzvme3kYLyhbqvj4+DB79my956j2LDd81TtGB8YbXs7zo3X+0Zf+msG3Bwx9sCgTCyyArGeTNNL2yhuaSX43w6hrUjFluCf/nRRJe3s7OTk5PLe5mbIeGN5f8+MI8HDi132neeGPE6YPdJbfH4zhlk9Ni3efz6QIH5beFoZCoeDgwYMEBQVRVVWFtbU1Li4uarF2lUGrOinU1NTQ1NTEqFGjsLe3/zeXzpeP8RYUFHD//ffz+++/6/1CVU/ZuDhlIbwxhteVAT/4SSrbKwy/1kmhTiy9Q7sYv6q+iTH/22v4QGdJfTwesYuDxmuPfZ/OPwWmLX8PvahMExw6dCin6zu4Y3m2SeMA3BzmxRu3R1Dd3MZVbxteS62Pn+ZGsWDVfop6qaXwzqcSkTY3UFtby6BBg8jJyWHkyJHs2bOHmJgYioqKsLa2RiwWk5mZiUKhwMvLC1dXV/r16/dvFiv0/VCRikGDBtHR0aFXSgfA1tYWS0tLWlqUP2pjwkChi1JZv163sX/6UDL7n00yeKx1uU06Hxxerk5Gh6YAkt9NZ/pyzT3vezPieee2IUaPBRD+2lYaLZwoKSkh1M+drOcN1wI7nz8OVxG6KBVbC+X3Pap/z5L7b1++H5EFLL6+Z0LyKka9vZslqSXU1NRgY2ODlZUVzc3Nao0rleKkjY0NFhYWuLi4EBAQQElJySXblKzPGa9IJGLatGmsXbu2y/NU8iYqclOSSTJQ5OrpTP2zta2trUkx4fKac7m9LS0tFBQU8OX19jxmpCxydolS9K7ziumGiAD2P29aTPihdUXc91MRgiBgb21Fbkoy/Ry6f58+Rr65m292nuDze5JYPqNnkrGFjbDwr2I+nxGOWy9EbP7Kq+W+v9s4XVSqjvWqjNbOzg4rKyuampqwsLBALpfj5OSEVCqlpaXlkoz59rllMyhVJm+88UbS0tL0Lp1lMhmZmZnEx8drndNby+gv/t7L27sND5V4WsP/XeukdjypnCWmyuR8dGcI44ZJNF57Y/1Bvs40PE7dmUMvjlV7VrceKWXeD103Je8KL0crtj11FTKZjEgTq6Y6ExPgRFSAI8v39LzaCeDRsf0ZYVNBUlISOTk5DB48mJaWFurq6qisrEQkEpGYmEhJSQkdHR0MGTLk33Jc/bvL5jlz5uDj40N4uO4nsSAIPPbYYwQFBREZGcm+fV3n0Xp5eSGRSDQ6NJyPlZUVjo6ONDRoJwAYu4weq8fY51wby5fXGy5nUS2Fuzc0ER0dja+vrzoBQCRSdjkwVmv94bVHuOadNI3XnpkQYbJKZPhrW1mxS+nkGxvia9LSXkVVs4zQRamU1rcqM7x8evZTyypq4os95Wx/IqFH46h4f+sZXt7TptEOxcfHh/Lycjw8PNSNuX19fSkvL0d2VmHlUuKiGO8999zDxo0b9R7fsGED+fn55Ofn89lnnzFv3rxux5wxYwZr1qzp8hzVfkYXuSnJfHOjYaLIleierUUiEWKxmL/vG4oxdhe6KJX049ozyJcPJrP+Yd2C9PoobVR2cGiTnnMX+7o6kpuSzDDjsjQBeP2fQoa/du5ec1OSeeIq0/edN3yYyVt/5rBwjDdvxfcsJ0sArnpnD4+PG0CkX8+TK4sa4YZPcnA729MIlA99Ozs7tUFbWVnh5OREXV3dlWm8Km1nfaxbt47Zs2cjEolISEigrq6uS4cUKNupbNq0qUu9XQ8PD7U+ry7i4uKMnoXPR9VG49cF0Xw03tHgse5dmatzvIFeribNeNFLd/D+ljyN135+NJkPTMhFliqU93qktA6A+5OH9siZ9WVWDVN/Kkfi7c7uBdH4mqi9peLdtNMU1nSw5BbTuit2plkGEa+msbvSWj3DqlrLnt+Y+1LTdr4kHFbFxcX063dOelWldt8V9vb2xMTEsHOn/lQ/CwsLdae4rjDWgFVG19LSQklJCQ0NDZw8eZLwQYM49OLYbkbQHk/XjyI3JZkHRxrYjOksH+8o0RK9uzpUwq6nDfeQd+b2z/ez8BdlP2SVM+uuEZ7dvEs3MgFm/lHP75kFbP5PMq9OCDRpHBW1bQpe+P04q+6JpDdqf5ZnN/JZWh7u7u60t7erBepKSkrw8PCgoaHhkitWuCSMV9cXYkhQfMaMGd16nU+0WHHNx4cIXZRK2KJUcs7oNuTclGSjjfjo0aM4ODgwZMgQXF1d1X2Ec1OSucq7+zFURCzexrt/HtB6fcH4CDL/Y9weTzh7beV15wKkbg5KD7nEBC/yLwdrNB4wL0+M5O9HTW9gtmRnPbO/zOD22IFkvzCWHji2AZj29QESg1y5bkjPezj8cEyGu7u72gMtFovVEQuJREJ5efklNfteEsYbEBBAYeE5ITZVp7fuuPrqq8nIyKCtTXe9akFlE0/9fkYtjyOg/GN35W02xoBn/lHPL78cVTs1Oj+EPp2XzL5nDO9y8FlWtc7rcrC3JzclmaFGLjWT38vg7q80E0G2PJ3Mx0b2ElYRsXgbWaeVDz5fV0dSHzS9BUtmUTOhi1LpkErJTElmWpRxK4zz2Xq8nrTjjfx8fzTOPQgpKQRlgo/KV2JpaYmrqyu1tbWXZFuUS8J4TdF2BqX86w033KDXGfbeFv2peqGLUvUG33NTkok2MBvugyqIXLIdBwcHLc+2nZ2dSTHhkuo6rdd/eXocLxnny2JvYZPWsnzsMF+jl/YqZn1zgJjFqaSnp1NbW8uux6KYMtz0RsAj39zFO5uO8NItEaw3sS2MCqkAkz/fx+uTw9h0XzAPxzgS5m2LpxF+raEellRUVKh7/TY0NKh7/KpSJJuami6ZlqAXJc57obSdAbKzs3n55ZdZtWqV1rEHvtvPjhN1Xb5/Vpg9z9+uXJqqxNY7a15NXtu146wzP9/pq1ck7+99x1nwh+Eyr0M9rfn1Ec3Ei1OnTikfWEbmbAN8PjWYUUM1VzMv/7yfHw7VGT0WwJ7/jMLFXrnbLKlt4tr3jU/3VCECDp6NMY9akkptDxOanG1EvDvWVqOkEqClQ86v+4r5bHsBFefVdQx0gvdvH0JbUz39+vVT97/qLFBXWVlJQ0MDISEhFzPme/nkNmtdgCAQFxfHunXrcHfXjI3kldZx2+eG9eP5fUZ/KisrcXJyQiKRaGheGZPU4Q/83YsFDp1n7ra2Ng4dOkRsbCwPLN/KjhLjZgB7K8h6QfPaqhpbGPOOaWWGt0Z4suS2cylid76fyqFak4YC4IsZ4SQO9u5xgoiKfx4IxU+iP4inUCjYv38/4eHhpKenM2TIEBwcHCguLtYQqDt9+jTW1tb4+fn9GwJ1l09u8/mIRCImT57ML7/8onUs2NfN4HFu+f4MLv0GEB4ernY8qchNSeZdA/1GxXRtoKYsow+cUWZM2dnZYWFhQUtLC5/NHcu2x41barae7YxY13KukMHLWSkB5KZxpmHP5F8PVjP81XP3uvbRZN67w7Q9NcB93x9i1JvbGBviS/YLV/X4x3ntZ7ks/FV/kb9IJMLX15dTp07h4OBAWVkZ7u7uNDQ0qFuBVlRUqB1YlpaWeHh4UFVVdUksnfv8zAvK5eTdd9/N+vXrdT4NjZnt4r3gq4d7Z+Y8/NI4vU/nmsYmRr9j3FIzNyWZ0tJSWlpaNDonRi1KNbq9yNWDnXl/eoxay7qhoYGDtfBupmlidV9ODyUhSDnLldY1cM17PSvp2//8VdhaWzH/+71sKeiZWoerrYhdz4zV+lsIgoBMJmP37t0MGjSIsrIyhg8fTnFxMXZ2dri7u6sF6rKysggODkYul3PixAmio6MvlkDd5TvzAgwYMACRSKThsVahUCj46z7DG2ClV/XezBn2apresTycja8sCl2UiqOjI5WVlRqe7f0pyXw4JdiosbYUNCqbjVVU4O/vT0JCAg9MSOTwS+OMGkfFnJW5XP9/yhxmXzcXpdNPYnoZXdTS7bzwyyE+mBHH6ntGmDwOQH27Mnf8g1RNxUqRSISlpSWCIGBjY6NRrHC+QJ3qmIuLC21tbbS2tv7rMd/LYuYF+PTTTykvL+fpp59GEATq6+spKyujtrYWDw8PJBIJif8zTnv4g0n9uXq4dm9gFb1V4PDaz7tZecjwGS/YAb6aG4Wbm5vG63K5nIjF2wweR8VnU8MZPVQzMH3bG6nkmagYm/XcaOxtlA6d/NJ6Jn1uuuaznSVkvTAOgKSladT30Jk1xMuOdQ8nIpPJqKioUGfy2draEhYWZpBA3Zkzyn7EgwcPvhgCdZevw0pFTU0N48ePJy4ujltuuUVtsJ0dT3K5nOe+2MafRqqs9JbCxuMSeOAB3WOZYnj6ruvuT1LZa4RoAECQuzW/PXrOu11UVMSxslrmrzetQmnJzUHcGn0ua276JzvIrjBdrmPvc1fhaGPFh//k8uGunlUWiYBPrrXDX6IUqrO1tSUzM5O4uDiOHDlCv379aGtr0ylQ5+npiZubm7pR2UUQqLu8l81ffvklt9xyC6Wlpbi7uxMbG6vT8WRpacm9o8VsezjCqPH1pTCC0oDuMDDF9t0y/cZuaWlp0jJ6/wnNH7JCoeCN20J592rjiuGP10qV93nWESMWi3GUt5i8jH7hj+MaqZorHxrNazcNMmksgLjXt/PIqv08cm0oWxbEd/+GLhCAB/9pY+WRNhwcHNQFJqWlpWr9ZpWzqnOihipt18bGBjs7O+rr6//VpfNlYbyDBg3il19+4d1336WxsREHB/1JdxKJhIqKCnJTkjFmRxaxeBvLftEt77JouvGplfrITUnmo9sN7+Q947tcohalUlNTQ25uLunp6TQ0NJA4PMKkAoeI17by5c4CdaJCY2MjuSnJPJPcfcbb+ahSNUtrlC1DJ8cMIDclGVPrElLz6whdlIqXkzL5ZWxgz1IiV2aVk/TGVrXxlpSUqL3NgiCoPcuqRA1nZ2d1Y26Vysa/mXF12SybQZnaFhMTw65du/TuRVRZXCNHjsTS0pLtR4p48AfjWm/01jLaG9jaizHhbQ9H4OnpqbWMW/LbPr7LNq5LnwjY+nAEVVVV6sSTwuombvjQtGSMYT72/PzQuXjbk6uy2Jhvuth6cpA7H04fwb5TVcxcYZqedWe+nT0ci9rTDBs2jPLycmxtbXF2dtYrUDdgwAB27txJUFAQ/fv3v5BL58t72azC1taWxMREtm/Xr9wgEonw8vKiqkq5l7sqJMCk5WpxjW5jyE1J5qsbDEu311cn3HksY0pyx3x0kNf/PKT1+gsTo9lvZEmfcHa804Wl6phmP0+lh9zVhN/p0YpWQhelIj2bkvp/02L47C7T86NTjytbxAzv587m+4dhoMKRXmatyOHrXJk6r76kpAQXFxdaW1uRSqXq0s/ONeIuLi6cOnXqX4v5XlbGC7p7Gp2PriL93JRk4v0MX9Bd98E+rtJjePHx8WydZ7h+U1cGvDMlmW2PGZ7U/O2+Kp3j2Z4t6TOWJ3bD3Z9rll3ufimZL2eaZnjDl2znzwPK/lKjg8UcfHEsPVn8RizexonSSv55Mo63JvWsvvfvgiYmrSqiAxEikYiWlhadAnWqvGeVusa/JVB3WS2bQem1jYqKYsuWLV3ufTMyMoiKitLKUW1qbWXkm3uM+kxdRqFQKEhPT+e+v42Lt6jGUigUVFVVUVZWRktLC2KxmJu+PWXUWHueTsDFQfuBtCmnkCfWGZ8frathebiJfYAdrCHzeeV4lZWVfL09j68Ome6NDhE78tODI6lubuWqt437++nivhhPbh3myMCBA8nKyiIxMVGdqNHW1kZpaSmNjY1YWloSERGBm5vbhVo6X/6hos48++yzhIWFcccdd+g9RyXM3lkEoDPG7je/nDqMhKGalVB5eXl4eHjg7e1t1Hg/TZFQV1eHp6cnvr6+ODk5qX8Yy//O4v92G75XlDgoSwHPx1TRux/mDCcsQLOS6I73U8k1Mad5139G4WJrRXp6OoEhEVz9nukFDp0blo9bmkoPIlMAuFrDnheuJicnh6CgIJqbm6mvr2fo0KGkpaUxYMAAbGxsaGtrY+jQoReqWOHS3vNu3LiR4OBggoKCeP3117WO19fXc8sttzB8+HDCwsL46quvuhxv5syZ3epbdS601kVuSjLPjDHcwzpn9VEtA+28PDdmyXr7D2UUFCh/EOc3U5t7XYxR4ZuyFv36W7kpyYwPMs73O+XLHCJe1Rzvx0eT2fK4aeGbpDd38vB3+/D09MSio4XclGT8TJR8bpEr7/XPnGLSnk9mxWwjdXXPo14KIf/dAg4e6vBRZWUlCoUCQRCwsrJSF+n/GwJ1//rMK5fLGTp0KH///TcBAQHExcWxatUqQkPP7amWLFlCfX09y5Yto7KykuDgYMrKyvTmlgqCQHx8PD/88APe3rolLcobWrn53T2o9CZGD3Dgs7u1f4A9SZ4QBIH09HRiYmLUT+Xnv01l3Unjx9KFsauDz6YMYXRIgNbrtc0tjHrb+Mqizu1DTb2mznx5vT0ODg60tLSQUSblkwOm7yVdbETseW4cAMlvpFJuen9xAO4aYsEr08eRl5eHnZ0dNTU1KBQK4uLiOHjwIL6+vvj6+l6IpmSX7sybkZFBUFAQgwYNwsbGhqlTp7Ju3TqNc0QiEY2NjQiCQFNTEx4eHl2mpYlEIqZMmcLPP/+s87ggCFzfyXABdpxu0Vmgr0qeMEb0O3RRKs9+t00dP6yoUKY7yWQyHr0mmO9u1t8kTddY+shNSeanOYYLzD3wQz5ROsZzdzS+uThA5OJtfJymGWbLTUnm2bGmqU3O+auVOpEz8fHxPHbrVex80nAlkvNp6FCqah4vbyT1mWSmjzBCl0gHa/IVnDwr1n769GkCAwMRBEHdmPvf6KzwrxuvIeJz8+fP58iRI/j5+REREcH//ve/bp9w06dP56efftK5lKlplqJvOxS5ZDvbjmqL3+16IZnvZhie/P/7CTmhi1Lx8fHhzJkzHDx4kMzMTNra2hg2bJjRSR0Reow4JEDC7gWGe6Pb0a8ikpuSzHNXG5eM8f62Iq0HzN1jh5L9gmlqkw/9copxbytDfe5nkzHmxJlueBM/zWTR77m8ODGcwy+N49EYe0xtxLL7yBns7e2RSqU4OjqqixXc3NxobGykvb39oi6d/3XjNUR8btOmTYwYMYKSkhKys7OZP3++TjH1zvj5+WFnZ8eJE9pSOPY2Xd/2Q2uPMWaxtrFED/YzeoaKezuDM+VtiMVi4uPjGTRokNoLnpuSzPe3uBk0jhz9s7CLiwtfXm9PgBHOzsgl23lipbZzaPboYNKN0N5SEboolZKqOvW/bayUoankQYbL4aqobFE++Jo7lA+Yp28M56/5MUaPo2L1/nKizv49Z4+L4IublfK6vz8UQ7iX4U4ma5HAmTNn8PT0pLS0FLFYTHm5Mj1V1V7nYmZc/evGa4j43FdffcXkyZMRiUQEBQURGBjI0aNdKy2IRCKmT5/ODz/8oHXMwab7NXDVWedHu54ZaoTY8HX0fw/AjM8O6wwlREVFGT0Ln/8DUS3Pv74vmIynDXccbTquuxGa81ntLYmRBTPXfrSf1/48V/wul8t56dr+vDnGtLrXuNe380uWsoInwENZZuhlYhFPu1xZolneJEUul9PW1sZgHxfWPjxarRz69axIHPTkzIb7OhIf7E9hYSHBwcEaAnU1NTUaAnUXa/b91x1WMpmMoUOHsnnzZvz9/YmLi2PlypWEhYWpz5k3bx5isZhXXnmF8vJyoqOjycnJwcura9XB+vp6xo0bx/bt27WW2Tmnypm2wrCO7k8nBzDnKu1OfFX1jYz5X6ZBY6joTQdU57FaW1vJzc0lJibGpLE2PjSc/j7aYnK7jpUxd/URo8YCWH2rJ83NzXh7e+Pr64ujo6PJzixL4GCne1255wSv/XXapLEA4n3tSLlWQmCgfu3oyvoWnlq7nxaZiOuHuXOVrwW+EjH79+/nmmuuIScnh8DAQBQKhboxd2ZmJsHBwRqVbL3ApR3nXb9+PY8//jhyuZw5c+awcOFCPvnkE0ApUFdSUsI999xDaWkpgiDw3HPPMXPmTIPGnjJlCo899pj6R90ZU41FJU5WVlaGVCrlno3GNZF99Xpfbk/QLxdjap1wVlYWoaGh2Nsrwz9vrUvnyxzDe/f62sHmZ3Q/XEwxvDdvDuKmaM04+lfb8nkzrcjosQD+eCiGQT7KRMh2mYyoHjYw2/tMEo52unfAgiDQ3NzMsWPH1Pnyrq6utLW1ERQUhEwmo7KyUi1QFxcXR3V1NXV1dYSGhvZmzPfSNt4Lybp16/jrr7948803tY6dqW5g/IfGybV8fL0jtijw9vZGIpHg6Kjc032Weox3t3fd5eF8emsWfqo/3HdPMiUlJbS3t2vNKD2Z0TvzzJp0/sgzrpG3rQXsf1FzvFapjJilphleiNienx48V+Aw67OtZJWZnlt8f6I/T1ynW2lFoVCwb98+IiMjOXXqFBUVFYSFhVFWVkZYWJiWQJ2/vz+7du0iMTFRLRXbC1y6oaILjaoVqFSq7V/u7+nCTaHG6Q7P+6uZlflyBg8erDZcgAeSh5pU4NBVnfB9Bqbqvn0GtWf7fPF31VjGXtfREu0i/Dfuiid1fpRRY7Wf7XvU2nHu+1e1Trl2sPHFgUfKW9WC7QDfPjCWtXOGGz2Ois93FzPu7R06j6kE6kpKSrCyssLCwgJPT0/q6+uRy+V4e3trCNSpjl8sgbrLfuYF5Z75hhtu4Prrr9d5PGXdYX7MMVJ6Av0Cc5GLUjEm4jdlqAX/napfCN2YmXPNbV70798fV1ftWHJ+cQWTvtCvpng+9iLIeilZnWddWlpKW1sbMpmM+zd3GDyOinmj/Hj0Gs1w28HTFdz1jeHX1JkHR/uz4GrlrCmVShm+VLcRGsq6B6MZIj73vQmCgFQqZd++fYhEymKFqKgotQi7m5sbR44c0RCoUygUHD9+nOjoaGxtTQ1KaXDlLpsBduzYwYcffsgXX3yh95yyikqu/kS7nK47dOU0A+w+Vsp9q43THu6tZbQFcKgXHWOfJlshFovVedZlZWU0Nzezs1zOss0lRo0Fuu/TVGdW50wq6Lkz644RXiyaeE5pRaFQcODAAeRyORKJBJlMhq+vLwcPHmTkyJFkZGQQERFBQ0MDdXV1BAcHs3v3bqKionBxcemNpfOVu2wGSEpKIjc3l6YmTQlRuVxOeXk52dnZnDl1kk+TrTA2EjFn9VFG6PjhJQ41vjl16KJU6pp17ylzU5JZc5thPX0UdJ+Zda0RSVAPpsr4/lCLOs9aleM7O2koBxeOMXygs4QuSmVHnnZJ5v0JEqPHUmVSFdUonYbTEwax59nR3bxLPz9mV2kksSgUCmQymXoJXVZWphZcb2lpUXcS9Pb2Vuc9X6y+RlfEzAvw4osvEhgYyNSpU6mtraW0tJSGhga148nJyYnCwkIUCgU1CgdmfGv8LLzvuVHY6ci3nvbpDnLKDS9xkQBbemnmTH8iFmdn3RWzLW1txL6hW9pHH6oHUm5urlrgD+Dej1NJrzRqKJytIf15zftsau9g5DL9bVu7IsrXke/vH6n+94xPt7O/3PSUxdkRdoyVgJeXF9XV1cTHx3P48GEGDBhAS0sLzc3NBAYG6hSoy8zMJCEhoTc6K1zZy2ZBEFi3bh2vvvoqDg4OvPvuu0gkEq0azI6ODrKzsxk5cqTJtaovXhvA9CTtmHBTWzsj39hl1FgXKiYM58JdpaWlzPzDOImc1XeH08/FipKSEo14fH55DZM+zTFqLIBDZ/sUdWbkolRMlVrv7IvYeayc+1cbFs/XhZejJWlPjOb06dPY2dlhb29PaWkpISEh7Nmzh6SkJA4fPoyfnx/W1tbk5+cTHR1NdnY2AwYMwMfHp6cx3yvXeKVSKWPGjCEgIID8/HxWrlzJwIED9Z6fnZ1NUFAQTk5OAFz/7jaKGoxf/ugzPGON7rvZoUQP1N9vx5jxDr04Vt0hobGxER8fH3W468eduaRsNlxS1cUK3r3anri4OCx7obJo4fWBzEgYqPHaqj35vPqXaTHhH+aOIMxPKSLUIZMzYonxetad+fP+4VQXn2DkyJHqENHRo0cRi8VYWlqqEzX27NlDVFQUjY2NlJWVERkZ2dOY75VrvKDMPrK3t+edd95BJpPxyCOP6D23vLycxsZGgoLOxWlyzlQz7Wvt5tfdsfXxOLxdnLRef/DTVLYbKT3cW7NwEvDWY1G4urr2uDUMwD9zg3X2Un529W5+P2ZcHZ6dBew7LyZ8NC+PyWuMd4oBRPk78f195/o5PfJtBqknjUuo6UysD3x2TyIlJSU4Ojri4ODA6dOniYyMVAvUqRQlBw4cyM6dO4mPj8fR0bEnS+cr22GlyjiaNm2a3kojFV5eXlrtRIb39yQ3JVn/t6iHse/uZcr/tI3hw/uu4kMji266c0DNNbC/1y4g6b39en9MxjrZrl2ex6rtx7ReXzY1kd3/STJqrLazMeGy2kb1a/5+fqyc6M6o/sbHhPcXK3O3a5raUCgU/PfGQawwsBBEF5kVEPvGbnU7FFdXV5qbm5HJZDoF6lTloBcq5ntFGK8K1T43P1+/1KulpSUuLi7U1dVpHTucksyDid03/dZ4T7224VlZWeHp6cmex2MwovczoYtS+WiT7palT96ZTM7zhntZu3sYbHnM8A6Er6YW6xzP1d6W3JRkBhlpd1e/n8nkT5SONGdnZzo6OvhwejQb54/s5p26Gf1/u3l0+VaqqqoYMmgQh18ah69hAp9aKIDoN3aTVdREa2urlkCdtbU1jo6O1NfXqz3RF6rO94pYNndmxYoVHDlyhJdeeknvOSXlVaxJP0Fhux0x/V2ZET9A43hDawcJbxrvEV1+ZyhJw5T718rKSmpqaggODua3rAKe+/OMUWN1FqqrrKykrKyMtjZl6eGEFYZLdZyf9H8+xi6j9z+bpDM5Yf/JSpM8+LkpyereQP379zfpmlTYWsL+hefuddXufF7927Q9NSg1s/Y8exX79+8nISFBQ6CuvLycsLAwMjIycHZ2JjIy0lTH1ZW95+1MY2Mjo0ePZufOnTq/zHaZnLs+38uxSs3W6WlPJuLjpDlPhi9KxdgFkRewLSVZrS4ZHx9vUhNvgC8nOCNIpXh5eSGRSDRCQmfOnGH81wUGj9XVcvmVb1NZa4R0z7hARz6apXuWNMXwUq7rR6CFMlSj4vU/slmxzzTVu43zYujvrSxwaGxtJ/5N46IA57P8RldGhA3TEKhTObWOHDlCU1MTSUlJpjquzMbbmRkzZnDfffeRkKDdMXtrfiXzVumeIcYF2vPRLM33/JJ1hoV/Gm4kKnKeH01BQQGenp4aOlszPtvD/rLWLt6pib8I/n6pd2bOrfPC9Wp+yWQyIo2s4tH3QJiwLJVTxjYURjtmXVhTzw0fmNaBcNRAZz6fHav+9/7TNcz4xvgwl4oPbxJzdWyYupNgQUEBDg4OnDp1CkEQSEpKMjXm27eNd+PGjSxYsAC5XM7cuXN57rnntM5JS0vj8ccfR3p2Jtq6VX+Mdv369fz666+8++67WseW7zjJ/2051eX1nJ/TLJXJGG5CedqMSHtuGeRAZKSmymFlQxNj3zW+8bYuamtrGfV+dq+MBSaEumYGEz1I2xtdVdfAGBMacG9/MhFLmVI3WSWP+58/isk3MSic/fxobM7OiGVlZdTV15NWIvDe9tJu3qmJlx1sf/Zq8vLy8PT0xMnJiX379uHt7Y1UKkUsFuPn52fK0rnvGq8h6pJ1dXUkJSWxceNG+vfvT0VFBT4+PnrHlEqlREVFsWPHDi0FysOl9UwxoJ/sZ9PDGB2k+Rmjl6VSY8KMcuCFqzQE9VS9Y69dnmfUOB/dOoBxkZqd+Dr3ZjJGBbMrA171Txav7jJcO7qrfbUpy+gRHvD+tEg8PDzUD9GDZ2q462vTZs7XJgxmcmx/5HI5GRkZJCQkIBKJSM8v595VhiV42FjCPw9GYG9vT0FBAcOHD2fz5s2MHDkSQRA4ffo0I0aM0Kt42gV9N1RkiLrkypUrmTx5stqh0ZXhAlhbW5OcnMw///yjdSzM1zBlxwdWHubu5Zqz7Y5nk3l9ovEtNyKXbKe1vZ2qqioNobrMp+O5ZZjhzUAe/vW0ljGocpFVnRENJXRRKq/oCHMBTLs2xqixVPpbuiaK3JRk5iXoT0LRRXYNXPXhAY3VT0R/D7KN8LZ35sX1BUxfno6lpSXu7u5UV1cDED9ETG5KMjkvXEViP+14fWfETrYUFRXh4uJCW1sbhYWFeHh4UFlZqQ4ptbW19apEziVvvIaoSx47doza2lrGjRtHTEwMK1as6HbcWbNm6e1ptGCstraxLvaWyJS1pZ0S0CeO6GdST9uYZbv4clsu/fv31xCqW3ZnLPufNS5e2p34+/361V80WKsjzNWZ3JRkAo2oegt7NY0P/tH2Jzx6fShZRt4jKK9ty5FzCRw21tZKIfnBxpfiZZcopX8d3D0pKdFMCrG2suKLe+PITUnmgylhWsUrNhbw5T0xSKVS2tvb8fX15cSJE4SGhl5QgbpL3ngNUZeUyWRkZWXx559/smnTJl599VWOHdNOHOhMbGwsBQUFOlUoHxyrnZvcFSMWb2PDwXMieqpuBLcP9zRqnK8PyUn83z6t+7O1VcZL7Y1oKBy6KJUV25VLPgcHB7XoGsATs3qvn/Cfzyaz8wnDY8If7arUGk8QBFqbm/nhdh8Mq5s6x/wf8rSquv5vRhK/P2ya2uT4Tw+RU1itU7wB4OoQH3JeGkf2c6N4OtqaT6ZF8tmNrnjYKNSxXpFIhCAI2NnZXVCBukveeA1RlwwICGD8+PE4Ojri5eXFmDFjyMnpev9jYWHBpEmTtJbgKu4YYZzCxlO/HGfsUs0f0auTItn4SKyed+gndFEqhTWNWq9nLUzm8VGGPxBeTy1XG4rqyd+Z3tKOdnd2Mqn88VRJCfn5+ezZs4eysjL8/PzY+tI4Vt0T0f0AnehAW5VksJeLSULyAEuz4GCe/tiYSCTCysqK0UN9CHFTMGhAf4qKitQrnKKiItzc3GhoaCAgIICioiJsbW2xtrZWN+7uDS55442LiyM/P5+TJ0/S0dHB6tWrmThxosY5kyZNYvv27chkMlpaWkhPTyckpPtOAjNmzGDt2rU6v8xFE42XVqmUKn9ETW1Kr5VcLsdW1sLKie5GC33f8EEm936epvX6A9dEGr0sD12UioOzs16JnH/mGiYm35V2NCjDXyOM+EVNWJ7HE78XER8fT2hoKO7u7ohEIob39zLJ8CIWb+O3fYUar+WmJPPMOOOE5AFm/lzcpZGp6nuLi4vx9vamqqoKS0tLRCIRTk5O9OvXTyOFUiqVEhAQQGlpaa91E7zkjdfKyooPPviAG264gZCQEO68807CwsL45JNP1AqTISEhjB8/nsjISEaOHMncuXMJD+++P25QUBDt7e1avXpV5CzUL03TFSPf2MXzK1LJyMhQFznsT0nm7jjjis3TS5WF5ufnxqqW5ZESwx8Jo/+3jyc3tWgJEoBSoN7YWVjXtsTa2ppnb/Fg12OG61ydaEJv6WVuSjIPxBu39Xjuj+NaD5h7xgRzwATRgK66KIpEIhwdHZFKpUilUjw9PamoqEAqlaobuNfW1qJQKNQzsre3N4GBgb1mvJd8qOhC88EHH9DY2MiCBQt0HpcrBIa/lmZ0JpWK842ivb2dqGXGZ/R8MiWIMSHa7UhLquu49kPd+c6GXlNnYhelYow+5PljqdI+hwwZYnQ99DezwogL1I4UtEilxJqgT5X17GjsbTWzmlJ+zebHA4ZnZrnZwq5ndX9fgiBQVFSEIAh4eXlx8OBBPDw8qK6uJiEhgYKCApydnfHw8GD//v2mZln13TjvhaayspIJEyaQlpam94moUCi49X9bOa69DTWIPx6MZpBYMwQ1cmkqTUb2j+1NXeWvpw1l5BD9WjimaEcLgkB9fT3Z2dnY2tri7u7O6doWnvjb8GJ/f0f4+6neiwl/MjWEMUM1VzxnahoZ/4HhYvnTo8W8eHOo1usqgbr9+/cTHx/Pli1bSExMpLS0FEdHR1xdXdUCddnZ2URERHTZ8F0PZuPtiptvvplXXnlFI/HjfI4ePYpga8/t32r3PjKEiSGuvD5FsyHYzoIK7v/eeOXEzGdH4WCrHexf+HMWvxwyPHkCei+bKgZ45Hp7HB0dkcvl+Pn5acTbjTU8XeoaAP/9OZM1h4x7ir40fgDTRg7Set2Ya3p3cjDXh2vvnRUKBYcPH8bR0ZG6ujq8vLwQi8UcPnyYuLg4tUCdo6Njr+c2X/J73ovBjBkzumzGrVAosLOzo6roDF9cZ0wR3zl+O1Kv9Ih22r+OGuzDoReN31fHLtvJp2namT+LJ8eQ/dwoo8bSlzwBSsN+/RrD9tVZKFt0RkREMHjwYC0/Qm5KMnHuhl9X+Gtb+WqH9r765cmxHDTyO3t142mOl9VpvZ6bkszDowxLEHn85zwOFNdovS4SiZBIJJw5c4bQ0FBKS0vV9eOdBeq6aklrKuaZF+WXHB8fz65du9SSLqolYFlZGbW1tXh4eFBVVcXIkSOxtrbm07QC/rfNuDI+FZ3bdqiY8vEuDlcan1upb+YctSyVWiOGuy/clqcm60+UMHYZnZ6eTnR0tNZs0y6VEmXk/lXfPcYsSsXwEg5YeG0/ZiRpZ8A1tLaTYGBlUfozo3C2O7fq6ejo4OTJk5SVlTF27FgOHjxIYGAgzc3NaoG60tJSDWUWIzEvm7vj7rvvZsaMGQwbNoz6+noqKytxcnLC19cXd3d3LCws1JUivr7Kgvza5lZGvb3HpM8b5m7Fz49qymkcLalh8nLj83O3LIhD4qqdvvfd9mMsSe1ZCxZBEKirq6O0tJQ7fjROu2fTvUEa2XGdMXYZveOJODycNe/xzJkznKxqZN4fhgvmfz51GKN06GwDTHwntVu/RrjYnlVz49Qi9B0dHfj4+CCVSnFxccHa2pqKigqGDRvGnj17GDVqFDY2NhdEBsdsvCidVq+88gobN24kKSmJV155BR8fHy1htZaWFo4ePUp0tObeNWZxKq0mZr1lP38VNtbnllStra3EvGn8AyHYBX55vHccPdnPjUIqlVJaWkplZSUuLi7qh1hjYyMJ7xpeDdTVnvr/ft3N8gOG61w9EOvO4xNGqP/dWe3TmHv8buYwogfpNuANOWd4ap3+Ek8XK3gn2Q4vLy+1CL0gCLS0tJCXl0dsbCw7d+4kMTGRI0eOEBAQgFhsXO72eZiNVx8nTpxg1qxZ3HHHHXz++efs2LEDOzv9+9qMjAyChoXjYG+LnfU54163v4jnf9cvr9MVDyeKmX/dOWdZZmYma46L+POYcc4n0N+CxdiKp3BreO/eELy9vbUeYmDcA+GnKZIuk2aMGSvQCf588twDIScnh0GDBuHs7MzXO/J4Y4thYnVLJgzg1lhtJxZAh1TKCD1L+yhfO76bm6D1HSsUCvbv3094eDhFRUU4Ozvj5eWFvb29zu/PCMzGawhPPvkkcXFxTJo0SefxNqmct9Yf4MecOjoAJytYc38sgd7Kyp/2jg6iXjdNMNwV2H12liouLqajowN3H1+S3jZOFB3g25lhxAzSjpfWNrUy6v+Mm9W7mjmvfTWVEiN+AV2NlbQolToTrquiooLa2lqCg5VZYk2trYw0cOXy/h3DuCZUewbu6OigrKyMZ9YV0Dkk7GwF6+YnIHHRFuUSBEHdy8nPz4/c3FySkpJ6o0+v2XgNYf/+/SxatIjvv/9e5/G1WUW88qf27BrvZ81Xc8+Vo419ayuVLaaldex4KglnGwuysrKIj49Xpgu+morUyG96iCOs66V46YYHwhkg0a2wYex4XRnw5px8Hl1nuKZUrh45IYDoRakYsiDPXjgGG0tL5HI5lZWVlJaWIpPJEIvFSCQSWuUCGw5WUnL6OHMmJODmqHtVJgiCukAmKSmJEydOEBISckF7FZmNtxOCIBAbG8tvv/2Gu7t2XGPa8j3klOj3b+5/4Spsz4YE/sgp4pl1pi2jF4z2JdGrg4EDB+LiovRKr9lziv/+ZYSQ1FnSn47H+WxiQGNjI6WlpVRXV/Pp/jYytLt4dklvxYTHe8D/ze+ldi5PjaSkqAgPDw8tCZ/Ve46z6K9CPe9U8v6tgUgsm9Wtb3x9fTVat6rIz8/H2dkZiUR/iqtCoeDo0aP4+fkhkUh6Y9YFs/EazuLFi3Fzc+Pee+/VOnbjBzs5XdN1a8vPp4cx6qzChqktU1T8cpefejkIyodLV/m2+kj2g9kRDtjZ2eHr64uXlxcWFhYUlZZy/ee918nwg+9S+ciIHJauxhqzKBVDny0vJXsQ6gLDh+suJunqYbB4nBvjIgZqtb45n+bmZvLy8rSclZ1ReebLysqIiTGtJFEH5iQNQ+mq0mhqdPet9e5feZgb/0/5Y7GwsCA3JZmxgwxT5zif29aUsPvYuTI+VUGCr5NxS7HUErh3UwvDhw/X6J3jL5HwxXV2RonJhy5K5Zu0bJ3H5s/svTrhbSnJZBpYoP9qag3vb6+ho0P3gzU3JZlBOoQw+rvaMD4xUl3N1BWqzDFVTbQuqqqqWL16da+qZXTFZWG8GzduJDg4mKCgIF5//XW95+3duxdLS0t+/PFHvecMGKDUaC4q0t57zUrob9D1nG5S/jBb2pXJyx/PjGbHU8YrRQDct/oIt7yj+SPf/OQ4Vt8Tqecd+gldlEpueZ3636rqly0PhvL2rYYLECzbVtutwoYx13Tz67rHcjgrQmAIu2pgxVbNrDO5XE5paSn79u3j5aucWHlnPyL9HPBztuLaoZ6svC8We2vDPcGqgvrOtLa28vPPPzNlyhSmTJmCQqFQ/4YuNH1+2WyIQJ3qvOuuuw47OzvmzJnDHXfcoXfMTz75hMrKSp566imtY/O/y2LLCcNDOLOiPHn+lnOGFrEoFVOFUPY/PxrbThlLpi6jbw6y4Y3pyjTKpqYm8vPziYpSlvEZ68w6XzyvM83NzcS9nWHwWF0Z6udbDvDOjuou328JZD4/hqYGZVJJU1OTeh9rQkGAFjKZjIceeoiPP/6YPXv2sHr1avbu3cuNN97IrFmziIiI6LVyv05cvnve3bt388orr7Bp0yYAli5dCsDzzz+vcd67776LtbU1e/fu5eabb+7SeGtqarjuuuvYtm1brzTjAs1E+293n2Lp38Y7nwDemxxC0iAXysrKqKiowNHRkS3FAp9mdP3D1oXKWDIyMjSUDSd9sJX8GsO95aEO8OPTF1472pCx3hxlxSA/ZQKFvmZqpiAIAnl5ecybN4/W1lZGjBjBrFmzuPrqq3sax+2Oy3fPa4hAXXFxMb/88gsPPfSQQWN6eHjQv39/Dh48qPO4KbUh4a9tpbhOWSk7K3Eg+58zTenwsZ+PMO6tdOzs7IiNjSUiIoIF4yNNKjYPXZTK9sNnkEgkGhI56+aP5ff7DV+W57Z0L1RnqAbm2I8P8Xmafidabkoy+uobxA4g8fYkJCSkWweUoZSXl/Phhx9yzTXXsHDhQq699lq1yOF11113oQ23S/q88RoiUPf444+zbNkyo77oGTNm6FWXTH3KuModFde9l84Hm5XhI1sbpdKhi9EyvtCEsjtfu+zc7GhlaUluSjK2Rv5FH/ypgNlfF2jpWw329TRJlyr9uO7QzK6FY/jyesM6jr2zrZQl64/oPCaTyfhhbjBLzmt24W4nYuGEUDqa63vcla+lpYW1a9cyefJkpk2bhqWlJb/++isbNmzg1VdfJTs7m+Zm01uF9hZXxLI5MDBQbeRVVVU4ODjw2Wefceutt+odt62tjdjYWHbv3q3T6E1JnFDhbw9//+ecYWw9Wsq8tcaFbFQsvG4AMxI10/y+33OSxX+dMnqs7Y+OwFNHfHv6Z7vJLjM8B9kZSNdh+KoO8oZ2cNj25Ci8nGxQKBRUV1dTWlpKS0sLPj4++Pr6Ym9vz6mqJsoapQz0tEfiYkdeXp7OmG93yOVyduzYwapVq9i/fz8TJkxg9uzZhIaGak0Gp06dol+/fhdr1r1897wymYyhQ4eyefNm/P39iYuLY+XKlYSFhek8/5577ul2z6ti7ty5TJ48mXHjxuk8bmq3OhUr7xnOiP5KlcqmpiZG/p9xLU46c/4s2SGTMcKEFixv3TKQCVHaws4tbW3EvmFcqub511RTU0NZWRmhoaGcPn2aG7/pOij8zqRBDLBtpba2Fnd3d3x9fXFxcelyOdzY2MiJEyf0xnw7IwgCR44cYfXq1WzatInY2Fhmz57NmDFj/tXl8HlcvnteQwTqTGXmzJl6l86g/HEO8zRJHQGA6V/nMO19pVBdfn4+qQ+G0t/NhHU0Z0XhOoWBbKysyE1Jpp+Lcdf39O+niNfxUHKwsyM3JRkbI34xoYtSefu3c3Iz7u7u1NfXI5fLGTBgQLfL8vqaSjw9PYmPj2fYsGEGOaCcnZ1pa2vTG/MFZU+i9957j+TkZF5++WViYmLIyMjgq6++Ijk5+VIy3C7p8zPvhUQulzNixAjS0tLU6gi6+PCXrXx4sGf7rM4hl8MldUxZbpyonIrxwS78312a2T2ZJ6uY/a1u51tX7H1mFI522g+TptY2Rr5p+Cwc6AJ/ni1XPH78OI6OjuqaaICdeaXcv0Zz2zB2sBtv3B6Bs53xChTn9/MFZdjq999/Z82aNTQ2NnLXXXcxbdq0blvjXAJcvsvmC80zzzxDZGQkkydP1ntOSUkJbW1t3Pzd6R591nd3RxA9QNkzwNQYrgpds5opy/yJwba8fpfuBJOnV2ex3sCyxbtjJTw7IUSdZjhixAh1QXtraytisZi9lSKOVLQyop8LycE+uNiZtqppb29n/vz5fPzxx2zbto3Vq1dz4MABbr75ZmbNmsWwYcMuRDz2QmE2XlM5ePAgzz//PGvXrtV7jkwmIzMzk/j4eL7dc5rXTYzhAgx1h18fPWd4z/1yiN8OVpo01qdTQ7jqPOXEmz7YzskamdFj9bQowVoE2S+Oo76+npycHKysrPDy8sLPzw8nJ6dejccePnyYBx54AJlMRkJCAjNnzmTMmDG9VShwsTEbr6kIgkB8fDw//vgjXl76O+kcOHCAgQMHYmNjw5GC08z4ybCicH38Pi+WwWfrhGua2xhtQl0vgK8NbH5O0/ByzlQz7esDRo/13cxQogfpVoU4UVbJzZ/pbkoOSufKp9fa4urigpWVFba2tgQGGtjxrBsEQaCsrIy1a9fy888/4+vrS3BwMFKpVGcP5j6G2Xh7whtvvIGNjQ0PPPCAzuNSqZSCggJ1xpNEIsHHx4fEN3fSIjP9K3p6bABzOjU9S3o9lbqui5r0krNwLNaWmjPP2Le2UdliXLJmV9rRALGLU9E1ZKiXNasfSMDKygqpVKpRr2wqTU1NrFu3jrVr19LS0sLUqVOZOnUq3t7edHR0EBsby759+y6IcuNFxGy8PaGoqIipU6eyadMm9Y/t/Nijp5c3/5d2hiMNVjhYWbD01jCG93fjSHEtt3+RbfJn2wH7Oomav/p9GqtNk47mxev6MT1RU8VQEARilqTRZoQNO1hA5ovaBqxQKKisrKS4uJj7/6yn9eyvI8rfiffvisTD6ZyM7IEDBxgwYACursZVXMlkMtLS0li1ahW5ublMnDiRWbNmMWTIEK0HQUFBAYGBgX11uazCbLw9QRAErr32Wt555x3s7e1pbGxUy8H6+vri7OzMXZ+lc6hcs1A/yteW7+9PorG1nXgDpUX18eN9kYT6e3LixAnkWDLpO9Ms2FUEu1/SNryWDimxrxsuyXpfQgBPXT9EQ12yvr4eLy8vJBIJzs5dJ0RWVVVRVVXFsGHDuv0sQRA4cOAAq1evZvPmzYwaNYqZM2cyatSovm6YhmA23p5w/PhxnnrqKbKzs5k2bRqPPPKIWg4WoLKxjbHv6N+TqqRWbnx/B6drjexx0on/JPdjaqw/ubm5xMTEcO272ylpMN75BLD9qSQ8HbUF1cvrW0j+X3q37w90t+F/N4rV6pISiQQPDw+Dl8H65GtUCIJASUkJa9as4ZdffqFfv37MnDmTW265BVtb45tn92HMxmsq27ZtY+nSpdx22218+OGH7Ny5U+vHll/eyKRPu+59887tw7ghzJfPth3n3bSupVm6Y8UEZ8LDw7Gzs+PbXSdY+o9pIarnrxnArFG6FRR/2ryDl3bqf9CEeYp4b3KwTolcQzl27Bhubm4asdbGxkZ+/fVX1qxZg1QqZdq0adx11114ehrXLfAywmy8vcGUKVNYsGCBlhSKQhAIfzWt2/eHeon48eFxyOQKIhebLo8D8P4Eb66JPdfGdM43mew5bXwntNEDHfls9kit11XL2sIOBxb8qqljbGMBK+6NJtLfNIUQFRUVFbz++ussW7aMLVu2sHr1avLy8pg0aRKzZs1i8ODBl1Q8duPGjSxYsAC5XM7cuXN57rnnNI4LgsCCBQtYv349Dg4OfP31113K5hiI2Xh7g19//ZV//vmHN954Q+vYhPe2c6rOsCVs+jOjcbaz5u6vsthbaLw2s4oAZ/jrCc39a2puGY/8qLsiRx+jB7vz8bRILM+uKARBoKamhgMHDmBnZ4e3tze5tZb8fqQWf3c7ZsT3Z4iPDl0ZI1AoFOTk5DBr1izs7e0ZM2YMs2bNIiEh4ZLcxxoi+rB+/Xref/991q9fT3p6OgsWLCA9vfstSDeYjbc36OjoIDo6mp07d+rowSMjaqnhhQDzr/Lj4eRg8ssbmPSp4R0IdJHx7Ciczusa+PfhUhb8ZFilkr0VvDoxlDEDHSktLaWqqgoXFxfkcjlisbiniv9qBEGguLiY1atXs27dOgYOHIi7uzthYWE88cQTvfIZFwpDqtcefPBBxo0bx7Rp0wAIDg4mLS1NIxXUBC7fwoSLiY2NDaNHjyYtLU3rmK21lVFJ+x9sL2H0klSCfJzJTUnu0R9i5LKdrNuvue+9LsyX3JRkfprbfXUNAuw+eIz8/HxcXFwYOXIkYWFhBAYGatX5mkJ9fT3ffPMNN998M/fffz8+Pj78/fff/PTTTyxdupQffvihx59xoTFU9KG7c3oTs/EaSVftQP96VHvv2BU1Mgh7NY3jlY0cSknm3vju1Sn18fzvJ3j2h2yt10P8PMhNSWbjvCi977W1tGDc8CFERUUhFovVDihVhY5UaryHXCqVsmHDBmbPns2ECROorKzkiy++IC0tjYceeggPD2UppLe3N8uXLzd6/IuNIaIPhpzTm1z2xtudsuT3339PZGQkkZGRJCUlkZPTdZe+UaNGcejQIZqamrSO+bhqi3UbwsSPM1m6Ppf/3DCU3JRkJkeYVuny+5Fabnx3O1K5dsZFPy9Xtj8Sycc3uDLASfmHtwR8nKy5a2QAySG6l8bnS+R0hUKhICsri6effppRo0axZcsW/vOf/5CVlUVKSgqDBg3SrQnWRVPzS4WAgAAKC89FCYqKivDz8zP6nN7kst7zGuJk2LVrFyEhIbi7u7NhwwZeeeWVbp0MCxcuZPDgwUydOlXr2Naj5cxbq9342hBsgKyF5wrBa5vbuOrt3RhbbDh5uJjXJinvsXOXBFdXV/z8/HB1dUUhQFlDG/bWFnjoiPeqaG9v58CBA8TFxek8LggCZ86cYfXq1fz2228EBQUxa9Ysxo8frxa0uxwwRPThzz//5IMPPlA7rB577DEyMgxXz9SD3qm7Tyd9dkdGRgZBQUEMGqSMZU6dOpV169ZpGG9S0rlyt4SEBJ16zeczc+ZMHn/8cZ3GO3aYGDDNeDuAiMXbWHlPJCP6e+LuaMehlGSkMhl3fZHJ0XLDWknvOVnLqVOnKC8vV3dJCAoK0vDiWorA3617TSlbW1ssLS1pbm7WaANSV1fHL7/8wtq1a7GwsGD69Omkpqbi5uZm7G33CTqLPsjlcubMmaMWfQB46KGHmDBhAuvXrycoKAgHBwe++uqrC3tNF3T0fxldDoSuZtUvvviCG2+8sdtxhw0bRkNDAxUVFTqLuV+5aSiv/HnMtIsGpn99gKlRPqTconyqW1tZ8fODSsW1rUdKmPdDXpfvtxE6sLKy0tmZ3hQEQeCzzz7jkUce4a+//mL16tWcOnWKyZMn88033zBgwIA+FY/9/vvvWbZsGQBOTk58/PHHBsnmTJgwgQkTJmi81lmRVCQS8eGHH/bCHRjGZW28xjgQUlNT+eKLL9ixo/v8XpFIxJQpU/jxxx95+OGHtY7fGeNPVVM7H2w1vTh/9f4KUvOq2PT4KGw6VcWMDfEjN8UPmVzO42sPsiW/VuN9dpaQMjGSgIDeyUhSKBQ0NzfzySefsHLlSq655hpeeOEFoqOjL9l47COPPKKxVZo4caLGaiswMJCtW7eqt0oPPPBAb8RjLzqXtfEa6kA4cOAAc+fOZcOGDQan4U2bNo3Jkyczb948nQ+Eh8cOYt6YQFauT2OxiWHc8hYF497ZydYnr1KX8wmCoN7HzhjYyvwoX3YdryC3xRVHexumxvkTInEx7QPPIggCp06dYs2aNfz2228EBwcTEhLCs88+yzXXXNOjsS80F2qrdClyWRtvXFwc+fn5nDx5En9/f1avXs3KlSs1zjlz5gyTJ0/m22+/ZejQoQaPrVIyzM/P1/s+kUhE3BB/0uLc8PLy4tZP0jleZbiEKkBdq4L9hfVESuwpLS2lvLwcBwcHfH19GTJkCBYWFjg5OXG1TMbAgQONGvt8ampq+Pnnn/nhhx+wsbFhxowZbN26FVdXV7Zt28Y333xzyRvvhdoqXYpc1sZriJNh0aJFVFdXq5e/VlZWZGZ2XWSgYvr06axdu5YXX3xR7zm+vr6cOHECHx8ffns4EYA2qZzFG/L4KbvcoM/ZdzAXiypbfH19iYmJ0drHisVi9u3bZ9Les729nU2bNrF69WqKioq4/fbbWblyJQEBARpjjR49GiennqVEXgwu1FbpUuSyDhVdaBobGxk9erTOSqPOpKen63UeldU1c9snGdTrUciwtYA1c4Yz1M+jy2s5cOAAgYGB3dbRwrlyvNWrV7N7926uvfZaZs+eTVRU1CXleDIFQ3tXHThwgNtuu40NGzYYteL6FzDnNl8opk+fzv333098fLzec06dOoW1tTX+/rozqARBoKGhgaxjp/hodw25dcrX3ewteGtyKEmDu1f/r6yspKamRqMZ9/mfUVBQwOrVq/nzzz8JDw9n5syZXHfddX1dJkYDQ+KxZ86c4eqrr2bFihUa+99LlCszznsxUPU06sp4JRIJhw8f1jLe1tZWSktL1dpXkYH+rI2JwMLCAkEQjJoFPT09OX78OAqFQmMVUF1dzU8//cSPP/6Ivb09M2bM4Pnnnzdohu6LXOit0qWEeebtIVKplKioKHbs2NFlRlFWVhahoaFYWVlRXl5OWVkZIpEIX19ffHx8emX2++GHH3B3d2f06NFs2LCB1atXU15ezh133MGMGTPw8/O7pJbF3cVjVezdu5eEhATWrFljUJuaywzzzHuhsLa2Zty4cWzevFmv11KhUODo6Mj+/fuxtLRELBarlTB6C4VCQV1dHe+99x7t7e3ccMMNvPbaa0RGRl5SBqvCkHis6rxnn32WG2644V+60ksXs/H2ArNmzeKtt97SMF5BEKivr6e0tJS6ujrc3NwQBIG4uLheS24QBIH8/HxWr17Nhg0biIiIoKGhgb17917yaYqGxGMB3n//fW6//Xb27jW9CdvlyqWXItMHUcWTGxoaaGpqoqCggD179lBUVISPjw8JCQmEhITg7OzcK31dKysr+fjjj7n22mt57rnnCAsLY+fOnXz33Xfcfffd/Pbbb71wVxeWC9EU/UrDPPP2ArW1tfj7+3P11VczadIk5s6dy8iRI7WE2Xx9fSktLTXJWdTa2sqGDRtYtWoVVVVVTJkyRd0doDNz586ltrZWzyiXDheqKfqVhNl4e8hvv/3GkiVLSE5OpqWlhRdffFHvHlPlETbUkyyXy9m1axerV68mMzOTG2+8kddff53w8HC975dIJEgkEp3HLiUMSV3NzMxUV25VVVWxfv16rKysumyKfiVh9jb3EFVoRhAERo8ezYoVK7rULDpy5Ag+Pj56c6gFQSAvL49Vq1axadMmoqKimDVrVp/qG2sIF7Ip+mWG2dt8oVA5n0QiEVOnTmXt2rUsWLBA7/m+vr4UFRVpGW95eTk//vgjP/30Ex4eHsyYMYP//ve/ODg4XNDr/7cwJB5rpmvMM28vUlFRwc0330xqaqreZa0gCLzxxhvMmzcPKysr/vjjD1avXk19fT133nkn06dP7zW1xt7CkHhsWloajz/+OFKpFC8vL7Zu7ZkutRk15vTIi8VNN93Ef//7X726THK5nAceeICmpiZOnz7NhAkTmDVrFqGhoZdsPLY7KaG6ujqSkpLYuHEj/fv31ytSYMYkzNKvFwtd6pKCIJCbm0tKSgpJSUlIpVIqKyvZv38/r7/+OmFhYZek4YJmPNbGxkYdj+3MypUrmTx5Mv379wcwG+5F4oo33u7UJQVB4LHHHiMoKIjIyEj27dvX5Xi33normzZtQqFQUFZWxnvvvUdycjKvvPIKsbGxZGRk8OOPP2JpaUl1dfWFuq1ew5B47LFjx6itrWXcuHHExMSwYsWKi32ZVyRXtMPKkBS9DRs2kJ+fT35+Punp6cybN6/L4m4HBwd1zNfKyoq77rqLP/74Q2s2WrZsGQqFsbqQFx9D4rEymYysrCw2b95Ma2sriYmJJCQkXOqldn2eK9p4DUnRW7duHbNnz0YkEpGQkKDuRdtVOOjJJ5/E2tqasWPH6l0Ojx49undv5gJhqF6xl5cXjo6OODo6MmbMGHJycszGe4G5opfNF6qFxXXXXce4ceMu2X2sMXSWEuro6GD16tVMnDhR45xJkyaxfft2ZDIZLS0tpKenExIS8i9d8ZXDFT3zXootLC41DInHhoSEMH78eCIjI7GwsGDu3LmEh4d3M7KZHiMIQk/+69Ps2rVLuP7669X/XrJkibBkyRKNcx544AFh5cqV6n8PHTpUKCkpuWjXaAwbNmwQhg4dKgwePFhYunSp1vG6ujrh5ptvFiIjI4XQ0FDhyy+//Beu0oyR6LW/K9p4pVKpEBgYKJw4cUJob28XIiMjhUOHDmmc88cffwjjx48XFAqFsHv3biEuLu5futqukclkwqBBg4SCggL1vRw+fFjjnMWLFwvPPPOMIAiCUFFRIbi7uwvt7e3/xuWaMRy99ndFL5svxRYWpmKI800kEtHY2IggCDQ1NeHh4XFZ6VddaVzxf7lLrYWFqRiiVzx//nwmTpyIn58fjY2NrFmz5pLsemDGMMx/ucsEwQDH2qZNmxgxYgQlJSVkZ2czf/58GhoaLtYlmullzMZ7mWBIPParr75i8uTJiEQigoKCCAwM5OjRoxf7Us30EmbjvUwwJB7bv39/Nm/eDChLEPPy8tR7ZDN9jyt+z3u5YIjz7aWXXuKee+4hIiICQRBYtmwZXl5e//KVmzEVc0ngRWbOnDnqXOdDhw5pHRcEgQULFrB+/XocHBz4+uuviY6O/heu1Mwlgrkk8FLhnnvuYePGjXqPdy6E+Oyzz5g3b95FvDozfQmz8V5kxowZg4eH/qZh+gohzJg5H7PxXmKYUghh5srEbLyXGIbEa82YAbPxXnIYEq81YwZ67m02YwIikWgg8IcgCFp1cyKR6CZgPjABiAfeEwRh5MW9QjN9AXOc9yIjEolWAeMAL5FIVAS8DFgDCILwCbAepeEeB1qAe/+dKzVzqWOeec2Y6aOY97xmzPRRzMZrxkwfxWy8Zsz0UczGa8ZMH8VsvGbM9FHMxmvGTB/FbLxmzPRRzMZrxkwf5f8BEsyv4h55NxMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the triples as points in R^3\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.scatter(xs[0],xs[1], xs[2])\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (18.0, 18.0) \n",
    "\n",
    "ax.view_init(-100,110)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wichmann-Hill (1982)\n",
    "\n",
    "Sum of 3 LCGs. Period is 6,953,607,871,644. \n",
    "\n",
    "    def WH(s1, s2, s3):\n",
    "        s1 = (171 * s1) % 30269\n",
    "        s2 = (172 * s2) % 30307\n",
    "        s3 = (170 * s3) % 30323\n",
    "        r = (s1/30269 + s2/30307 + s3/30323) %  1\n",
    "        return [r, s1, s2, s3]\n",
    "\n",
    "#### The right way, the wrong way, and the Microsoft way.\n",
    "WH generally not considered adequate for statistics, but was (nominally) the PRNG in Excel for several\n",
    "generations. Excel did not allow the seed to be set, so analyses were not reproducible. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![mcCullough](./RandPics/notWichmannHill08.png)\n",
    "McCullough, B.D., 2008. Microsoft Excel's 'Not The Wichmann–Hill' random number generators\n",
    "_Computational Statistics & Data Analysis_, _52_, 4587&ndash;4593\n",
    "doi:10.1016/j.csda.2008.03.006"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mersenne Twister (MT)  Matsumoto & Nishimura (1997)\n",
    "\n",
    "+ example of \"twisted generalized feedback shift register\"\n",
    "+ period $2^{19937}-1$, a Mersenne Prime\n",
    "+ $k$-distributed to 32-bit accuracy for all $k \\in \\{1, \\ldots, 623\\}$. \n",
    "+ passes DIEHARD and most of TestU01 (see below)\n",
    "+ standard in many packages:\n",
    "    - GNU Octave, Maple, MATLAB, Mathematica, Python, R, Stata\n",
    "    - Apache, CMU Common Lisp, Embeddable Common Lisp, Free Pascal, GLib, PHP, GAUSS, IDL, Julia, Ruby, SageMath, Steel Bank Common Lisp, Scilab, Stata, GNU Scientific Library, GNU Multiple Precision Arithmetic Library, Microsoft Visual C++.\n",
    "    - SPSS and SAS offer MT, as does C++ (v11 and up)\n",
    "+ generally considered adequate for statistics (but not for cryptography); however, will trouble that in this work, esp. for \"big data\"\n",
    "+ usual implementation has 624-dimensional state space, but TinyMT uses only 127 bits\n",
    "+ seeding complicated, since state is an array\n",
    "+ can take a while to \"burn in,\" especially for seeds with many zeros\n",
    "+ output for close seed states can be close\n",
    "+ 2002 update improves seeding\n",
    "+ completely predictable from 624 successive outputs\n",
    "+ problems discovered in 2007 (see TestU01, below)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Python implementation of MT19937 from Wikipedia \n",
    "# https://en.wikipedia.org/wiki/Mersenne_Twister#Python_implementation\n",
    "\n",
    "def _int32(x):\n",
    "    # Get the 32 least significant bits.\n",
    "    return int(0xFFFFFFFF & x)\n",
    "\n",
    "class MT19937:\n",
    "\n",
    "    def __init__(self, seed):\n",
    "        # Initialize the index to 0\n",
    "        self.index = 624\n",
    "        self.mt = [0] * 624\n",
    "        self.mt[0] = seed  # Initialize the initial state to the seed\n",
    "        for i in range(1, 624):\n",
    "            self.mt[i] = _int32(\n",
    "                1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)\n",
    "\n",
    "    def extract_number(self):\n",
    "        if self.index >= 624:\n",
    "            self.twist()\n",
    "\n",
    "        y = self.mt[self.index]\n",
    "\n",
    "        # Right shift by 11 bits\n",
    "        y = y ^ y >> 11\n",
    "        # Shift y left by 7 and take the bitwise and of 2636928640\n",
    "        y = y ^ y << 7 & 2636928640\n",
    "        # Shift y left by 15 and take the bitwise and of y and 4022730752\n",
    "        y = y ^ y << 15 & 4022730752\n",
    "        # Right shift by 18 bits\n",
    "        y = y ^ y >> 18\n",
    "\n",
    "        self.index = self.index + 1\n",
    "\n",
    "        return _int32(y)\n",
    "\n",
    "    def twist(self):\n",
    "        for i in range(624):\n",
    "            # Get the most significant bit and add it to the less significant\n",
    "            # bits of the next number\n",
    "            y = _int32((self.mt[i] & 0x80000000) +\n",
    "                       (self.mt[(i + 1) % 624] & 0x7fffffff))\n",
    "            self.mt[i] = self.mt[(i + 397) % 624] ^ y >> 1\n",
    "\n",
    "            if y % 2 != 0:\n",
    "                self.mt[i] = self.mt[i] ^ 0x9908b0df\n",
    "        self.index = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### xorshift family\n",
    "\n",
    "Originated by Marsaglia, 2003.\n",
    "\n",
    "Vigna, S., 2014. Further scramblings of Marsaglia's xorshift generators. https://arxiv.org/abs/1404.0390\n",
    "\n",
    "128-bit xorshift+  Implemented in Python package randomstate https://pypi.python.org/pypi/randomstate/1.10.1\n",
    "\n",
    "    uint64_t s[2];\n",
    "\n",
    "    uint64_t xorshift128plus(void) {\n",
    "\t    uint64_t x = s[0];\n",
    "\t    uint64_t const y = s[1];\n",
    "\t    s[0] = y;\n",
    "\t    x ^= x << 23; // a\n",
    "\t    s[1] = x ^ y ^ (x >> 17) ^ (y >> 26); // b, c\n",
    "\t    return s[1] + y;\n",
    "    }\n",
    "    \n",
    "    \n",
    "1024-bit xorshift+\n",
    "\n",
    "    uint64_t s[16];\n",
    "    int p;\n",
    "    \n",
    "    uint64_t next(void) {\n",
    "        const uint64_t s0 = s[p];\n",
    "        uint64_t s1 = s[p = (p + 1) & 15];\n",
    "        const uint64_t result = s0 + s1;\n",
    "        s1 ^= s1 << 31; // a\n",
    "        s[p] = s1 ^ s0 ^ (s1 >> 11) ^ (s0 >> 30); // b, c\n",
    "    return result;\n",
    "}\n",
    "\n",
    "xorshift+ passes all the tests in BigCrush, has 128-bit state space and period $2^{128}-1$, but is \n",
    "only $(k-1)$-dimensionally equidistributed, where $k$ is the dimension of the distribution of the xorshift\n",
    "generator from which it's derived. E.g., for the 128-bit version, xorshift+ is only 1-dimensionally equidistributed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other non-cryptographic PRNGs\n",
    "\n",
    "See http://www.pcg-random.org/ and the talk http://www.pcg-random.org/posts/stanford-colloquium-talk.html\n",
    "\n",
    "PCG family permutes the output of a LCG; good statistical properties and very fast and compact. Related to Rivest's RC5 cipher.\n",
    "\n",
    "Seems better than MT, xorshift+, et al. \n",
    "\n",
    "    // *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org\n",
    "    // Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)\n",
    "\n",
    "    typedef struct { uint64_t state;  uint64_t inc; } pcg32_random_t;\n",
    "\n",
    "    uint32_t pcg32_random_r(pcg32_random_t* rng)\n",
    "    {\n",
    "        uint64_t oldstate = rng->state;\n",
    "        // Advance internal state\n",
    "        rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);\n",
    "        // Calculate output function (XSH RR), uses old state for max ILP\n",
    "        uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;\n",
    "        uint32_t rot = oldstate >> 59u;\n",
    "        return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PRNGs based on cryptographic hash functions\n",
    "\n",
    "Cryptographic hash functions have several basic properties:\n",
    "\n",
    "1. produce fixed-length \"digest\" of an arbitrarily long \"message\": $H:\\{0, 1\\}^* \\rightarrow \\{0, 1\\}^L$.\n",
    "1. inexpensive to compute\n",
    "1. non-invertible (\"one-way,\" hard to find pre-image of any hash except by exhaustive enumeration)\n",
    "1. collision-resistant (hard to find $M_1 \\ne M_2$ such that $H(M_1) = H(M_2)$)\n",
    "1. small change to input produces big change to output (\"unpredictable,\" input and output effectively independent)\n",
    "1. equidistributed: bits of the hash are essentially random \n",
    "\n",
    "Summary: _as if_ $H(M)$ is random $L$-bit string is assigned to $M$ in a way that's essentially unique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 step of SHA-256\n",
    "By User:kockmeyer - Own work, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=1823488\n",
    "\n",
    "![SHA-2](./RandPics/SHA-2.svg)\n",
    "\n",
    "\\begin{equation*}  \\mbox{Ch} (E,F,G) \\equiv (E\\land F)\\oplus (\\neg E\\land G) \\end{equation*}\n",
    "\\begin{equation*}  \\mbox{Ma} (A,B,C) \\equiv (A\\land B)\\oplus (A\\land C)\\oplus (B\\land C) \\end{equation*}\n",
    "\\begin{equation*}  \\Sigma _0 (A) \\equiv (A\\!\\ggg \\!2)\\oplus (A\\!\\ggg \\!13)\\oplus (A\\!\\ggg \\!22) \\end{equation*}\n",
    "\\begin{equation*}  \\Sigma _1 (E) \\equiv (E\\!\\ggg \\!6)\\oplus (E\\!\\ggg \\!11)\\oplus (E\\!\\ggg \\!25) \\end{equation*}\n",
    "\\begin{equation*} \\boxplus  \\mbox{ is addition mod } 2^{32}\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple, hash-based PRNG\n",
    "\n",
    "Generate a random string $S$ of reasonable length, e.g., 20 digits.\n",
    "\n",
    "\\begin{equation*}  \n",
    "   X_i = {\\mbox{Hash}}(S+i),\n",
    "\\end{equation*}\n",
    "\n",
    "where $+$ denotes string concatenation, and the resulting string is interpreted as a (long) hexadecimal number.\n",
    "\n",
    "**\"Counter mode.\"  Hash-based generators of this type have unbounded state spaces.**\n",
    "\n",
    "Implementation in Python by Ron Rivest: http://people.csail.mit.edu/rivest/sampler.py\n",
    "\n",
    "Implementation in angular-js by Chris Jerdonek: https://github.com/cjerdonek/quick-sampler\n",
    "\n",
    "Implementation in JavaScript by Philip Stark: https://www.stat.berkeley.edu/~stark/Java/Html/sha256Rand.htm\n",
    "\n",
    "#### Variation to take advantage of incremental hashing\n",
    "\n",
    "\\begin{equation*}  X_i = {\\mbox{Hash}}(S+0*i),\n",
    "\\end{equation*}\n",
    "\n",
    "where $+$ denotes string concatenation and $0*i$ denotes a string of $i$ zeros.\n",
    "\n",
    "This is the method used in the [`cryptorandom`](https://github.com/statlab/cryptorandom) python package.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some trainwrecks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Microsoft Visual Studio system.random, 2011\n",
    "\n",
    "![Microsoft](./RandPics/visualStudiosysRand.png)\n",
    "\n",
    "https://connect.microsoft.com/VisualStudio/feedback/details/634761/system-random-serious-bug"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![XP2](./RandPics/xp2k.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aristocrat Leisure of Australia Slot Machines\n",
    "\n",
    "<img src=\"./RandPics/wiredRussianSlot.png\" alt=\"Wired Slots\" height=\"600px\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./RandPics/russiaLeisure.png\" alt=\"Russian hack\" height=\"600px\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dual Elliptic Curve\n",
    "\n",
    "![dual EC](./RandPics/dualEC-1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "![dual EC](./RandPics/dualEC-2.png)\n",
    "Bernstein, D.J., T. Lange, and R. Niederhagen, 2016. Dual EC: A Standardized Backdoor, in _The New Codebreakers, Essays Dedicated to David Kahn on the Occasion of his 85th Birthday_, Ryan, P.Y.A., D. Naccache, and J-J Quisquater, eds., Springer, Berlin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "## GnuPG RNG bug (18 years, 1998--2016) \n",
    "\n",
    "\"An attacker who obtains 4640 bits from the RNG can trivially predict the next 160 bits of output\" \n",
    "\n",
    "https://threatpost.com/gpg-patches-18-year-old-libgcrypt-rng-bug/119984/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## RSA\n",
    "\n",
    "https://www.schneier.com/blog/archives/2012/02/lousy_random_nu.html\n",
    "\n",
    "![Lenstra](./RandPics/lenstraEtal.png)\n",
    "\n",
    "\"An analysis comparing millions of RSA public keys gathered from the Internet was announced in 2012 by Lenstra, Hughes, Augier, Bos, Kleinjung, and Wachter. They were able to factor 0.2% of the keys using only Euclid's algorithm. They exploited a weakness unique to cryptosystems based on integer factorization. If n = pq is one public key and n′ = p′q′ is another, then if by chance p = p′, then a simple computation of gcd(n,n′) = p factors both n and n′, totally compromising both keys. Nadia Heninger, part of a group that did a similar experiment, said that the bad keys occurred almost entirely in embedded applications, and explains that the one-shared-prime problem uncovered by the two groups results from situations where the pseudorandom number generator is poorly seeded initially and then reseeded between the generation of the first and second primes.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## PHP\n",
    "\n",
    "![php](./RandPics/argyrosKiayias.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Bitcoin on Android: Java nonce collision \n",
    "\n",
    "\"In August 2013, it was revealed that bugs in the Java class SecureRandom could generate collisions in the k nonce values used for ECDSA in implementations of Bitcoin on Android. When this occurred the private key could be recovered, in turn allowing stealing Bitcoins from the containing wallet.\"\n",
    "\n",
    "https://en.wikipedia.org/wiki/Random_number_generator_attack#Java_nonce_collision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Debian OpenSSL\n",
    "\n",
    "![openSSL](./RandPics/openSSL.png)\n",
    "\n",
    "Valgrind and Purify warned about uninitialized data. \n",
    "Only remaining entropy in seed was the process ID. \n",
    "\n",
    "Default maximum process ID is 32,768 in Linux.\n",
    "\n",
    "Took 2 years (2006--2008) to notice the bug."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "XKCD \n",
    "![xkcd](https://imgs.xkcd.com/comics/security_holes.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Generating a random integer uniformly distributed on $\\{1, \\ldots, m\\}$\n",
    "\n",
    "\n",
    "#### Naive method\n",
    "A standard way to generate a random integer is to start with $X \\sim U[0,1)$ and define $Y \\equiv 1 + \\lfloor mX \\rfloor$. \n",
    "\n",
    "In theory, that's fine. But in practice, $X$ is not really $U[0,1)$ but instead is derived by normalizing a PRN\n",
    "that's uniform on $w$-bit integers. Then, unless $m$ is a power of 2, the distribution of $Y$ isn't uniform on $\\{1, \\ldots, m\\}$. For $m < 2^w$, the ratio of the largest to smallest selection probability \n",
    "is, to first order,  $1+ m 2^{-w}$. (See, e.g., Knuth v2 3.4.1.A.)\n",
    "\n",
    "For $m = 10^9$ and $w=32$, $1 + m 2^{-w} \\approx 1.233$. That could easily matter.\n",
    "\n",
    "For $m > 2^{w}$, at least $m-2^w$ values will have probability 0 instead of probability $1/m$. \n",
    "\n",
    "If $w=32$, then for $m>2^{32}=4.24e9$, some values will have probability 0. Until relatively recently, R did not support 64-bit integers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## R random integer generator / sample()\n",
    "\n",
    "From version 3.5.x and earlier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "if (dn > INT_MAX || k > INT_MAX) {\n",
    "\t    PROTECT(y = allocVector(REALSXP, k));\n",
    "\t    if (replace) {\n",
    "\t\tdouble *ry = REAL(y);\n",
    "\t\tfor (R_xlen_t i = 0; i < k; i++) ry[i] = floor(dn * ru() + 1);\n",
    " \t    } else {\n",
    "#ifdef LONG_VECTOR_SUPPORT\n",
    "\t\tR_xlen_t n = (R_xlen_t) dn;\n",
    "\t\tdouble *x = (double *)R_alloc(n, sizeof(double));\n",
    "\t\tdouble *ry = REAL(y);\n",
    "\t\tfor (R_xlen_t i = 0; i < n; i++) x[i] = (double) i;\n",
    "\t\tfor (R_xlen_t i = 0; i < k; i++) {\n",
    "\t\t    R_xlen_t j = (R_xlen_t)floor(n * ru());\n",
    "\t\t    ry[i] = x[j] + 1;\n",
    "\t\t    x[j] = x[--n];\n",
    "\t\t}\n",
    "#else\n",
    "\t\terror(_(\"n >= 2^31, replace = FALSE is only supported on 64-bit platforms\"));\n",
    "#endif\n",
    "\t    }\n",
    "\t} else {\n",
    "\t    int n = (int) dn;\n",
    "\t    PROTECT(y = allocVector(INTSXP, k));\n",
    "\t    int *iy = INTEGER(y);\n",
    "\t    /* avoid allocation for a single sample */\n",
    "\t    if (replace || k < 2) {\n",
    "\t\tfor (int i = 0; i < k; i++) iy[i] = (int)(dn * unif_rand() + 1);\n",
    "\t    } else {\n",
    "\t\tint *x = (int *)R_alloc(n, sizeof(int));\n",
    "\t\tfor (int i = 0; i < n; i++) x[i] = i;\n",
    "\t\tfor (int i = 0; i < k; i++) {\n",
    " \t\t    int j = (int)(n * unif_rand());\n",
    " \t\t    iy[i] = x[j] + 1;\n",
    " \t\t    x[j] = x[--n];\n",
    " \t\t}\n",
    " \t}\n",
    "```\n",
    "```\n",
    "/* Our PRNGs have at most 32 bit of precision, and all have at least 25 */\n",
    "\n",
    "static R_INLINE double ru()\n",
    "{\n",
    "    double U = 33554432.0;\n",
    "    return (floor(U*unif_rand()) + unif_rand())/U;\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a numerical example in R 3.5.1 (h/t [Duncan Murdoch](https://stat.ethz.ch/pipermail/r-devel/2018-September/076827.html)):\n",
    "\n",
    " > m <- (2/5)*2^32  \n",
    " > x <- sample(m, 1000000, replace = TRUE)  \n",
    " > table(x %% 2)  \n",
    " >     0      1  \n",
    " > 399850 600150\n",
    " \n",
    "See Ottoboni and Stark, [Random Problems with R](https://arxiv.org/abs/1809.06520) and\n",
    "Stark and Ottoboni, [Simple Random Sampling: Practice makes imperfect](https://arxiv.org/abs/1810.10985).\n",
    "\n",
    "R changed how it generates random integers starting with version 3.6.\n",
    "Here's the corresponding output for R version 4.0:\n",
    "\n",
    " > m <- (2/5)*2^32  \n",
    " > x <- sample(m, 1000000, replace = TRUE)  \n",
    " > table(x %% 2)  \n",
    " >     0      1  \n",
    " > 499958 500042  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### More accurate method\n",
    "A better way to generate a (pseudo-)random integer on $\\{1, \\ldots m\\}$ from a (pseudo-random) $w$-bit integer in practice is as follows:\n",
    "\n",
    "1. Set $\\mu = \\log_2(m-1)$.\n",
    "1. Define a $w$-bit _mask_ consisting of $\\mu$ bits set to 1 and $(w-\\mu)$ bits set to zero.\n",
    "1. Generate a random $w$-bit integer $Y$. \n",
    "1. Define $y$ to be the bitwise `and` of $Y$ and the mask.\n",
    "1. If $y \\le m-1$, output $x = y+1$; otherwise, return to step 3.\n",
    "\n",
    "This is how random integers are generated in numpy by `numpy.random.randint()`.\n",
    "However, `numpy.random.choice()` does something else that's biased: it finds the closest integer to $mX$.\n",
    "\n",
    "In `R`, one would generally use the function `sample(1:m, k, replace=FALSE)` to draw pseudo-random\n",
    "integers. \n",
    "Before version 3.6, `sample()` used a variant of the faulty `1 + floor(m*X)` approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How can you tell whether a sequence is random?\n",
    "\n",
    "http://dilbert.com/strip/2001-10-25"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tests of PRNGS\n",
    "\n",
    "+ Theoretical analyses, e.g., Knuth (1969), Marsaglia (1968)\n",
    "\n",
    "+ Statistical tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Knuth (1969) _The Art of Computer Programming_, v.2 \n",
    "\n",
    "+ 11 types of behavior: equidistribution, series, gaps, poker, coupon collector, permutation frequency, runs, max of $t$, collisions, birthday spacings, serial correlation\n",
    "+ tests on subsequences, spectral test\n",
    "+ Many $\\chi^2$-based tests\n",
    "+ Kolmogorov-Smirnov test for uniformity\n",
    "+ Sphere-packing\n",
    "+ **MORE**\n",
    "\n",
    "#### Marsaglia (1996) DIEHARD tests\n",
    "\n",
    "+ Birthday spacings\n",
    "+ Overlapping permutations of 5 random numbers\n",
    "+ Ranks of binary matrices of various dimensions\n",
    "+ Monkeys at typewriters: count overlapping \"words\" in strings of bits\n",
    "+ Count the 1s in bytes; translate to \"words.\"\n",
    "+ Parking lot test, 100 &times; 100 square lot. Count non-collisions.\n",
    "+ Minimum distance test: Min distance between 8,000 random points in a 10,000 &times; 10,000 square. \n",
    "+ Sphere-packing in a cube at random; diameter of smallest sphere.\n",
    "+ Squeeze test: Multiply 231 by random floats on (0,1) until hitting 1.\n",
    "+ Overlapping sums of 100 random (0,1) floats.\n",
    "+ Runs test for random floats\n",
    "+ #wins and #rolls in 200,000 games of craps\n",
    "\n",
    "#### L’Ecuyer and Simard (2007) TestU01 http://dl.acm.org/citation.cfm?doid=1268776.1268777\n",
    "\n",
    "+ Kolmogorov-Smirnov, Cr&aacute;mer-von Mises, Anderson-Darling, clustering, runs, gaps, hits in partition of a hypercube (collisions, empty cells, time between visits, ...), birthday spacings, close pairs, coupon collector, sum collector, complexity of bit strings (__linear complexity__, jump complexity, jump size complexity, Lempel-Ziv complexity), spectral tests on bit strings, autocorrelation of bits, runs and gaps in bits, ..., ranks of binary matrices, longest runs, Hamming weights, random walks, close pairs of binary sequences, \n",
    "\n",
    "### NIST Tests\n",
    "\n",
    "+ http://csrc.nist.gov/groups/ST/toolkit/rng/stats_tests.html\n",
    "+ http://csrc.nist.gov/groups/ST/toolkit/rng/documents/SP800-22rev1a.pdf\n",
    "\n",
    "![NIST](./RandPics/NISTTests.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pigeonhole principle\n",
    "\n",
    "If you put $N>n$ pigeons in $n$ pigeonholes, then at least one\n",
    "pigeonhole must contain more than one pigeon.\n",
    "\n",
    "## Corollary\n",
    "At most $n$ pigeons can be put in $n$ pigeonholes if at most\n",
    "one pigeon is put in each hole."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sizes of some coops & flocks\n",
    "\n",
    "| Expression     | full      |  scientific notation | \n",
    "|:---------------:|----------:|----------:|\n",
    "|$2^{32}$        | 4,294,967,296 | 4.29e9 |\n",
    "|$2^{64}$        |   18,446,744,073,709,551,616 | 1.84e19|\n",
    "|$2^{128}$       |            | 3.40e38 |\n",
    "|$2^{32 \\times 624}$ | | 9.27e6010|\n",
    "|$13!$      |                           6,227,020,800 | 6.23e9|\n",
    "|$21!$      |     51,090,942,171,709,440,000 | 5.11e19|                        \n",
    "|$35!$      |                                    | 1.03e40|\n",
    "|$2084!$     |                                                 | 3.73e6013 |\n",
    "|${50 \\choose 10}$   |        10,272,278,170 | 1.03e10|\n",
    "|${100 \\choose 10}$ | 17,310,309,456,440 | 1.73e13|\n",
    "|${500 \\choose 10}$ |         |     2.4581e20 |\n",
    "|$\\frac{2^{32}}{{50 \\choose 10}}$   | 0.418 | |\n",
    "|$\\frac{2^{64}}{{500 \\choose 10}}$ | 0.075 | |\n",
    "|$\\frac{2^{32}}{7000!}$ | | $<$ 1e-54,958 |\n",
    "|$\\frac{2}{52!}$ | | 2.48e-68 |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "5566f748-1f2c-452b-89cd-99a857c0b0fc"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example: frequency of samples using PIKK\n",
    "\n",
    "Use `PIKK` to generate samples using RANDU and MT.\n",
    "\n",
    "Tally empirical probability of each sample; test for uniformity by range.\n",
    "\n",
    "A (true) simple random sample of size $k$ from a population of size $n$ has chance $1/{n \\choose k}$ of\n",
    "resulting in each of the ${n \\choose k}$ possible subsets.\n",
    "\n",
    "The joint distribution of the number of times each _possible_ sample is selected in $N$\n",
    "independent random samples is multinomial with ${n \\choose k}$ categories,\n",
    "equal category probabilities $1/{n \\choose k}$, and $N$ draws.\n",
    "\n",
    "The pidgeonhold arguments prove that the actual distribution cannot be exactly multinomial, but\n",
    "how bad is the approximation?\n",
    "\n",
    "We can test the hypothesis using as the test statistic the range of category counts.\n",
    "\n",
    "Since we can't trust simulations to give an accurate result (that's the problem we are studying!), we need to rely on asymptotics.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "nbpresent": {
     "id": "bedd360c-6dfc-4070-b7e4-b3bcde5b550e"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def PIKK(n: int, k: int, gen: callable=np.random) -> set:\n",
    "    try: \n",
    "        x = gen.random(n)\n",
    "    except TypeError:\n",
    "        x = np.array([gen.random() for i in np.arange(n)])\n",
    "    return set(np.argsort(x)[0:k])\n",
    "\n",
    "# utility function\n",
    "\n",
    "def tallySamples(n, k, prng=Random.random, reps=int(10**7)):\n",
    "    uniqueSampleCounts = dict()\n",
    "    for i in np.arange(reps):\n",
    "        sam = frozenset(PIKK(n, k, prng))\n",
    "        if sam in uniqueSampleCounts:\n",
    "            uniqueSampleCounts[sam] += 1\n",
    "        else:\n",
    "            uniqueSampleCounts[sam] = 1\n",
    "    return uniqueSampleCounts\n",
    "\n",
    "def distrNormalRange(w, n):\n",
    "    '''\n",
    "    The CDF of the range of n IID standard normals evaluated at w\n",
    "    '''\n",
    "    innerInt = lambda x: norm.pdf(x)*(norm.cdf(x+w) - norm.cdf(x))**(n-1)\n",
    "    tmp = scipy.integrate.quad(innerInt, -np.inf, np.inf)\n",
    "    return n*tmp[0]\n",
    "\n",
    "def distrMultinomialRange(w, n, k):\n",
    "    '''\n",
    "    CDF of the range of multinomial variables, evaluated at w\n",
    "    n draws, k categories each having probability 1/k\n",
    "    '''\n",
    "    cutoff = (w - 1/(2*n))*np.sqrt(k/n)\n",
    "    return distrNormalRange(cutoff, k)\n",
    "\n",
    "def aveSamRange(sam): # empirical mean of the sample range\n",
    "    m = 0\n",
    "    totCnt = 0\n",
    "    for s, cnt in sam.items():\n",
    "        m += np.ptp([x[i] for i in s])*cnt\n",
    "        totCnt += cnt\n",
    "    return m/totCnt\n",
    "\n",
    "def aveSamVar(sam):  # empirical mean of the sample variance\n",
    "    m = 0\n",
    "    totCnt = 0\n",
    "    for s, cnt in sam.items():\n",
    "        m += np.var([x[i] for i in s], ddof=1)*cnt\n",
    "        totCnt += cnt\n",
    "    return m/totCnt\n",
    "\n",
    "def summarizeSams(counts, n, k, verbose=False):\n",
    "    reps = np.sum(counts)\n",
    "    vals = sp.special.binom(n, k)\n",
    "    obsVals = len(counts)\n",
    "    ect = reps/vals\n",
    "    maxct = np.amax(counts)\n",
    "    minct = np.amin(counts)\n",
    "    rangect = maxct-minct\n",
    "    pvalueRange = 1-distrMultinomialRange(rangect, reps, vals)\n",
    "    pvalueChi2 = chisquare(counts)\n",
    "    if verbose:\n",
    "        print(f'reps: {reps} possible values: {vals} observed values: {obsVals}' \\\n",
    "              f'expected count: {ect} min count: {minct} max count: {maxct}' \\\n",
    "              f'range: {rangect} p value for range: {pvalueRange}'\\\n",
    "              f'p value for chi-square: {pvalueChi2}')\n",
    "    return [reps, vals, obsVals, ect, maxct, minct, rangect, pvalueRange, pvalueChi2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "nbpresent": {
     "id": "18d48f30-b1be-4de6-af1b-fc93d71853bc"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# set up the simulation\n",
    "reps = int(10**7)\n",
    "n = 30\n",
    "k = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "nbpresent": {
     "id": "35399f1c-a4a3-4879-9313-b7bb93d61221"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# RANDU\n",
    "lcg = lcgRandom(seed=100) # set seed of RANDU to 100\n",
    "uniqueSamplesRU = tallySamples(n, k, lcg, reps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "nbpresent": {
     "id": "e3c029c4-4def-4342-8aa4-9ae6b4d4bb32"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reps: 10000000 possible values: 435.0 observed values: 435expected count: 22988.505747126437 min count: 21235 max count: 23533range: 2298 p value for range: -4.440892098500626e-16p value for chi-square: Power_divergenceResult(statistic=2780.769287, pvalue=0.0)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABB4AAAPxCAYAAABdE27qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqQklEQVR4nO3df8zu913X8dfb1uJPKLgznG1nixyIlSxjHrrGSMwgC+1CdohI0gXTZltsCisR1EjnEowaE2RGYmWuGaFhNYM5FeRoumzLNOAfFHaGY6xuZYcq9NDKCtOKqdJUPv5xrurN3fvc93Wfc17r3dPHI7lzruv7/Xy+1/d7n+vT6+y567rvWWsFAAAAoOH3vdAnAAAAAFy+hAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqufKFP4DBe9rKXreuvv/6FPg0AAABgh49//OO/udY6tte+F1V4uP7663P69OkX+jQAAACAHWbmV8+3z0ctAAAAgJqtwsPM3DIzj8zMmZm5Z4/9MzP3bvZ/cmZes2Pf/TPzuZn51K4575yZz2zG/+TMXH3RVwMAAAAcKQeGh5m5Ism7ktya5MYkb5qZG3cNuzXJ8c3XnUnevWPfjya5ZY9DfyTJ16y1XpXkl5O8/bAnDwAAABxt27zj4aYkZ9Zaj661nkny/iQnd405meSBdc5DSa6emVckyVrrZ5J8fvdB11ofXms9u7n7UJJrL/QiAAAAgKNpm/BwTZLHdtw/u9l22DH7eUuSD+61Y2bunJnTM3P6ySefPMQhAQAAgBfaNuFh9ti2LmDM3gefeUeSZ5O8b6/9a633rLVOrLVOHDu252/mAAAAAI6obX6d5tkk1+24f22Sxy9gzPPMzB1JvjnJN661tgoVAAAAwIvHNu94+FiS4zNzw8xcleS2JKd2jTmV5PbNb7e4OclTa60n9jvozNyS5HuTvHGt9fQFnDsAAABwxB0YHjY/APLuJB9K8ukkH1hrPTwzd83MXZthDyZ5NMmZJD+c5Dufmz8zP57kZ5N89cycnZm3bnb9UJI/muQjM/OJmbnvUl0UAAAAcDTMi+kTDidOnFinT59+oU8DAAAA2GFmPr7WOrHXvm0+agEAAABwQYQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGq2Cg8zc8vMPDIzZ2bmnj32z8zcu9n/yZl5zY5998/M52bmU7vmfNnMfGRmPrv580sv/nIAAACAo+TA8DAzVyR5V5Jbk9yY5E0zc+OuYbcmOb75ujPJu3fs+9Ekt+xx6HuSfHStdTzJRzf3AQAAgMvINu94uCnJmbXWo2utZ5K8P8nJXWNOJnlgnfNQkqtn5hVJstb6mSSf3+O4J5O8d3P7vUm+5QLOHwAAADjCtgkP1yR5bMf9s5tthx2z25evtZ5Iks2fL99r0MzcOTOnZ+b0k08+ucXpAgAAAEfFNuFh9ti2LmDMBVlrvWetdWKtdeLYsWOX4pAAAADAF8g24eFskut23L82yeMXMGa333ju4xibPz+3xbkAAAAALyLbhIePJTk+MzfMzFVJbktyateYU0lu3/x2i5uTPPXcxyj2cSrJHZvbdyT5qUOcNwAAAPAicGB4WGs9m+TuJB9K8ukkH1hrPTwzd83MXZthDyZ5NMmZJD+c5Dufmz8zP57kZ5N89cycnZm3bnZ9f5LXz8xnk7x+cx8AAAC4jMxal+RHMXxBnDhxYp0+ffqFPg0AAABgh5n5+FrrxF77tvmoBQAAAMAFER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqNkqPMzMLTPzyMycmZl79tg/M3PvZv8nZ+Y1B82dmVfPzEMz84mZOT0zN12aSwIAAACOigPDw8xckeRdSW5NcmOSN83MjbuG3Zrk+ObrziTv3mLuDyT5O2utVyf5vs19AAAA4DKyzTsebkpyZq316FrrmSTvT3Jy15iTSR5Y5zyU5OqZecUBc1eSL97c/pIkj1/ktQAAAABHzJVbjLkmyWM77p9N8totxlxzwNzvTvKhmfmHORdA/txeDz4zd+bcuyjyyle+covTBQAAAI6Kbd7xMHtsW1uO2W/udyT5nrXWdUm+J8mP7PXga633rLVOrLVOHDt2bIvTBQAAAI6KbcLD2STX7bh/bZ7/sYjzjdlv7h1JfmJz+1/k3McyAAAAgMvINuHhY0mOz8wNM3NVktuSnNo15lSS2ze/3eLmJE+ttZ44YO7jSf7C5vY3JPnsRV4LAAAAcMQc+DMe1lrPzszdST6U5Iok96+1Hp6Zuzb770vyYJI3JDmT5Okkb95v7ubQfyXJP56ZK5P872x+jgMAAABw+Zi1dv+4hqPrxIkT6/Tp0y/0aQAAAAA7zMzH11on9tq3zUctAAAAAC6I8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQI3wAAAAANQIDwAAwPO87nV7324/1mH2XeyxL9VjAPsTHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACo2So8zMwtM/PIzJyZmXv22D8zc+9m/ydn5jXbzJ2Z79rse3hmfuDiLwcAAAA4Sq48aMDMXJHkXUlen+Rsko/NzKm11n/aMezWJMc3X69N8u4kr91v7sy8LsnJJK9aa/3OzLz8Ul4YAAAA8MLb5h0PNyU5s9Z6dK31TJL351ww2OlkkgfWOQ8luXpmXnHA3O9I8v1rrd9JkrXW5y7B9QAAAABHyDbh4Zokj+24f3azbZsx+839qiRfPzM/NzM/PTNfd5gTBwAAAI6+Az9qkWT22La2HLPf3CuTfGmSm5N8XZIPzMxXrLV+z7Fn5s4kdybJK1/5yi1OFwAAADgqtnnHw9kk1+24f22Sx7ccs9/cs0l+YvPxjJ9P8rtJXrb7wdda71lrnVhrnTh27NgWpwsAAAAcFduEh48lOT4zN8zMVUluS3Jq15hTSW7f/HaLm5M8tdZ64oC5/zrJNyTJzHxVkquS/ObFXhAAAABwdBz4UYu11rMzc3eSDyW5Isn9a62HZ+auzf77kjyY5A1JziR5Osmb95u7OfT9Se6fmU8leSbJHbs/ZgEAAAC8uG3zMx6y1now5+LCzm337bi9krxt27mb7c8k+cuHOVkAAADgxWWbj1oAAAAAXBDhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEegBed173uhT4DADh6Dnp9/EK8fu71GM9t22/fYbcfdv/u7Qcd56B9+83Z5tiHeaz235t/V/GFIDwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADVbhYeZuWVmHpmZMzNzzx77Z2bu3ez/5My85hBz/8bMrJl52cVdCgAAAHDUHBgeZuaKJO9KcmuSG5O8aWZu3DXs1iTHN193Jnn3NnNn5rokr0/yaxd9JQAAAMCRs807Hm5Kcmat9eha65kk709ycteYk0keWOc8lOTqmXnFFnN/MMnfTLIu9kIAAACAo2eb8HBNksd23D+72bbNmPPOnZk3Jvn1tdYv7vfgM3PnzJyemdNPPvnkFqcLAAAAHBXbhIfZY9vudyicb8ye22fmDyV5R5LvO+jB11rvWWudWGudOHbs2IEnCwAAABwd24SHs0mu23H/2iSPbznmfNv/VJIbkvzizPyXzfZfmJk/fpiTBwAAAI62bcLDx5Icn5kbZuaqJLclObVrzKkkt29+u8XNSZ5aaz1xvrlrrV9aa718rXX9Wuv6nAsUr1lr/ddLdWEAAADAC+/KgwastZ6dmbuTfCjJFUnuX2s9PDN3bfbfl+TBJG9IcibJ00nevN/cypUAAAAAR86B4SFJ1loP5lxc2Lntvh23V5K3bTt3jzHXb3MeAAAAwIvLNh+1AAAAALggwgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAwAAAFAjPAAAAAA1wgMAAABQIzwAAAAANcIDAAAAUCM8AAAAADXCAy8pr3vddtsOc7yLmb/NcS/l8bc91qV+zIOOt3t/43va9GI73xcb31+4cBe6fs43r/X6dNBjHXRehx1zvvHnu31Yu79PF/P6u/s4ex33fI+x1/UcdLxtz3e/58Fex9p9Luc7t/2eY3tdw2Fs83za7zH2uobzfc+2+Xvbff+g7+n59u13Tdsc80IeZ69jbPMY2zxXz3cu+z3Xtj3H8x3nYjSOfzn+20d4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqhAcAAACgRngAAAAAaoQHAAAAoEZ4AAAAAGqEBwAAAKBGeAAAAABqtgoPM3PLzDwyM2dm5p499s/M3LvZ/8mZec1Bc2fmnTPzmc34n5yZqy/JFQEAAABHxoHhYWauSPKuJLcmuTHJm2bmxl3Dbk1yfPN1Z5J3bzH3I0m+Zq31qiS/nOTtF301AAAAwJGyzTsebkpyZq316FrrmSTvT3Jy15iTSR5Y5zyU5OqZecV+c9daH15rPbuZ/1CSay/B9QAAAABHyDbh4Zokj+24f3azbZsx28xNkrck+eBeDz4zd87M6Zk5/eSTT25xugAAAMBRsU14mD22rS3HHDh3Zt6R5Nkk79vrwdda71lrnVhrnTh27NgWpwsAAAAcFVduMeZskut23L82yeNbjrlqv7kzc0eSb07yjWut3TEDAAAAeJHb5h0PH0tyfGZumJmrktyW5NSuMaeS3L757RY3J3lqrfXEfnNn5pYk35vkjWutpy/R9QAAAABHyIHveFhrPTszdyf5UJIrkty/1np4Zu7a7L8vyYNJ3pDkTJKnk7x5v7mbQ/9Qki9K8pGZSZKH1lp3XcqLAwAAAF5Y23zUImutB3MuLuzcdt+O2yvJ27adu9n+lYc6UwAAAOBFZ5uPWgAAAABcEOEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgBrhAQAAAKgRHgAAAIAa4QEAAACoER4AAACAGuEBAAAAqBEeAAAAgJqtwsPM3DIzj8zMmZm5Z4/9MzP3bvZ/cmZec9DcmfmymfnIzHx28+eXXppLAgAAAI6KA8PDzFyR5F1Jbk1yY5I3zcyNu4bdmuT45uvOJO/eYu49ST661jqe5KOb+wAAAMBlZJt3PNyU5Mxa69G11jNJ3p/k5K4xJ5M8sM55KMnVM/OKA+aeTPLeze33JvmWi7sUAAAA4Ki5cosx1yR5bMf9s0leu8WYaw6Y++VrrSeSZK31xMy8fK8Hn5k7c+5dFEnyP2fmkS3OGc5rZs9tL0vym5fymJfCc8e9lMff9liX+poOOt7u/Ycd/0I7audzCV3U2rhULuPvLy9OR2JdbOtC189+8xqvT9uexzaPedjz2jn+fLcvxIUc6zz/Ttn6uNt8vw4zZpvH2rH9wLVx2O/JQc+1bc99vzm7xx/m+33QuRz2OPud+7brbou/p333HXZ9H3Rdl2LbhTwP9tvf/nfuHsc/9OvGi/TfPn/yfDu2CQ97XfLacsw2c/e11npPkvccZg4c1sycXmudeKHPA44aawOez7qAvVkbsDdrY7uPWpxNct2O+9cmeXzLMfvN/Y3NxzGy+fNz2582AAAA8GKwTXj4WJLjM3PDzFyV5LYkp3aNOZXk9s1vt7g5yVObj1HsN/dUkjs2t+9I8lMXeS0AAADAEXPgRy3WWs/OzN1JPpTkiiT3r7Uenpm7NvvvS/JgkjckOZPk6SRv3m/u5tDfn+QDM/PWJL+W5Nsu6ZXB4fg4D+zN2oDnsy5gb9YG7O0lvzZmrUP9yAUAAACArW3zUQsAAACACyI8AAAAADXCA5eNmbluZv79zHx6Zh6emb+62f5tm/u/OzMndoy/fmb+18x8YvN13459f3ZmfmlmzszMvTPnfpPuzHzRzPzzzfafm5nrv+AXCoewz7p458x8ZmY+OTM/OTNX75jz9s1z/JGZ+aYd260LLhuHXRteM3ip2Gdt/L3NuvjEzHx4Zv7EjjleN7isHXZdeM14PuGBy8mzSf76WutPJ7k5ydtm5sYkn0ryF5P8zB5zfmWt9erN1107tr87yZ1Jjm++btlsf2uS/7bW+sokP5jkH3QuBS6Z862LjyT5mrXWq5L8cpK3J8lm321J/kzOPe//6cxcsTmWdcHl5FBrY8NrBi8F51sb71xrvWqt9eok/zbJ9yVeN3jJONS62PCasYPwwGVjrfXEWusXNrd/O8mnk1yz1vr0WuuRbY8zM69I8sVrrZ9d53766gNJvmWz+2SS925u/8sk3/hcpYSjaJ918eG11rObYQ8luXZz+2SS96+1fmet9Z9z7rcV3WRdcLm5gLWxJ2uDy80+a+N/7Bj2h5M89xPqvW5w2buAdbGnl/K6EB64LG3emvS1SX7ugKE3zMx/nJmfnpmv32y7JsnZHWPObrY9t++x5Nyvi03yVJI/dqnOG5r2WRdvSfLBze3/9xzfeO75b11w2dpybSReM3iJ2b02Zubvz8xjSb49////2fW6wUvKlusi8ZrxewgPXHZm5o8k+VdJvntXhdztiSSvXGt9bZK/luTHZuaLk+xVFp+rl/vtgyPrfOtiZt6Rc28ffN9zm/aYvvbZvt8cOPIOsTa8ZvCSstfaWGu9Y611Xc6ti7ufG7rHdK8bXJYOsS68ZuwiPHBZmZnfn3P/MXjfWusn9hu7eUvgb21ufzzJryT5qpwrjzvfWnttksc3t88muW7zWFcm+ZIkn7+U1wCX2vnWxczckeSbk3z75u1+yY7n+MZzz3/rgsvOYdaG1wxeSrb499SPJfnWzW2vG7wkHGZdeM14PuGBy8bmM1A/kuTTa61/tMX4Y8/98KOZ+Yqc++Euj661nkjy2zNz8+aYtyf5qc20U0nu2Nz+S0n+3Y7/wQZHzvnWxczckuR7k7xxrfX0jimnkty2+cnKN+Tcuvh564LLzWHXhtcMXir2WRvHdwx7Y5LPbG573eCyd9h14TXj+eYyuhZe4mbmzyf5D0l+Kcnvbjb/rSRflOSfJDmW5L8n+cRa65tm5luT/N2ceyvt/0nyt9da/2ZzrBNJfjTJH8y5z/d+11przcwfSPLPcu5zXZ9Pctta69EvyAXCBdhnXdybc2vjtzbbHnruJy5v3mL+lpxbG9+91vrgZrt1wWXjsGvDawYvFfusjbcm+erNtl9Nctda69c3c7xucFk77LrwmvF8wgMAAABQ46MWAAAAQI3wAAAAANQIDwAAAECN8AAAAADUCA8AAABAjfAAAAAA1AgPAAAAQM3/Bb0qzS1K2uEQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1296x1296 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "[reps2, vals, obsVals, ect, maxct, minct, rangect, pvalueRange, pvalueChi2] =\\\n",
    "    summarizeSams(list(uniqueSamplesRU.values()), n, k, verbose=True)\n",
    "\n",
    "ns, bins, patches = plt.hist(list(uniqueSamplesRU.values()), maxct-minct+1, density=True, facecolor='blue', alpha=0.75)\n",
    "plt.rcParams['figure.figsize'] = (18.0, 18.0) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "bd8190f0-8cf4-494e-b112-88d8f150b711"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "RANDU with PIKK hits all 435 samples, but not equally often. Based on the chi-square test and the range test, the $p$-value for the hypothesis that these $10^7$ samples are simple random samples is essentially 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "43989d52-bcb4-4dc7-91fc-171dc10ee1d5"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Hypothesis tests for MT and RANDU\n",
    "\n",
    "$5 \\times 10^7$ replications, $n=30$, $k=2$.\n",
    "\n",
    "Sample variance of a population\n",
    "with 28 zeros, one 1, one -1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "nbpresent": {
     "id": "3b8762e7-ed54-45a8-9ca7-35a00df0478d"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.set_option('display.float_format', '{:.4g}'.format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "nbpresent": {
     "id": "7feffce4-fd19-4787-a799-954b33f4945a"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PRNG</th>\n",
       "      <th>seed</th>\n",
       "      <th>Chi-squared</th>\n",
       "      <th>Df</th>\n",
       "      <th>P-value</th>\n",
       "      <th>Range</th>\n",
       "      <th>Range P-value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>RANDU</td>\n",
       "      <td>100</td>\n",
       "      <td>[ 12155.0519878]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.]</td>\n",
       "      <td>8709</td>\n",
       "      <td>-4.441e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>100</td>\n",
       "      <td>[ 407.4165142]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.81558425]</td>\n",
       "      <td>2286</td>\n",
       "      <td>0.08578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>100</td>\n",
       "      <td>[ 409.6081138]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.79423828]</td>\n",
       "      <td>2196</td>\n",
       "      <td>0.171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>SD</td>\n",
       "      <td>100</td>\n",
       "      <td>[ 355.2852094]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.99767171]</td>\n",
       "      <td>2189</td>\n",
       "      <td>0.1798</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>MT</td>\n",
       "      <td>100</td>\n",
       "      <td>[ 454.6745488]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.23778712]</td>\n",
       "      <td>2174</td>\n",
       "      <td>0.1996</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>MT</td>\n",
       "      <td>429496729</td>\n",
       "      <td>[ 493.676701]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.02488693]</td>\n",
       "      <td>2089</td>\n",
       "      <td>0.341</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>233424280</td>\n",
       "      <td>[ 461.0437666]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.17825415]</td>\n",
       "      <td>2026</td>\n",
       "      <td>0.4741</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>MT</td>\n",
       "      <td>233424280</td>\n",
       "      <td>[ 429.9387088]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.54602838]</td>\n",
       "      <td>1965</td>\n",
       "      <td>0.6143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>429496729</td>\n",
       "      <td>[ 482.0878312]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.05502463]</td>\n",
       "      <td>1964</td>\n",
       "      <td>0.6166</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>233424280</td>\n",
       "      <td>[ 381.5704324]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.96666455]</td>\n",
       "      <td>1846</td>\n",
       "      <td>0.8545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>429496729</td>\n",
       "      <td>[ 451.5718156]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.27056103]</td>\n",
       "      <td>1792</td>\n",
       "      <td>0.9248</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SD</td>\n",
       "      <td>233424280</td>\n",
       "      <td>[ 289.7548432]</td>\n",
       "      <td>434</td>\n",
       "      <td>[ 0.99999999]</td>\n",
       "      <td>1422</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         PRNG       seed       Chi-squared   Df        P-value  Range  \\\n",
       "0       RANDU        100  [ 12155.0519878]  434          [ 0.]   8709   \n",
       "9   MT_choice        100    [ 407.4165142]  434  [ 0.81558425]   2286   \n",
       "6      SHA256        100    [ 409.6081138]  434  [ 0.79423828]   2196   \n",
       "1          SD        100    [ 355.2852094]  434  [ 0.99767171]   2189   \n",
       "3          MT        100    [ 454.6745488]  434  [ 0.23778712]   2174   \n",
       "5          MT  429496729     [ 493.676701]  434  [ 0.02488693]   2089   \n",
       "10  MT_choice  233424280    [ 461.0437666]  434  [ 0.17825415]   2026   \n",
       "4          MT  233424280    [ 429.9387088]  434  [ 0.54602838]   1965   \n",
       "8      SHA256  429496729    [ 482.0878312]  434  [ 0.05502463]   1964   \n",
       "7      SHA256  233424280    [ 381.5704324]  434  [ 0.96666455]   1846   \n",
       "11  MT_choice  429496729    [ 451.5718156]  434  [ 0.27056103]   1792   \n",
       "2          SD  233424280    [ 289.7548432]  434  [ 0.99999999]   1422   \n",
       "\n",
       "    Range P-value  \n",
       "0      -4.441e-16  \n",
       "9         0.08578  \n",
       "6           0.171  \n",
       "1          0.1798  \n",
       "3          0.1996  \n",
       "5           0.341  \n",
       "10         0.4741  \n",
       "4          0.6143  \n",
       "8          0.6166  \n",
       "7          0.8545  \n",
       "11         0.9248  \n",
       "2               1  "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "US = pd.read_csv('uniqueSampleSummary.csv')\n",
    "US_cols = ['PRNG', 'seed', 'Chi-squared', 'Df', 'P-value', 'Range', 'Range P-value']\n",
    "US[US_cols].sort_values(['Range P-value', 'P-value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "nbpresent": {
     "id": "5f5265f4-a888-4ec3-83ff-48545d4394c9"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PRNG</th>\n",
       "      <th>seed</th>\n",
       "      <th>Avg Sample Var</th>\n",
       "      <th>Var Bias</th>\n",
       "      <th>Var Bias/SE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>RANDU</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06879</td>\n",
       "      <td>-0.0001749</td>\n",
       "      <td>-6.462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06888</td>\n",
       "      <td>-8.172e-05</td>\n",
       "      <td>-3.019</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>100</td>\n",
       "      <td>0.0689</td>\n",
       "      <td>-6.341e-05</td>\n",
       "      <td>-2.343</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06892</td>\n",
       "      <td>-4.549e-05</td>\n",
       "      <td>-1.681</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SD</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06892</td>\n",
       "      <td>-4.337e-05</td>\n",
       "      <td>-1.602</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>MT</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06892</td>\n",
       "      <td>-4.176e-05</td>\n",
       "      <td>-1.543</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>429496729</td>\n",
       "      <td>0.06893</td>\n",
       "      <td>-3.774e-05</td>\n",
       "      <td>-1.394</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>MT</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06893</td>\n",
       "      <td>-3.407e-05</td>\n",
       "      <td>-1.259</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>SD</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06894</td>\n",
       "      <td>-2.605e-05</td>\n",
       "      <td>-0.9624</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06895</td>\n",
       "      <td>-2.036e-05</td>\n",
       "      <td>-0.7522</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>MT</td>\n",
       "      <td>429496729</td>\n",
       "      <td>0.06895</td>\n",
       "      <td>-2.012e-05</td>\n",
       "      <td>-0.7433</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>MT</td>\n",
       "      <td>429496729</td>\n",
       "      <td>0.06895</td>\n",
       "      <td>-1.746e-05</td>\n",
       "      <td>-0.645</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>429496729</td>\n",
       "      <td>0.06895</td>\n",
       "      <td>-1.397e-05</td>\n",
       "      <td>-0.5161</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06896</td>\n",
       "      <td>-8.697e-06</td>\n",
       "      <td>-0.3214</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>429496729</td>\n",
       "      <td>0.06896</td>\n",
       "      <td>-6.167e-06</td>\n",
       "      <td>-0.2279</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>SD</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06896</td>\n",
       "      <td>-3.447e-06</td>\n",
       "      <td>-0.1274</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SD</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06897</td>\n",
       "      <td>-4.472e-07</td>\n",
       "      <td>-0.01653</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>MT_choice</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06897</td>\n",
       "      <td>-1.372e-07</td>\n",
       "      <td>-0.005071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06897</td>\n",
       "      <td>1.453e-06</td>\n",
       "      <td>0.05368</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>MT</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06897</td>\n",
       "      <td>2.723e-06</td>\n",
       "      <td>0.1006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>429496729</td>\n",
       "      <td>0.06897</td>\n",
       "      <td>4.303e-06</td>\n",
       "      <td>0.159</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>RANDU</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06898</td>\n",
       "      <td>1.452e-05</td>\n",
       "      <td>0.5366</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>MT</td>\n",
       "      <td>100</td>\n",
       "      <td>0.06898</td>\n",
       "      <td>1.768e-05</td>\n",
       "      <td>0.6534</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>SHA256</td>\n",
       "      <td>233424280</td>\n",
       "      <td>0.06898</td>\n",
       "      <td>1.849e-05</td>\n",
       "      <td>0.6833</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         PRNG       seed  Avg Sample Var   Var Bias  Var Bias/SE\n",
       "1       RANDU        100         0.06879 -0.0001749       -6.462\n",
       "19  MT_choice        100         0.06888 -8.172e-05       -3.019\n",
       "13     SHA256        100          0.0689 -6.341e-05       -2.343\n",
       "15     SHA256  233424280         0.06892 -4.549e-05       -1.681\n",
       "3          SD        100         0.06892 -4.337e-05       -1.602\n",
       "7          MT        100         0.06892 -4.176e-05       -1.543\n",
       "23  MT_choice  429496729         0.06893 -3.774e-05       -1.394\n",
       "9          MT  233424280         0.06893 -3.407e-05       -1.259\n",
       "5          SD  233424280         0.06894 -2.605e-05      -0.9624\n",
       "21  MT_choice  233424280         0.06895 -2.036e-05      -0.7522\n",
       "11         MT  429496729         0.06895 -2.012e-05      -0.7433\n",
       "10         MT  429496729         0.06895 -1.746e-05       -0.645\n",
       "17     SHA256  429496729         0.06895 -1.397e-05      -0.5161\n",
       "20  MT_choice  233424280         0.06896 -8.697e-06      -0.3214\n",
       "22  MT_choice  429496729         0.06896 -6.167e-06      -0.2279\n",
       "4          SD  233424280         0.06896 -3.447e-06      -0.1274\n",
       "2          SD        100         0.06897 -4.472e-07     -0.01653\n",
       "18  MT_choice        100         0.06897 -1.372e-07    -0.005071\n",
       "12     SHA256        100         0.06897  1.453e-06      0.05368\n",
       "8          MT  233424280         0.06897  2.723e-06       0.1006\n",
       "16     SHA256  429496729         0.06897  4.303e-06        0.159\n",
       "0       RANDU        100         0.06898  1.452e-05       0.5366\n",
       "6          MT        100         0.06898  1.768e-05       0.6534\n",
       "14     SHA256  233424280         0.06898  1.849e-05       0.6833"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "biases = pd.read_csv('statBias.csv')\n",
    "biases_lfs = biases[biases['method']=='least freq sample']\n",
    "cols = ['PRNG', 'seed', 'Avg Sample Var', 'Var Bias', 'Var Bias/SE']\n",
    "biases[cols].sort_values(['Var Bias/SE', 'PRNG', 'seed'], ascending = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "9db040dc-5ddc-43b9-905a-bdfebc480f88"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## PRNGs in Common packages\n",
    "\n",
    "| Package/Lang   | default    | other      | SRS algorithm  |\n",
    "|:---------------|:----------:|:----------:|:--------------:|\n",
    "| SAS 9.2        | MT         | 32-bit LCG | Floyd's ordered hash or Fan et al. 1962 |\n",
    "| SPSS 20.0      | 32-bit LCG | MT1997ar   | trunc + rand indices |\n",
    "| SPSS &le; 12.0 | 32-bit LCG |            |                |\n",
    "| STATA 13       | KISS 32    |            | PIKK           |\n",
    "| STATA 14       | MT         |            | PIKK           |\n",
    "| R              | MT         |            | trunc + rand indices |\n",
    "| python         | MT         |            | mask + rand indices  |\n",
    "| MATLAB         | MT         |            | trunc + PIKK         |\n",
    "\n",
    "**Key.** MT = Mersenne Twister. LCG = linear congruential generator. PIKK = assign a number to each of the $n$ items and sort. The KISS generator combines 4 generators of three types: two multiply-with-carry\n",
    "generators, the 3-shift register SHR3 and the congruential generator CONG.\n",
    "\n",
    "Prior to April 2011, in Stata 10 exhibited predictable behavior:\n",
    "95.1% of the $2^{31}$ possible seed values resulted in the first and second draws from rnormal() having the same sign."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Is MT adequate for statistics?\n",
    "\n",
    "+ Know from pigeonhole argument that $L_1$ distance between true and desired is big for modest sampling & permutation problems.\n",
    "\n",
    "+ Know from equidistribution of MT that large ensemble frequencies will be right, but expect dependence issues\n",
    "\n",
    "+ Looking for problems that occur across seeds, large enough to be visible in $O(10^5)$ replications\n",
    "\n",
    "+ Examined simple random sample frequencies, derangements, partial derangements, Spearman correlation, etc.\n",
    "    \n",
    "+ _Must_ be problems, but where are they?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "8066555a-8af2-40b5-b264-a0fac9da86b7"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Best Practices\n",
    "\n",
    "+ Use a source of real randomness to set the seed with a substantial amount of entropy, e.g., 20 rolls of 10-sided dice.\n",
    "+ Record the seed so your analysis is reproducible.\n",
    "+ Use a PRNG at least as good as the Mersenne Twister, and preferably a cryptographically secure PRNG. Consider the PCG family. \n",
    "+ Avoid standard linear congruential generators and the Wichmann-Hill generator.\n",
    "+ Use open-source software, and record the version of the software.\n",
    "+ Use a sampling algorithm that does not \"waste randomness.\" Avoid permuting the entire population.\n",
    "+ Be aware of discretization issues in the sampling algorithm; many methods assume the PRNG produces $U[0,1]$ or $U[0,1)$ random numbers, rather than (an approximation to) numbers that are uniform on $w$-bit binary integers.\n",
    "+ Consider the size of the problem: are your PRNG and sampling algorithm adequate?\n",
    "+ Avoid \"tests of representativeness\" and procedures that reject some samples. They alter the distribution of the sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "e155dd7f-847c-4154-9b55-5ace34c45a67"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Recommendations\n",
    "\n",
    "+ Replace the standard PRNGs in R and Python with PRNGs with unbounded state spaces, and cryptographic or near-cryptographic quality\n",
    "    - Consider using AES in counter mode, since Intel chips have hardware support for AES\n",
    "+ Replace R's sample() with the mask algorithm **Done in 2019!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "8c651f97-ced8-491b-b960-1971e68f20cd"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## References\n",
    "\n",
    "1. Argyros, G. and A. Kiayias, 2012. PRNG: Pwning Random Number Generators.\n",
    "https://media.blackhat.com/bh-us-12/Briefings/Argyros/BH_US_12_Argyros_PRNG_WP.pdf\n",
    "1. Bernstein, D.J., T. Lange, and R. Niederhagen, 2016. Dual EC: A Standardized Back Door, in _The New Codebreakers, Essays Dedicated to David Kahn on the Occasion of his 85th Birthday_, Ryan, P.Y.A., D. Naccache, and J-J Quisquater, eds., Springer, Berlin.\n",
    "1. Cormen, T.H., C.E. Leiserson, R.L. Rivest and C. Stein, 2009. _Introduction to Algorithms, 3rd edition_, MIT Press.\n",
    "1. Fishman, G.S., and L.R. Moore, 1981. In Search of Correlation in Multiplicative Congruential Generators with Modulus 2\\*\\*31-1, _Computer Science and Statistics: Proceedings of the 13 Symposium on the Interface_, William F. Eddy, ed., Springer Verlag, New York.\n",
    "1. Knuth, D., 1997 _The Art of Computer Programming, V.II: Seminumerical methods_, 3rd edition, Addison-Wesley, Boston.\n",
    "1. L'Ecuyer, P. and R. Simard, 2007. TestU01: A C Library for Empirical Testing of Random Number Generators, _ACM Trans. Math. Softw._, _33_, http://doi.acm.org/10.1145/1268776.1268777\n",
    "1. Marsaglia, G., 1968. Random Numbers Fall Mainly in the Planes, _PNAS_, _61_, 25&ndash;28.\n",
    "1. Marsaglia, G., 2003. Xorshift RNGs. _Journal of Statistical Software_, _8_, 1&ndash;6.\n",
    "1. Matsumoto, M., and T. Nishimura, 1998. 8). Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator, _ACM Transactions on Modeling and Computer Simulation_, _8_, 3&ndash;30. doi:10.1145/272991.272995\n",
    "1. McCullough, B.D., 2008. Microsoft's 'Not the Wichmann-Hill' random number generator. _Computational Statistics and Data Analysis_, _52_ (10), 4587–4593. http://dx.doi.org/10.1016/j.csda.2008.03.006\n",
    "1. NIST Computer Security Division, _Random Number Generation_ http://csrc.nist.gov/groups/ST/toolkit/rng/\n",
    "1. O'Neill, M.E., 2015. PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation, submitted to _ACM Transactions on Mathematical Software_. http://www.pcg-random.org/pdf/toms-oneill-pcg-family-v1.02.pdf\n",
    "1. http://www.pcg-random.org/\n",
    "1. Shannon, C.E., 1948. A Mathematical Theory of Communication, _Bell System Technical Journal, 27_, 379&ndash;423, 623&ndash;656. \n",
    "1. Vitter, J.S., 1985. Random Sampling with a Reservoir, _ACM Transactions on Mathematical Software, 11_, 37&ndash;57.\n",
    "1. Wikipedia articles, including https://en.wikipedia.org/wiki/Mersenne_Twister, https://en.wikipedia.org/wiki/Linear_congruential_generator, https://en.wikipedia.org/wiki/Comparison_of_hardware_random_number_generators, https://en.wikipedia.org/wiki/Pseudorandom_number_generator, https://en.wikipedia.org/wiki/List_of_random_number_generators, https://en.wikipedia.org/wiki/Random_number_generator_attack "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}