1_theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Warp Simulation\n",
    "\n",
    "## Pseudorandom number generator (PRNG)\n",
    "PRNG is an algorithm to generate a sequence of numbers whose properties approximate the properties of sequences of random numbers. In fact, the sequence is not truly random, and it is completely determined by an initial value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# package importing\n",
    "import random\n",
    "from matplotlib import pyplot as plt\n",
    "import time\n",
    "import statistics\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "# for speed up\n",
    "from numba import jit\n",
    "\n",
    "# settings\n",
    "np.set_printoptions(threshold=sys.maxsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# operating system settings\n",
    "import os\n",
    "os.environ['NUMBAPRO_LIBDEVICE'] = '/usr/local/cuda-10.1/nvvm/libdevice'\n",
    "os.environ['NUMBAPRO_NVVM'] = '/usr/local/cuda-10.1/nvvm/lib64/libnvvm.so'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "521657541\n"
     ]
    }
   ],
   "source": [
    "import collections\n",
    "\n",
    "\"\"\"\"\n",
    "Implements the Subtractive PRNG used by C#: https://rosettacode.org/wiki/Subtractive_generator#C.2B.2B.\n",
    "seed is used to define a new generation for the PRNG. A new call of the srand function is always either a new generation based on that seed alone.\n",
    "\n",
    ":param seed: the seed used for generation. It has to be a value within 1-10^9. Keep in mind this parameter is not persisted between calls. One function call has its own instance. \n",
    "\"\"\"\n",
    "def srand(seed):\n",
    "    r = collections.deque(maxlen=55)\n",
    "    s = collections.deque(maxlen=55)\n",
    "    m = 10**9\n",
    "    s.extend([seed, 1])\n",
    "    s.extend((s[n-2] - s[n-1]) % m for n in range(2, 55))\n",
    "    r.extend(s[(34 * (n+1)) % 55] for n in range(55))\n",
    "    for _ in range(219-54):\n",
    "        r.append((r[0] - r[31]) % m)\n",
    "    while True:\n",
    "        r.append((r[0] - r[31]) % m)\n",
    "        yield r[54]\n",
    "\n",
    "# Verify our SPRNG works\n",
    "x = srand(1029202)\n",
    "print(next(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Calculates the probability between N times of iterations.\n",
    ":param seed: the seed used by the PRNG to calculate the rate of roll. This could be any positive integers from 1 to 10^9.\n",
    ":param rate: the rate of the roll. By default it uses 0.02.\n",
    ":param iterations: the N amount of iterations to perform to check for a probability\n",
    ":param total_samples: the amount of total samples to check against.\n",
    "\"\"\"\n",
    "# We're using Numba to speed up our increments, since the LLVM opcodes is going to help us alot.\n",
    "@jit(parallel=True, forceobj=True)\n",
    "def get_probability(seed, rate, rate_increment, iterations, total_samples, debug=False):\n",
    "\n",
    "    tl_mean_list = []\n",
    "\n",
    "    prng = srand(seed)\n",
    "    s = 0\n",
    "    start = time.time()\n",
    "\n",
    "    while s < iterations:\n",
    "        trial_list = []\n",
    "        a = 0\n",
    "\n",
    "        while a < total_samples:\n",
    "            trial = 1\n",
    "            x = next(prng) / 10**9\n",
    "            rate_current = rate\n",
    "            while x > rate_current:\n",
    "                trial += 1\n",
    "                x = next(prng) / 10**9\n",
    "                if trial > 50:\n",
    "                    rate_current += rate_increment\n",
    "                \n",
    "            a += 1\n",
    "            trial_list.append(trial)\n",
    "        \n",
    "        eee = sum(trial_list) / iterations\n",
    "        tl_mean_list.append(eee)\n",
    "        s += 1\n",
    "\n",
    "        if s % 100 == 0 and debug is True:\n",
    "            endf = time.time()\n",
    "            print('progress =','%.2f' % (s*100/iterations),'% (time =','%.2f' % (endf-start),'s)')\n",
    "  \n",
    "    return tl_mean_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean = 34.512800999999975\n",
      "variance = 0.45533237577477464\n",
      "sd = 0.6747832065002616\n",
      "[33.879, 35.863, 34.849, 33.938, 34.554, 33.218, 33.842, 35.805, 33.283, 34.776, 34.117, 35.658, 35.205, 34.833, 35.082, 34.177, 33.998, 34.016, 34.647, 35.536, 34.753, 34.902, 34.593, 34.625, 34.114, 35.035, 34.108, 35.364, 33.95, 35.957, 34.975, 34.421, 34.749, 34.302, 34.077, 35.049, 34.373, 33.882, 35.882, 34.999, 33.89, 34.813, 36.01, 35.178, 34.354, 33.225, 34.441, 34.493, 33.975, 33.271, 34.846, 34.696, 36.106, 35.175, 33.729, 34.596, 33.679, 34.348, 33.859, 34.088, 33.097, 34.36, 34.428, 34.716, 33.779, 35.93, 34.121, 34.667, 34.15, 35.432, 33.712, 34.624, 34.924, 34.994, 34.486, 34.515, 34.523, 33.951, 35.388, 34.937, 33.896, 34.288, 33.346, 34.73, 33.947, 34.885, 34.638, 35.431, 34.938, 33.755, 34.06, 34.441, 34.294, 33.639, 33.803, 35.232, 35.035, 34.822, 33.283, 35.15, 34.179, 34.226, 34.827, 34.323, 34.242, 34.502, 34.611, 34.107, 34.312, 34.09, 34.998, 34.126, 34.479, 34.785, 34.079, 35.079, 34.986, 34.189, 34.221, 34.917, 34.862, 34.154, 34.374, 34.817, 34.216, 34.591, 36.132, 34.692, 35.169, 34.898, 33.946, 34.916, 34.833, 34.532, 33.287, 33.654, 35.433, 34.471, 34.499, 33.818, 34.156, 33.886, 34.521, 34.275, 33.977, 35.287, 34.145, 34.535, 34.437, 35.003, 34.39, 34.456, 35.489, 34.444, 34.701, 34.875, 35.939, 33.477, 33.144, 36.406, 34.821, 33.791, 34.804, 34.016, 35.437, 34.274, 34.075, 34.652, 33.399, 35.188, 35.17, 34.859, 33.533, 34.806, 34.519, 33.092, 34.533, 34.631, 34.878, 33.269, 34.343, 35.651, 35.075, 34.278, 34.711, 34.201, 33.285, 35.016, 34.597, 35.149, 34.317, 35.496, 35.123, 33.077, 32.44, 34.993, 35.213, 34.133, 34.6, 34.307, 33.835, 34.394, 35.039, 33.591, 34.935, 35.96, 33.699, 33.382, 33.868, 34.02, 34.896, 34.603, 35.224, 34.852, 34.746, 33.504, 33.583, 35.04, 34.106, 34.073, 34.633, 35.162, 35.068, 34.094, 34.169, 34.53, 34.34, 33.539, 35.925, 35.062, 32.996, 34.51, 35.044, 35.017, 34.115, 34.387, 35.008, 35.775, 34.136, 34.18, 35.599, 34.147, 35.146, 33.988, 34.416, 34.974, 35.607, 33.784, 34.968, 34.696, 34.898, 34.305, 34.304, 34.542, 34.722, 34.738, 34.524, 34.163, 34.127, 35.933, 34.902, 33.947, 34.316, 34.363, 33.81, 34.033, 34.33, 34.817, 34.25, 33.709, 35.49, 35.157, 34.915, 34.246, 35.39, 35.312, 34.455, 33.423, 34.992, 35.712, 34.714, 34.389, 34.324, 34.386, 35.321, 34.719, 34.333, 32.981, 34.685, 33.869, 34.444, 35.503, 35.645, 35.718, 34.771, 35.302, 36.131, 35.344, 34.974, 35.262, 34.503, 33.6, 35.506, 34.015, 34.1, 33.766, 34.836, 34.621, 33.308, 35.134, 35.532, 33.681, 35.035, 34.616, 34.552, 34.701, 34.981, 34.685, 33.746, 33.96, 35.13, 34.286, 34.546, 34.285, 34.612, 35.37, 35.299, 34.659, 34.83, 34.916, 34.957, 34.012, 33.499, 33.525, 34.802, 34.742, 33.659, 33.694, 33.76, 33.72, 34.864, 34.965, 33.607, 33.965, 33.228, 34.01, 34.843, 34.514, 35.476, 35.706, 33.488, 34.729, 34.988, 33.637, 34.618, 35.571, 33.107, 33.759, 35.306, 34.831, 34.803, 35.577, 35.917, 34.903, 34.746, 34.635, 34.85, 34.713, 34.189, 35.051, 34.507, 34.071, 34.377, 34.005, 34.085, 33.412, 33.955, 33.802, 33.955, 35.161, 34.442, 34.253, 34.381, 33.718, 33.729, 33.435, 34.863, 33.637, 34.648, 35.238, 34.82, 35.686, 34.243, 34.901, 34.953, 35.421, 33.353, 36.604, 35.237, 33.268, 34.255, 35.394, 34.478, 33.339, 35.067, 35.721, 34.547, 34.732, 34.836, 33.835, 35.483, 33.82, 34.936, 34.984, 35.316, 34.071, 34.729, 34.324, 34.241, 32.523, 35.034, 34.781, 34.451, 35.013, 34.758, 34.44, 34.491, 34.278, 34.059, 33.88, 35.066, 35.216, 33.261, 34.897, 34.187, 34.218, 34.363, 35.627, 34.713, 34.915, 34.406, 33.951, 33.703, 35.908, 34.142, 35.307, 34.251, 34.416, 35.176, 34.784, 34.317, 34.354, 34.162, 34.878, 34.13, 34.843, 34.634, 34.71, 34.846, 35.939, 34.176, 34.928, 35.226, 34.123, 33.621, 35.388, 36.287, 34.735, 34.559, 34.179, 33.65, 34.28, 36.027, 34.607, 33.732, 34.736, 35.057, 33.944, 34.645, 35.728, 34.578, 34.462, 34.013, 34.351, 34.245, 34.363, 33.878, 35.515, 33.613, 34.273, 34.164, 34.278, 34.664, 35.02, 34.277, 35.15, 33.22, 35.148, 34.33, 33.887, 35.038, 33.095, 34.669, 32.674, 34.182, 35.716, 33.319, 34.871, 34.126, 34.858, 34.744, 33.662, 34.05, 35.129, 34.921, 35.084, 33.337, 34.51, 33.928, 35.095, 34.658, 34.521, 35.019, 34.26, 33.751, 34.908, 34.046, 34.839, 34.212, 34.661, 35.771, 35.223, 34.553, 35.099, 34.045, 34.821, 33.979, 34.438, 34.201, 35.125, 34.676, 34.551, 34.303, 35.293, 35.335, 34.049, 34.544, 34.242, 33.685, 33.767, 35.306, 33.306, 33.627, 34.179, 34.147, 34.764, 34.732, 34.366, 33.777, 34.132, 34.091, 34.641, 34.202, 34.435, 33.786, 34.206, 34.399, 35.169, 34.696, 34.505, 35.407, 35.39, 33.325, 35.557, 33.831, 35.37, 33.615, 34.102, 34.461, 35.224, 34.121, 33.609, 36.129, 34.397, 34.84, 33.68, 34.962, 33.516, 34.542, 34.657, 35.796, 34.673, 34.143, 34.663, 34.821, 34.754, 34.739, 34.111, 35.58, 34.148, 35.372, 35.284, 35.172, 34.212, 35.058, 34.506, 33.47, 35.317, 34.913, 34.698, 34.569, 34.605, 33.013, 34.209, 34.626, 35.584, 33.317, 35.326, 35.006, 34.262, 33.795, 35.19, 34.675, 34.651, 34.514, 34.711, 34.806, 33.62, 34.96, 35.434, 34.271, 34.682, 34.652, 35.143, 34.545, 34.91, 33.917, 34.668, 35.139, 34.073, 34.296, 34.916, 35.608, 35.845, 34.628, 34.193, 35.208, 35.01, 35.345, 33.435, 33.131, 34.164, 33.8, 34.264, 34.971, 34.624, 34.766, 33.239, 35.057, 35.797, 34.406, 33.659, 34.971, 35.042, 33.902, 34.786, 35.68, 35.073, 34.148, 35.077, 33.342, 34.843, 34.161, 34.189, 34.24, 33.66, 34.558, 34.672, 34.419, 34.354, 34.598, 34.79, 34.789, 34.686, 34.713, 33.434, 34.018, 35.183, 34.363, 33.281, 34.286, 33.515, 34.767, 35.25, 36.686, 34.507, 34.486, 34.964, 35.455, 33.361, 34.963, 34.169, 33.711, 34.701, 34.334, 34.864, 33.948, 33.938, 34.574, 34.798, 34.677, 32.788, 33.811, 34.889, 34.776, 33.73, 34.425, 34.69, 33.647, 35.115, 34.274, 33.957, 34.859, 34.677, 34.549, 34.217, 34.817, 33.98, 34.183, 33.746, 33.716, 33.458, 33.611, 34.744, 35.181, 34.129, 35.475, 34.061, 34.417, 36.089, 34.867, 34.525, 35.049, 34.85, 34.37, 34.802, 33.035, 35.162, 34.145, 34.156, 34.011, 34.761, 35.357, 34.514, 33.935, 34.657, 33.278, 34.143, 34.986, 35.644, 35.11, 35.709, 34.826, 35.178, 34.843, 33.882, 34.302, 34.448, 34.899, 33.219, 34.798, 34.681, 33.685, 33.736, 34.592, 33.006, 35.278, 34.231, 33.96, 34.458, 34.565, 34.258, 34.141, 35.056, 33.782, 34.677, 33.711, 34.946, 34.649, 34.629, 34.232, 34.756, 34.408, 33.927, 34.741, 34.994, 33.856, 35.016, 33.909, 34.668, 34.804, 33.762, 34.005, 33.428, 33.936, 35.743, 34.926, 34.273, 35.563, 33.67, 34.05, 33.084, 34.675, 34.757, 33.002, 34.012, 34.791, 34.425, 34.364, 34.387, 34.052, 33.849, 34.175, 35.377, 34.925, 34.214, 35.116, 34.708, 35.279, 35.18, 34.837, 34.182, 33.761, 33.456, 33.443, 34.59, 33.311, 35.525, 34.236, 33.316, 35.08, 34.865, 35.536, 33.675, 34.399, 35.154, 34.485, 34.696, 33.757, 33.722, 34.869, 34.535, 34.88, 34.114, 34.138, 33.81, 33.75, 33.691, 34.259, 34.736, 35.189, 34.915, 35.835, 35.093, 33.717, 35.996, 34.689, 34.104, 33.876, 33.951, 34.163, 33.531, 33.299, 34.132, 34.232, 34.885, 34.965, 35.882, 35.139, 35.167, 35.445, 34.591, 34.668, 34.868, 34.448, 34.046, 33.622, 35.044, 34.968, 35.543, 33.185, 33.465, 33.427, 34.872, 33.422, 34.891, 34.853, 34.619, 34.609, 34.064, 36.052, 34.95, 34.539, 33.94, 34.526, 34.42, 34.213, 35.032, 34.906, 34.366, 34.492, 33.937, 34.646, 34.217, 33.648, 34.818, 36.029, 35.428, 34.22, 34.641, 35.092, 34.21, 33.873, 34.314, 34.398, 35.445, 35.053, 33.452, 34.341, 34.609, 35.419, 35.643, 34.962, 35.474, 34.212, 35.846, 34.293, 35.777, 33.352, 34.313, 35.042, 35.309, 34.305, 34.498, 34.525, 34.755, 34.838, 34.273, 34.534, 34.122, 34.083, 34.142, 33.846, 33.806, 34.603, 33.293, 34.415, 34.378, 33.809, 34.039, 34.81, 33.694, 33.042, 34.345, 34.886, 34.168, 34.894, 34.126, 34.506, 34.315, 34.758, 33.774, 34.002, 34.005, 33.766, 34.277, 35.236, 34.349, 35.86, 33.968, 35.276, 35.048, 36.122, 35.274, 33.839, 35.136, 34.149, 34.976, 35.166, 34.114, 35.276, 34.058, 35.518, 33.907, 31.862, 36.038, 35.983, 34.26, 33.614, 34.039]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzYAAAJGCAYAAACJGhA0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAj/0lEQVR4nO3df3TV9X348VcwEKiQIBQSOAah1onOgi3tIGqdwwij6MFjulmrp87DdO1BOmCbJedYLbRHmHWFeQ5ofzisWzl09Ayn5QjH0SNdV0Ch86jtpNrKgRYS1lkSjCNwIN8/OvI1CsgNN9y85PE4555jPp/PvXnB+RjyPO/kvss6Ojo6AgAAILE+pR4AAADgVAkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrlpR7g7Y4cORK7d++OQYMGRVlZWanHAQAASqSjoyP2798fI0eOjD59Trwm0+vCZvfu3VFbW1vqMQAAgF5i165dce65557wml4XNoMGDYqI3w1fWVlZ4mkAAIBSaW1tjdra2s5GOJFeFzZHf/yssrJS2AAAACf1KyrePAAAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSKy/1AABwpho9f223nrdj8fQiTwKQnxUbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASK+81AMAkMfo+Wu79bwdi6cXeRIA6MqKDQAAkJ6wAQAA0hM2AABAegWHza9//eu45ZZbYujQoTFgwID40Ic+FFu3bu0839HREffcc0+MGDEiBgwYEPX19fHKK68UdWgAAIC3Kihsfvvb38bll18effv2jaeeeip+9rOfxd/93d/FOeec03nN/fffHw8++GA8/PDDsWXLljj77LNj6tSpceDAgaIPDwAAEFHgu6L97d/+bdTW1saKFSs6j40ZM6bzvzs6OmLp0qVx9913x4wZMyIi4rHHHovq6up4/PHH41Of+lSRxgYAAPj/ClqxeeKJJ+KjH/1o/Mmf/EkMHz48PvzhD8c3v/nNzvOvvfZaNDU1RX19feexqqqqmDhxYmzatOmYr9ne3h6tra1dHgAAAIUoKGx++ctfxkMPPRQXXHBBrF+/Pj73uc/F5z//+fj2t78dERFNTU0REVFdXd3ledXV1Z3n3m7RokVRVVXV+aitre3OnwMAADiDFRQ2R44ciY985CNx3333xYc//OG444474vbbb4+HH3642wM0NjZGS0tL52PXrl3dfi0AAODMVFDYjBgxIi6++OIuxy666KLYuXNnRETU1NRERERzc3OXa5qbmzvPvV1FRUVUVlZ2eQAAABSioLC5/PLLY/v27V2O/fznP4/zzjsvIn73RgI1NTWxYcOGzvOtra2xZcuWqKurK8K4AAAA71TQu6LNnTs3LrvssrjvvvviT//0T+PZZ5+Nb3zjG/GNb3wjIiLKyspizpw58ZWvfCUuuOCCGDNmTHzxi1+MkSNHxvXXX98T8wMAABQWNh/72MdizZo10djYGAsXLowxY8bE0qVL4+abb+685q677oq2tra44447Yt++fXHFFVfEunXron///kUfHgAAIKLAsImIuPbaa+Paa6897vmysrJYuHBhLFy48JQGAwAAOFkF/Y4NAABAb1Twig0AUFqj56/t1vN2LJ5e5EkAeg8rNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJBeeakHAIDjGT1/bbeet2Px9CJPAkBvZ8UGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0rOPDQA9rrv70QDAybJiAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOnZxwYAzhDd3U9ox+LpRZ4EoPis2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB65aUeAAB6i9Hz13breTsWTy/yJAAUyooNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApGcfGwDghLq7v0+EPX6A08eKDQAAkJ6wAQAA0hM2AABAegWFzZe+9KUoKyvr8hg7dmzn+QMHDsSsWbNi6NChMXDgwGhoaIjm5uaiDw0AAPBWBa/Y/P7v/37s2bOn8/GjH/2o89zcuXPjySefjNWrV8fGjRtj9+7dccMNNxR1YAAAgLcr+F3RysvLo6am5h3HW1pa4pFHHomVK1fG5MmTIyJixYoVcdFFF8XmzZtj0qRJpz4tAADAMRS8YvPKK6/EyJEj4wMf+EDcfPPNsXPnzoiI2LZtWxw6dCjq6+s7rx07dmyMGjUqNm3adNzXa29vj9bW1i4PAACAQhS0YjNx4sR49NFH48ILL4w9e/bEggUL4uMf/3i89NJL0dTUFP369YvBgwd3eU51dXU0NTUd9zUXLVoUCxYs6NbwANAbnMo+LwAUR0FhM23atM7/HjduXEycODHOO++8+Od//ucYMGBAtwZobGyMefPmdX7c2toatbW13XotAADgzHRKb/c8ePDg+L3f+7149dVXo6amJg4ePBj79u3rck1zc/MxfyfnqIqKiqisrOzyAAAAKMQphc0bb7wRv/jFL2LEiBExYcKE6Nu3b2zYsKHz/Pbt22Pnzp1RV1d3yoMCAAAcT0E/ivbXf/3Xcd1118V5550Xu3fvjnvvvTfOOuusuOmmm6KqqipmzpwZ8+bNiyFDhkRlZWXMnj076urqvCMaAADQowoKm1/96ldx0003xf/8z//EsGHD4oorrojNmzfHsGHDIiJiyZIl0adPn2hoaIj29vaYOnVqLF++vEcGBwAAOKqgsFm1atUJz/fv3z+WLVsWy5YtO6WhAAAACnFKv2MDAADQGwgbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6ZWXegAAKLbR89eWegQATjMrNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgvVMKm8WLF0dZWVnMmTOn89iBAwdi1qxZMXTo0Bg4cGA0NDREc3Pzqc4JAABwXN0Om+eeey6+/vWvx7hx47ocnzt3bjz55JOxevXq2LhxY+zevTtuuOGGUx4UAADgeLoVNm+88UbcfPPN8c1vfjPOOeeczuMtLS3xyCOPxNe+9rWYPHlyTJgwIVasWBE//vGPY/Pmzcd8rfb29mhtbe3yAAAAKES3wmbWrFkxffr0qK+v73J827ZtcejQoS7Hx44dG6NGjYpNmzYd87UWLVoUVVVVnY/a2trujAQAAJzBCg6bVatWxU9+8pNYtGjRO841NTVFv379YvDgwV2OV1dXR1NT0zFfr7GxMVpaWjofu3btKnQkAADgDFdeyMW7du2Kv/zLv4ynn346+vfvX5QBKioqoqKioiivBQAAnJkKWrHZtm1b7N27Nz7ykY9EeXl5lJeXx8aNG+PBBx+M8vLyqK6ujoMHD8a+ffu6PK+5uTlqamqKOTcAAECnglZsrr766njxxRe7HLvtttti7Nix8YUvfCFqa2ujb9++sWHDhmhoaIiIiO3bt8fOnTujrq6ueFMDAAC8RUFhM2jQoLjkkku6HDv77LNj6NChncdnzpwZ8+bNiyFDhkRlZWXMnj076urqYtKkScWbGgAA4C0KCpuTsWTJkujTp080NDREe3t7TJ06NZYvX17sTwMAANCprKOjo6PUQ7xVa2trVFVVRUtLS1RWVpZ6HADeYvT8taUegWR2LJ5e6hGAxAppg6Kv2AAAnKruRrSQgjNXtzboBAAA6E2EDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9+9gAnIFstAnAe40VGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOmVl3oAgPeK0fPXdut5OxZPL/IkAHDmsWIDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6dnHBgDoMd3d3wmgUFZsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID372AAkZX8QAPj/rNgAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAeuWlHgCgtxk9f22pRwAACmTFBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID0hA0AAJCesAEAANKzjw1Aidk3BwBOnRUbAAAgPWEDAACkJ2wAAID0Cgqbhx56KMaNGxeVlZVRWVkZdXV18dRTT3WeP3DgQMyaNSuGDh0aAwcOjIaGhmhubi760AAAAG9VUNice+65sXjx4ti2bVts3bo1Jk+eHDNmzIif/vSnERExd+7cePLJJ2P16tWxcePG2L17d9xwww09MjgAAMBRZR0dHR2n8gJDhgyJr371q/HJT34yhg0bFitXroxPfvKTERHx8ssvx0UXXRSbNm2KSZMmndTrtba2RlVVVbS0tERlZeWpjAbQLd6lDPLasXh6qUcAiqiQNuj279gcPnw4Vq1aFW1tbVFXVxfbtm2LQ4cORX19fec1Y8eOjVGjRsWmTZuO+zrt7e3R2tra5QEAAFCIgsPmxRdfjIEDB0ZFRUV89rOfjTVr1sTFF18cTU1N0a9fvxg8eHCX66urq6Opqem4r7do0aKoqqrqfNTW1hb8hwAAAM5sBYfNhRdeGM8//3xs2bIlPve5z8Wtt94aP/vZz7o9QGNjY7S0tHQ+du3a1e3XAgAAzkzlhT6hX79+8cEPfjAiIiZMmBDPPfdc/P3f/33ceOONcfDgwdi3b1+XVZvm5uaoqak57utVVFRERUVF4ZMDAAD8n1Pex+bIkSPR3t4eEyZMiL59+8aGDRs6z23fvj127twZdXV1p/ppAAAAjqugFZvGxsaYNm1ajBo1Kvbv3x8rV66MZ555JtavXx9VVVUxc+bMmDdvXgwZMiQqKytj9uzZUVdXd9LviAYAANAdBYXN3r174zOf+Uzs2bMnqqqqYty4cbF+/fq45pprIiJiyZIl0adPn2hoaIj29vaYOnVqLF++vEcGBwAAOOqU97EpNvvYAKVmHxvIyz428N5yWvaxAQAA6C2EDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPTKSz0AAECxjJ6/tlvP27F4epEnAU43KzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkVFDaLFi2Kj33sYzFo0KAYPnx4XH/99bF9+/Yu1xw4cCBmzZoVQ4cOjYEDB0ZDQ0M0NzcXdWgAAIC3KihsNm7cGLNmzYrNmzfH008/HYcOHYopU6ZEW1tb5zVz586NJ598MlavXh0bN26M3bt3xw033FD0wQEAAI4qL+TidevWdfn40UcfjeHDh8e2bdviyiuvjJaWlnjkkUdi5cqVMXny5IiIWLFiRVx00UWxefPmmDRpUvEmBwAA+D+n9Ds2LS0tERExZMiQiIjYtm1bHDp0KOrr6zuvGTt2bIwaNSo2bdp0zNdob2+P1tbWLg8AAIBCFLRi81ZHjhyJOXPmxOWXXx6XXHJJREQ0NTVFv379YvDgwV2ura6ujqampmO+zqJFi2LBggXdHQMogdHz13breTsWTy/yJAAAv9PtFZtZs2bFSy+9FKtWrTqlARobG6OlpaXzsWvXrlN6PQAA4MzTrRWbO++8M77//e/HD3/4wzj33HM7j9fU1MTBgwdj3759XVZtmpubo6am5pivVVFRERUVFd0ZAwAAICIKXLHp6OiIO++8M9asWRM/+MEPYsyYMV3OT5gwIfr27RsbNmzoPLZ9+/bYuXNn1NXVFWdiAACAtyloxWbWrFmxcuXK+Nd//dcYNGhQ5+/NVFVVxYABA6KqqipmzpwZ8+bNiyFDhkRlZWXMnj076urqvCMaAADQYwoKm4ceeigiIq666qoux1esWBF/9md/FhERS5YsiT59+kRDQ0O0t7fH1KlTY/ny5UUZFgAA4FgKCpuOjo53vaZ///6xbNmyWLZsWbeHAgAAKMQp7WMDAADQG3R7HxsAgDOdfb2g97BiAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOnZxwYAOON1dz8aoPewYgMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADp2ccG6PW6u7/EjsXTizwJANBbWbEBAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9MpLPQBw5hg9f22pRwBIrbtfR3csnl7kSaD3sWIDAACkJ2wAAID0hA0AAJCesAEAANITNgAAQHrCBgAASE/YAAAA6dnHBgDgNLOvFxSfFRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIr7zUAwAA8N4yev7abj1vx+LpRZ6EM4kVGwAAID1hAwAApCdsAACA9AoOmx/+8Idx3XXXxciRI6OsrCwef/zxLuc7OjrinnvuiREjRsSAAQOivr4+XnnllWLNCwAA8A4Fh01bW1uMHz8+li1bdszz999/fzz44IPx8MMPx5YtW+Lss8+OqVOnxoEDB055WAAAgGMp+F3Rpk2bFtOmTTvmuY6Ojli6dGncfffdMWPGjIiIeOyxx6K6ujoef/zx+NSnPnVq0wIAABxDUX/H5rXXXoumpqaor6/vPFZVVRUTJ06MTZs2HfM57e3t0dra2uUBAABQiKKGTVNTU0REVFdXdzleXV3dee7tFi1aFFVVVZ2P2traYo4EAACcAUr+rmiNjY3R0tLS+di1a1epRwIAAJIpatjU1NRERERzc3OX483NzZ3n3q6ioiIqKyu7PAAAAApR1LAZM2ZM1NTUxIYNGzqPtba2xpYtW6Kurq6YnwoAAKBTwe+K9sYbb8Srr77a+fFrr70Wzz//fAwZMiRGjRoVc+bMia985StxwQUXxJgxY+KLX/xijBw5Mq6//vpizg0AANCp4LDZunVr/NEf/VHnx/PmzYuIiFtvvTUeffTRuOuuu6KtrS3uuOOO2LdvX1xxxRWxbt266N+/f/GmBgAAeIuCw+aqq66Kjo6O454vKyuLhQsXxsKFC09pMAAAgJNV8ndFAwAAOFUFr9gAZDF6/tpSjwAAnCZWbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9+9jAe0B392vZsXh6kScBoDeyrxdnAis2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ6wAQAA0hM2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkF55qQcAAICIiNHz13b7uTsWTy/iJGRkxQYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQnrABAADSs48N9CKn8v79AABnMis2AABAesIGAABIT9gAAADpCRsAACA9YQMAAKQnbAAAgPSEDQAAkJ59bAAASK+7e8HtWDy9yJNQKlZsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID372MAZrLvv+Q8A0NtYsQEAANITNgAAQHrCBgAASE/YAAAA6QkbAAAgPWEDAACkJ2wAAID07GMDAADvUd3ds27H4ulFnqTnWbEBAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9OxjcxLOpPf/LkR3/166q7t/n1nmBABOv9P9fZ7vK3uOFRsAACA9YQMAAKQnbAAAgPR6LGyWLVsWo0ePjv79+8fEiRPj2Wef7alPBQAAnOF6JGy++93vxrx58+Lee++Nn/zkJzF+/PiYOnVq7N27tyc+HQAAcIbrkbD52te+FrfffnvcdtttcfHFF8fDDz8c73vf++If/uEfeuLTAQAAZ7iiv93zwYMHY9u2bdHY2Nh5rE+fPlFfXx+bNm16x/Xt7e3R3t7e+XFLS0tERLS2thZ7tG470v5mt57Xm/4MPaG7fy/d1d2/zyxzRpz+WQGA7jnd35dk+XzFdnSOjo6Od7226GHzm9/8Jg4fPhzV1dVdjldXV8fLL7/8jusXLVoUCxYseMfx2traYo922lUtLfUE7y1Z/j6zzAkAdN/p/vf+vf753s3+/fujqqrqhNeUfIPOxsbGmDdvXufHR44ciddffz2GDh0aZWVlJZysZ7S2tkZtbW3s2rUrKisrSz0OSbmPKAb3EcXgPqIY3EccT0dHR+zfvz9Gjhz5rtcWPWze//73x1lnnRXNzc1djjc3N0dNTc07rq+oqIiKioouxwYPHlzssXqdyspK/+NyytxHFIP7iGJwH1EM7iOO5d1Wao4q+psH9OvXLyZMmBAbNmzoPHbkyJHYsGFD1NXVFfvTAQAA9MyPos2bNy9uvfXW+OhHPxp/8Ad/EEuXLo22tra47bbbeuLTAQAAZ7geCZsbb7wx/vu//zvuueeeaGpqiksvvTTWrVv3jjcUOBNVVFTEvffe+44fv4NCuI8oBvcRxeA+ohjcRxRDWcfJvHcaAABAL9YjG3QCAACcTsIGAABIT9gAAADpCRsAACA9YQMAAKQnbHrAQw89FOPGjevcPbeuri6eeuqpiIh4/fXXY/bs2XHhhRfGgAEDYtSoUfH5z38+WlpaSjw1vc2J7qOIiL/4i7+I888/PwYMGBDDhg2LGTNmxMsvv1zCiemN3u0+OqqjoyOmTZsWZWVl8fjjj5/+QenV3u0+uuqqq6KsrKzL47Of/WwJJ6Y3OpmvR5s2bYrJkyfH2WefHZWVlXHllVfG//7v/5ZoYrIRNj3g3HPPjcWLF8e2bdti69atMXny5JgxY0b89Kc/jd27d8fu3bvjgQceiJdeeikeffTRWLduXcycObPUY9PLnOg+ioiYMGFCrFixIv7rv/4r1q9fHx0dHTFlypQ4fPhwiSenN3m3++iopUuXRllZWYmmpLc7mfvo9ttvjz179nQ+7r///hJOTG/0bvfRpk2b4o//+I9jypQp8eyzz8Zzzz0Xd955Z/Tp49tVTo59bE6TIUOGxFe/+tVjBszq1avjlltuiba2tigv75E9U3mPONF99MILL8T48ePj1VdfjfPPP78E05HF2++j559/Pq699trYunVrjBgxItasWRPXX399aYek13vrfXTVVVfFpZdeGkuXLi31WCTz1vto0qRJcc0118SXv/zlUo9FUhK4hx0+fDhWrVoVbW1tUVdXd8xrWlpaorKyUtRwXO92H7W1tcWKFStizJgxUVtbW4IJyeBY99Gbb74Zn/70p2PZsmVRU1NT4gnJ4Hhfj77zne/E+9///rjkkkuisbEx3nzzzRJOSW/39vto7969sWXLlhg+fHhcdtllUV1dHX/4h38YP/rRj0o9Kon4TrqHvPjii1FXVxcHDhyIgQMHxpo1a+Liiy9+x3W/+c1v4stf/nLccccdJZiS3u7d7qPly5fHXXfdFW1tbXHhhRfG008/Hf369SvhxPRGJ7qP5s6dG5dddlnMmDGjxFPS253oPvr0pz8d5513XowcOTJeeOGF+MIXvhDbt2+Pf/mXfynx1PQ2x7uPNm/eHBERX/rSl+KBBx6ISy+9NB577LG4+uqr46WXXooLLrigxJOTgR9F6yEHDx6MnTt3RktLS3zve9+Lb33rW7Fx48Yu35S2trbGNddcE0OGDIknnngi+vbtW8KJ6Y3e7T5qaWmJvXv3xp49e+KBBx6IX//61/Ef//Ef0b9//xJPTm9yvPvo1Vdfjb/6q7+K//zP/4yBAwdGRERZWZkfReOYTubftaN+8IMfxNVXX+1HY3mH491H+/bti8svvzwaGxvjvvvu67x+3LhxMX369Fi0aFEJpyYLYXOa1NfXx/nnnx9f//rXIyJi//79MXXq1Hjf+94X3//+930jykl5+330VgcPHoxzzjknvvWtb8VNN91UgunI4uh9NGDAgHjwwQe7/GLu4cOHo0+fPvHxj388nnnmmdINSa93oq9HbW1tMXDgwFi3bl1MnTq1BNORxdH7aP78+fGBD3wg/vEf/zFuueWWzvM33nhjlJeXx3e+850STkkWfhTtNDly5Ei0t7dHxO9WaqZOnRoVFRXxxBNPiBpO2lvvo7fr6OiIjo6O456Ho47eRwsWLIg///M/73LuQx/6UCxZsiSuu+66Ek1HFif6evT8889HRMSIESNO40RkdPQ+Gj16dIwcOTK2b9/e5fzPf/7zmDZtWommIxth0wMaGxtj2rRpMWrUqNi/f3+sXLkynnnmmVi/fn20trbGlClT4s0334x/+qd/itbW1mhtbY2IiGHDhsVZZ51V4unpLU50H/3yl7+M7373uzFlypQYNmxY/OpXv4rFixfHgAED4hOf+ESpR6cXOdF9VFNTc8w3DBg1alSMGTOmBNPSW53oPvrFL34RK1eujE984hMxdOjQeOGFF2Lu3Llx5ZVXxrhx40o9Or3Iie6jsrKy+Ju/+Zu49957Y/z48XHppZfGt7/97Xj55Zfje9/7XqlHJwlh0wP27t0bn/nMZ2LPnj1RVVUV48aNi/Xr18c111wTzzzzTGzZsiUiIj74wQ92ed5rr70Wo0ePLsHE9EYnuo92794d//7v/x5Lly6N3/72t1FdXR1XXnll/PjHP47hw4eXenR6kRPdR3CyTnQf7dq1K/7t3/4tli5dGm1tbVFbWxsNDQ1x9913l3psepl3+3o0Z86cOHDgQMydOzdef/31GD9+fDz99NN+T4uT5ndsAACA9OxjAwAApCdsAACA9IQNAACQnrABAADSEzYAAEB6wgYAAEhP2AAAAOkJGwAAID1hAwAApCdsAACA9IQNAACQ3v8DYeKwPwPBz2MAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Number of iterations\n",
    "N = 1000\n",
    "# Sample Size\n",
    "sN = 1000\n",
    "\n",
    "# Execute. Keep in mind it may take a while for the function to complete since it's going through a large number of iterations.\n",
    "tl_mean_list = get_probability(102927162, 0.02, 0.02, N, sN)\n",
    "\n",
    "var = statistics.variance(tl_mean_list)\n",
    "\n",
    "print('mean =',sum(tl_mean_list)/sN)\n",
    "print('variance =',var)\n",
    "print('sd =',np.sqrt(var))\n",
    "print(tl_mean_list)\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (10,7))\n",
    "ax.hist(tl_mean_list, bins = 50)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "31.862\n",
      "36.686\n"
     ]
    }
   ],
   "source": [
    "print(min(tl_mean_list))\n",
    "print(max(tl_mean_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean = 34.43346938775508\n",
      "variance = 1.9873416565656568\n",
      "sd = 1.4097310582397116\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy0AAAJGCAYAAABWVGorAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhlElEQVR4nO3dfZBV9X348c/KwgWTZRWE4BYWjDYSUEiCxBJTjQlBGELROBmTMSkam2nNEqO0tG5bY2nGLEwSi2MdtI4PZBI0YSaoiRNspALJKI1gaKR18CE+rMpDx4ZdweTKsPf3R+v+3MjTvbt397Ps6zVzZryHczif3S8XeXPuvdSUSqVSAAAAJHVcXw8AAABwOKIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkFptb1+wo6MjXn311airq4uamprevjwAAJBEqVSK119/PRoaGuK44w59P6XXo+XVV1+NcePG9fZlAQCApFpbW2Ps2LGH/PFej5a6urqI+N/Bhg8f3tuXBwAAkmhvb49x48Z1NsKh9Hq0vPWSsOHDh4sWAADgiG8b8UZ8AAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDUyoqWAwcOxHXXXRennHJKDBs2LE499dT4+te/HqVSqVrzAQAAA1xtOQcvW7YsVqxYEStXrozJkyfH5s2b4/LLL4/6+vq46qqrqjUjAAAwgJUVLY8++mjMnz8/5s6dGxEREyZMiHvuuSd+8YtfVGU4AACAsl4e9pGPfCTWrVsXTz/9dERE/Md//Ef8/Oc/jzlz5hzynGKxGO3t7V02AACAo1XWnZZrr7022tvbY+LEiTFo0KA4cOBA3HDDDXHppZce8pyWlpZYsmRJtwcF4Ngy4doHKzrvhaVze3gSALIr607LD37wg/je974Xq1atiieeeCJWrlwZ3/rWt2LlypWHPKe5uTna2to6t9bW1m4PDQAADBxl3WlZvHhxXHvttfHZz342IiLOPPPMePHFF6OlpSUWLFhw0HMKhUIUCoXuTwoAAAxIZd1peeONN+K447qeMmjQoOjo6OjRoQAAAN5S1p2WefPmxQ033BCNjY0xefLk+OUvfxk33nhjfPGLX6zWfAAAwABXVrTcfPPNcd1118WXv/zl2L17dzQ0NMSf//mfx9e+9rVqzQcAAAxwZUVLXV1dLF++PJYvX16lcQAAALoq6z0tAAAAvU20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDUyoqWCRMmRE1NzTu2pqamas0HAAAMcLXlHPz444/HgQMHOh9v27YtPvnJT8ZnPvOZHh8MAAAgosxoGTVqVJfHS5cujVNPPTXOO++8Hh0KAADgLWVFy9u9+eab8d3vfjcWLVoUNTU1hzyuWCxGsVjsfNze3l7pJQEAgAGo4jfi33fffbFnz5647LLLDntcS0tL1NfXd27jxo2r9JIAAMAAVHG03HHHHTFnzpxoaGg47HHNzc3R1tbWubW2tlZ6SQAAYACq6OVhL774Yjz88MPxwx/+8IjHFgqFKBQKlVwGAACgsjstd911V4wePTrmzp3b0/MAAAB0UXa0dHR0xF133RULFiyI2tqK38cPAABwVMqOlocffjheeuml+OIXv1iNeQAAALoo+1bJrFmzolQqVWMWAACAd6j408MAAAB6g2gBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAILWyo+WVV16Jz3/+8zFy5MgYNmxYnHnmmbF58+ZqzAYAABC15Rz8m9/8Js4555w4//zz4yc/+UmMGjUqnnnmmTjxxBOrNR8AADDAlRUty5Yti3HjxsVdd93Vue+UU0457DnFYjGKxWLn4/b29jJHBAAABrKyXh72wAMPxFlnnRWf+cxnYvTo0fHBD34wbr/99sOe09LSEvX19Z3buHHjujUwAAAwsJQVLb/+9a9jxYoV8Yd/+Ifx0EMPxZVXXhlXXXVVrFy58pDnNDc3R1tbW+fW2tra7aEBAICBo6yXh3V0dMRZZ50V3/jGNyIi4oMf/GBs27Ytbr311liwYMFBzykUClEoFLo/KQAAMCCVdafl5JNPjkmTJnXZ9/73vz9eeumlHh0KAADgLWVFyznnnBPbt2/vsu/pp5+O8ePH9+hQAAAAbykrWq655prYtGlTfOMb34hnn302Vq1aFf/yL/8STU1N1ZoPAAAY4MqKlunTp8eaNWvinnvuiTPOOCO+/vWvx/Lly+PSSy+t1nwAAMAAV9Yb8SMiPvWpT8WnPvWpaswCAADwDmXdaQEAAOhtogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApFZWtPzDP/xD1NTUdNkmTpxYrdkAAACittwTJk+eHA8//PD//wlqy/4pAAAAjlrZxVFbWxtjxow56uOLxWIUi8XOx+3t7eVeEgAAGMDKjpZnnnkmGhoaYujQoTFjxoxoaWmJxsbGQx7f0tISS5Ys6daQABy9Cdc+2KvXe2Hp3F69Hv1fpb9G/Vo7uO48531P6S/Kek/L2WefHXfffXesXbs2VqxYEc8//3z88R//cbz++uuHPKe5uTna2to6t9bW1m4PDQAADBxl3WmZM2dO539PmTIlzj777Bg/fnz84Ac/iCuuuOKg5xQKhSgUCt2bEgAAGLC69ZHHJ5xwQrzvfe+LZ599tqfmAQAA6KJb0bJ379547rnn4uSTT+6peQAAALooK1r+6q/+KjZs2BAvvPBCPProo3HRRRfFoEGD4nOf+1y15gMAAAa4st7T8vLLL8fnPve5eO2112LUqFHx0Y9+NDZt2hSjRo2q1nwAAMAAV1a03HvvvdWaAwAA4KC69Z4WAACAahMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACC1bkXL0qVLo6amJq6++uoeGgcAAKCriqPl8ccfj9tuuy2mTJnSk/MAAAB0UVG07N27Ny699NK4/fbb48QTT+zpmQAAADpVFC1NTU0xd+7cmDlz5hGPLRaL0d7e3mUDAAA4WrXlnnDvvffGE088EY8//vhRHd/S0hJLliwpezDIasK1D1Z03gtL5/bwJEA5Kn3uVmogPOd7+3sKDFxl3WlpbW2Nr371q/G9730vhg4delTnNDc3R1tbW+fW2tpa0aAAAMDAVNadli1btsTu3bvjQx/6UOe+AwcOxMaNG+Of//mfo1gsxqBBg7qcUygUolAo9My0AADAgFNWtHziE5+IJ598ssu+yy+/PCZOnBh/8zd/845gAQAA6K6yoqWuri7OOOOMLvve9a53xciRI9+xHwAAoCd06x+XBAAAqLayPz3s961fv74HxgAAADg4d1oAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASK2saFmxYkVMmTIlhg8fHsOHD48ZM2bET37yk2rNBgAAUF60jB07NpYuXRpbtmyJzZs3x8c//vGYP39+/Od//me15gMAAAa42nIOnjdvXpfHN9xwQ6xYsSI2bdoUkydP7tHBAAAAIsqMlrc7cOBArF69Ovbt2xczZsw45HHFYjGKxWLn4/b29kovCQAADEBlR8uTTz4ZM2bMiN/97nfx7ne/O9asWROTJk065PEtLS2xZMmSbg0J5Dfh2gcrOu+FpXN7eJLq6e2vsdLr9bb+MudAYC0Orre/L/3p9zXoL8r+9LDTTz89tm7dGv/+7/8eV155ZSxYsCD+67/+65DHNzc3R1tbW+fW2trarYEBAICBpew7LUOGDInTTjstIiKmTZsWjz/+eNx0001x2223HfT4QqEQhUKhe1MCAAADVrf/nZaOjo4u71kBAADoSWXdaWlubo45c+ZEY2NjvP7667Fq1apYv359PPTQQ9WaDwAAGODKipbdu3fHn/7pn8aOHTuivr4+pkyZEg899FB88pOfrNZ8AADAAFdWtNxxxx3VmgMAAOCguv2eFgAAgGoSLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgtbKipaWlJaZPnx51dXUxevTouPDCC2P79u3Vmg0AAKC8aNmwYUM0NTXFpk2b4qc//Wns378/Zs2aFfv27avWfAAAwABXW87Ba9eu7fL47rvvjtGjR8eWLVvi3HPP7dHBAAAAIsqMlt/X1tYWEREjRow45DHFYjGKxWLn4/b29u5cEgAAGGAqjpaOjo64+uqr45xzzokzzjjjkMe1tLTEkiVLKr1MWhOufbCi815YOreHJ6HStTjW+b5wrOovv7b7y5z0PH9G6Hm9/T3tL2vYX+bsCRV/elhTU1Ns27Yt7r333sMe19zcHG1tbZ1ba2trpZcEAAAGoIrutCxcuDB+/OMfx8aNG2Ps2LGHPbZQKEShUKhoOAAAgLKipVQqxVe+8pVYs2ZNrF+/Pk455ZRqzQUAABARZUZLU1NTrFq1Ku6///6oq6uLnTt3RkREfX19DBs2rCoDAgAAA1tZ72lZsWJFtLW1xcc+9rE4+eSTO7fvf//71ZoPAAAY4Mp+eRgAAEBvqvjTwwAAAHqDaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqYkWAAAgNdECAACkJloAAIDURAsAAJCaaAEAAFITLQAAQGqiBQAASE20AAAAqZUdLRs3box58+ZFQ0ND1NTUxH333VeFsQAAAP5X2dGyb9++mDp1atxyyy3VmAcAAKCL2nJPmDNnTsyZM6caswAAALxD2dFSrmKxGMVisfNxe3t7tS8JAAAcQ6oeLS0tLbFkyZJqX6ZiE659sK9HqKre/vpeWDq3ovOO9XUgD7/WoP861p+/ffH1VXrNSv9/D5Wq+qeHNTc3R1tbW+fW2tpa7UsCAADHkKrfaSkUClEoFKp9GQAA4Bjl32kBAABSK/tOy969e+PZZ5/tfPz888/H1q1bY8SIEdHY2NijwwEAAJQdLZs3b47zzz+/8/GiRYsiImLBggVx991399hgAAAAERVEy8c+9rEolUrVmAUAAOAdvKcFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNRECwAAkJpoAQAAUhMtAABAaqIFAABITbQAAACpiRYAACA10QIAAKQmWgAAgNQqipZbbrklJkyYEEOHDo2zzz47fvGLX/T0XAAAABFRQbR8//vfj0WLFsX1118fTzzxREydOjUuuOCC2L17dzXmAwAABriyo+XGG2+ML33pS3H55ZfHpEmT4tZbb43jjz8+7rzzzmrMBwAADHC15Rz85ptvxpYtW6K5ublz33HHHRczZ86Mxx577KDnFIvFKBaLnY/b2toiIqK9vb2SeXtcR/GNXr1eb3/d/eXr6+05+8KxvvaVGgjfF88L4FjTX37v7u3ff4/170s1vDVLqVQ6/IGlMrzyyiuliCg9+uijXfYvXry49OEPf/ig51x//fWliLDZbDabzWaz2Wy2g26tra2H7ZCy7rRUorm5ORYtWtT5uKOjI/7nf/4nRo4cGTU1NdW+/DGrvb09xo0bF62trTF8+PC+HocKWMNjg3Xs/6xh/2cN+z9r2P9VuoalUilef/31aGhoOOxxZUXLSSedFIMGDYpdu3Z12b9r164YM2bMQc8pFApRKBS67DvhhBPKuSyHMXz4cE/ufs4aHhusY/9nDfs/a9j/WcP+r5I1rK+vP+IxZb0Rf8iQITFt2rRYt25d576Ojo5Yt25dzJgxo6zhAAAAjkbZLw9btGhRLFiwIM4666z48Ic/HMuXL499+/bF5ZdfXo35AACAAa7saLnkkkviv//7v+NrX/ta7Ny5Mz7wgQ/E2rVr4z3veU815uMQCoVCXH/99e946R39hzU8NljH/s8a9n/WsP+zhv1ftdewpnTEzxcDAADoO2X/45IAAAC9SbQAAACpiRYAACA10QIAAKQmWgAAgNRES3IbN26MefPmRUNDQ9TU1MR99913yGP/4i/+ImpqamL58uW9Nh9HdqQ1vOyyy6KmpqbLNnv27L4ZloM6mufhU089FX/yJ38S9fX18a53vSumT58eL730Uu8Py0EdaQ1//zn41vbNb36zbwbmHY60hnv37o2FCxfG2LFjY9iwYTFp0qS49dZb+2ZYDupIa7hr16647LLLoqGhIY4//viYPXt2PPPMM30zLAfV0tIS06dPj7q6uhg9enRceOGFsX379i7H/O53v4umpqYYOXJkvPvd746LL744du3a1e1ri5bk9u3bF1OnTo1bbrnlsMetWbMmNm3aFA0NDb00GUfraNZw9uzZsWPHjs7tnnvu6cUJOZIjreFzzz0XH/3oR2PixImxfv36+NWvfhXXXXddDB06tJcn5VCOtIZvf/7t2LEj7rzzzqipqYmLL764lyflUI60hosWLYq1a9fGd7/73Xjqqafi6quvjoULF8YDDzzQy5NyKIdbw1KpFBdeeGH8+te/jvvvvz9++ctfxvjx42PmzJmxb9++PpiWg9mwYUM0NTXFpk2b4qc//Wns378/Zs2a1WWNrrnmmvjRj34Uq1evjg0bNsSrr74an/70p7t/8RL9RkSU1qxZ8479L7/8cukP/uAPStu2bSuNHz++9E//9E+9PhtH52BruGDBgtL8+fP7ZB7Kd7A1vOSSS0qf//zn+2Ygynao30vfbv78+aWPf/zjvTMQZTvYGk6ePLn0j//4j132fehDHyr93d/9XS9OxtH6/TXcvn17KSJK27Zt69x34MCB0qhRo0q33357H0zI0di9e3cpIkobNmwolUql0p49e0qDBw8urV69uvOYp556qhQRpccee6xb13KnpZ/r6OiIL3zhC7F48eKYPHlyX49DhdavXx+jR4+O008/Pa688sp47bXX+nokjlJHR0c8+OCD8b73vS8uuOCCGD16dJx99tmHfSknue3atSsefPDBuOKKK/p6FMrwkY98JB544IF45ZVXolQqxSOPPBJPP/10zJo1q69H4ygUi8WIiC53qI877rgoFArx85//vK/G4gja2toiImLEiBEREbFly5bYv39/zJw5s/OYiRMnRmNjYzz22GPdupZo6eeWLVsWtbW1cdVVV/X1KFRo9uzZ8Z3vfCfWrVsXy5Ytiw0bNsScOXPiwIEDfT0aR2H37t2xd+/eWLp0acyePTv+9V//NS666KL49Kc/HRs2bOjr8ajAypUro66urmdezkCvufnmm2PSpEkxduzYGDJkSMyePTtuueWWOPfcc/t6NI7CW3+wbW5ujt/85jfx5ptvxrJly+Lll1+OHTt29PV4HERHR0dcffXVcc4558QZZ5wRERE7d+6MIUOGxAknnNDl2Pe85z2xc+fObl2vtltn06e2bNkSN910UzzxxBNRU1PT1+NQoc9+9rOd/33mmWfGlClT4tRTT43169fHJz7xiT6cjKPR0dERERHz58+Pa665JiIiPvCBD8Sjjz4at956a5x33nl9OR4VuPPOO+PSSy/1nqR+5uabb45NmzbFAw88EOPHj4+NGzdGU1NTNDQ0dPlbX3IaPHhw/PCHP4wrrrgiRowYEYMGDYqZM2fGnDlzolQq9fV4HERTU1Ns27at1+6EudPSj/3sZz+L3bt3R2NjY9TW1kZtbW28+OKL8Zd/+ZcxYcKEvh6PCr33ve+Nk046KZ599tm+HoWjcNJJJ0VtbW1MmjSpy/73v//9Pj2sH/rZz34W27dvjz/7sz/r61Eow29/+9v427/927jxxhtj3rx5MWXKlFi4cGFccskl8a1vfauvx+MoTZs2LbZu3Rp79uyJHTt2xNq1a+O1116L9773vX09Gr9n4cKF8eMf/zgeeeSRGDt2bOf+MWPGxJtvvhl79uzpcvyuXbtizJgx3bqmaOnHvvCFL8SvfvWr2Lp1a+fW0NAQixcvjoceeqivx6NCL7/8crz22mtx8skn9/UoHIUhQ4bE9OnT3/GRj08//XSMHz++j6aiUnfccUdMmzYtpk6d2tejUIb9+/fH/v3747jjuv6xZtCgQZ13Q+k/6uvrY9SoUfHMM8/E5s2bY/78+X09Ev+nVCrFwoULY82aNfFv//Zvccopp3T58WnTpsXgwYNj3bp1nfu2b98eL730UsyYMaNb1/bysOT27t3b5W/cn3/++di6dWuMGDEiGhsbY+TIkV2OHzx4cIwZMyZOP/303h6VQzjcGo4YMSKWLFkSF198cYwZMyaee+65+Ou//us47bTT4oILLujDqXm7Iz0PFy9eHJdcckmce+65cf7558fatWvjRz/6Uaxfv77vhqaLI61hRER7e3usXr06vv3tb/fVmBzGkdbwvPPOi8WLF8ewYcNi/PjxsWHDhvjOd74TN954Yx9OzdsdaQ1Xr14do0aNisbGxnjyySfjq1/9alx44YU+TCGRpqamWLVqVdx///1RV1fX+T6V+vr6GDZsWNTX18cVV1wRixYtihEjRsTw4cPjK1/5SsyYMSP+6I/+qHsX79Znj1F1jzzySCki3rEtWLDgoMf7yON8DreGb7zxRmnWrFmlUaNGlQYPHlwaP3586Utf+lJp586dfT02b3M0z8M77rijdNppp5WGDh1amjp1aum+++7ru4F5h6NZw9tuu600bNiw0p49e/puUA7pSGu4Y8eO0mWXXVZqaGgoDR06tHT66aeXvv3tb5c6Ojr6dnA6HWkNb7rpptLYsWNLgwcPLjU2Npb+/u//vlQsFvt2aLo42PpFROmuu+7qPOa3v/1t6ctf/nLpxBNPLB1//PGliy66qLRjx45uX7vm/wYAAABIyXtaAACA1EQLAACQmmgBAABSEy0AAEBqogUAAEhNtAAAAKmJFgAAIDXRAgAApCZaAACA1EQLAACQmmgBAABS+39HC3cpoYfd/AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Number of iterations\n",
    "N = 100\n",
    "# Sample Size\n",
    "sN = 49\n",
    "\n",
    "# Execute. Keep in mind it may take a while for the function to complete since it's going through a large number of iterations.\n",
    "tl_mean_list = get_probability(int(time.time()), 0.02, 0.02, N, sN)\n",
    "\n",
    "var = statistics.variance(tl_mean_list)\n",
    "\n",
    "print('mean =',sum(tl_mean_list)/sN)\n",
    "print('variance =',var)\n",
    "print('sd =',np.sqrt(var))\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (10,7))\n",
    "ax.hist(tl_mean_list, bins = 50)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[467478574, 512932792, 539453717, 20349702, 615542081]\n"
     ]
    }
   ],
   "source": [
    "import collections\n",
    "\n",
    "_ten2nine = 10**9\n",
    "\n",
    "class Subtractive_generator():\n",
    "    \n",
    "    def __init__(self, seed=292929):\n",
    "        self.r = collections.deque(maxlen=55)\n",
    "        s = collections.deque(maxlen=55)\n",
    "        s.extend([seed, 1])\n",
    "        s.extend((s[n-2] - s[n-1]) % _ten2nine for n in range(2, 55))\n",
    "        self.r.extend(s[(34 * (n+1)) % 55] for n in range(55)) \n",
    "        for n in range(219 - 54):\n",
    "            self()\n",
    "     \n",
    "    def __call__(self):\n",
    "        r = self.r\n",
    "        r.append((r[0] - r[31]) % _ten2nine)\n",
    "        return r[54]\n",
    "     \n",
    "if __name__ == '__main__':\n",
    "    srand = Subtractive_generator()\n",
    "    print([srand() for i in range(5)])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}