foxqstm/StochasticFactorizatoin.ipynb

## StochasticFactorizatoin.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "6bad5e5e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count= 9958\n",
      "99.58 %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOVElEQVR4nO3dXbBddXnH8e+vifiugDk6mJcmzKS0uRDFU4S+Uqma0E4zneEiWEUZmQxT07HtRQnjTDsdb7S2HYcBSTM0VfsCtspoSqNph744HQskTBEIEDyCJcfQJpQO7egFRp5e7BXcbk7OXjnseDj/fD8ze85a//Xfaz/PTvJjnbX22qSqkCS14UcWuwBJ0uQY6pLUEENdkhpiqEtSQwx1SWrI8sV64RUrVtTatWsX6+UlaUm65557nqyqqRNtX7RQX7t2Lfv371+sl5ekJSnJf8y33dMvktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSFjQz3JriRHkjxwgu1Jcn2SmST3Jblg8mVKkvroc6T+KWDjPNs3Aeu7x1bgphdeliRpIcaGelV9BXhqnimbgc/UwJ3AmUnOmVSBkqT+JnFH6Urg0ND6bDf2xOjEJFsZHM2zZs2aBb/g2u1/9wPr3/zoLz1v7GTHT3buC61jMWr2PXrxvR74Ho2bC+29R8f7OhUmcaE0c4zN+b9TqqqdVTVdVdNTUyf86gJJ0gJNItRngdVD66uAwxPYryTpJE0i1HcDV3afgrkIeLqqnnfqRZJ06o09p57kFuASYEWSWeD3gJcAVNUOYA9wGTADfAe46lQVK0ma39hQr6orxmwv4IMTq0iStGDeUSpJDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqSK9QT7IxycEkM0m2z7H9tUn+NsnXkhxIctXkS5UkjTM21JMsA24ENgEbgCuSbBiZ9kHgwao6H7gE+KMkZ0y4VknSGH2O1C8EZqrq0ap6BrgV2Dwyp4BXJwnwKuAp4NhEK5UkjdUn1FcCh4bWZ7uxYTcAPwEcBu4HPlRVz47uKMnWJPuT7D969OgCS5YknUifUM8cYzWy/i7gXuCNwJuBG5K85nlPqtpZVdNVNT01NXWSpUqSxukT6rPA6qH1VQyOyIddBdxWAzPAY8CPT6ZESVJffUJ9H7A+ybru4ucWYPfInMeBSwGSvAE4D3h0koVKksZbPm5CVR1Lsg3YCywDdlXVgSTXdNt3AB8BPpXkfgana66tqidPYd2SpDmMDXWAqtoD7BkZ2zG0fBh452RLkySdLO8olaSGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDWkV6gn2ZjkYJKZJNtPMOeSJPcmOZDkXyZbpiSpj+XjJiRZBtwIvAOYBfYl2V1VDw7NORP4JLCxqh5P8vpTVK8kaR59jtQvBGaq6tGqega4Fdg8MufdwG1V9ThAVR2ZbJmSpD76hPpK4NDQ+mw3NuzHgLOS/HOSe5JcOakCJUn9jT39AmSOsZpjP28FLgVeDvxbkjur6pEf2FGyFdgKsGbNmpOvVpI0rz5H6rPA6qH1VcDhOeZ8uaq+XVVPAl8Bzh/dUVXtrKrpqpqemppaaM2SpBPoE+r7gPVJ1iU5A9gC7B6Z80XgZ5MsT/IK4G3AQ5MtVZI0ztjTL1V1LMk2YC+wDNhVVQeSXNNt31FVDyX5MnAf8Cxwc1U9cCoLlyQ9X59z6lTVHmDPyNiOkfWPAx+fXGmSpJPlHaWS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktSQXqGeZGOSg0lmkmyfZ95PJvlekssnV6Ikqa+xoZ5kGXAjsAnYAFyRZMMJ5n0M2DvpIiVJ/fQ5Ur8QmKmqR6vqGeBWYPMc834D+DxwZIL1SZJOQp9QXwkcGlqf7caek2Ql8KvAjvl2lGRrkv1J9h89evRka5UkjdEn1DPHWI2sfwK4tqq+N9+OqmpnVU1X1fTU1FTPEiVJfS3vMWcWWD20vgo4PDJnGrg1CcAK4LIkx6rqC5MoUpLUT59Q3wesT7IO+BawBXj38ISqWnd8OcmngNsNdEn64Rsb6lV1LMk2Bp9qWQbsqqoDSa7pts97Hl2S9MPT50idqtoD7BkZmzPMq+r9L7wsSdJCeEepJDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSG9Qj3JxiQHk8wk2T7H9l9Lcl/3+GqS8ydfqiRpnLGhnmQZcCOwCdgAXJFkw8i0x4Cfr6o3AR8Bdk66UEnSeH2O1C8EZqrq0ap6BrgV2Dw8oaq+WlX/063eCayabJmSpD76hPpK4NDQ+mw3diIfAL4014YkW5PsT7L/6NGj/auUJPXSJ9Qzx1jNOTH5BQahfu1c26tqZ1VNV9X01NRU/yolSb0s7zFnFlg9tL4KODw6KcmbgJuBTVX135MpT5J0Mvocqe8D1idZl+QMYAuwe3hCkjXAbcB7q+qRyZcpSepj7JF6VR1Lsg3YCywDdlXVgSTXdNt3AL8LvA74ZBKAY1U1ferKliTNpc/pF6pqD7BnZGzH0PLVwNWTLU2SdLK8o1SSGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhrSK9STbExyMMlMku1zbE+S67vt9yW5YPKlSpLGGRvqSZYBNwKbgA3AFUk2jEzbBKzvHluBmyZcpySphz5H6hcCM1X1aFU9A9wKbB6Zsxn4TA3cCZyZ5JwJ1ypJGiNVNf+E5HJgY1Vd3a2/F3hbVW0bmnM78NGq+tdu/Q7g2qraP7KvrQyO5AHOAw4usO4VwJMLfO5SZc+nB3s+PbyQnn+0qqZOtHF5jx1kjrHR/xL0mUNV7QR29njN+QtK9lfV9Avdz1Jiz6cHez49nMqe+5x+mQVWD62vAg4vYI4k6RTrE+r7gPVJ1iU5A9gC7B6Zsxu4svsUzEXA01X1xIRrlSSNMfb0S1UdS7IN2AssA3ZV1YEk13TbdwB7gMuAGeA7wFWnrmRgAqdwliB7Pj3Y8+nhlPU89kKpJGnp8I5SSWqIoS5JDVlyoT7uKwuWiiSrk/xTkoeSHEjyoW787CT/kOTr3c+zhp5zXdf3wSTvGhp/a5L7u23XJ5nrI6YvGkmWJfn37v6G5ntOcmaSzyV5uPvzvvg06Pm3ur/XDyS5JcnLWus5ya4kR5I8MDQ2sR6TvDTJZ7vxu5Ks7VVYVS2ZB4MLtd8AzgXOAL4GbFjsuhbYyznABd3yq4FHGHwNwx8A27vx7cDHuuUNXb8vBdZ178OybtvdwMUM7hf4ErBpsfsb0/tvA38F3N6tN90z8Gng6m75DODMlnsGVgKPAS/v1v8aeH9rPQM/B1wAPDA0NrEegV8HdnTLW4DP9qprsd+Yk3wTLwb2Dq1fB1y32HVNqLcvAu9gcJftOd3YOcDBuXpl8Gmki7s5Dw+NXwH8yWL3M0+fq4A7gLfz/VBvtmfgNV3AZWS85Z5XAoeAsxl8wu524J0t9gysHQn1ifV4fE63vJzBHagZV9NSO/1y/C/LcbPd2JLW/Vr1FuAu4A3Vfca/+/n6btqJel/ZLY+Ov1h9Avgd4NmhsZZ7Phc4CvxZd8rp5iSvpOGeq+pbwB8CjwNPMLhv5e9puOchk+zxuedU1THgaeB14wpYaqHe6+sIlpIkrwI+D/xmVf3vfFPnGKt5xl90kvwycKSq7un7lDnGllTPDI6wLgBuqqq3AN9m8Gv5iSz5nrvzyJsZnGZ4I/DKJO+Z7ylzjC2pnntYSI8L6n+phXpTX0eQ5CUMAv0vq+q2bvi/0n3DZffzSDd+ot5nu+XR8RejnwZ+Jck3GXzb59uT/AVt9zwLzFbVXd365xiEfMs9/yLwWFUdrarvArcBP0XbPR83yR6fe06S5cBrgafGFbDUQr3PVxYsCd0V7j8FHqqqPx7atBt4X7f8Pgbn2o+Pb+muiK9j8N31d3e/4v1fkou6fV459JwXlaq6rqpWVdVaBn92/1hV76Htnv8TOJTkvG7oUuBBGu6ZwWmXi5K8oqv1UuAh2u75uEn2OLyvyxn8exn/m8piX2hYwIWJyxh8UuQbwIcXu54X0MfPMPhV6j7g3u5xGYNzZncAX+9+nj30nA93fR9k6FMAwDTwQLftBnpcTFnsB3AJ379Q2nTPwJuB/d2f9ReAs06Dnn8feLir988ZfOqjqZ6BWxhcM/gug6PqD0yyR+BlwN8w+PqVu4Fz+9Tl1wRIUkOW2ukXSdI8DHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUkP8HUg9mCGiVVBQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD6CAYAAABApefCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzk0lEQVR4nO3dd5hU5dn48e9NlyIdRNqCKIqKoCsg2BBBihHzan4vtqjRFzWaaExiQJO89hDzRo2RqCS22FADVkAURBGRsvTe27LAUneXBbY+vz9mZnfKmZkzM2f6/bmuvXbmzCnP2XLPM/fTxBiDUkqpzFYn2QVQSikVfxrslVIqC2iwV0qpLKDBXimlsoAGe6WUygIa7JVSKguEDfYi0llE5ojIOhFZIyL3u7e3EpGvRGST+3vLIMcPF5ENIrJZRMY5fQNKKaXCk3D97EWkA9DBGLNURJoBS4BrgduAQ8aYCe4g3tIY8zu/Y+sCG4GhQD6wGLjBGLM21DXbtGljcnJyorohpZTKRkuWLDlgjGkb7PV64U5gjNkD7HE/LhGRdUBHYDRwuXu3N4FvgN/5Hd4P2GyM2QogIpPdx4UM9jk5OeTl5YUrmlJKKTcR2RHq9Yhy9iKSA/QFFgLt3W8EnjeEdhaHdAR2eT3Pd29TSimVQLaDvYg0BaYADxhjiu0eZrHNMm8kImNFJE9E8vbv32+3WEoppWywFexFpD6uQP+OMWaqe/M+dz7fk9cvtDg0H+js9bwTUGB1DWPMJGNMrjEmt23boGknpZRSUbDTG0eAV4F1xphnvV76FLjV/fhW4BOLwxcDp4tINxFpAIxxH6eUUiqB7NTsBwG3AFeIyHL310hgAjBURDbh6m0zAUBEThWR6QDGmErgPmAmsA74wBizJg73oZRSKgQ7vXHmYZ17BxhisX8BMNLr+XRgerQFVEopFTsdQauUUllAg71SSjmkvLKaD/J2kYqLQoVN4yillLJn4pzN/G32JhrWq8PoPqk1pEhr9kop5ZADR8sAKD5RmeSSBNJgr5RSSXS4tJwzHpnBwq0H43odDfZKKWWhsqqabzfGfzT/sl2HKa+q5uVvt8T1OhrslVLKwotzNnPra4uYm4CAnwga7JVSysKOg8eA2jx8OEXHKtiwtyTkPvmHj/HJ8t0xly0aGuyVUipCby/YQcGR4z7bfvLKfPJ2HA553HUvzef+ycvjWLLgNNgrpbLCzoPH2Hag1Gfb+KmrmLPeag7H4A4cLeP3H6/mttcX+WzfuO9ozeNgUw7sK7b3KSEeNNgrpbLCpX+Zw+D/+8Zn23uLdnL7G4sjOk91tWvA1OFjFU4VLSE02CulVBbQYK+UUhGwMxHC9FV74l6OSGmwV0pllS/X7GW7X+7e39fr97Hb3QAbbJqbYHl5gPlb4jtAKhoa7JVSGeHA0TJ+8+EKTlRUhdxv7FtLuNwvd+/vZ2/ksWjbIQdLl3wa7JVSaeGDxbsCetN4e3r6Ov6zJJ/PVgSufLpxX+j+75FIwQktbdFgr5RKCw9NWcnVL3wXdj/XSqq1CktOMOy5ufEqVtrQYK+USgmFJSfYU3Q85D6l5dYpmvsnL2PqUuuRqSUOz0ApoZL1IRz3KrtnvvsdB0uZuXqfE8UKK+x89iLyGnA1UGiMOce97X2gp3uXFsARY0wfi2O3AyVAFVBpjMl1pNRKqYzT76nZAGyfMCriYz9ZHpi6iZdo0zhj38oL2Db02bmUV1UDgZ9InGanZv8GMNx7gzHmv40xfdwBfgowNcTxg937aqBXSsXskY9WpcRKUJHG5u82HQjY5gn0iWBnwfG5IpJj9Zq43or+H3CFw+VSSik2F5bww9ZD7PBqmH1n4U6eGH1O1OmUbBXrsoSXAPuMMZuCvG6AL0XEAK8YYybFeD2lVBa58lnrhtVkBnpja1hV6ok12N8AvBfi9UHGmAIRaQd8JSLrjTGWvz0RGQuMBejSpUuMxVJKJYMxhvKqahrWq5vsosSdhBxWFbl4p6ai7o0jIvWA/wLeD7aPMabA/b0Q+AjoF2LfScaYXGNMbtu2baMtllIqSXYfOU638dPp+fsvKC2z3wNmwNOzeW3etjiWTEFsXS+vBNYbY/KtXhSRJiLSzPMYGAasjuF6SqkU9ub87TWPI+nuuLf4BI9/vjYOJVLewgZ7EXkP+AHoKSL5InKH+6Ux+KVwRORUEZnuftoemCciK4BFwDRjzBfOFV0plQwrdh2hMoG9SJQz7PTGuSHI9tssthUAI92PtwLnxVg+pVQKWb27iNETv+fnl5/GQ8PPTHZx0pYxiW9k1hG0SinbCktOALBuT3GSS2JfrO2e/oenQBf/qGiwV0pllERVmNOtn78Ge6VUWDf+cwGjJ35f89wAd76Zx8Q5mwE4UVHFR8us56ZJtDSteMedBnulVFjztxxkxa4jPn3LZ63bx19mbgBgwoz17C9J3mLaoaRLDTwV5sZRSqmQPLn8REnXvHkyabBXSqksoMFeKaWygAZ7pVTEUiGNcrSskookDO5KgVuPigZ7pTLUB3m7uPPNwAUzYpJCjZ3n/O9M7nD6/iKQQj8KW2Kd9VIplaIe+s/KZBch7uZu3J8W5/SXjE8HWrNXKsOVVVqv26qsvZ+3K+TrqbBKVjQ02CuV4dI0Njkm1vsvPl6RlLYBp2mwV0qpEB7/fC0PTF5e89zu4Kfj5an1iUqDvVIqYomYHOy7TfHPnVuxCuXTVu2peWw3jfOPbzZHdN2UXalKKZVYu48cj2pKAidjSCJ7oNzy6qKojkvcRGihr1RemVqpH+2No1SaGDThawC2TxiV5JIESpf5Z1KZzo2jlIqLqmrDS99sSbnc8pz1hawpKAq5T6gPK/6vBYuhD76/3HJ7sKC7r9g1/4/tT0op9gaowV6pLPXRst38+Yv1PD9rY7KL4uP2NxYz6oV5cb/O1AinZO7/9GyqqtO3a5OdNWhfE5FCEVntte1REdktIsvdXyODHDtcRDaIyGYRGedkwZVS9pgg9eDjFa4a/dEy+4uDZ7vqCBpAJETVfuv+o+wtSuxMoXZq9m8Awy22P2eM6eP+mu7/oojUBSYCI4BewA0i0iuWwiqVjZbuPEx1itUoi46VJ7sIjpq/+UDN41DZl8oqQ8mJwDdHqzfMUCn4oc/NZcCfZkdSxJiFDfbGmLnAoSjO3Q/YbIzZaowpByYDo6M4j1JZa8HWg/zXP+bzytytzp88im46nnz2ivzQOfV0c9vri23td887Sxj5wncB22/85wKni+S4WHL294nISneap6XF6x0B73HH+e5tSimbCo4cB2DjvpKozxEqnQCR9aSJtC/495sPxD1NtGFv6J+Nk11Pv9lQ2/d/95HjNYu2rEyDN79og/1LwGlAH2AP8FeLfaz+hIL+2EVkrIjkiUje/v3JGUyhVLxt2lfC2X/8IqHD74Pl7B29hsUlCotPcNO/FvLA5GVxvfZVz8+N/SRR9py5+V8LnT5l3EQV7I0x+4wxVcaYauCfuFI2/vKBzl7POwEFIc45yRiTa4zJbdu2bTTFUirlDX1uLqXlVbwwe1OyixKVSPqCexqAN+47Gq/iJN2Og8eCvpZqYw+iCvYi0sHr6Y+B1Ra7LQZOF5FuItIAGAN8Gs31lMo0hyNs4IzHUHqnzlh0rIKZa/aGuE5qNS5b8cTlJz9fy6crgtZJI3KiIrVG0Nrpevke8APQU0TyReQO4BkRWSUiK4HBwK/c+54qItMBjDGVwH3ATGAd8IExZk2c7kOptGI3dieidhgupx/Ove8uxaqzkOe86TTr5r/mbXPsXK86eC4nhJ0uwRhzg8XmV4PsWwCM9Ho+HQjolqlUplm68zBV1YYLc1o5et6iYxUAVFRFHzHDBdtIat5Wbws7D1mnMlItjRFKOpU1WjqCVikH/Nc/5vOTl39gc6Gz+em/zNwA+M66mGrCvVl43myqqw1b9sc/f5+ucTve5dZgr1QYS3Yctj3a8cpnv2XO+sKw+9mtS1ckYDBVrGkcy3NanPLFOZsZ8tdvw3aVtOOpaetiPoe3aH8GZSk2s2UoGuyVCuO6l+Yz5K/f2N4/lj7xmSxvx2EA9hQdj/lcr31vPx8e7xTN3W8tCfn66Bfn0fvRmfEthA06xbHKWkeOlSMiND+pfth9Sx2eGdJ2A62jV40Pq3sxJjDIJmrtVv+rOPFJIpQvQvREgtQZbaw1e5V1yiqreHfhTvo8/hXnPfZlsosTF0vctehQoom9TtSS4z1v+4GjZez06v/uafcIJRUaaOP9Vqg1e5U1jDE8+9VG8g8f56MIp7dNJ8fKK7nupfm29w8W6NbvLea0tk2pX9eZOqHnzSXecXXCjPVMmLGevN9fyca9JWw7UBrnK6YHDfYqa6wpKObvXweuC7p6dxGntjiJFbuOMPjMdkkombMqKmOvI+44WMrw57/jZ4O68ccfxTZZrSd94+m1k6hadO6Ts2zve6y8ipxx0+JYmuTTYK+yRrC0xdV/n0efzi1YvusIqx+7iqYNI/u3OHDUd11YJ4OZk+cqOHKc7QdKGdijDRA6h37gqGuE77Jd4dNBwfina2pr9imQM0mQu97KS3YRamjOXilci0kAVEUxeOmq56KZiKv2OsYY/vHNZnYFGZwUK8+VrnpuLjdaTNwVa+gNP2jL73rZE+uZuWZfsotQQ4O9UjE6WOo7z02kDZ/7ist45osN3Pr6ooDXnOzAUhJkquF4NQx6YrrnHhKVs1fWNNgrlQTeQdyz1F2iF/6OpldMLCmYmpG2Gu2TQoO9SluHSst58etNtvtvRxrbjDGs31sccbkivo6D54roujZ+bv67RDSPjtQe0/vRmSzY6rvg3dSl+Qya8LXt86nYaAOtSlvjpqzky7X7OL9rSwae1sbx8786bxtPOjwsPxSruO5EY2a4oB6P9xMRwNSeudhr3VbPPT30n5VUptjausmkc+MoFcQxd9qjKk4BY9XuxIx8DBWMkz0XvP8nCyfefDznrJNNLbUpQIO9Sgv/+GYz5/6v7/wikQbCdJpX3UnBcvPx/HGE+ll7SqOxPrE0jaPSwjNfBB/yno79tlPpjSdUQ60T5fQ/h+d6GuwTS2v2KmtEGlyijUXRvvlYBd1w51pbUMzAP83mcGnwZQ6jmYAs2M/Kanuw84c7h6ZxfPkPznOaBnuVMsorq1m+60iyi5FwVrHymhfnMWVJfthjJ36zmYKiE8zbfMCRsuwpOs79k5dR5rV+amVV7eNoavolJ3z794vfd+VyKMJ1iSNlZw3a10SkUERWe237i4isF5GVIvKRiLQIcux291q1y0UkdcYNq5T01LS1XDvxe9urPaVSKsSbnbaEcPuszC/i1x+uiNv1g3ni87V8sryAWetqR34u3HYoxBHBc/+eYH68wnf8gKdCH+/ZL9PNrkOxz/Mfip2a/RvAcL9tXwHnGGN6AxuB8SGOH2yM6WOMyY2uiCqTFRaf4J63l1BaVsnqAlef9iMR1nCyNWbsKz7BtJWu5Qrj9b63fNcRtnrNGmn3Z731QOg37FX5RRwNMqJXxUfYYG+MmQsc8tv2pTHG85taAHSKQ9lUFvi/LzcwY/VePl9ZEPU5UrWGH0osZf7brE3kjJvGX7+0brT2z/OHu9Qb87fzl5nrfcrlfYY/fLw68CAv+0sCc803/nNhiD70wo9enBemVMppTuTsfwbMCPKaAb4UkSUiMtaBa6kMk46BOhwnewdZ1aQ/XLILgPUxrsDk/bOfOGdL2OuCdZ49WFCfsdp6Bads/SSWbDEFexF5BKgE3gmyyyBjzPnACOBeEbk0xLnGikieiOTt378/lmKpFFRaVsm6PcGnHvAOkJHG/3gFDydyypVV1UxZkk+1jYFfdi/nGUTmXb5wPW68p30oLD7BU9PW1szJ482zqSKK2T8Dz2V9jkoHzq0iF3U/exG5FbgaGGKC/FaNMQXu74Ui8hHQD7CcD9YYMwmYBJCbm6t/DRlkZf4RrnnxewA2PTXCZ+Uj7190pKE1HT4VTPpuK898sQEDXH+BdbYz1H1YveapSdeN4AfmnWrp9/TsoPvtKT4BuFI78XIwzl0M01XdOvH9yBNVzV5EhgO/A64xxlhOwi0iTUSkmecxMAwInfxTGenhj1bVPPavTdY8zdCP9gdKXI3NdhqdrWr2VtusavahPoUUHavgllcDp0+2ssLBrq/Pz9pkuf3+95c7do1MUlVtLNs/nGKn6+V7wA9ATxHJF5E7gBeBZsBX7m6VL7v3PVVEprsPbQ/ME5EVwCJgmjHmi7jchUope4tOhEzZWEl2rLeqXPuXadhz3/Leop1xu56VYxbTHnuCvXdF0PvDtX/XyzcdrKXviGCBlWC9bcorqy23K1hdEL/5mMKmcYwxN1hsfjXIvgXASPfjrcB5MZVOpZUT7v7UA/7kShNsnzAqYJ9IpszddqCUuRv3c+vAnJDXTdQbxcZ9Rxk/dRU39OsS87mMgQ/zdnFNn1MtXgv9VuDJ/1vl1RdsPRgQTP81b1sMJfU1fuqq8DuplKRz4yjHXPjkrKCrIYVjlYa47qX5HCot58b+1sHVyRkh4/WGEayMM9fsZcrSfLYeKKVXh5MjOqcnFWY12njMpAURl1FlB50uQTkm2kAfTPHxCkfPFw8j//adrf3838w8KY6DR8v4xXvLXPu433LCNTzrHPAqGhrsVXLZiFvBgp93d83Xv99GzrhpEY++jdVam20TdicjO1pWyTmPzgy5T7zm7w81mZpKjHimJDXYq6Sy6oxTM4ozzF++d4pk8iLXQKO97q6DViLuNh/lf1603fNF4J63l1g2ynqLV80+llHMKvVpsFcpQcRekKyuNvz5i/XsPuI1aZTN4JrMfvl2B2gt33kkvgUJ4Q+frEnatVX8aQOtSipbi1571eDXFBTz0jdbWLj1oM/grNrzOVo8x9i5z8oqE/WnifsnL+fXH0Q3S6bKDhrsVVLVpHHEYkUjBP+kvqcnSmW1oX5dr32T3VHfT9GxCr7ZWBh2P+92B59PK1HQhlsVigZ7lRKsJg8L1bUyWEU5VWr2d/57MYu3H+bU5o0AVxpn077YJi5TKhaas1cJFXIhahu1c8spBWzmPhJZ+1+8/TAABUWuBmNjDEOfs5wWSqmE0GCvkipk8A8TxC2nOHAwoKfjQuZKBaPBXqWEUEHa+w0hXJomEWmceZsiX+81WG+cssrQ3SxVdonnUo0a7FXchZzC13J/i+WSLDj5b1HstSj20bJKcsZN47MV1v3Ob351IRVV1REN4Jq61Hrx8I+Xa992lRjaQJvFlu08TI92TWnWqH7SyuDdJTGStIkJ0nzrxHw5+YddMzuWVwWfnbH/07M5FGTE6dKdhwO2rSmIbBZQlZ3sjrSOhtbss9Tx8ip+/I/53PXWkrhfy17Dq71AHyxPH8+Pv1aCBXqA6ausl+NTKpk02GepimpXrXVVfvzmz7biX/N2qh4TzxqRUomiOXsVN0kPkfbS80qpGGmwz1LpHlyNweedKlyNqLKqmtHudXDD0S6XKhNpsM9wc9YXkjNuWlzXtoxEJNkWuyHXf8bMExVVvL1gR82KTgCHj1XYmk7geEUVB0tT42ellJPsrEH7mogUishqr22tROQrEdnk/t4yyLHDRWSDiGwWkXFOFlzZ84Z7/dFga1smItcd7BKlZZVMW7UH8G3EDVWiYOX1DvjPzdrI7z9ezfTVeyIrKK4lFW/858KIj1PKCcmez/4NYLjftnHAbGPM6cBs93MfIlIXmAiMAHoBN4hIr5hKqxyT6N4rVp6ctrbmsSBB/9Kt4rt3V8byquqa7pJlldUcPOrqKVPqtXJWCtyuUkkVNtgbY+YCh/w2jwbedD9+E7jW4tB+wGZjzFZjTDkw2X2ciqOJczbzlFcQrZH0lthAnqAcTEWIfu7efvX+8ppBUde9NJ//LHENYBKEdxbuoN9Ts2IrqFIJEs9/02gHVbU3xuwBMMbsEZF2Fvt0BHZ5Pc8H+kd5PWXTX2ZuAOCRUa4PUalQow238hK4y+n3l+5JsQcbKOXZfiDIm8bK3Ud4e8FO174p+GanVCLFs4HWKswE/ZcTkbEikiciefv3749jsVSibTtQGnafaSv3UBVhRPbMLBmMJ9ADPPjB8ojOrVQyxLNuFm3Nfp+IdHDX6jsAVqs05AOdvZ53AoJOBGKMmQRMAsjNzdV6mA1V1YaqakODetG/Zyf6Bx3sep6GWoAxkxYw5sLOlvtFW97vopi8TKlMEm2U+BS41f34VuATi30WA6eLSDcRaQCMcR+nHHLaw9M54/czojo2Udmd+VuiC7KTF9dmADUFo1Ts7HS9fA/4AegpIvkicgcwARgqIpuAoe7niMipIjIdwBhTCdwHzATWAR8YY3RF4xQT70DqRDfGeZu1Vq5UrMKmcYwxNwR5aYjFvgXASK/n04HpUZdOBfhq7T5+/s4Slv9xWETHHSwt5+0FO7h5QFcgMQ23VpOF3f/eMs7v2pJep57Ml2v32TqP92RtX66xd4xS6Sie/5c6xXGa+euXG6ioMuw4eCyi437z4QoABnRvRY92zSKu0RtjMAbq1LH/17jMYqrf2esLmb0+/ELcweg6rkpFR6dLSAP7S8oY++889rrXM42Ef2gur4wub/PSt1vo/vB0Sk5U2D7mjjfzorpWKLG8USiV6hZt8x/S5Byt2aeBC92DguwOMvLmH9o9fdMj/bg4eZGrwfTg0fKkLnaiVCabva6QXw/rGZdza80+TTmV27O7slMqDM5SKtPF8/9Mg32G8//bqV3eNfRfVWlZpeUnCe0FqVR60mCf5bwbal//fhs546ZRUVXN2f87k5u8uk2Gq3AUHbOfy1dKJZ4Ge1Xjsc9cE6jtPnwcgEXbAxuLgk0xXBxBw61Sypp2vVSAa1riSLtMBpvKOFSuvrI6MH0TakrkzYVH2XU4sq6gSqlA8RzkqME+jXjXqu3WAKJZnOTZrzbWPM4ZN43NT42oPZ/F/lc++23E11BKBYpnsNc0jgqwaJvvYKgTldU1M1cO+eu3vDB7EznjppEzbloyiqeUioIG+wznn345XuE7t7xVReLAUd81WP0/HXjX/K1GySqloqNdL5Vj/jR9naPne+aLDY6eTykVHxrs04gT68YGW9UplKv/Pi/oaz9sPRhLcZRSCaLBPsE+XVFAzrhpHLaYEdIOuyNePY56LbptebyN00U66ZpSKjqaxskQH+Tt4pfvLQNgq42l+vx5586/2WBv6cZgEyvpgiBKZRcN9gn0yrdbYj6HJ0hPmLE+5nMppbKHBvsMUFVt2F9SFn5HLyf8euXM8FoDVil/53Q82ef5mz/rxwNXnh6wX8/2zWK6TqsmDZhx/yUB29s2axj22ME928Z07VSg/ewV4GqgtcrpPfD+ci58alZAl0kruw8fJ2fcNC54clbNtqJjFdzzzlIni6oSoPlJiZtq+s6Lu/s8v+yMtjxw5RkB+/n/fV7du0PQc57bsbnP8y6tGvPJvYOotoh4dW0ks7u2bhJ2n2wWdbAXkZ4istzrq1hEHvDb53IRKfLa548xlziNOdGbxuqd/7MVBQAcORa+0bfa7/jyqmqufE5HwKab5ifV5/XbL7R8bfNTI5j14GW8fPMFjl3v2r4dbe977+DTah6/eOP5/Ofuiyz3e+b63tSvW/s/8aPzOtC5VWPLv/FsmWI7JRtojTEbjDF9jDF9gAuAY8BHFrt+59nPGPN4tNdT4X20bHdUx0WaAlLJ4x0Mzu/S0nKfenXr0KNdU87rXFtzfuWW8IF/3IgzAXjsmrN5aHjwBTT+NqZPyPP8uG8nn+e5Oa1Y9MgQnv9v3+MiSVlY7duice0nm1sGdOWGfl0C9rGziuZ4932ngnRI4wwBthhjdjh0PmUh3Dw3E+fE3gCs4u+Ukxvx/3I7hd/RQjd3qsIT9C87I3ieumG9ugD07tScoWe1r9neoJ7vv/3sX1/GXZd1565Lu7N9wihuHZjDf3kF7Gm/vNhn/9F97NXyu7etTau0a9bIcv1iq3UVOrdsbOv8n9w7iJv6d2HL0yN54tpz6HlKYHvBfVcEtiv4a9wwO6YIcyrYjwHeC/LaRSKyQkRmiMjZDl1PqaT73fDoaoQLHh7CM9efF7D9w7svYlSIHDdQk8/2hMg3f9aP7RNG8eqtuQH7tmrSgMljB/DOnf19Aq1/rv+0tk0ZP+KsoGnGs09tbrkdYMb9l/Dn684NWeZINW9cn+0TRvHu//QPus8pJzeia+smPPXjc6kbpPp+60Vd6du5Re15E9jGkYpiDvYi0gC4BvjQ4uWlQFdjzHnA34GPQ5xnrIjkiUje/v32+pCnmgVbD5Izbhpfr99nY+/QtXRjTMAc8ZHk/I+XV0U9cEvZ43R+9cKcVky88fyQ+3j+avz/FoZ41dy9DejeumbN4L/f0Bdw9ZjZPmFUyOucVL+ujRLDWR1OpmXjBrb2DfXJtEPzRgFpmIGntal57Lndr351KbMevNSyx46/Ln4Ntt/9bjD//Gngm2IqNQekZM7eywhgqTEmIMIZY4qNMUfdj6cD9UWkjf9+7tcnGWNyjTG5bdumZxeqMZMWAPDOgp0xn2vinM30fvRLCktO+Gy3m9IbPXEefZ/4KuZyKJdBPVoHbBvWqzbAdmjeyOe1G/t3oVeH2u6K153vSotc2+dUy/M/PLL2U8KT155T8/jvN/Rlyj21DZwdW5wEwM0DukZSfABGntuBWy/qyjPX9w67b/PGgbXgS063/NeljleE+tuYvrbL07ZZQ+68pBsA8353BZ1CpG+auFMtTRvVo0e7ZrRsEvoN5q07+nH7wJyA7fXqWqSSvDZ5fr7+vH8/8ZTq89nfQJAUjoicAuwzxhgR6YfrzSXjJ1Nx4vc1bdVeAAqLo2s83bjvqAOlUB6N6gXWdLu3bVrzOKd1E/YUnQjYp3bf4N0C5/zmcnJa1wa6mwd05fcfrwZg1LkdqFNH6Nq6MTsOHuPFG8+npUUgBnj99gtpFaKWXbeO8Njoc4K+Hs5rt11ouS7x4DPbcdel3Rl7aXdaN23I5sISy+OvOvsURp3bgYdHnUXThvVoflJ9Hhp+Jg/ZSIe9ftuFTF+1hw7NrYOxv0tOd1cY/WL7Zae35baBOVyY04r73lsaEFzP69yc3UeO1zx/cOgZvPztlppPR+ksppq9iDQGhgJTvbbdLSJ3u59eD6wWkRXAC8AYE81qGhlCQjzz5/kx1cmWPmcpwrsXyue/qG2Y9KRN/FMNP7nAVWPv3ak5C8YP4dEf9QKgQd3af62xl3avqfl7p1+m3HMRX/7qUrq1aRKQljmrw8lccWa7mlz757+4mB/GX0GrJg3c4y0C/y4G92zHeV456nCG9mpPj3ZNw+/oVr9uHRo3CKwf1q0jjB95Fq2b+g588i9ho/p1mXjT+XRscVLE+fPOrRpz12Wnhd/R/7iWvm8OdeoIj15zdti2EY/rLujE2seHB20XcFrKLktojDkGtPbb9rLX4xeBF2O5RrbyNMS9u8i3g1MWv1fGVfe2Tdi6v5RLerRlQ58SzmjfjHO8Bv14atM92jXlo58PpIW7Bv3n63rTu3MLfnJBJxrVr8uN/buyp+gE913Rg0c/XcvaPcXc3L8rjRrUoVH9Otw+KKfmnBd0bRW0PP456WaN6jteu7TKX3u7b3APGje0l7v35ql933N5j6jK5aQe7YKP6LXzr5RJVa3s6HOUInzftYP/pVVVm5o0zNsO5P+Vr+Fnn8IXa/b6bPOkaUSs885NG7n+VYwx9PXq316njnCLV/68Qb06jB95FgBP/fgcxvTrTBd3imb9EyNIJ7+5Knhf+1CaNKwXtgE4Xl68sS+n+wX4f/+sH6/M3UJTi08l4NvV88GhPZm+aq/lfh4z7r+EU05ulHZtYjpdQpS+27SfS575OmCOmVDs1CTKK6s5Vl4Z9HUnRuFmknfuDN49L5ibBgQOvvGkP/zTC4seHsLiR6607A8eTqP6dbkwJ3jtXTnv6t6nBvS3v/SMtrxz5wDLfv4Ag3q0oWOLk7hlQNegaS3vIxvUqxO2gTgVac0+So99tpZdh46z89AxzvCb/CmWVMvACV8HneNmU2EJJSeCvxFkk/VPDKdhvTq23/y6tGrMzkOuefmt5ll59Jpe3NivC51b+fYIaXeyK9c+8LTWvPb9NvpEkBNXsTuvU3NW5BfF9RoCfD/uirheIxVosI+S8RvcYodvjLE+MtRkZrsOHQ/6WrZpZNEP/Md9O1JtDI/+6Gw+X7WHP7h7tAD8etgZHC4t59HP1tKldWO6t2lCm6YNWbTdNd9/w3p1ObdT8MFDV/Zqz6pHh2VEr4x0MnnsRQEL8MTbu//Tnxv/udBnW6I+UUfzCdIuTeNEqXZwS+znOlxaTmHJCXLGTYv9ZGmscYPaAL59wig+uMt6Ai3vni7e7rn8NP42pi8tmzTgJot5Um4dmMPqx66iU8vGfP2by/ng7otqetPYoYE+8U5qUNfW9Max8P8fHnhaG045uZHPa/W8UkBdW9mbziEaka5EFwkN9jGLNtrX/lL7PvEV/Z6a7Uxx0tgL7oZRT6qkX7dWAYOVfhh/BaseG+az7X/cA3PaeHX9s3oTFhGa+s2D8sz1vdny9MhYi67SmFWt3T/oev4mX7ihL/WCVDZSnaZxohXiDVg7R0bu619fVtM4+pthtb1Afhg/hHV7ihnxt+94645+loNqxo84i18OOd2n5m33Y7eIYDGoUmWBk+rX5XiYDhaetEpOmyZseXpkwvrbx0N6vkUlWFllVcBc8U6mcTLV67cFzrfu+Xg89ecDfbZ3b9uU1k0bsn3CKC72G5Z/VoeT2T5hVO2oSD916ohliuX9sQNqpvbt3y1wugOV3T77xaCaQXB2pHOgB63Z23LVc3PZfvCYT9/hqBpovfbeUlhK/uHjtqeLTTd9u7Tg0jPa8vjos1mxq4gpS/Pp3rYJsx+8jMKSMssGVqf17+4K8Mnq861SW492zYIOuurfrTWfriiwPSFcOtBgH8bMNXvZfvBY0NejbaV/aMpKgJCTP6WLHu2asrnQdy6eF8b0pW4d4acX5bChWwlTlubzj5vOR0Rof3Ijio5VBDmbUsn3zPW9ue+KHpYTwsXT6t3FcTu3pnHCWJl/xHJ7JHl5Ywx/mbme7QdLA1677qX50RUshXz+i4t5129wk3d/9Z6nuKbUPfOU2lkgvRvAZj14WfwLqVQEGtWvGzB+Jpyrzm7PXZd1D79jkmiwj5Gden3+4eNMnLOFssrAGQPTwZe/utTn+be/vZxvf3s579zZn21/Gkmj+nUZ2KMNc3872PY5PWvhtmhcP6LJuJRKFR/6ra37yi25jB9xVpJKE56mccIINsghm+Yj867heOe/u/otDtGldWOm/nwgB4+GXzTFM4rV02CrVLppWM+6rtygXh3KU7BipzV7C/tLynh7gb3ldK1S9sbAJ8t388OW9Ju6f/uEUXx876Ca554pff/109yQy8R5nN+lJUN7Wa+a5K154/o8/999+PfP+kVfWKWSyFPha39yQ591ej+77+IgRySX1uwtXPjULCD4yjwQfqTb/ZOXA7Dl6ZEBywumkok3ns+97y712dancwu2TxjFnqLjNTXvK20E8Ehd2zczeyKp7NL+5EY+6/RaLXwOMOWegUlto9Ng72exe64UgIoqE7Qfvedd3SrN8+3G2jV0T3t4uqPlc9qo3h3o0+UKluw4TDe/tIzdVYGUykYd3Quj/Nii0tK7U3NW+k3gdkHXlrRr1pDCkuhWn4uVBns/dhfprgn2aTDO4vHRZ/PHT9YAruD++1FnUUekpuwdW5wUdO1NpZS1Nk0bsvmpEZaDrT6972Kfua488/v0796az1YUJKyM3jTY+8mk+eL/c/dFlJyoZPCZ7RjdpyPNGtYLOqe3UipydubJWfqHoUEbcxMppmAvItuBEqAKqDTG5Pq9LsDfgJHAMeA2Y8xS//OkrvTqcvP+2AF8vb6Qy85oy9Rlu7mga8uaN69I1/xUSjmjlddCJ8msajlRsx9sjDkQ5LURwOnur/7AS+7vKWdNQRG9OpwcsD3YL6dmugT3Dr/9cEWcSmbPxidH0KBenZopAgb2CN64rJTKPvH+bDEa+LdxWQC0EBF7y7on0GvztjHqhXlMWbrbZ3uovvS1E6G5ov2HS/LjVLrgPPO/3zv4NBqkwMdEpZSvU5unzjiSWGv2BvhSRAzwijFmkt/rHYFdXs/z3dv2xHhdRz3++VoANuwtpl+EsyMOmvB1wibaWjB+CAP+NJtLTm/DQ1edybmdmrO/pMznY6JSKjUs+8NQGtZPnUpYrMF+kDGmQETaAV+JyHpjzFyv162yIJb1ZREZC4wF6NIlcJWhRAkocJAG2z1FJ2oe3/P2kvgVyG3Wg5dxSvNGAW8s8V7FRykVHatFyZPZ/yOmYG+MKXB/LxSRj4B+gHewzwc6ez3vBFj2O3J/KpgEkJubmxIto3YLMWP1XkevO+mWCyivqqbawC/fW0azhvV0/hilVEyiDvYi0gSoY4wpcT8eBjzut9unwH0iMhlXw2yRMSalUjje/HP0/s/fX7yT301ZFdcyvHH7hVzes13N8345rTJqTm2lstmtA3P4ZHn69bNvD3zkbqCsB7xrjPlCRO4GMMa8DEzH1e1yM66ul7fHVtz48/+Y5f00HoG+ReP6TP/lJZRXVpPTpknA66ekUAOPUio253dpmbRrRx3sjTFbgfMstr/s9dgA90Z7jVSwfm/8FhPQFZSUUomiI2i9+Ofor3p+ruV+sRjWqz1PXHsO7XVqX6VUAmmw9xOv1vLfjzqL3p1a0K9bq/hcQCmlQtBg7yVeC5Kc3q4pd16SusuVKaUyX+r0+M9Qtw3M4StdY1UplWRas/dSVV3Nd5uCTfNj38s3n8+Qs9rz/KyN3HN5DwdKppRSscmqYH+4tJwN+0oY0N16SoQ3f7C3FGEwDw3vyT2XnVYzX85vrzozpvMppZRTsirY3/SvhazdU8zWp0c6Oq/7Z/ddzLmdmoffUSmlkiSrgv3aPc72mV/+x6G0aKyTkCmlUl9WNdB6ulU61elGA71SKl1kVc3ewzWwN7o0ztSfD+RYWRUFR447WyillIqjrAr2gqtW76nZHy4tZ+G2g7aPX/G/w3R5P6VUWsquYC8CxmAMHDxaxgVPzrJ1nM5ho5RKd9mVs3d/rzaGvcUnQu7r8fkvLo5fgZRSKkGyqmbvLf9w6Jz7Tf27cO/gHpza4qQElUgppeInK4O9MXDXW8GXEvz2t5fTtXXg3PJKKZWusirYe7peHi2rDLqPNsIqpTJRVgV7jwufsm6YXfXoMJo10kCvlMo8WRXsK6qCD6fSHjdKqUwWdW8cEeksInNEZJ2IrBGR+y32uVxEikRkufvrj7EVN3pV1daBvmOLkzTQK6UyXiw1+0rg18aYpSLSDFgiIl8ZY9b67fedMebqGK7jiMrq6oBtbZs15PtxVyShNEoplVhR1+yNMXuMMUvdj0uAdUBHpwrmtCc+930PatqwHosfuTJJpVFKqcRyZFCViOQAfYGFFi9fJCIrRGSGiJztxPWi8faCnT7PVz06LEklUUqpxIs52ItIU2AK8IAxxn8O4aVAV2PMecDfgY9DnGesiOSJSN7+/ftjLVaNL1bvIWfcNJ9tW58eWbPAiFJKZYOYgr2I1McV6N8xxkz1f90YU2yMOep+PB2oLyJtrM5ljJlkjMk1xuS2bds2lmLVOHi0jLvfXuqzbf0Twx1duEQppdJBLL1xBHgVWGeMeTbIPqe490NE+rmvZ3+ayRh9u9H3E8Lc3w6mUf26ibq8UkqljFh64wwCbgFWichy97aHgS4AxpiXgeuBe0SkEjgOjDGuyeQTYsrSfJ/nXVo3TtSllVIqpUQd7I0x8wizAogx5kXgxWivEakJM9Zz5Fg5E67rDcD3m2s/RGhfeqVUNsuoKY5f/nYLkxfvAuBERVXN9n//rF+yiqSUUikho4K9t6enrwNgaK/2XHqGMw2+SimVrjIy2JeWVfLvH3YA8Jfreye5NEoplXwZGeyve2k+AHUEWjRukOTSKKVU8mVksD9wtAyAj+8dlOSSKKVUasjQYF/O7YNy6N2pRbKLopRSKSEjgz1Atza6rKBSSnlkbLAffV7KTsCplFIJl5HB/ldXnkHzxrq8oFJKeWRksO+q0yIopZSPjAz2OgeOUkr5yshgf86pzZNdBKWUsvSbYWfQv1sr2jVrmNDrSgInobQtNzfX5OXlJbsYSimVNkRkiTEmN9jrGVmzV0op5UuDvVJKZQEN9koplQU02CulVBbQYK+UUllAg71SSmUBDfZKKZUFNNgrpVQWSMlBVSKyH9gR5eFtgAMOFicd6D1nvmy7X9B7jlRXY0zQBbdTMtjHQkTyQo0iy0R6z5kv2+4X9J6dpmkcpZTKAhrslVIqC2RisJ+U7AIkgd5z5su2+wW9Z0dlXM5eKaVUoEys2SullPKTMcFeRIaLyAYR2Swi45JdnliISGcRmSMi60RkjYjc797eSkS+EpFN7u8tvY4Z7773DSJyldf2C0Rklfu1F0REknFPdohIXRFZJiKfu59n+v22EJH/iMh69+/6oiy451+5/6ZXi8h7ItIo0+5ZRF4TkUIRWe21zbF7FJGGIvK+e/tCEcmxVTBjTNp/AXWBLUB3oAGwAuiV7HLFcD8dgPPdj5sBG4FewDPAOPf2ccCf3Y97ue+5IdDN/bOo635tEXARIMAMYESy7y/EfT8IvAt87n6e6ff7JnCn+3EDoEUm3zPQEdgGnOR+/gFwW6bdM3ApcD6w2mubY/cI/Bx42f14DPC+rXIl+wfj0A/3ImCm1/PxwPhkl8vB+/sEGApsADq4t3UANljdLzDT/TPpAKz32n4D8Eqy7yfIPXYCZgNXUBvsM/l+T3YHPvHbnsn33BHYBbQC6gGfA8My8Z6BHL9g79g9evZxP66HaxCWhCtTpqRxPH9EHvnubWnP/RGtL7AQaG+M2QPg/t7OvVuw++/ofuy/PRU9DzwEVHtty+T77Q7sB153p67+JSJNyOB7NsbsBv4P2AnsAYqMMV+Swffsxcl7rDnGGFMJFAGtwxUgU4K9Vb4u7bsZiUhTYArwgDGmONSuFttMiO0pRUSuBgqNMUvsHmKxLW3u160ero/6Lxlj+gKluD7eB5P29+zOU4/Gla44FWgiIjeHOsRiW1rdsw3R3GNU958pwT4f6Oz1vBNQkKSyOEJE6uMK9O8YY6a6N+8TkQ7u1zsAhe7twe4/3/3Yf3uqGQRcIyLbgcnAFSLyNpl7v+Aqa74xZqH7+X9wBf9MvucrgW3GmP3GmApgKjCQzL5nDyfvseYYEakHNAcOhStApgT7xcDpItJNRBrgarT4NMllipq71f1VYJ0x5lmvlz4FbnU/vhVXLt+zfYy7lb4bcDqwyP1xsUREBrjP+VOvY1KGMWa8MaaTMSYH1+/ua2PMzWTo/QIYY/YCu0Skp3vTEGAtGXzPuNI3A0SksbusQ4B1ZPY9ezh5j97nuh7X/0v4TzbJbshwsEFkJK5eK1uAR5Jdnhjv5WJcH8tWAsvdXyNx5eVmA5vc31t5HfOI+9434NUzAcgFVrtfexEbDTlJvvfLqW2gzej7BfoAee7f88dAyyy458eA9e7yvoWrF0pG3TPwHq42iQpctfA7nLxHoBHwIbAZV4+d7nbKpSNolVIqC2RKGkcppVQIGuyVUioLaLBXSqksoMFeKaWygAZ7pZTKAhrslVIqC2iwV0qpLKDBXimlssD/Bxb6QDi+KPiBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count= 1225\n",
      "100.0 %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOVElEQVR4nO3dXbBddXnH8e+vifiugDk6mJcmzKS0uRDFU4S+Uqma0E4zneEiWEUZmQxT07HtRQnjTDsdb7S2HYcBSTM0VfsCtspoSqNph744HQskTBEIEDyCJcfQJpQO7egFRp5e7BXcbk7OXjnseDj/fD8ze85a//Xfaz/PTvJjnbX22qSqkCS14UcWuwBJ0uQY6pLUEENdkhpiqEtSQwx1SWrI8sV64RUrVtTatWsX6+UlaUm65557nqyqqRNtX7RQX7t2Lfv371+sl5ekJSnJf8y33dMvktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSFjQz3JriRHkjxwgu1Jcn2SmST3Jblg8mVKkvroc6T+KWDjPNs3Aeu7x1bgphdeliRpIcaGelV9BXhqnimbgc/UwJ3AmUnOmVSBkqT+JnFH6Urg0ND6bDf2xOjEJFsZHM2zZs2aBb/g2u1/9wPr3/zoLz1v7GTHT3buC61jMWr2PXrxvR74Ho2bC+29R8f7OhUmcaE0c4zN+b9TqqqdVTVdVdNTUyf86gJJ0gJNItRngdVD66uAwxPYryTpJE0i1HcDV3afgrkIeLqqnnfqRZJ06o09p57kFuASYEWSWeD3gJcAVNUOYA9wGTADfAe46lQVK0ma39hQr6orxmwv4IMTq0iStGDeUSpJDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqSK9QT7IxycEkM0m2z7H9tUn+NsnXkhxIctXkS5UkjTM21JMsA24ENgEbgCuSbBiZ9kHgwao6H7gE+KMkZ0y4VknSGH2O1C8EZqrq0ap6BrgV2Dwyp4BXJwnwKuAp4NhEK5UkjdUn1FcCh4bWZ7uxYTcAPwEcBu4HPlRVz47uKMnWJPuT7D969OgCS5YknUifUM8cYzWy/i7gXuCNwJuBG5K85nlPqtpZVdNVNT01NXWSpUqSxukT6rPA6qH1VQyOyIddBdxWAzPAY8CPT6ZESVJffUJ9H7A+ybru4ucWYPfInMeBSwGSvAE4D3h0koVKksZbPm5CVR1Lsg3YCywDdlXVgSTXdNt3AB8BPpXkfgana66tqidPYd2SpDmMDXWAqtoD7BkZ2zG0fBh452RLkySdLO8olaSGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDWkV6gn2ZjkYJKZJNtPMOeSJPcmOZDkXyZbpiSpj+XjJiRZBtwIvAOYBfYl2V1VDw7NORP4JLCxqh5P8vpTVK8kaR59jtQvBGaq6tGqega4Fdg8MufdwG1V9ThAVR2ZbJmSpD76hPpK4NDQ+mw3NuzHgLOS/HOSe5JcOakCJUn9jT39AmSOsZpjP28FLgVeDvxbkjur6pEf2FGyFdgKsGbNmpOvVpI0rz5H6rPA6qH1VcDhOeZ8uaq+XVVPAl8Bzh/dUVXtrKrpqpqemppaaM2SpBPoE+r7gPVJ1iU5A9gC7B6Z80XgZ5MsT/IK4G3AQ5MtVZI0ztjTL1V1LMk2YC+wDNhVVQeSXNNt31FVDyX5MnAf8Cxwc1U9cCoLlyQ9X59z6lTVHmDPyNiOkfWPAx+fXGmSpJPlHaWS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktSQXqGeZGOSg0lmkmyfZ95PJvlekssnV6Ikqa+xoZ5kGXAjsAnYAFyRZMMJ5n0M2DvpIiVJ/fQ5Ur8QmKmqR6vqGeBWYPMc834D+DxwZIL1SZJOQp9QXwkcGlqf7caek2Ql8KvAjvl2lGRrkv1J9h89evRka5UkjdEn1DPHWI2sfwK4tqq+N9+OqmpnVU1X1fTU1FTPEiVJfS3vMWcWWD20vgo4PDJnGrg1CcAK4LIkx6rqC5MoUpLUT59Q3wesT7IO+BawBXj38ISqWnd8OcmngNsNdEn64Rsb6lV1LMk2Bp9qWQbsqqoDSa7pts97Hl2S9MPT50idqtoD7BkZmzPMq+r9L7wsSdJCeEepJDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSG9Qj3JxiQHk8wk2T7H9l9Lcl/3+GqS8ydfqiRpnLGhnmQZcCOwCdgAXJFkw8i0x4Cfr6o3AR8Bdk66UEnSeH2O1C8EZqrq0ap6BrgV2Dw8oaq+WlX/063eCayabJmSpD76hPpK4NDQ+mw3diIfAL4014YkW5PsT7L/6NGj/auUJPXSJ9Qzx1jNOTH5BQahfu1c26tqZ1VNV9X01NRU/yolSb0s7zFnFlg9tL4KODw6KcmbgJuBTVX135MpT5J0Mvocqe8D1idZl+QMYAuwe3hCkjXAbcB7q+qRyZcpSepj7JF6VR1Lsg3YCywDdlXVgSTXdNt3AL8LvA74ZBKAY1U1ferKliTNpc/pF6pqD7BnZGzH0PLVwNWTLU2SdLK8o1SSGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhpiqEtSQwx1SWqIoS5JDTHUJakhhrokNcRQl6SGGOqS1BBDXZIaYqhLUkMMdUlqiKEuSQ0x1CWpIYa6JDXEUJekhhjqktQQQ12SGmKoS1JDDHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUEENdkhrSK9STbExyMMlMku1zbE+S67vt9yW5YPKlSpLGGRvqSZYBNwKbgA3AFUk2jEzbBKzvHluBmyZcpySphz5H6hcCM1X1aFU9A9wKbB6Zsxn4TA3cCZyZ5JwJ1ypJGiNVNf+E5HJgY1Vd3a2/F3hbVW0bmnM78NGq+tdu/Q7g2qraP7KvrQyO5AHOAw4usO4VwJMLfO5SZc+nB3s+PbyQnn+0qqZOtHF5jx1kjrHR/xL0mUNV7QR29njN+QtK9lfV9Avdz1Jiz6cHez49nMqe+5x+mQVWD62vAg4vYI4k6RTrE+r7gPVJ1iU5A9gC7B6Zsxu4svsUzEXA01X1xIRrlSSNMfb0S1UdS7IN2AssA3ZV1YEk13TbdwB7gMuAGeA7wFWnrmRgAqdwliB7Pj3Y8+nhlPU89kKpJGnp8I5SSWqIoS5JDVlyoT7uKwuWiiSrk/xTkoeSHEjyoW787CT/kOTr3c+zhp5zXdf3wSTvGhp/a5L7u23XJ5nrI6YvGkmWJfn37v6G5ntOcmaSzyV5uPvzvvg06Pm3ur/XDyS5JcnLWus5ya4kR5I8MDQ2sR6TvDTJZ7vxu5Ks7VVYVS2ZB4MLtd8AzgXOAL4GbFjsuhbYyznABd3yq4FHGHwNwx8A27vx7cDHuuUNXb8vBdZ178OybtvdwMUM7hf4ErBpsfsb0/tvA38F3N6tN90z8Gng6m75DODMlnsGVgKPAS/v1v8aeH9rPQM/B1wAPDA0NrEegV8HdnTLW4DP9qprsd+Yk3wTLwb2Dq1fB1y32HVNqLcvAu9gcJftOd3YOcDBuXpl8Gmki7s5Dw+NXwH8yWL3M0+fq4A7gLfz/VBvtmfgNV3AZWS85Z5XAoeAsxl8wu524J0t9gysHQn1ifV4fE63vJzBHagZV9NSO/1y/C/LcbPd2JLW/Vr1FuAu4A3Vfca/+/n6btqJel/ZLY+Ov1h9Avgd4NmhsZZ7Phc4CvxZd8rp5iSvpOGeq+pbwB8CjwNPMLhv5e9puOchk+zxuedU1THgaeB14wpYaqHe6+sIlpIkrwI+D/xmVf3vfFPnGKt5xl90kvwycKSq7un7lDnGllTPDI6wLgBuqqq3AN9m8Gv5iSz5nrvzyJsZnGZ4I/DKJO+Z7ylzjC2pnntYSI8L6n+phXpTX0eQ5CUMAv0vq+q2bvi/0n3DZffzSDd+ot5nu+XR8RejnwZ+Jck3GXzb59uT/AVt9zwLzFbVXd365xiEfMs9/yLwWFUdrarvArcBP0XbPR83yR6fe06S5cBrgafGFbDUQr3PVxYsCd0V7j8FHqqqPx7atBt4X7f8Pgbn2o+Pb+muiK9j8N31d3e/4v1fkou6fV459JwXlaq6rqpWVdVaBn92/1hV76Htnv8TOJTkvG7oUuBBGu6ZwWmXi5K8oqv1UuAh2u75uEn2OLyvyxn8exn/m8piX2hYwIWJyxh8UuQbwIcXu54X0MfPMPhV6j7g3u5xGYNzZncAX+9+nj30nA93fR9k6FMAwDTwQLftBnpcTFnsB3AJ379Q2nTPwJuB/d2f9ReAs06Dnn8feLir988ZfOqjqZ6BWxhcM/gug6PqD0yyR+BlwN8w+PqVu4Fz+9Tl1wRIUkOW2ukXSdI8DHVJaoihLkkNMdQlqSGGuiQ1xFCXpIYY6pLUkP8HUg9mCGiVVBQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fin\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy\n",
    "import math\n",
    "from random import randint\n",
    "\n",
    "def Exp(N,expN):\n",
    "    if N>10:\n",
    "        pi=2*N**0.5/math.log(N)\n",
    "        r=randint(5,N)\n",
    "        if r==N:\n",
    "            r=N-2\n",
    "        s=N-r\n",
    "        if r<=10:\n",
    "            return expN[r],expN\n",
    "        else:\n",
    "            expN[N]=max(pi-expN[r]-expN[s],expN[N])\n",
    "        if s<=10:\n",
    "            return expN[s],expN\n",
    "        else:\n",
    "            expN[N]=max(pi-expN[r]-expN[s],expN[N])\n",
    "        return expN[N],expN\n",
    "    else:\n",
    "        return expN[N],expN\n",
    "\n",
    "n_iter=100\n",
    "Nmax=10000\n",
    "primes=list(sympy.primerange(2,Nmax+1))\n",
    "prime_train=[2,3,5,7]\n",
    "expN=[1]*(Nmax+2)\n",
    "expN[2]=1\n",
    "expN[3]=1\n",
    "expN[4]=1\n",
    "expN[5]=1\n",
    "expN[6]=2\n",
    "expN[7]=1\n",
    "expN[8]=1\n",
    "expN[9]=1\n",
    "expN[10]=2\n",
    "\n",
    "PrimeTrue=[]\n",
    "PrimeCheck=[]\n",
    "for n in range(10,Nmax+1):\n",
    "    BOOL=True\n",
    "    if n in primes:\n",
    "        PrimeTrue.append(1)\n",
    "    else:\n",
    "        PrimeTrue.append(0)\n",
    "    pi=2*n**0.5/math.log(n)\n",
    "    \n",
    "    BOOL=True\n",
    "    for m in range(n_iter):\n",
    "        r=randint(5,n)\n",
    "        if r==n:\n",
    "            r=n-2\n",
    "        s=n-r\n",
    "        g=math.gcd(s,r)\n",
    "\n",
    "        if BOOL:\n",
    "            if g!=1:\n",
    "                PrimeCheck.append(0)\n",
    "                prob=0\n",
    "                BOOL=False\n",
    "                continue\n",
    "            if (n%2==0 or n%3==0) or (n%5==0 or n%7==0):\n",
    "                PrimeCheck.append(0)\n",
    "                prob=0\n",
    "                BOOL=False\n",
    "                continue\n",
    "        if m==0:\n",
    "            ermax,expN=Exp(r,expN)\n",
    "            esmax,expN=Exp(s,expN)\n",
    "        else:\n",
    "            er,expN=Exp(r,expN)\n",
    "            es,expN=Exp(s,expN)\n",
    "            ermax=max(er,ermax)\n",
    "            esmax=max(es,esmax)\n",
    "        prob=(ermax+esmax)/pi\n",
    "    if BOOL:\n",
    "        if prob >0.5:\n",
    "            PrimeCheck.append(1)\n",
    "        else:\n",
    "            PrimeCheck.append(0)\n",
    "\n",
    "x=[]\n",
    "Hit=[]\n",
    "count=0\n",
    "for n in range(Nmax-10):\n",
    "    C=PrimeCheck[n]\n",
    "    T=PrimeTrue[n]\n",
    "    if C==T:\n",
    "        count+=1\n",
    "        x.append(n)\n",
    "        Hit.append(1)\n",
    "    else:\n",
    "        x.append(n)\n",
    "        Hit.append(0)\n",
    "\n",
    "print('count=',count)\n",
    "print(round(count/Nmax*100,2),'%')\n",
    "\n",
    "plt.bar(x,Hit)\n",
    "plt.show()\n",
    "plt.plot(expN)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "count=0\n",
    "for n in range(Nmax-10):\n",
    "    C=PrimeCheck[n]\n",
    "    T=PrimeTrue[n]\n",
    "    if C==T and T==1:\n",
    "        count+=1\n",
    "        x.append(n)\n",
    "        Hit.append(1)\n",
    "    else:\n",
    "        x.append(n)\n",
    "        Hit.append(0)\n",
    "\n",
    "print('count=',count)\n",
    "print(round(count/(len(primes)-4)*100,2),'%')\n",
    "\n",
    "plt.bar(x,Hit)\n",
    "plt.show()\n",
    "\n",
    "print('Fin')"
   ]
  }
 ],
 "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}