2_examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6fe8025a-c5aa-4dec-8d47-1ee2961daddb",
   "metadata": {},
   "source": [
    "# Chapter two worked examples\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c0e5f067-c446-45c9-9e47-1a5b3549abd3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import ant\n",
    "from math import cos, sin, pi, fabs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "11378ebc-68ac-4963-a5ed-3a3b6c85d4aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# To list available functions in the ant module: type ant., Tab\n",
    "# Docstring for the module: type ant. + Shift + Tab\n",
    "# ant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5cc9f83-e73c-4455-b5fb-b006d4967f42",
   "metadata": {},
   "source": [
    "## Problem 0: Bonus Fourier transform check!\n",
    "\n",
    "\n",
    "The aim here is to show empirically that that the totient at $n$ is equal to the first term in the discrete \n",
    "Fourier transform of the gcd function evaluated from $1$ to $n$. This claim can be found in the Wikipedia \n",
    "page on the [Euler totient function](https://en.wikipedia.org/wiki/Euler%27s_totient_function#Fourier_transform).\n",
    "\n",
    "\n",
    "\n",
    "The discrete Fourier transform operates on a vector $x$ of $N$ (real / complex) values, producing a new vector $X$ of $N$ complex values.\n",
    "The transform is defined for the $k$th element of the resultant transform vector as:\n",
    "\n",
    "\n",
    "$\\begin{align}\\large{X_k = \\sum_{n=0}^{N-1}x_n \\cdot e^{-2 \\pi i \\frac{k}{N}n}}\\end{align}$\n",
    "\n",
    "\n",
    "Note the indexing runs $0 \\dots N-1$. The DFT in the context of the totient function will use an index $k$ running $1 \\dots n$ \n",
    "to follow the notation of the article cited above. An element of the vector $x$ is defined $x_k = gcd(k, n) \\textrm{ for } k \\in \\{1, \\dots, n\\}$.\n",
    "For example taking $n=10$ results in $x = \\{ 1, 2, 1, 2, 5, 2, 1, 2, 1, 10 \\}$. Now we have the DFT defined for the $m$th element of $X$ as\n",
    "\n",
    "\n",
    "$\\begin{align}\\large{X_m = \\sum_{k=1}^{n} x_k \\cdot e^{-2 \\pi i \\frac{m}{n}k}}\\end{align}$\n",
    "\n",
    "\n",
    "The important detail to verify here is that the last term in this sum is equivalent to the first term of the previous sum.\n",
    "\n",
    "\n",
    "The claim is that $\\varphi(n)$ equals the real part of $X_1$. Writing this in terms of the\n",
    "real cosine component of the above equation: \n",
    "\n",
    "\n",
    "$\\begin{align}\\textrm{Claim: }\\varphi(n) = \\sum_{k=1}^{n} gcd(k, n) \\cdot cos \\left( -2 \\pi \\frac{k}{n} \\right) \\end{align}$\n",
    "\n",
    "\n",
    "Also the imaginary part should be zero. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3eeda12e-3838-4740-93a1-1ddb3f72ad6c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 Totient is 1    ...DFT is 1\n",
      "2 Totient is 1    ...DFT is 1\n",
      "3 Totient is 2    ...DFT is 2\n",
      "4 Totient is 2    ...DFT is 2\n",
      "5 Totient is 4    ...DFT is 4\n",
      "6 Totient is 2    ...DFT is 2\n",
      "7 Totient is 6    ...DFT is 6\n",
      "8 Totient is 4    ...DFT is 4\n",
      "9 Totient is 6    ...DFT is 6\n",
      "10 Totient is 4   ...DFT is 4\n",
      "11 Totient is 10  ...DFT is 10\n",
      "12 Totient is 4   ...DFT is 4\n",
      "13 Totient is 12  ...DFT is 12\n",
      "14 Totient is 6   ...DFT is 6\n",
      "15 Totient is 8   ...DFT is 8\n",
      "16 Totient is 8   ...DFT is 8\n",
      "17 Totient is 16  ...DFT is 16\n",
      "18 Totient is 6   ...DFT is 6\n",
      "19 Totient is 18  ...DFT is 18\n",
      "20 Totient is 8   ...DFT is 8\n",
      "21 Totient is 12  ...DFT is 12\n",
      "22 Totient is 10  ...DFT is 10\n",
      "23 Totient is 22  ...DFT is 22\n",
      "24 Totient is 8   ...DFT is 8\n",
      "25 Totient is 20  ...DFT is 20\n",
      "26 Totient is 12  ...DFT is 12\n",
      "27 Totient is 18  ...DFT is 18\n",
      "28 Totient is 12  ...DFT is 12\n",
      "29 Totient is 28  ...DFT is 28\n",
      "30 Totient is 8   ...DFT is 8\n",
      "\n",
      "  The DFT matched the totient in all cases...\n",
      "     ...and meanwhile the imaginary part was zero as expected.\n"
     ]
    }
   ],
   "source": [
    "space_buffer = 4\n",
    "real_matches, imag_is_zero = True, True\n",
    "for n in range(1, 31):\n",
    "    DFT_real, DFT_imag, totient = 0., 0., ant.Totient(n)\n",
    "    for k in range(1, n+1):\n",
    "        DFT_real += ant.gcd(k, n) * cos(2 * pi * k / n)\n",
    "        DFT_imag += ant.gcd(k, n) * -1 * sin(2 * pi * k / n)\n",
    "    DFT_real    = int(round(DFT_real, 8))\n",
    "    this_buffer = space_buffer - len(str(n)) - len(str(totient))\n",
    "    print(n, 'Totient is', totient, ' '*this_buffer, '...DFT is', DFT_real)\n",
    "    if not DFT_real == totient: real_matches = False\n",
    "    if fabs(DFT_imag) > 1e-8:   imag_is_zero = False\n",
    "\n",
    "if not real_matches: print('\\n  Oh no: There was at least one mismatch between the DFT and the totient...')\n",
    "else:                print('\\n  The DFT matched the totient in all cases...')\n",
    "if not imag_is_zero: print('     ...and meanwhile the imaginary part of the DFT was somewhere non-zero.')\n",
    "else:                print('     ...and meanwhile the imaginary part was zero as expected.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "241d1e6d-39f9-4800-9582-b7125b767a68",
   "metadata": {},
   "source": [
    "## Problem 2.1c) Find n for which $\\varphi(n)$ = 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3b845bb6-1a0f-49cc-beaa-529e44f6d07e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13 12\n",
      "21 12\n",
      "26 12\n",
      "28 12\n",
      "36 12\n",
      "42 12\n"
     ]
    }
   ],
   "source": [
    "for n in range(1, 625):\n",
    "    if ant.Totient(n) == 12: print(n, ant.Totient(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c192fe9-bc6f-42aa-b73c-cd44146a6f79",
   "metadata": {},
   "source": [
    "## Does totient $\\varphi(n)$ have a lower bound?\n",
    "\n",
    "Making a scatter plot could give some indication..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "a915b751-3152-45e4-9a76-83a1d6a91d3f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "200\n",
      "400\n",
      "600\n",
      "800\n",
      "1000\n",
      "1200\n",
      "1400\n",
      "1600\n",
      "1800\n",
      "2000\n",
      "2200\n",
      "2400\n",
      "2600\n",
      "2800\n",
      "3000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fe1fa809970>"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAGdCAYAAADjWSL8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABoMUlEQVR4nO2dfZBVxZn/n3svl4FLhhsvCsOLjNaK2SQIlThZhU1pNAqkMmoqf0BmCCG1G3UMYzYVt3ZJqAruVo3cbFXcSpXx3mxe3GyVBq1Fd63EkGCC5GUGkwhWUFZjGUHCiygDwygMyMzz+4NfH/r29Dmnz/s5d76fqimYe/v06dPnzOlvf/vp7hwzMwEAAAAAZIx80gUAAAAAAPADRAwAAAAAMglEDAAAAAAyCUQMAAAAADIJRAwAAAAAMglEDAAAAAAyCUQMAAAAADIJRAwAAAAAMsmkpAsQFWNjY3To0CFqbW2lXC6XdHEAAAAAYAAz0/DwMM2ZM4fyeWevpWlFzKFDh+jSSy9NuhgAAAAA8MGBAwdo3rx5jmmaVsS0trYS0flKmD59esKlAQAAAIAJJ0+epEsvvdRqx51oWhEjhpCmT58OEQMAAABkDJNQEAT2AgAAACCTQMQAAAAAIJNAxAAAAAAgk0DEAAAAACCTQMQAAAAAIJNAxAAAAAAgk0DEAAAAACCTQMQAAAAAIJNAxAAAAAAgk3gSMbVajRYtWmStgrtkyRL66U9/an3PzHTvvffSnDlzaOrUqfSxj32MXnzxxYY8zpw5Q3fffTddfPHFNG3aNLr11lvpL3/5S0Oa48eP05o1a6hcLlO5XKY1a9bQiRMn/F8lAAAAAJoOTyJm3rx5VK1W6Q9/+AP94Q9/oBtvvJFuu+02S6j827/9G91///30wAMP0O9//3tqa2ujm2++mYaHh608vvzlL9MTTzxBmzdvpt/85jf09ttvU2dnJ42Ojlppuru76fnnn6etW7fS1q1b6fnnn6c1a9aEdMkAAAAAaAo4IBdddBF/73vf47GxMW5ra+NqtWp9NzIywuVymev1OjMznzhxgovFIm/evNlKc/DgQc7n87x161ZmZt67dy8TEe/cudNKMzAwwETEL730knG5hoaGmIh4aGgo6CUCAAAAICa8tN++Y2JGR0dp8+bN9M4779CSJUvotddeoyNHjtCyZcusNC0tLXT99ddTf38/ERE999xz9O677zakmTNnDi1cuNBKMzAwQOVyma655horzbXXXkvlctlKo+PMmTN08uTJhh8AAAAAREO9XqfLLruM6vV6YmXwLGL27NlD73nPe6ilpYV6enroiSeeoA984AN05MgRIiKaNWtWQ/pZs2ZZ3x05coQmT55MF110kWOamTNnjjvvzJkzrTQ6Nm3aZMXQlMtluvTSS71eGgAAAAAMqVartH//fqpWq4mVwbOIed/73kfPP/887dy5k+666y5au3Yt7d271/pe3TqbmV2301bT6NK75fPVr36VhoaGrJ8DBw6YXhIAAAAAPLJ+/Xpqb2+n9evXJ1aGSV4PmDx5Ml1xxRVERNTR0UG///3v6Vvf+hb98z//MxGdd1Jmz55tpT969KjlzrS1tdHZs2fp+PHjDW7M0aNHaenSpVaaN954Y9x533zzzXEuj0xLSwu1tLR4vRwAAAAA+KCnp4d6enoSLUPgdWKYmc6cOUOXX345tbW10bZt26zvzp49Szt27LAEytVXX03FYrEhzeHDh+mFF16w0ixZsoSGhobod7/7nZXm2WefpaGhISsNAAAAAOKju7ubJk2aRN3d3UkXpQFPTszXvvY1+sQnPkGXXnopDQ8P0+bNm+mZZ56hrVu3Ui6Xoy9/+ct033330YIFC2jBggV03333UalUsi66XC7T3//939M999xDM2bMoEqlQv/4j/9IV111Fd10001ERPT+97+fVqxYQbfffjt95zvfISKiO+64gzo7O+l973tfyJcPAAAAADvq9Tpt2LCBBgcHiYjoscceo0ceeSThUl3Ak4h54403aM2aNXT48GEql8u0aNEi2rp1K918881ERPRP//RPdPr0afriF79Ix48fp2uuuYZ+/vOfU2trq5XHv//7v9OkSZNo5cqVdPr0afr4xz9O//mf/0mFQsFK8/DDD9OXvvQlaxbTrbfeSg888EAY1wsAAAAAQ2QBQ0S0cuXKBEsznhwzc9KFiIKTJ09SuVymoaEhmj59etLFAQAAADJDvV6narVKb775Jp06dYry+Tx9+9vfjiUGxkv7jb2TAAAAANCAmD49ZcoUam9vj03AeAUiBgAAAABEdGEBu6VLl1J7ezv19fXRvn37UilgiHxMsQYAAABAcyGGj4aHh60YmH379iVbKAPgxAAAAAATlHq9TjNmzKB169bR/v37iYgSX8DOC3BiAAAAgAlIvV6n3t5eGh0dJSKiQqFAfX19qR060gEnBgAAAJhgyAIml8tRpVKhBx54IFMChggiBgAAAJgwiMDdDRs20OjoKBUKBXrwwQfp2LFjngRMGnawJsI6MQAAAMCE4bLLLqP9+/dTpVKh1tZWWr9+vS/3ReTT3t4eegAw1okBAAAAABFdCN6dMWNGKFOn6/U6DQ8PU6VSSTwAGE4MAAAA0KSowbthOCdRujBEcGIAAACACY8ueDcM52T9+vWpmYYNJwYAAABoMrq7u+lHP/oREZ2fOp2lmUdwYgAAAIAJzGOPPWb9P0sCxisQMQAAAEATIE97XrlyJRUKBerq6mpaAUMEEQMAAABkGnntl/3791O1WqVHHnmEzp07R4888ojnfJzWfknL+jACiBgAAAAgo9TrdWvfo5GREV8Bt6oI6u3ttRUp1WrVEkppACIGAAAAyCDd3d1011130djYGBERTZkyxdfaL0KYEJ0PAh4dHbUVKWmamUQEEQMAAABkCuGcPProo9ZnlUqF+vr6fOUnhElfXx898MADjiKlp6fH9yJ5UQARAwAAAGQEsfbL/v37acqUKVbwrpe9j9S4FiFMiM67Mn63IkgCrBMDAAAAZASxWq6ftV/q9TpVq1UaHh6mwcHBcSvuRr0SrylYJwYAAABoIoR7IvY+UgWMyawhOfZFN2Qkx7ukbRaSHXBiAAAAgJRSr9dpw4YNdOLECRobG7N1SUxcFOHEmAwXJenKwIkBAAAAmoBqtUqDg4M0NjZGhULBNuDWZNaQl6BcXX5pdGfgxAAAAAApo7u7mx577DH60Ic+RH/+85+JiKivry/RgNu43Bkv7fekyEoBAAAAAE+IIZ8DBw7Q2NgY7d69m86dO5d0sYjovDsjhqPSAoaTAAAAgBQgFq8Ts48KhQKtXLkysfLYTcVO0/RriBgAAAAgQXSL17W2tnre+8j0PKYxLbotBtIWFwMRAwAAACSE6r7kcjkqlUq+Vt91Exhe9z3SBfdi7yQAAAAAUL1epx/96EfW762trTQ2NkbvvPOOryEbN4Hhdd8jeSVfIY7StncScZMyNDTERMRDQ0NJFwUAAACwqNVqXKlUOJ/PMxExEXGpVOJarWZ8fHt7+7j0dp8Hpb29nYmI29vbQ83XDi/tN5wYAAAAICbE3kfy2i+1Ws2T+2LnuOicE7eymKRLnfsigXViAAAAgBjo7u62ho9yuRxddNFFvtZ+cVt513R/pbTslaSCFXsBAACAFKHGvzz44IOedp6WcZvqvH79eioUCjQ6OuoYgJtmh8UUiBgAAAAgYmQx0dXV5Xn36RkzZtCMGTOMpjb39PTQAw88EOo2BGkFw0kAAABABMjDPkRkvPmiihj2IaLUDf1EAYaTAAAAgIQQzsm6deusANwgrsf69eupUqlQpVIJdegnbQvX+QFODAAAABASYvbR6OgoEZFrcG2SILAXAAAAABbVapVGR0cpl8tRpVJJrYAhQmAvAAAAMOHp7u6mSZMmUXd3tyUMvM4+Mh3aCXMICIG9KQbDSQAAAOJg0qRJNDo6SoVCgc6dO+crD9OhnbQOAYUJhpMAAACAiBGuyIc+9CEqFAq0cuVK33mZDu00wxBQmMCJAQAAAHwwEVyRJIATAwAAAESAHJOSVlekGaZOmwInBgAAADBA3vsoze5L1h0iODEAAABASIjF6+S9j9Lmvsik1SGKAk8iZtOmTfSRj3yEWltbaebMmfSpT32KXn755YY0n//85ymXyzX8XHvttQ1pzpw5Q3fffTddfPHFNG3aNLr11lvpL3/5S0Oa48eP05o1a6hcLlO5XKY1a9bQiRMn/F0lAAAA4AOxeN3g4KD1mde9j+KmGaZOm+JJxOzYsYPWrVtHO3fupG3bttG5c+do2bJl9M477zSkW7FiBR0+fNj6eeqppxq+//KXv0xPPPEEbd68mX7zm9/Q22+/TZ2dndYKh0Tnbbvnn3+etm7dSlu3bqXnn3+e1qxZE+BSAQAAAHPk1XfF4nW1Wo0eeeQRo2O9bNrolE9S8S12505VzA0H4OjRo0xEvGPHDuuztWvX8m233WZ7zIkTJ7hYLPLmzZutzw4ePMj5fJ63bt3KzMx79+5lIuKdO3daaQYGBpiI+KWXXjIq29DQEBMRDw0NebwqAAAAgLm9vZ2JiAuFAtdqNV/HEhG3t7cHLkOQPLxSq9W4vb2dK5WK9txRl8lL+x0oJmZoaIiIiCqVSsPnzzzzDM2cOZOuvPJKuv322+no0aPWd8899xy9++67tGzZMuuzOXPm0MKFC6m/v5+IiAYGBqhcLtM111xjpbn22mupXC5baVTOnDlDJ0+ebPgBAAAAvKCbfeRn64CwNm2U41vickCq1WrDrtlq+VMVc+NXKY2NjfEtt9zCH/3oRxs+37x5M//4xz/mPXv28JNPPsmLFy/mD37wgzwyMsLMzA8//DBPnjx5XH4333wz33HHHczM3NfXxwsWLBiXZsGCBXzfffdpy7Nx40ZL9co/cGIAAACYEofzIZwOv+5OkLKZnNtv+cIiFiemt7eX/vjHPzZEaxMRrVq1ij75yU/SwoUL6ZZbbqGf/vSn9Kc//Yl+8pOfuIkpyuVy1u/y/+3SyHz1q1+loaEh6+fAgQM+rgoAAMBEJgqXQXVQhNNRrVZjL5vJubMUGOxLxNx999305JNP0vbt22nevHmOaWfPnk3t7e30yiuvEBFRW1sbnT17lo4fP96Q7ujRozRr1iwrzRtvvDEurzfffNNKo9LS0kLTp09v+AEAAADs0AXfRtGAq8LBrxgJo2ypGgoKAy8Wz9jYGK9bt47nzJnDf/rTn4yOeeutt7ilpYV/+MMfMvOFwN5HH33USnPo0CFtYO+zzz5rpdm5cycCewEAAISGCFylAEM0Yuilq6vLdggm6eGZrOGl/fYkYu666y4ul8v8zDPP8OHDh62fU6dOMTPz8PAw33PPPdzf38+vvfYab9++nZcsWcJz587lkydPWvn09PTwvHnz+Omnn+Zdu3bxjTfeyIsXL+Zz585ZaVasWMGLFi3igYEBHhgY4Kuuuoo7OzuNywoRAwAAQIcQFaVSiYmI8/m8Z4GhzuApFAqxzyJqViITMaQJnCUifuihh5iZ+dSpU7xs2TK+5JJLuFgs8vz583nt2rX8+uuvN+Rz+vRp7u3t5UqlwlOnTuXOzs5xaY4dO8arV6/m1tZWbm1t5dWrV/Px48eNywoRAwAAQKVWq1mCo1Kp+HZIRJCtyKOrq4srlQpXKpVUOi5ZcoMiEzFZAiIGAACAoFarcaVS4Xw+73ntF50A0H2WxJoupqhlS7OogYhhiBgAAAAXkBef87p4nak4SbMwUMuWZsEV22J3AAAAQJoR05uXLl1qLT7ndfE60xk9cU5N9rrwnVq2ZpmllGNmTroQUeBlK28AAADNRb1epw0bNtCJEydobGyM2tvbad++fUkXKzQuu+wy2r9/f9NdF5G39htODAAAgKajWq3S4OAgjY2NUaFQyLzjoNIsTkpQIGIAAAA0BfLidV6Hj8TwTHd3t9EwjdfhnLD3PVKHh+LeWTotO1ljOAkAAEBTMGPGDBocHCQi8jzMIoZnCoUCjY6Ouh7vdTgn6uGfuIeXojwfhpMAAABMGOr1Ok2bNs0SMPl83niYRQ78bW9vp5UrV1KlUqHh4WFHl8Fkd2ndjthhDP/ozhf38FJqhrMinimVGJhiDQAAzY+8eJ34CTp9Wv7MZNq03XTlqKYxxz09Ou6p41gnhiFiAACgmVGX/af/v31AV1eXr3zsFrIzEQx2jXxUjX/coiJu0QQRwxAxAADQrIS1dYDpuUxW621m4MQkAEQMAAA0F7Vazdq0kXysvCvyCOqapHm12zBIWqRBxDBEDAAANBNBY18E6saNfpbhT7qRj5qkRRq2HQAAANBUVKtVGh0dJaLzs4+6urqM137RzeIhItq/fz9Vq9WGzxOfbUPJr8GSprpwA+vEAAAASD31ep2q1SqtX7/eeG8ip7VM/ORnkm8YNPOWAiZgnRgAAACZp7u7myZNmkTd3d2+Nld0chSc8nNzQqJ2KqLOP2mnJ0zgxAAAAEgVwiU5cOCAtffRuXPnfOWRRqcladJ+fXBiAAAAZA6x99G6deto//79NGXKFCoUCrRy5UrPeVWr1YaYFy/EFROSlCOSpZgXVyIOMk4MzE4CAIBsIS9c52f6tIyXGURJrQUT1yygrM2mwuwkAAAAmUHEvgwPDxPR+dlHJjtPh4XOtQni5NihOi9xxb5s2LAh9GtJDTGIqkSAEwMAAOmnVqtZ7ksulwvNMfCz7ktXV5d1/ijci7jXX7FbEyftYLE7hogBAIA0o9v7yOu+Ryb5B90MMkziHtbJ2jCSACKGIWIAACCthLn3UZgNdVYbfUGtVuNKpcKVSiWz18Dsrf3GFGsAAACx0d3dTT/60Y+IiKhQKASOfUnbdOEgU7uDIuqCiFJTH37AFGsAAACp5LHHHrP+H0bwbtqmC0cREKzDbkuFSqVClUolNfURNXBiAAAARI5wKC655BLavXs3rVy5kh555JHEyhGVUxKXE5M2BypM4MQAAABIHLF43YwZM6xpvm+++SadO3cuFAHjZ7G4qJ0SP9sj+CFtDlRSQMQAAACIhA0bNtDg4CANDg4SEfludO3Eih9BElbjn/T+Q37FUtLlDhsMJwEAAIiEGTNm0ODgIOXzefr2t7/t252wGzpJQxBt1oZzslBuDCcBAACIHbWX39fXR+3t7UYCRhzb3d2tDVjVuSdxDd3oUMsUhcMRRZ7NNgwFJwYAAEAoBOnli2MLhQKNjo6mynUxIQqHIwuuSRTAiQEAABAbwjFYunTpuF6+k5sgB/6KY1euXKl1CuT4lzTGdQR1OGQnSq0TL3mqdeO3rtJYx1oiXHQvUbBiLwAARIu6dYBuuX6npfzFd3bf684l/jU5JkuIaxIrGfu9PrVu/NZVknWMXawBAABESr1ep97e3oYVYnWOgZND4WVxNjn+Ja7dn+N0IcQ1rVy5MtCCdWrd+K2rrMTOICYGAACAJ+r1Oq1bt47GxsYCzzwyPV+c8TATNRYlLSAmBgAAQCR0d3fTXXfdRWNjY0RE9N73vjdyYRHXUv6CrLgQACIGAACAB+S9jyqVCvX19UV+zrhFRRxTtzMTOJtyIGIAAAA40t3dTZMmTaLu7m5auXIlFQoF6urqomPHjnlq6P023FGJiiSFhIm7BKFjQORhxgmB2UkAABAMMRMon89bM2f8HF+r1ZjZ/4wXNZ+wSHIGjsk1NeMsLBO8tN8QMQAAAMZRq9Ws6b6lUokLhQJ3dXUZHWc3FdqvGAmrMVfPH5U4CivfqMqXdiBiGCIGAAD80tHRYa1VUigUjBpR3ZoxciMcpEEOqzGP2tkwWTcnraRJMEHEMEQMAAD4RQgYIjJu1IRAqFQq2sYwDUM3XV1dkTbUbnWQZtI0dAURwxAxAADghVqtxpVKhSuVCs+fP5+JiDs6Ojwd79RwxzWUoyOuBjpNboZX0lR2L+03FrsDAABgLfBGRLEs8hblgnLq4nhp3zwSNILF7gAAABghb94YZLl7r0S59os6fTnIFG1Mc043cGIAAGACItyJ4eFhGhwcbKol9sN0XpwcIzg80RCZE7Np0yb6yEc+Qq2trTRz5kz61Kc+RS+//HJDGmame++9l+bMmUNTp06lj33sY/Tiiy82pDlz5gzdfffddPHFF9O0adPo1ltvpb/85S8NaY4fP05r1qyhcrlM5XKZ1qxZQydOnPBSXAAAAAr1ep1mzJhB69atc928Meh5knIwwlwcTzhGS5cuHXc9cW+HADR4CbZZvnw5P/TQQ/zCCy/w888/z5/85Cd5/vz5/Pbbb1tpqtUqt7a28pYtW3jPnj28atUqnj17Np88edJK09PTw3PnzuVt27bxrl27+IYbbuDFixfzuXPnrDQrVqzghQsXcn9/P/f39/PChQu5s7PTuKwI7AUAgPGIIFfyMH06yHnSMNslDHTXYxoMm6ag2SwQ2+yko0ePMhHxjh07mJl5bGyM29rauFqtWmlGRka4XC5zvV5nZuYTJ05wsVjkzZs3W2kOHjzI+Xyet27dyszMe/fuZSLinTt3WmkGBgaYiPill14yKhtEDAAAXEDMPiqVSlwqlbhSqXhqVO1mF9lNW4664Y47f7/nkxcNDCLo7M7fjAIpNhHzyiuvMBHxnj17mJn51VdfZSLiXbt2NaS79dZb+XOf+xwzM//iF79gIuLBwcGGNIsWLeKvf/3rzMz8/e9/n8vl8rjzlctl/sEPfqAty8jICA8NDVk/Bw4cgIgBAABubEj9NqaqEyF+D6OB9kPUTo98fUEEQtj5qNfbbI4XszcR43t2EjPTV77yFfroRz9KCxcuJCKiI0eOEBHRrFmzGtLOmjXL+u7IkSM0efJkuuiiixzTzJw5c9w5Z86caaVR2bRpkxU/Uy6X6dJLL/V7aQAA0BTU63WaNm0a3XXXXTQ6Okq5XM737CN1NpH4feXKlbYxI1ESxewmOY5n/fr1VCgUaHR0NFDMiyjnAw88EChGx+56497hmyhlM7b8KqUvfvGL3N7ezgcOHLA+++1vf8tExIcOHWpI+4UvfIGXL1/OzMwPP/wwT548eVx+N910E995553MzNzX18dXXnnluDRXXHEFb9q0SVseODEAAHAB1X0hD6vv+qEZHAH1GppxqCYMor7XkTsxd999Nz355JO0fft2mjdvnvV5W1sbEdE4t+To0aOWO9PW1kZnz56l48ePO6Z54403xp33zTffHOfyCFpaWmj69OkNPwAAMBGp1+vU29tLo6OjRESUz+epq6vLyAnw28uWHYGgPXX5+Dh7/aqrEeYsp2ZAzGx78803Y1tPyBUv6mhsbIzXrVvHc+bM4T/96U/a79va2vgb3/iG9dmZM2e0gb2PPvqolebQoUPawN5nn33WSrNz504E9gIAgAuyA+MnDiOMXnbQPOTjm8HhSTNObpP6nTyzLcr7EVlg71133cXlcpmfeeYZPnz4sPVz6tQpK021WuVyucyPP/4479mzh7u6urRTrOfNm8dPP/0079q1i2+88UbtFOtFixbxwMAADwwM8FVXXYUp1gAA4IKfQNKgu02HvS9SWLtfm5Z3IuMkEnXDa2J/rSjrLjIRQ9LYqvzz0EMPWWnGxsZ448aN3NbWxi0tLXzddddZs5cEp0+f5t7eXq5UKjx16lTu7Ozk119/vSHNsWPHePXq1dza2sqtra28evVqPn78uHFZIWIAABOFrq4uLhQK3NXV5auBDtM5yQJZK2+UOAnQpMQedrFmiBgAwMRBHj7yg9uaL3Iar9+lhTQ0zmlG1EmlUmEi4kql4lhHUdYhRAxDxAAAmhvZfZH/74TbsI+XoYWoCbuRhPuix068iN/t6ivK+oSIYYgYAEBz48d9URseXcyDk9sSRzyEXVn9Iq9UHFfZ1fOHGS8UNqKeVedFV8643CyIGIaIAQA0H3LDYeq+2B2v+92NON2MIMv869ylpFwYN+GYNF7qOa6yQ8QwRAwAoLkIaw8ekVcYAiGN6NylOB0klTQ4MW6uSpB8ogAihiFiAADNQ9C1X1S89KjDaLjibLizILTiQtRFqVSyhowEaXOEZCBiGCIGANAchC1gRJ5xDiGkucFsZkS953K5cSImzWIPIoYhYgAA2UbEvIhedFgCxitZc2ImOrrgW6ep82kEIoYhYgAA2aRWq1nCRfSiTRogr0IhSmEB0RINJrPH8vl85l0viBiGiAEAZA956Ej8mM4+8jpkE8UQj7rmSNgN6UQXR7p7pta5iWuX9nqEiGGIGABA9pCnA+fzeSMB43fIIIqGTJTfbbVXL8jlnOixNbqZTkL0ihlYJrOw/NRjnMIHIoYhYgAA2cGrEElrwx5mQ6dzddLoICRZJnHvvcZL+SlznM8ZRAxDxAAA0o/f4Re5QUkyeDPKBjwKVycK4mrc1fVu4l7/Bk5MzEDEAADSjDoUEHQoKAlHJuxzxrWsfZjEVU55qLG9vT1VDlzYQMQwRAwAIN34HQqwIy0rwQahmRtmL+jcNZ0TkwWR5weIGIaIAQCkE7vhH5NGqZkbLubmvz5TZIE7EUUdRAxDxAAA0oVb/IudC6EL4rVzb9KwT08SZHVRNzuE65LUzttJAxHDEDEAgHSgLkJmF/9iJzjUIF7dEvK6tLrf00JUw1BenIu0C7y03rs4gIhhiBgAQPLUajVLvPiNf1EbW+HkiB666rxkIW4iqoBgL05M2kVCWu9dHEDEMEQMACB55Bkl+Xw+lsXf0tg4p3GYK6kypOHa0w5EDEPEAACSQZ0mHOVaHroG0U8j6XRMGI1u2DOxsob8HES1JYN6viD3LGmhBRHDEDEAgPiR135JkxPihFuZ/To7qpjLWr2EgRoPFdfifUHduKTdPIgYhogBAMSL3FBnyXHwOuPJa76iIUy6d58E8nBiLpfLzOq6Sd8riBiGiAEAxIPa2/YiYJJuLKIsQxquLWni3hqgWYCIYYgYAEA8yL1trw6MONbLDsRpJIhgiULsRCmgTOKQ0iTgsriIIkQMQ8QAAKJFXpDM66Jk6pRgEezpFIeQtoZGEDTeJYr4iyhjOuS87RYwTDqmRMakLGkqLzNEDDNDxAAAokNuuP28/HXxIm5OTNoaGkHQmUdZdWK6urpsN/BMk+CEE5NRIGIAAGEjGq5SqRQoWFPk09XVZXxM2hoaQVrLFRWq+5KlIO6sABHDEDEAgHCp1WoNM02CNNxRuyoTTVjEhezAxTFVeqLipf3OEwAAAFeq1ar1/8985jO0b98+6unp0aat1+s0Y8YMmjFjBtXr9XHfr1+/ntrb22n9+vWhlrFer9Nll11GGzZsoP379zeUWXynK49Jnl6PizqvMPNzykf+rlqt0ujoKBUKBerr63N8BrzQ3d1NkyZNou7u7kBljfLY1BKDqEoEODEAgDDx4m7IM5bijGGRZzupZfXj/oS9SF3YDlRY+eny0QXtRuVwyesL+SmrKWmNq1LBcBJDxAAAgtHV1cW5XI5LpZLnRsvL+iBhNoxhbx8QNGg3jDLEkZ8uHydBGDZenrW0TWePAogYhogBAPhHjn8x6bmaNg5OjaWf3nHUjVJWGr2g1Go1LpVKnM/nrWDruOs2aZckbAEcBIgYhogBAPhHHg4y6R2bNkBOwxZ+GoikG76sow4ZmQ7phIF675IWjE7PUtzPGUQMQ8QAALwhNyJ2/zc51vQcYZc5LJJuTONCnWmkOjFxnD+s5yrq8sCJSQCIGACACXarrgq89ELTvPy8wE2gTRR3J+x4H+bw73cU9yKNz6QKRAxDxAAA3Onq6rKGEewCOHUvfbuGQG10wgoODbPhkcsY9vBWWtEFzpo4D2JbCNO6CCJ4/abxShZEKkQMQ8QAAJyRBQwR+Zq1ozYEdk6MncvjRtDjnfI0HSrLOl6DtJkbXRov9e5nGn7cYiIL9xwihiFiAAB6xPRnuWHzGgfhZQq1SO+n4Uibk5OFmB4VUYemQdpyubw6MV7IgphICogYhogBAOgJKmAEcfSkw2rowipr2Nccdn5ehv5AeoGIYYgYAEAjojETmzfm8/lAvfK0N45Bh438CgIv5wqrXIKwnCs/5/aSJkj6pImjvBAxDBEDALiAOp3WT7BmmMvvq2WLolEI6nL4PT5qh0rOXx3WiyKGyO7cQdIESZ80cZQXIoYhYgAAF2JX8vm87+m0UcdHRNUoBBVHfo+Puqcu5y/Hu8j1F1UZJqITk8SyARAxDBEDwERHdl9MBYzXRiqMF3oUs4XS3jCGhdcA6ziIou6TvJ9JOEUQMQwRA8BERp4+ncvljBs5r8MFYb/g/eSna+CyNkThRNYEWRR1n+T9TKL+IWIYIgaAiYq6LojX4F23nr1f5ySKoQjm8PdjSgvqUGBWBFmWnJi0PieRipgdO3ZwZ2cnz549m4mIn3jiiYbv165d2/ACISK+5pprGtKMjIxwb28vz5gxg0ulEt9yyy184MCBhjSDg4P82c9+lqdPn87Tp0/nz372s3z8+HHjckLEADCxUIM6vUyfjtJd0Z0jCEnEKCSBHO8S5tYA4AJpdewiFTFPPfUUb9iwgbds2WIrYlasWMGHDx+2fo4dO9aQpqenh+fOncvbtm3jXbt28Q033MCLFy/mc+fOWWlWrFjBCxcu5P7+fu7v7+eFCxdyZ2encTkhYgCYOASZfSSODzsuxekcQUhrwxMG6n1IW7xLlvE7ZT4JYhtOshMxt912m+0xJ06c4GKxyJs3b7Y+O3jwIOfzed66dSszM+/du5eJiHfu3GmlGRgYYCLil156yahsEDEAND+q++LWY49qSMcLQddGSWvDEwSvU6OD1EEa6s+uDEGfDSeyJH4TFzHlcpkvueQSXrBgAX/hC1/gN954w/r+F7/4BRMRDw4ONhy3aNEi/vrXv87MzN///ve5XC6PO1+5XOYf/OAH2rKMjIzw0NCQ9XPgwAGIGACaHPFiNnVfTF7kIo3fIQy3RsVPY5KlBsgU2WkR4iXM+xjFsWFhV4Yon400iDdTEhUxmzdv5h//+Me8Z88efvLJJ3nx4sX8wQ9+kEdGRpiZ+eGHH+bJkyePy+vmm2/mO+64g5mZ+/r6eMGCBePSLFiwgO+77z5tWTZu3DguFgciBoDmQ278vK7dYurEBFnYzq1RibK3nSXkmBevQ4BwYoIdE5Soh/oSFTEqhw4d4mKxyFu2bGFmexFz00038Z133snM50XMlVdeOS7NFVdcwZs2bdKeB04MAM1PR0eH5x2J/ZBkI5mGRjYqEPOSXdRlC6L8G0yViGE+Lz6q1SozRzecpIKYGACaC3XqtJfGLw2xMKakYbgjbNSYl6j2NgLhI8SmOsrhZf0lr6RKxLz11lvc0tLCP/zhD5n5QmDvo48+aqU5dOiQNrD32WeftdLs3LkTgb0ATFDkIR4i4o6ODk/He4mFiVs8NOuUafk61NilqPY2Us8LvCNES6lUaohXkn9KpVKk9RupiBkeHubdu3fz7t27mYj4/vvv5927d/P+/ft5eHiY77nnHu7v7+fXXnuNt2/fzkuWLOG5c+fyyZMnrTx6enp43rx5/PTTT/OuXbv4xhtv1E6xXrRoEQ8MDPDAwABfddVVmGINwAQjrNlHaXZims150c00ilOoNVt9xoVunSUhPOMe9otUxGzfvl0bQLt27Vo+deoUL1u2jC+55BIuFos8f/58Xrt2Lb/++usNeZw+fZp7e3u5Uqnw1KlTubOzc1yaY8eO8erVq7m1tZVbW1t59erVWOwOgAlCrVbjUqnkOfBTHq7ICs3mHIh7UCqVErmuZqvPKBGuS7FYHCdahBOTRD1i2wGGiAEgq6hDR25bB8iNVlQiJoqGMalZJXJgbZjBy+IzIT693oMgbplTeaKu37DPY+JahXXvdENFaQi0hohhiBgAsoZ4qYq9coiI8/m869YB8vBBVA1KGDEcatmSGPaQzxnk/HKMi9+1XpzK5jWN7vO46jfs86j5hXFt4tmTZ/ep8S75fN4SoEkPxUHEMEQMAFlDHj7ystBcHPEVYcymURueZnBi5EYwaB1NRCfGpNxBr61WqzV0DFTHxe1cSQARwxAxAGQNee2JpF+igjBf6mlpIIKga/Sw1ot/onSL7AJ145hdFBSIGIaIASALyI1iV1cXFwoFo52nxXFeV+yNi6RdlqgQjS52lQ6HKIY/1eG9UqnEuVyOi8Vi6GsrRQVEDEPEAJB25ABerz1RuTFNwxi+ShKNfVS9enUIKuk6bwZHK0x0wiWMNXniiinSARHDEDEApJ0gDX3YTkwYcQ1yWZJo7IM27naxF+p1JB1rkmTjGpSw6k5+5uSZfLo1XfyeE05MwkDEAJA+vAaWeg1g9PvSDdIw2rlCWXMMdKIyyeEju3uS1no1KVdQAdbV1cW5XM6KHxPPXBhbAKSpXiFiGCIGgDShTp82fYl7eel7SRvmrAwvrlCaGgoZu/uTZHnTWld2mDx/Xq9JDZyWXRcRPxb0uVUFaxocLogYhogBIC3oXr5hBBcGESJJvbDT1FAI5PuDjRn9E/ZMNnXNJCGUc7lcKLOL1GcxTaIRIoYhYgBIC+JlGfaut34EQdKzmtLSUOhiKjDjKHl006Kj2i06Lc+iDogYhogBIC1EEXDqN1+/TkhSL/wopuDKjWRUAibO+grbAYn7PutmF6kBunEHUye9dAFEDEPEAJAmkgi6tZtt46ccWR9+UsWLGDaKqrGKs77CPFcS91mcUzecp963qMslypL0NHqIGIaIASAJopgWmwbhkXUnRtRFXDEvWXZiolyBWOd0OJ0zqfsGJyYFQMQAEC9dXV0NQYgyTkNCUTUafl/IclnDmLWU1MyeKIcjvORnmjYtMRp24tdv+eS/C3mzRROB7aee1ec9i88wRAxDxAAQF+JFJ8+kMH3pyVZ6VNa1V0dGTh/EzYl7aEJucKKuVy/XZpo2iaEcHUHdRKdp0fLeRWGLA7uhoCw9wwKIGIaIASAO5F5msVg02vtIdToqlQqXSqVUuAa68rkdG2bgcRDkoQe1ToOiXksYDoHXdGHUZ5DhErfnolarNezELk+Llv9GwrgnuvOLmWYdHR3c3t7OHR0dDb/7iQ2DE5MgEDEARIcQH/ILu1KpOKZXXQK5dxeW+5EESZZX17BGEQQaxjUGdWTCLEPQwFWdYJSdSDFk5LfxdxMPbn9DzNww88zk+DQBEcMQMQBERa1WG/fCdlt8S35pus0aSktshClJlFeISF1sRVjlCfueBO39J+3EyHkUi8VxQ0Tyui6lUinQOdxEhnodXV1d49w3p13hg9Zl1M88RAxDxAAQFbIDY9Lb1AXvZk2opA3RyEW5SF2UvfUs33/VgRSiRR0m8lt/ur8Xu3Sy2xLn1GjZ0YriHkLEMEQMAFEhXuJ2AkbtJeqGN/y8BIP20LPacKqNmi6OKKrzRlVfaR/OEOicG1nEODmQXp5Xt+FWHSKdLKKimBptV94oBRNEDEPEABAmXoYW1LgD3RoXfl6Cdi93ry/9tDecKrLr0p7BmCEdWRCUOqdDFhp+y667f/JnXkV51Avh2T1vUd5DiBiGiAEgDPy8KMUxYnaE3Wwlry9BP05M1mJt7Hq9cSxBn4X6kVHLG1bMjBqom8vluFgscj6fd51557XcstPjNoQUZayQSXnjXAAPIoYhYgAIitwT9bNiaJKOQVy91DCp1Wqcy+VcZ3r5ydfk3mXN4VHL66f8tVrNEixCdMvxLoVCoWH4KKy60T2fbuVP+v7EeX6IGIaIAcAv6swXP8F7psGJYaKLK4hrufYgqA1a2CLGtPEJ0qNPwsXx68QI4ZLL5SzR6BSoG8WzLO6JvEO1W/nDrGM/ecV5jyFiGCIGAD+oS6T7nX2QRK9RPmeWhkZkwRWF8IujLpJ2CUxQxbm6potuUbgoy5LkJotpv18QMQwRA4BXarVaw4vdb2OahAsjzpsV4ZK1WB030nwNduJFOC5xlzmJGBO7MqTxfjFDxDAzRAwApuiGM7wG4/qZIhomcQ6FhHGuoLE6SQ0HRNH4meTp9bzyEvxynIs8fOM1b79lUP+WTB3DrA3vhQlEDEPEAGCC1+Bd0fiqMRteh3LCfskGEU5ej/VzLtUN0DWmclq3uvFThjDEZRQC1SRPkzSyw6GLc7Grb9Nrsnv27bBb9t9U8Hu5ZvW60j5c5AZEDEPEAOCE3+Bduxe5V1ES9ks2rU6MzuVym+kSpPEKq9xx5hFE9IrnWATqivgWUccmcS6m1+RVxDgt+29y7iBiNkonJg6XByKGIWIAsEN2X9wEjJ8Gx+3corccZr5xY9L4iEZMDdo13c05K3XhlaACVo650onDUqkUiaBN6r44xdC4lSmKMsfh8kDEMEQMACqqK+BksQvCfmHZ5ed0njQ26ibDAPJwQprKnjRBXSR1KrpwYsLYhiGNwzDq8+SlbH6uJwlhpAIRwxAxAMh4jX0Rx5i6B055mDguTi9G+SXuZQgn7EBNu2PlehK/2/Wcg5zHT5o0OgsmyMNEsuMiD4GKvbvCiMHS3TP1vrrlHfR+u7mUfvL3+rcgu4fiuLCeZS9AxDBEDAAybps2CuSXntyL89tDDaNnKwswk3xMzxlWr1vkE0UP3qSMJq6Q1++SRq5TWXiL51gWtF7qyE4I6/Iwva9y3kHqM6x8wji3qKekygQRwxAxAMiYBiXKL3PVbfDSEwu7B+fl/FE7MapLYLLvjV8mkhOjOn9yoK6Tk2daR04NsV2+UToxdi5lEuvHROECBQEihiFiwMTGdBjH7Ti/eOnlp60xdUN1CdLoZGQFWSioexSF/Vyk7TmL0glL27V6BSKGIWLAxEVeJ8PpBWnitARxK0xFU5xWdRhxMMJ5CSOQtNkwdUUqlQoXi8Vxa7n4rdOgz2/cIj/svNT8hEAyiX9Lo+CBiGGIGDDxEC8jeZ0MpxeT3BO06xVG2VuU849rBo/X67GLEQJ67OpIFzgq/4hn1m/dBn1+w06XBHLZRH2brAydxmuCiGGIGDCxkB2NUqnkusiWOCYqJyboeaPC7Xzq97qGIaq6SBNBYk/kNOL3Uqk0znUpFoucz+e5q6srcD2YPr9xOTZ+03vJ02StpaAxVUkBEcMQMWDiIA8fuc0+SgNp7Pkx64e2on7Bp7UudOUyLat4HsXwkPi/yeKKURN3fUdxPpFnkrtgRw1EDEPEgIlBrTZ+5+m0k9beodw4pMUZigsTx8LUxZKHM3O5nOUa+JnBlfXg3jDOZ3dvwlqXJulnTwdEDEPEgImBaHjFMJLuheT20ovSSvc79TjqYRxdOcNa9VXkGXbjFWVav/UtP1tCvMgBu6VSyTUPp3x168I4XYfT50GI8hl0G16Nyj2yyzcN4gYihiFiQPPiNabEzX52e0l6CRLU5evnBRxHQK16XWFa9GGU2UseQdN6rW+dyJDdwKANoShDpVLx/NwGrfs4hYSat9u9CRO7fKO8VlMgYhgiBjQnptOnBbIb0tHRYQX8ehFCcoPi5WUaxIlR85HPG9ZLXb6uoNsruJXZz/Fe6i5M18bpezG7SI51CWvBP7tn0m43aN3nar35qRedWI/K9QlS3qhIw3AvRAxDxIDmQrzs5PgXkxeIXS/PS28rDS81mbB6imloMOxIQ2+Y+XwdlUolaxaRcEVyuVzogaVenRXd5+pnXurRj1h3y9+0g5D0ffZK1OWGiGGIGNA81Gq1cWtrOE2ftuvReh2G8oJft8YPXnvXcfZ0w8rfbz5hxeKIGCF1VpHsfgQ5l+5YN7fDz5TiMB0qr8fIf7cijek1p51MOzE7duzgzs5Onj17NhMRP/HEEw3fj42N8caNG3n27Nk8ZcoUvv766/mFF15oSDMyMsK9vb08Y8YMLpVKfMstt/CBAwca0gwODvJnP/tZnj59Ok+fPp0/+9nP8vHjx43LCREDmgF5+EhY924vDr+9pDAaJa9xM1Giij+vDpQfku5ZBzm/OlSk/ritOxR1OZN4roMgyqtuppiGv420E6mIeeqpp3jDhg28ZcsWrYipVqvc2trKW7Zs4T179vCqVat49uzZfPLkSStNT08Pz507l7dt28a7du3iG264gRcvXsznzp2z0qxYsYIXLlzI/f393N/fzwsXLuTOzk7jckLEgCwjGhS5J+wlFkXXc3XDy0s2qhiBMBHXI4u/rDgxcZxfdlxk4SLqy89sLdPz+3EF/dZtFOLBJGYpbleomYhtOEkVMWNjY9zW1sbVatX6bGRkhMvlMtfrdWZmPnHiBBeLRd68ebOV5uDBg5zP53nr1q3MzLx3714mIt65c6eVZmBggImIX3rpJaOyQcSArFKrNa79Yrf6btjj7V5enGnrVYY9tNDMiHrw6ri4DZ34deKifpaiuO+izHHMWJqIJCZiXn31VSYi3rVrV0O6W2+9lT/3uc8xM/MvfvELJiIeHBxsSLNo0SL++te/zszM3//+97lcLo87X7lc5h/84AfasoyMjPDQ0JD1c+DAAYgYkDnUIRCnBkUXWBlF3EtYY/lhNyZybzjMPWLiEj9u+erKofb+vZStq6uLc7lcg7snnBcTx0WtP/nc4js1JkoXQ2MS2xIGQeKKdC6L/HkYs7HcyjCRBXhiIua3v/0tExEfPHiwId3tt9/Oy5YtY2bmhx9+mCdPnjwur5tvvpnvuOMOZmbu6+vjBQsWjEuzYMECvu+++7Rl2bhxo7ZXAREDsoIsTNxWjRWNhprOraH28zIMq1cYZu+yVqs1rAwb5m69ajmj6hW75WtXDi9lkxteub5kl8/vUJF8bru6lZ9n+Zn1W59xuIV2LkvU7otJmdLkzEQprBIXMYcOHWpI94UvfIGXL1/OzPYi5qabbuI777yTmc+LmCuvvHJcmiuuuII3bdqkLQucGJBF5N6pqYARx9n1Fp1eLH5ehlE6On6RG5Sw94vy2wv2en1RODGqEye7emK2UalUCmU9HDtnRcbEifFyLi/DVVE6MXE7ImlzYuyc4LCYkMNJKoiJAVlAvJRFL9lLzzhJQRI3Jg160iTZW1bjXIQz5WVWm1/iuG5xjjim8IML2L0vZFetqZwYEdj7jW98w/rszJkz2sDeRx991Epz6NAhbWDvs88+a6XZuXMnAntB0yAaYRGfUCqVfPdOwxhCSSuinoTIc2soo75ep/zjrmvZ2RD1I56nMJb+91qOMM5jl1fWn+MsoMb8ODlfUd+PSEXM8PAw7969m3fv3s1ExPfffz/v3r2b9+/fz8znp1iXy2V+/PHHec+ePdzV1aWdYj1v3jx++umnedeuXXzjjTdqp1gvWrSIBwYGeGBggK+66ipMsQZNgzwc4tabCRLkZ9pLDvJSCrsRkx0WL/XEbLapnZ/y+hnO8IKXMqnCTh4yKhaLoVynk4gQQ0QdHR2hN2RhukdufzdBhrfczhlmnl7O6/YeEVPndVtsyH9r8hCkPDsyLjEZqYjZvn27NoB27dq1zHxhsbu2tjZuaWnh6667jvfs2dOQx+nTp7m3t5crlQpPnTqVOzs7+fXXX29Ic+zYMV69ejW3trZya2srr169Govdgcwjv+RMh0PUhtmLfW8adxGkcQ5zOEF+kYpyh7GHkFxGP+UVx0Q1nGFSJlW8CAEjOzGVSsVTnqblkD9XV48OU9CpcTxe85bvv9vfjfg9zNiOKPL0cl6n86kiRS2jThyr74Uw/9adwLYDDBED0oXXYRH1WKceZRDCaJyjdGLCIiwnJu5hKrk+5L2zVKdCd3wWnRhxLr/PgNzIwokZn0YMN4q9sNzudZTvHicgYhgiBqQHtXfpZJX7bWy9vlySetm6lSEtAiPJGAw74RJ0t20TgRL2dUT9jCXVyKYNE4EWVQchCiBiGCIGpANZwAjxolq08gtHtni9WLdqWreevR9HKCzUPXr8XrNAd0zQxsxPOYIgxyvIi9HJ06Lla3Ern65RsxvicHIvgiDydRpa8SvU0/AcJ4GdOFHrWJ316PdvKykgYhgiBiSP/HKRg1LVF7f6shdpgzgxdg2T+DzK6ZFuZZRn0kTlxAR9YcfRo5eFi24hunw+73t2iHr98jPm5MSE2dCZODFeBKicNsnnOG50zpxO+Mv1INLJsx6z5FJBxDBEDEgOtZdoOvso7LF53ZL8SVrK4mWrzniIgjS/sMU9sFv+X94CwO91+B1mibvedM++nZDSOTFZGBrxg13HQ8SuuQn/ND//JkDEMEQMSI4oe4leXk6iHEktEKY2UEFFWtIv5qDnV8WtLkYqjGvMUoyQyFcV2n4EWJaR/1bUYTI7wRZXvSRR/xAxDBED4iXKXqIXu18tR9wvH/llLA+PhTFEEeZQRxznV2fxqLOLdJsuhnGNYeThNDwZBLt8nYR20vc9TOzeE+owkUkHKK56SaL+IWIYIgbEQ61WswJUvQRaun0uI79EvMZChInbudUYIPFvWMNkSffITWx70TAVi8VxMS4ms4u8uhA60epW316fOdNjTDDNNywxngZny64TIv6vDhOZdoDCuidenrG4gIhhiBgQD/KLiIg8iwsT0eHlJRLlC8etrOL7MIVL2lHrRH0eiEi7noqp2DStc7lRDJonc/TDR275hiXGw8gnaB7y8WmM6wmrrsMEIoYhYkC0yEMFYnaJU6Bq0uPaXjFxHbq6uhqm/5pci+n1+qmXsHvsJj1U4bhMnjyZu7q6rCnR+Xze8ZnwWg+qKNR9budeOLk7Xp0ftzRxuQNh5RPGMyvcWN3fgs5ZkXf29lL+sB0TU+fOT95BgYhhiBgQHV1dXQ3j1qakscdjh0lZ/SwPH6ZbEMYxTsfa5ScEqRqgG2S5fK/lc7tW+XuntF7zdUuTpWecOVynRr734jNdPJg81OqlfF7uuZdym6SP+75CxDBEDIgOubFSe1NBe7lOhBUnYHIeE9dIdWK8XoPT5356+36G3ZxmTsn5yb1ntdESdeAWz+DXwdDdD9NeudcymdStXB41/7Q6izJe3Q27upefHScnRs0/LicmTOcmrnePACKGIWJA+Ig/3o6ODtuXkJhlIAL1wuy9eOmVhXUeQa1mv+JrVOcMI61bHibXJDtvYjq0WONF3nDRyznVBsDpmoJcb1jPiVp2tf6a0X2Rn3n5b1q34FwaieodEeW7RwARwxAxIBzEeHc+n7dmIdn98dZqtYbdhKMaw47DiVF7dE7uUxTnDCutWx66WBOxiq64h/K1i2fAy9o78rnsGj6/Lk3QY01RRayp05NGTOtTiBXhsMh/23YuS5qI6h0BJyYmIGJAGKhrezj98cq9U6cXnPqSTFtPTrXQxXVlpcft9yUrOy7ytcr2v5+85Z5rGEJUPS6MnrFJWbIwVGSKTqwKRH3KQtXJRQPhAxHDEDEgGOIlJy8Nr5sqqx6jvvR0DYvuJZkmZNGi9rrTWF4VL426LFCcdhoP0oDbHetXfKjHhSEuwhBCWcJJmOvqs5kEXBaAiGGIGOAP1SXRxQGYvOjDsPjjenHqho/CGiKI+hq8NjiqIBMiVUyHtgtUDquRj8KJCXrOMO93kHKGiek165yYpMiiUIqqzBAxDBEDvCOP+8urq6q98aC7LtudW80njt6xGusQJB8/boPdcab1alJH8n0TcQ3q9OhSqWR8fXZlM5l1EpYQNkW9v7rzR/Wcxe3uyGJMDr41GbaNojH2mmcY9RVE4PohqnsMEcMQMcA74g/SdM+SSqXiK1BTxU5IRNkzUx2noGP9di8zt2uwO8705WhSRyIvVbSogtQUXdlqtVqDc2daXqfy+20gdMOa4v7GOVQShwtnF7sl/22aDNtG0Rh7zTOJYcGg1w0nJkIgYoAJ8ovQtEEzEQBuLwenhibItbiVXy17WHE5YQ+N+B3iqNUurJ5aLBat+yoaNDXWxW/Z1M/lBlQNAHa6FqfnxE3w2OUr5+nnvoTdMHkRbrrj5PVYnESLOjzk5f7pXNa4rjvM+jY9h65uoxScXoGIYYgY4I7sgATt8Zp+p4oIvw2NDrsGUSeY0hpULGPaSxSNuTwNVg3a9FPHXnqpav7ysXJj63RvTMujxmiF7eB5uW4/+bnlL+6nuJfqpqLyPVWHj0zK7HT+sK/dJO8oz2l3Dq/3JG4gYhgiBjgjCxi73rnXISETQROliNCVQZ42rIvzCZs4epVqGp1wkZ0Yv/cvyLXIPV2vAaROAlhdqyULAbleHAl1OwcRYyTq0K4eo/hbDROnexXVOZ3O4dcdiwuIGIaIAXrEH6tYtCyfz9v+4XrpnZj27uJ+Yakrjsbd4wsb2XWRe+Beh4pUoiq3n3ydjklbY+MHu2sQ1627l2lzCryS9fLHDUQMQ8QAPfIqnKJhV/EzPp5E784O+YUZdKE2r0R5DlWUifvodZps2K6L13OZHpO2OIWwsBP1SbsVURJn+bNeV8wQMcwMEQMaES9IIV7E7BTdH3oWe01qY5D1lxjzhWsSe1UJ90xdBt7rfcrK/c1KOZ1wE4zNcI1poxnqFCKGIWLABdQevG4WkNz7FQ2k3z2C/AoKr26OLrDR7wwnL73/MGIxnBq3jo4OzuVyDasly0MMpjsmq9ekznIRDo7JFhFurk0Y99kpTiGoMA3iCLk9CyKNuk6O/HdnN8U7LMHtp0795Kd+r9t7y+/1+CmbybMYdR1HAUQMQ8SAC3908hCSXcMrei9hzFaSY0+89Iqc0uq+E5+p62L46YGp1++Uh3xev7099Xp0Q0WiPPKu4X7qU53dIt8Xt2t2qne3z2RM3Aev998Lfo43OUaeqi/Xs3y8/KwEvQ6vZVU/Nz2/Wzq7v5cg1+enbCbHhFXnUd07HRAxDBEz0ZEbRZOZQKIX2dHR4dtlEH/k8vnCdmLU1YK9OBNO5wnTiTHtwcvr8+gEjLz8v584EbGNgLxejK7uonRidMLWrqfstRxeCMOJ0TkA8q7eOiemvb29QYSG0Zv3U/9pc2L8unl+/ubhxGQUiJiJi87GdkqrNjJezuPnBeYV+cUVxgwj0x6V3+txyr9Wq1nDdZMnT7ZmiKnui3pOkzKr5RXHBL2vQdLrhK1dmiD3NEzE86YGTItyykv5J7Eabtj5xdk4C/xeQ9qelaiAiGGImImIKkhM4kNMGhm3Y6N+ociNsV+XR8b0OL/X59TLlKdEy46L2qv1IxDV8ro5RmFdt5toi6unHBby82Y35Kf+vZi4iGFdX9j5JSEMov7bzToQMQwRM9FQF3XT/aGbWLimtrRT2jCuRd5N2a4xdnv5xjUEYZdOvg4hXuQZRh0dHYEcDJN752eoxuQYL+niwsv90gWvy9dTqVS4WCzapvHjlAUtd1THhP130ixCI8nrgIhhiJiJhhxTYfdHZ/KitUsj2+dRYxpc7PaSiauHqXNA2tvbG4aJhBgxHQ4zGebzcj+9Bu+qx2ZhGMipPPKzItLJQbi6471cX5AGz089xnWMn3zT9lz4JcnrgIhhiJiJgNyjnD9/fkNQoZxGvFzVwEO7PHUvY9lJiKJ3opZTdmJMj/PynZ9y2SHqdf78+eOGidTr8FpeJ7fDy1CNn+BdeQl8WQCktcetK58ulsrNiQlyfVk+JqwhnrQ9F36BE5MwEDHNTa02fu0X9Xu1J6/2LLy81EycAb/XIefrFMejK28YvSWnht4pf3GcGCJSf0xefjqHQHcut+sMe0hAnE8ezkprD1sdBtLt8uy0vECYjbvfOgqzwXQSw075p/X+TkQgYhgippmR41/E7Ba1R6kLQFQbazdBonuphdED1DXcJuu8hFUeu3x157fr3csr5uoWpXNbKFAVQHLDa9LQqq6alwbISQzaLcomX7tboHAcDpjuGVJFvRA1bisbu9WdWjdBym2XLkwBoRv6VfP3K3TCcmuAMxAxDBHTjIgXgRxroXspuL10db1sp/MFffHIL1BZYMmNkF3DGUV57PJ1G3JRdxiWG0uvC+2JOhH5eY01UhdWC+oWqPfFy7F+0vgtoypc5N/tREvQ58q0boJcX5jPtk7EhCWa4j5uogIRwxAxzYb8Ii2VSraxLSa9ZbsXZlQ9Md0sHb/bA7jhpXwmjZeYpaI6LeoaIqZ1qgomrzN7xPHyAmomx5gsEGZaH1E7MbLbJW+NIG9a6sW1CqOMYbtLUQ4fpdFRgRPjDYgYhohpJuThI6d9j9Qequ57J8LuZevcCz/r0XjBS/ns0ooGUx0mkq/BCaf7YVI+p/vlp0crjonCHQkLWSDJQld+7k3r34QknIGozgmXo/mAiGGImGZAvNjdgkXFS0zeX0fu7eliX3S9QtMl9916fep53QIrTc5jip9rUTdElEWXWLY/n887rusio7sfTs6LF6vftNet3t+gG1Y6nctPGrU8stASQlcO1FWDdoPi5TkJCz+OiZ98QfaBiGGImKxTq9XGxV/YDR3YjdmLhkG3uqhuM0A/PTndsfJ5g+7D5Ld3aep4VCqVcY6Lbmq0F3T3w6k86ncmwzrie7t8k+qd685rJ9LUaw5TpAQp70Q4N0g3EDEMEZNl5EZQBH+6vdidHBH1WNkpCDo+r2t8woofCNKoOZVBBOKqU6O9ukVezu90PV7rS278TJyYOHFyWZIWLDqScGTUcyddByB9QMQwREwWUYdh1EbVLTDQ9KUY5svTzgUyOc6tDE49Vbe6sEuvTk+32+XZpLymjbEaG+SlnuSp1HYNbpKNoW7YSl0gT1dH8ud+F/Lzks4tvR9XxOl++BH2YaUJk7Q8WxMNiBiGiMkaJnsfyS9a3Us3DnvabmjA62wjk7I6vcTc6kIcL9ZzEa6L7HCJz8RxTmVyqm/TYGg/9SRPpbYrXxz3XWB3/+X7IF+nW5nle+LnWfZ67U7PitcGU372dc+R14DusNKESdznS8u5kwYihiFisoLotbkF74q04kXr1YkJa4hHfWH7yTeMoSK7KcPi/x0dHeMcl3al16wb8nGrP/V4L06M12nvsgiTnZgkho5Ul1A3FOnFgfDrxHi5Z07XEfTvQH4W4MQ037mTBiKGIWKygCwK7IJ3/YoE9ZigvRq5rEHXeFHL4vUa7Y7v6uoaFwztFhSt4laWsHryTo6G27niQh7KEmXJ5XK+hsXCwqROom78kr4voPmBiGGImLQjiwKn2BevL0ydW6LmaZqP3JsUPXA5yNhvr9ikATcpm1gMTcwukgVMsVhsmGFkWla3svg9Tj3WRMhF3RjbuUvi/7LAllcl1s12cypnmNdhklfUImMiOwQgHhIVMRs3bhzXC5w1a5b1/djYGG/cuJFnz57NU6ZM4euvv55feOGFhjxGRka4t7eXZ8yYwaVSiW+55RY+cOCAp3JAxKQb8aLVuRryS1jX0PjN15RarXF6d7s0jOC0H4sujen5dI2o+p2MOgQnYlzshm3Usorf3Rpkv8Mabm5Y0ja9PLVcdlbkepFnbzmV2U00hPFM2l2Hyb1x+9xrGgCiJnER88EPfpAPHz5s/Rw9etT6vlqtcmtrK2/ZsoX37NnDq1at4tmzZ/PJkyetND09PTx37lzetm0b79q1i2+44QZevHgxnzt3zrgcEDHpQzfTxClWwG04KOzeuzhWFgjCeTE5l18RI2MnNoRjJZwVeQaXugWA07WpcQx+NsHUfW7iLiXVQKrPlFynQvypIlK+l0Ebfjt3MChu4slPeq95AhAFiYuYxYsXa78bGxvjtrY2rlar1mcjIyNcLpe5Xq8zM/OJEye4WCzy5s2brTQHDx7kfD7PW7duNS4HREz6UDft02H3Eo0izkXNW264/ATehtFIy3nYLf8vuwJe1/dwEx8m16QTmk75xile1N2t5bLJdSevUOzFAfNLFENKfvedcrpGODEgDSQuYkqlEs+ePZsvu+wyXrVqFb/66qvMzPzqq68yEfGuXbsajrn11lv5c5/7HDMz/+IXv2Ai4sHBwYY0ixYt4q9//evG5YCISQdyg9fR0aGN03CaMaNi4uZ4RTRyYSz2FrRMcgOrOjDymi6qeJAdA6cyBBlqYG6cCm/qssTRuxfnV3c41wkWHWGXMWg9OxFWWeG6gLSSqIh56qmn+L//+7/5j3/8I2/bto2vv/56njVrFr/11lv829/+lomIDx482HDM7bffzsuWLWNm5ocffpgnT548Lt+bb76Z77jjDtvzjoyM8NDQkPVz4MABiJgUIA/NtEsxB/L/ZavdrQEWeTm5OW6oTkKYdn+QhqFWG7/VgjqzyM7xkIeF7MrgV2jIx8nBrn4dgKDoxK8ovxrHYno/wi6j3XnDEA5hlRWuC0grqZqd9Pbbb/OsWbP4m9/8piViDh061JDmC1/4Ai9fvpyZ7UXMTTfdxHfeeafteXQBxRAxySEck2KxaBtb4tQY6V7y4ju1YTdBPq+cjyoG4mwYnGI17Fwhkx6+XRpToaIiH6cbrokSp2FEO/Gb1FCWW7mdPgcAXCBVIob5vADp6emJdDgJTkx6kB0TsciaPA7vFG/iZyjECTWuRBwfdO8aP7Ek4nN1uEiOw3Ea9vAbCyHqQB2Kciq3uoNymPfE6VghksSO2fLicupu216X64/CEfLyHDnFoDjt8m1aFvXvrBncmqTPD5IhVSJmZGSE586dy//yL/9iBfZ+4xvfsL4/c+aMNrD30UcftdIcOnQIgb0ZQB3a0DkmqgsS9UtKHp4Jc4qr27CA/L2dC6SKGNOZI3ZDX3Z1KQ8BOeWvO0eUs1l0daQOp1UqlYaAcLcyOZUnSFmdzmWap+78an37HdZU8w7rWsOus6ydHyRDoiLmnnvu4WeeeYb//Oc/886dO7mzs5NbW1t53759zHx+inW5XObHH3+c9+zZw11dXdop1vPmzeOnn36ad+3axTfeeCOmWKeYWq3W4CzopqXqevhuw0d+yyL3aMVwVi6XC1UoufV8xbL/hULBcoGEgyAWqSsWi56W1HdzYuzqUjgxIqBaV15d/ZluKWAqQtW08tCUEHHFYrHBiVGv28/2BX7KaoJXVxFOTPbOD5IhUREj1n0pFos8Z84c/vSnP80vvvii9b1Y7K6trY1bWlr4uuuu4z179jTkcfr0ae7t7eVKpcJTp07lzs5Ofv311z2VAyImHmq18VsH6F6kdg2s3yESt7LYra4aFeL6dENF4sduuXrTeBQ33Boyk/sRNnJDL7tN8r1qb9cvJpgldPcQLgIA/kjVcFJSQMREjzzdVjTSdjNp3HpUIq065GPSE1OHseT9bsJyEZwQTpRuTRe1blR3QScodO6LnzKq9e82CyzMXq9uaFEWlOr9DlKeNPTWdYIlDeUCIItAxDBETJSI3rXOfRHfm1jkamOtiwdw6s2K44XzIa+86lR2tSF1cmvsGld5WEg3LVo3TCTXi3yt6owfuYE3CV41GYqyE4kyXp0Du2EeWdzaBSybOHBO4k4ug9e9suxEYpChsSgFC8QQmGhAxDBETJSIxsUteFc0KnaNkRrM6rU3LvIV7ofJUITqTujK4XQtzPp9i0yW/5eFhCwodPXjFOipDr2YiA8vjpZpI24XcKtbS8aujCYi1SkPO3HmlK/8nd3/3fAq+IIQ57kASAMQMQwRExVu03XlGAg7J0YWH36W9nfL3+14J6Hk1MNWh89M9y2S604NWlWvyc2pqNUuLIgnnKcw4olMkB0oUQY5MFd2YtS1ZEzcIpPzm4pcp8/VhQ7t7r3XskQFnBgw0YCIYYiYsBEvUrnx0mHXszUZOjJB5BlVz1TniohGTx028tqomLoDAl3jJTsPJtOy/eJ0bvET9q7McQBXA4D0AxHDEDFhIRpxufcvuwlqsKqul6vGgIh8nXqXTnkFXajOCXlqtDyjhogapmsLl8FLfIrbdFw1dsRN2Ph1EdTz6upVFkuqY+TFfYqKoNcb1rAaACB8IGIYIiYMZMdEbdDkRs4tpkDXIHo5r1PeYSEaaDXWRRZwupgbtVxudeBUftVl8tKYeq0fNR5Ire8gblkcRPk8wK0BIFkgYhgiJihy/IcudkUXsxHUJVBdH+H8OOURpNesxrjIP3LMj5OLojoZdi6FWl92s4505/HqHujcIXk2lewwyc6aLtbIq+sVxB2yiwNyO4/dZ2G5NXBmAIgXiBiGiAlCrVZraND9vLz99GZlJ8I07sXpPG5DOzrxYtdgu12P6qLYlUV1Pkzqx2tdqs6Xrl7FtHJdXI18Pr/n9nJskLqxO7fud7/AmQEgXiBiGCLGDzonRKx1YjejQ5eHSW9a5zrIn7mtj+LkbKiNohrHYzdFWp1NYzfjxs6FcnItREMoXK2gM6rs6lu4J/LqwDonxqnMaj37XVrf1FkSdSOEFZbeB2BiAxHDEDFe0cW/iMZEbvTVHrZdg+TUa9XFvJjg1MvXfWe3/L88pGLX4MnuhXo+8X8vK+CGGWMiiwNdTIssCJwaXlOHwSmdWx7q97p6gEgAAMhAxDBEjBfU2BARhyL3kOV4D7mnrGuknFwD2SXwuk6MW9yH+F3MMFJ/5P2UVHTXIbtCdjOl7OrBLs+gDoMqAGX3Ql1Ez0u+ftO55WEnDrM4PRsAEA8QMQwRY4JoYNT1T+xsfV1D7aURk92NoA2t3bCFU7CuGFbS7RrttECbLm7EZOjCRCQ4CTPxvSxORFmcgq2jFCV+0wY5JshxAIDsARHDEDFuyI2jXcyIOkSga+id8pfzsttDR02vfm43XCELIp3bIjtKOjGi5uvkpMjix61cXqjVauPW37EbevEyTCTn76VOvV5XGHVgSpznAgAkC0QMQ8TYoTboOkfErfHzEm/hNAtGPp9dvIgaHCoP5ag7RrttRSAf6xQTo6ZTg2TV/P26BDohZrcWj5fpzk4OknqtTtfuxUXS/R4mcGIAmDhAxDBEjA55qMVux2eTBs6pcdTNLlKFiIrcWOu+t1uETv5c3Y/IrufuJpjsHBu7zRqdzqWrG93sIHkrBychYepGeBGbdvn6cT7glgAAwgAihiFidKgBocz2gZeiIdI1THbxKOosJoFb42YXCyIPQelEjC7eRhUCdm6STjDpGn8vsS9OLpDqutjVp1N9hRHzYpo+zngXAACQgYhhiBgZ0bh0dHSM24FajlnROSZuDZPcONvt6uwlkLRWq40TLSLfYrFo/V+3u7PT7BwRzzN//nzb4F4/jb/sPInzi2EueU0Y3TotbnmaihWnoS0IFABA1oCIYYgY5gsNtxAFomevBtvaBZaaIAen+hlKkEVQPp+3XdtFdS50ZVU/lxteNUC2vX18cK+Xa25vb1w/RxZQ6jVEtTqv7Cr5CVj2ct1+6woAALwCEcMTW8TIDos6/MJ8QcTk8/lxC9rZxaXoYjrU77yu+aKuDqz+6BwX5vEBr3bOg9zwCkEnx86YuAu6oSJ5vRmnGCDdujom9eI1LZwYAEAzARHDE1fEqM6ILGLkhk4derFrpHSug0kQq11eonFXZxapTozqpMjIU71lQWNXdj+xIXbX7XVBOa91BAAAEx2IGJ6YIkYdNpGdFiEO5IBep2EO1SnRrfPiFjwr8tXFuajiRazBIucpD3vJqHE8sjPiVbTI16SblZTL5RpmD7nt6eT1PtmJQYgcAMBEBSKGJ56IkRtG3WquOkHg1FCKRtzJdbCLk1AdDTUeRXZidFsAOJVZ97nqmpjGbcjX6ORWuYkjr/EiTu6R3zwBAKBZgIjhiSVi5PVf1IZRnplkNytHRo6zKJVKDS6JLo0afyHPwikWi+McFxGTonM01BWB1fzU9VWcZkHZxYLI1yR+FzOe1ONM6sntGKftDLzErcCZAQBMFCBieGKIGN2qtXZDO3azctTGUf7Ozg1QnQS7DRd15XJyGNR1bOTgXDUWR5ePLK50sSuqu6TLx9QBsXNT1OPlcrjl5XROODMAgIkCRAw3v4ip1WrjhEKxWBw3e8guiFc4HGq8h90S/+rnQjw5xbqINPPnz2/Ix26tFHXzRiGudMcIASfPYBJlUQOb5XyKxaKjI2R3Pl39203zdnNidHm5uSxwYgAAEwWIGG5eEaObNSMaXDW+Qz1GJ2x0M2/snBc1nkUVDGKNFFkY6f5VXRX1+pyma6sBx3K+us0U7bYLcHI25Ot1cj4gLAAAIHwgYrg5RYxTsKlo3GWHRPe9PAtHXlFWdiXUdUeEcFFjXNxm66jOizrM42WTQXU9G6fy62JT1LgVE6HkZd0bAAAA4QARw80lYtQpyk4CQHVqhNiRG3/dbBu5kVfjUNThISGS3HaDZm50PFRxIL4T5dHlp5ZFng7tFAgrf26Xp9MwkN19gPMCAADRAhHDzSNi7KYoy+u3ODXeIp3dvkYCOVBVF6ibz+ctJ0bM5JHLJcqgDknpxJH6nchXFkmysNFtXGlXfqehMF2Z3I43OQcAAIDwgIjh5hAxtVrj6rvq3kKFQkG7NooubsYtBsVu+X95+rHq7qhDOm6r2apDQ2JoRx0eUheeMxne0eUtX5/b0BCcGAAASAcQMZxtEaMTFqIRlmfY6NY6Uff20TXs4hxOM4vU7QjkoR6nYR9xnDx7SI1D0a2MS0Tj1ofxE5MCtwQAALINRAxnV8SowzTyMJA6hVgIE90MIzs3RIgcNUhX/MyfP18bt6IO86hDWaoTogowOR95RWGnoR2TbQScAngBAABkD4gYzqaIUddJEfEnuinEqqshB9uqzsn8+fNtHRdZVIj/q2vH6Bahk10cMRNKTit/L1wj3V5MOjGjDocFjWUBAACQHSBiOFsiRo43kYdW5IZcnpGkGzJSnQenoSKi8wvjqSJCt3aMLiZFdX9068KoeapCRL5eeVhJF9sTdFYRAACA7AARw9kRMTr3Rawmq3MoZDdDHTISs4pUMSIP7XR0dIzbP0jNv1KpjNtrSV7vxcn90QXu2k2JtnNi7IBgAQCA5sdL+51jZqYm5OTJk1Qul2loaIimT5+edHFsmTRpEo2OjhIRUT6fpw9/+MO0e/duGh0dpUqlYqXr6+ujarVK+/fvJyKiXC5HU6dOpSlTptDy5cvpZz/7GQ0ODmrPUalUqK+vj371q1/RY489RmNjY8TMVCgUaN68ebR//34qFAr0wAMPEBHRhg0b6MSJEzQ2NkaVSoWGhoasMor8Wltbaf369UREVK1WaenSpdTf30/r16+3ytne3k779u2jyy67rOF3AAAAwA5P7XfUiiop0uzEyM5GsVjULuOfy+UaPpf3NpLXidH9iGPF7s/0/4edVIemo6OjwQnK5/MNLov8u906M/KsJJKGnxBsCwAAwA8YTuL0ihjd8JEsSIRQ0IkUeYaS3Y8YyhGCQl7bRUx7FuJIneIs0snBt24CxG59FgAAAMAPXtrvvB+rB/ijXq/Tj370o3GfFwoFKpVKlM/n6QMf+AANDQ3R2NgY5XI5qlQq1NHRQUREg4ODdOrUKW3euVyOiIh+9rOf0fDwMFUqFbrkkktocHCQcrkcLV++nK677jrK5XLEzJTP52np0qU0PDxMpVLJ+iEiWrVqFbW3t9Py5cupWq3S+vXrqaenx7qGyy67jOr1OhER9fT00L59++iRRx6hffv2WekAAACAqEFMTAzU6/WGOBMiomKxSOfOnbMExXvf+95xMS25XM6KmRHH2SHiZ+RYFjm/QqFA5XLZ+kzEtTjFruhiWRDfAgAAIEq8tN9wYiKmXq9Tb28vDQ4O0tjYGBUKBarVajRnzhwS+lEIlEKh0HAsM9O77747TsDI7ow4bvny5ZaDQ0Q0MjJifZ/L5Wh0dJSOHz9OpVKJKpUKLV++3HJsRIDu+vXrqb293fZ3u88AAACAJIATEyHd3d3jho9KpRLddttt9L//+7+2Q0NulEolOn36tCWCisXiOLcml8vRRRddRH19fURE1NvbS6Ojo5gxBAAAINV4ab8nxVSmCYUYPtJNeT516pQ2LsYLIyMjJGvPd999d1waZqbW1taGGBUR30JE1lRoOCoAAACyCpyYkKnX67Ru3TrXGBa/iPVhVMR6Mf39/Q1rtiDQFgAAQJaAE5MQuuGjsBDDQ0TnZylhGAgAAMBEJ/WBvQ8++CBdfvnlNGXKFLr66qvp17/+ddJFGkd3dzfl8/nQBEwul7MCcLu6uqi9vZ0efPBBOnbsGPX19SGwFgAAAKCUDyc9+uijtGbNGnrwwQfpb//2b+k73/kOfe9736O9e/fS/PnzHY+NazgpLPcln8/TqlWrMAwEAABgQuOl/U61iLnmmmvowx/+MNVqNeuz97///fSpT32KNm3a5Hhs1CKmXq/Tl770JW1QrRO6mJYpU6ZQX18fhAsAAIAJT1PExJw9e5aee+65ccMmy5Yto/7+/oRKdZ5KpULHjx83Tl8oFKilpQViBQAAAAiR1IqYt956i0ZHR2nWrFkNn8+aNYuOHDkyLv2ZM2fozJkz1u8nT56MrGxuAqZYLFKxWIRoAQAAACIk9YG9Yk8gATOP+4yIaNOmTVQul62fSy+9NLIyiVlCMqVSiWq1GjEznT17lt555x06duwYBAwAAAAQEal1Yi6++GIqFArjXJejR4+Oc2eIiL761a/SV77yFev3kydPRiZkdIvYAQAAACBeUuvETJ48ma6++mratm1bw+fbtm2jpUuXjkvf0tJC06dPb/gBAAAAQPOSWieGiOgrX/kKrVmzhjo6OmjJkiX0H//xH/T6669jiAYAAAAA6RYxq1atomPHjtG//uu/0uHDh2nhwoX01FNPUXt7e9JFAwAAAEDCpHqdmCCkYRdrAAAAAHjDS/ud2pgYAAAAAAAnIGIAAAAAkEkgYgAAAACQSSBiAAAAAJBJIGIAAAAAkEkgYgAAAACQSSBiAAAAAJBJIGIAAAAAkEkgYgAAAACQSVK97UAQxELEJ0+eTLgkAAAAADBFtNsmGwo0rYgZHh4mIqJLL7004ZIAAAAAwCvDw8NULpcd0zTt3kljY2N06NAham1tpVwuF2reJ0+epEsvvZQOHDiAfZlcQF2Zg7oyB3VlDurKG6gvc6KqK2am4eFhmjNnDuXzzlEvTevE5PN5mjdvXqTnmD59Oh5yQ1BX5qCuzEFdmYO68gbqy5wo6srNgREgsBcAAAAAmQQiBgAAAACZBCLGBy0tLbRx40ZqaWlJuiipB3VlDurKHNSVOagrb6C+zElDXTVtYC8AAAAAmhs4MQAAAADIJBAxAAAAAMgkEDEAAAAAyCQQMQAAAADIJBAxHnnwwQfp8ssvpylTptDVV19Nv/71r5MuUuzce++9lMvlGn7a2tqs75mZ7r33XpozZw5NnTqVPvaxj9GLL77YkMeZM2fo7rvvposvvpimTZtGt956K/3lL3+J+1JC51e/+hXdcsstNGfOHMrlcvQ///M/Dd+HVTfHjx+nNWvWULlcpnK5TGvWrKETJ05EfHXh4lZXn//858c9Z9dee21DmolSV5s2baKPfOQj1NraSjNnzqRPfepT9PLLLzekwbN1HpO6wrN1nlqtRosWLbIWq1uyZAn99Kc/tb7PxDPFwJjNmzdzsVjk7373u7x3717+h3/4B542bRrv378/6aLFysaNG/mDH/wgHz582Po5evSo9X21WuXW1lbesmUL79mzh1etWsWzZ8/mkydPWml6enp47ty5vG3bNt61axffcMMNvHjxYj537lwSlxQaTz31FG/YsIG3bNnCRMRPPPFEw/dh1c2KFSt44cKF3N/fz/39/bxw4ULu7OyM6zJDwa2u1q5dyytWrGh4zo4dO9aQZqLU1fLly/mhhx7iF154gZ9//nn+5Cc/yfPnz+e3337bSoNn6zwmdYVn6zxPPvkk/+QnP+GXX36ZX375Zf7a177GxWKRX3jhBWbOxjMFEeOBv/mbv+Genp6Gz/76r/+a169fn1CJkmHjxo28ePFi7XdjY2Pc1tbG1WrV+mxkZITL5TLX63VmZj5x4gQXi0XevHmzlebgwYOcz+d569atkZY9TtSGOay62bt3LxMR79y500ozMDDARMQvvfRSxFcVDXYi5rbbbrM9ZqLWFTPz0aNHmYh4x44dzIxnywm1rpjxbDlx0UUX8fe+973MPFMYTjLk7Nmz9Nxzz9GyZcsaPl+2bBn19/cnVKrkeOWVV2jOnDl0+eWX02c+8xn685//TEREr732Gh05cqShnlpaWuj666+36um5556jd999tyHNnDlzaOHChU1dl2HVzcDAAJXLZbrmmmusNNdeey2Vy+Wmq79nnnmGZs6cSVdeeSXdfvvtdPToUeu7iVxXQ0NDRERUqVSICM+WE2pdCfBsNTI6OkqbN2+md955h5YsWZKZZwoixpC33nqLRkdHadasWQ2fz5o1i44cOZJQqZLhmmuuof/6r/+in/3sZ/Td736Xjhw5QkuXLqVjx45ZdeFUT0eOHKHJkyfTRRddZJumGQmrbo4cOUIzZ84cl//MmTObqv4+8YlP0MMPP0y//OUv6Zvf/Cb9/ve/pxtvvJHOnDlDRBO3rpiZvvKVr9BHP/pRWrhwIRHh2bJDV1dEeLZk9uzZQ+95z3uopaWFenp66IknnqAPfOADmXmmmnYX66jI5XINvzPzuM+anU984hPW/6+66ipasmQJ/dVf/RX98Ic/tILj/NTTRKnLMOpGl77Z6m/VqlXW/xcuXEgdHR3U3t5OP/nJT+jTn/607XHNXle9vb30xz/+kX7zm9+M+w7PViN2dYVn6wLve9/76Pnnn6cTJ07Qli1baO3atbRjxw7r+7Q/U3BiDLn44oupUCiMU45Hjx4dp1QnGtOmTaOrrrqKXnnlFWuWklM9tbW10dmzZ+n48eO2aZqRsOqmra2N3njjjXH5v/nmm01df7Nnz6b29nZ65ZVXiGhi1tXdd99NTz75JG3fvp3mzZtnfY5nazx2daVjIj9bkydPpiuuuII6Ojpo06ZNtHjxYvrWt76VmWcKIsaQyZMn09VXX03btm1r+Hzbtm20dOnShEqVDs6cOUP/93//R7Nnz6bLL7+c2traGurp7NmztGPHDquerr76aioWiw1pDh8+TC+88EJT12VYdbNkyRIaGhqi3/3ud1aaZ599loaGhpq6/o4dO0YHDhyg2bNnE9HEqitmpt7eXnr88cfpl7/8JV1++eUN3+PZuoBbXemYyM+WCjPTmTNnsvNMBQ4NnkCIKdbf//73ee/evfzlL3+Zp02bxvv27Uu6aLFyzz338DPPPMN//vOfeefOndzZ2cmtra1WPVSrVS6Xy/z444/znj17uKurSzstb968efz000/zrl27+MYbb2yKKdbDw8O8e/du3r17NxMR33///bx7925rGn5YdbNixQpetGgRDwwM8MDAAF911VWZmtrJ7FxXw8PDfM8993B/fz+/9tprvH37dl6yZAnPnTt3QtbVXXfdxeVymZ955pmGacGnTp2y0uDZOo9bXeHZusBXv/pV/tWvfsWvvfYa//GPf+Svfe1rnM/n+ec//zkzZ+OZgojxyLe//W1ub2/nyZMn84c//OGGaXsTBbFWQLFY5Dlz5vCnP/1pfvHFF63vx8bGeOPGjdzW1sYtLS183XXX8Z49exryOH36NPf29nKlUuGpU6dyZ2cnv/7663FfSuhs376diWjcz9q1a5k5vLo5duwYr169mltbW7m1tZVXr17Nx48fj+kqw8Gprk6dOsXLli3jSy65hIvFIs+fP5/Xrl07rh4mSl3p6omI+KGHHrLS4Nk6j1td4dm6wN/93d9Z7dkll1zCH//4xy0Bw5yNZyrHzBzczwEAAAAAiBfExAAAAAAgk0DEAAAAACCTQMQAAAAAIJNAxAAAAAAgk0DEAAAAACCTQMQAAAAAIJNAxAAAAAAgk0DEAAAAACCTQMQAAAAAIJNAxAAAAAAgk0DEAAAAACCTQMQAAAAAIJP8P9S8wO2oF3l/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "x, y = [], []\n",
    "for i in range(1, 3001):\n",
    "    if not i % 200: print(i)\n",
    "    x.append(i)\n",
    "    y.append(ant.Totient(i))\n",
    "plt.scatter(x, y, c='k', s=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b310ade1-8e08-40b3-8b14-8976c7f7599a",
   "metadata": {},
   "source": [
    "## Problem 2.5\n",
    "\n",
    "Show the Dirichlet product of the Mobius function and $\\nu(n)$ (the number of unique prime factors of $n$) is $1$ or $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "85fe4f76-8389-47cf-a5f4-2a560df5763d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 0\n",
      "2 1\n",
      "3 1\n",
      "4 0\n",
      "5 1\n",
      "6 0\n",
      "7 1\n",
      "8 0\n",
      "9 0\n",
      "10 0\n"
     ]
    }
   ],
   "source": [
    "# empirical check establishes the pattern: n = prime gives 1, 0 otherwise\n",
    "for n in range(1, 100000, 1):\n",
    "    if n < 11:                                                           print(n, ant.Dirichlet(ant.Mobius, ant.Nu, n))\n",
    "    elif ant.prime(n) and not ant.Dirichlet(ant.Mobius, ant.Nu, n) == 1: print(\"Prime anomaly at\", n)\n",
    "    elif not ant.prime(n) and ant.Dirichlet(ant.Mobius, ant.Nu, n):      print(\"Composite anomaly at\", n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "3254248c-1657-4af3-9add-4365a258e383",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5217 zeros, 783 ones, 0 something-elses\n"
     ]
    }
   ],
   "source": [
    "total0, total1, totalOther, n = 0, 0, 0, 6000\n",
    "for i in range(1, n+1): \n",
    "    product     = ant.Dirichlet(ant.Mobius, ant.Nu, i)\n",
    "    if product == 0:   total0 += 1\n",
    "    elif product == 1: total1 += 1\n",
    "    else:              totalOther += 1\n",
    "        \n",
    "print(total0, 'zeros,', total1, 'ones,', totalOther, 'something-elses')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de1225c4-f755-4c71-880e-6269248a693d",
   "metadata": {},
   "source": [
    "#### 2.9 Extending the totient $\\varphi$ for $x \\in \\mathbb{R}$\n",
    "\n",
    "This problem involves two arguments to the extended totient: A real $x$\n",
    "and our usual integer $n$. So let's try and establish the pattern of\n",
    "behavior by observation.\n",
    "\n",
    "\n",
    "Problem: Suppose $x \\in R$, $x \\ge 1$: Let $\\varphi(x, n)$ denote the number of positive integers $\\le x$ that are relatively prime to $n$.\n",
    "(Notice no qualification of $n$ in relation to $x$.)\n",
    "\n",
    "\n",
    "a) Show $\\varphi(x, n) = \\sum_{d|n} \\mu(d) \\lfloor \\frac{x}{d} \\rfloor$.\n",
    "\n",
    "\n",
    "b) Show $\\sum_{d|n} \\varphi(\\frac{x}{d}, \\frac{n}{d}) = \\lfloor x \\rfloor$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "6fffa191-e2b2-4b10-8708-1bf6834219fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Part a) test passed\n",
      "\n",
      "when part b) hits a case of x/d < 1: Do a diagnostic print with these columns:\n",
      "n, x, [x], d, x/d, n//d, x-totient(x/d, n/d)\n",
      "\n",
      "6 1.2 1 2 0.6 3 0\n",
      "6 1.2 1 3 0.4 2 0\n",
      "6 1.2 1 6 0.2 1 0\n",
      "6 2.2 2 3 0.733 2 0\n",
      "6 2.2 2 6 0.367 1 0\n",
      "6 4.2 4 6 0.7 1 0\n",
      "16 1.2 1 2 0.6 8 0\n",
      "16 1.2 1 4 0.3 4 0\n",
      "16 1.2 1 8 0.15 2 0\n",
      "16 1.2 1 16 0.075 1 0\n",
      "16 2.2 2 4 0.55 4 0\n",
      "16 2.2 2 8 0.275 2 0\n",
      "16 2.2 2 16 0.138 1 0\n",
      "16 4.2 4 8 0.525 2 0\n",
      "16 4.2 4 16 0.263 1 0\n",
      "16 6.2 6 8 0.775 2 0\n",
      "16 6.2 6 16 0.388 1 0\n",
      "21 1.2 1 3 0.4 7 0\n",
      "21 1.2 1 7 0.171 3 0\n",
      "21 1.2 1 21 0.057 1 0\n",
      "21 2.2 2 3 0.733 7 0\n",
      "21 2.2 2 7 0.314 3 0\n",
      "21 2.2 2 21 0.105 1 0\n",
      "21 4.2 4 7 0.6 3 0\n",
      "21 4.2 4 21 0.2 1 0\n",
      "21 6.2 6 7 0.886 3 0\n",
      "21 6.2 6 21 0.295 1 0\n",
      "21 17.2 17 21 0.819 1 0\n",
      "\n",
      "Part b) test passed\n",
      "\n"
     ]
    }
   ],
   "source": [
    "test_x_values = [1.2, 2.2, 4.2, 6.2, 17.2]\n",
    "test_n_values = [1, 6, 16, 21]\n",
    "\n",
    "# print(\"\\nx, n, ex-tot\\n\")\n",
    "\n",
    "# part a) Extended totient\n",
    "part_a_exception = False\n",
    "for n in test_n_values:\n",
    "    for x in test_x_values:\n",
    "        extended_totient = ant.ExtendedTotient(x, n)\n",
    "        to_prove_sum = 0\n",
    "        for d in ant.divisors(n):\n",
    "            to_prove_sum += ant.Mobius(d) * ant.floor(x/d)\n",
    "        # print(x, n, extended_totient)\n",
    "        if not extended_totient == to_prove_sum:\n",
    "            part_a_exception = True\n",
    "            print(\"Exception, part a:\", n, x, extended_totient, to_prove_sum)\n",
    "            \n",
    "if not part_a_exception: print(\"\\nPart a) test passed\\n\")\n",
    "else: print(\"\\nPart a) test failed\\n\")\n",
    "                  \n",
    "# part b)\n",
    "print('when part b) hits a case of x/d < 1: Do a diagnostic print with these columns:')\n",
    "print('n, x, [x], d, x/d, n//d, x-totient(x/d, n/d)\\n')\n",
    "\n",
    "part_b_exception = False\n",
    "for n in test_n_values:\n",
    "    for x in test_x_values:\n",
    "        extended_totient_sum = 0\n",
    "        for d in ant.divisors(n):\n",
    "            extended_totient_sum += ant.ExtendedTotient(x/d, n/d)\n",
    "            if (x/d < 1): print(n, x, ant.floor(x), d, round(x/d, 3), n//d, ant.ExtendedTotient(x/d, n/d))\n",
    "        if not extended_totient_sum == ant.floor(x):\n",
    "            part_b_exception = True\n",
    "            print(\"Exception, part b:\", n, x, extended_totient_sum, ant.floor(x))\n",
    "        \n",
    "if not part_b_exception: print(\"\\nPart b) test passed\\n\")\n",
    "else: print(\"\\nPart b) test failed\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc1b39f9-c171-4780-9f7c-e2e9ec24b0fb",
   "metadata": {},
   "source": [
    "### Problem 2.9 continued\n",
    "\n",
    "\n",
    "From looking into part b) by hand: The sum of extended totients works as advertised; but the \n",
    "extended totient is passed real argument $x/d < 1$, a consequence of allowing $x < n$. \n",
    "Example: $x = 17$, $n = 21$: $d = \\{ 1, 3, 7, 21 \\}$ and extended totient argument pairs\n",
    "are $(17, 21), \\ (5.67, 7), \\ (2.43, 3), \\ (.81, 1)$.\n",
    "The extended totient returns zero in such cases; and this does not alter the\n",
    "sum. So define the extended totient to be zero for $x < 1$ and everything is fine.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d649c7b7-a321-4851-b908-7486848c12f7",
   "metadata": {},
   "source": [
    "## Problem 2.12\n",
    "\n",
    "\n",
    "If $f(n) > 0$ for all $n$ and if $a(n)$ is real with $a(1) \\ne 0$, prove that\n",
    "\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "g(n) =  \\prod_{d|n}f(d)^{a(\\frac{n}{d})} \\;\\; \\textrm{ if and only if } \\;\\; f(n)=\\prod_{d|n}g(d)^{a^{-1}(\\frac{n}{d})},\n",
    "\\end{align}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "475d0a13-2e13-4f2e-a08f-8cbe3bb0f265",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 3, 5, 6, 10, 15, 30]\n"
     ]
    }
   ],
   "source": [
    "print(ant.divisors(30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "b519ec09-4610-454d-8a60-d63aa7ac95c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "1\n",
      "2\n",
      "2\n",
      "8\n",
      "8\n",
      "16\n",
      "640\n"
     ]
    }
   ],
   "source": [
    "print(ant.Totient(1))\n",
    "print(ant.Totient(2))\n",
    "print(ant.Totient(3))\n",
    "print(ant.Totient(4))\n",
    "print(ant.Totient(15))\n",
    "print(ant.Totient(30))\n",
    "print(ant.Totient(60))\n",
    "print(ant.Totient(2400))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "988f2407-eff0-4411-b963-a3cd9518d374",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(GeneralizedTotient(1, 1))\n",
    "print(GeneralizedTotient(1, 2))\n",
    "print(GeneralizedTotient(1, 3))\n",
    "print(GeneralizedTotient(1, 4))\n",
    "print(GeneralizedTotient(1, 15))\n",
    "print(GeneralizedTotient(1, 30))\n",
    "print(GeneralizedTotient(1, 60))\n",
    "print(GeneralizedTotient(1, 2400))\n",
    "\n",
    "print(GeneralizedTotient(0, 1))\n",
    "print(GeneralizedTotient(0, 2))\n",
    "print(GeneralizedTotient(0, 3))\n",
    "print(GeneralizedTotient(0, 4))\n",
    "print(GeneralizedTotient(0, 15))\n",
    "print(GeneralizedTotient(0, 30))\n",
    "print(GeneralizedTotient(0, 60))\n",
    "print(GeneralizedTotient(0, 2400))\n",
    "\n",
    "print(GeneralizedTotient(2, 1))\n",
    "print(GeneralizedTotient(2, 2))\n",
    "print(GeneralizedTotient(2, 3))\n",
    "print(GeneralizedTotient(2, 4))\n",
    "print(GeneralizedTotient(2, 15))\n",
    "print(GeneralizedTotient(2, 30))\n",
    "print(GeneralizedTotient(2, 60))\n",
    "print(GeneralizedTotient(2, 2400))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9d2cc69-590e-4918-b352-6e84c1888507",
   "metadata": {},
   "source": [
    "## 2.15 Prove an equality involving a sum of generalized Euler totient functions. \n",
    "\n",
    "\n",
    "Totient $\\varphi (n)$ is the number of relative primes less than $n$. \n",
    "\n",
    "\n",
    "Generalized totient $\\varphi_k(n) = \\sum_{m<n \\textrm{ and } (m,n)=1}m^k$.\n",
    "\n",
    "\n",
    "When $k=0$ we have $\\varphi_{0}(n) = \\varphi (n)$.\n",
    "\n",
    "\n",
    "For $k=2$ and $n=7$ we have relative primes $1, 2, 3, 4, 5, 6$ and $\\varphi_2(7) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 + 6^2$. \n",
    "\n",
    "\n",
    "Problem: Show \n",
    "\n",
    "\n",
    "$\\begin{align}\\sum_{d|n} \\frac{\\varphi_k(d)}{d^k} = \\frac{1^k + 2^k + \\cdots + n^k}{n^k}\\end{align}$\n",
    "\n",
    "\n",
    "The following code shows that equality holds pretty reliably... maybe with the occasional Python rounding error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5e46657-c498-42ed-89c9-0cfbac9e6480",
   "metadata": {},
   "outputs": [],
   "source": [
    "first_n = 1\n",
    "last_n = 60\n",
    "first_k = 0\n",
    "last_k = 5\n",
    "precision = 3\n",
    "for n in range(first_n, last_n + 1,1):\n",
    "    for k in range(first_k, last_k + 1, 1):\n",
    "        lhs_sum = 0\n",
    "        for d in ant.divisors(n): lhs_sum += ant.GeneralizedTotient(k, d) / d**k\n",
    "        rhs_sum = 0\n",
    "        for i in range(1, n+1, 1): rhs_sum += i**k\n",
    "        rhs = rhs_sum / (n**k)\n",
    "        \n",
    "        ratio_sum = round(lhs_sum, precision)\n",
    "        sum_ratio = round(rhs, precision)\n",
    "        if not ratio_sum == sum_ratio: \n",
    "            print(\"ERROR...........\")\n",
    "            print(n, k, ratio_sum, sum_ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b91993f1-86fb-4281-b043-878aa70a2630",
   "metadata": {},
   "source": [
    "## Problem 2.13\n",
    "\n",
    "This problem is thematic to chapters 1 and 2 as it brings together several important ideas. The problem \n",
    "is easily solved by Chip Hurst in a couple of lines but I belabor matters.\n",
    "\n",
    "\n",
    "The problem is advertised as \"*The product form of the Mobius inversion formula*\": \n",
    "Start with a function $f(n) > 0 \\; \\forall \\; n$ and a second real-valued function $a(n)$ \n",
    "where $a(1) \\ne 0$. Show that \n",
    "\n",
    "\n",
    "$\\begin{align}\n",
    "\\large{g(n)=\\prod_{d|n}{f(d)^{a(n/d)}} \\; \\iff \\; f(n)=\\prod_{d|n}{g(d)^{a^{-1}(n/d)}}.}\n",
    "\\end{align}$\n",
    "\n",
    "\n",
    "The exponentiation is not a typo! Here $a$ and $a^{-1}$ are Dirichlet inverse functions. \n",
    "(The problem statement in the text uses $b$ in place of $a^{-1}$.)\n",
    "\n",
    "\n",
    "***Note: Mobius inversion sets up a relationship between two functions f and g\n",
    "via the Mobius function. This is not the same thing as Dirichlet inversion \n",
    "where $f * f^{-1} = f^{-1} * f = I.$***\n",
    "\n",
    "\n",
    "Zeroth belabor: Confirm the recursive formula for the Dirichlet inverse for $f(n) = n$.\n",
    "\n",
    "\n",
    "First belabor: Demonstrate the unit function u(n) = 1 is the Dirichlet inverse of \n",
    "the Mobius function.\n",
    "\n",
    "Second belabor: The Mobius Inversion Formula presents a function's inverse\n",
    "in terms of a sum involving the Mobius function. This can be recursively substituted \n",
    "to give the function in terms of itself. Check this\n",
    "empirically using $f(n) = n$.\n",
    "\n",
    "\n",
    "Third belabor: With $f$ as above and defining the second function $a$ the same as\n",
    "$f$: $a(n)=n$: Does the product form of the Mobius Inversion Formula yield $f$ via $g$?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8649ca7a-15db-4d9a-90f3-202c88543637",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Zeroth Belabor\n",
    "\n",
    "print(\"Problem 2.13 on the product form of the Mobius Inversion Formula...            \")\n",
    "print(\"   Part 0: Show the Dirichlet inverse by recursion works properly (gives I).   \")\n",
    "print(\"                                                                               \")\n",
    "print(\"Let's use the function f(n) = N(n) = n                                         \")\n",
    "print(\"  Recall that the identity function is [1, 0, 0, 0, ...] = I                   \")\n",
    "print(\"  The Dirichlet product of a function and its inverse gives I (commutative).   \")\n",
    "print(\"                                                                               \")\n",
    "print(\"  f-1(1) = 1/f(1); and then for n > 1:                                         \")\n",
    "print(\"  f-1(n) = -1/f(1) multiplied by S where                                       \")\n",
    "print(\"                                                                               \")\n",
    "print(\" S = sum over d divides n, d less than n: Of E where                           \")\n",
    "print(\"                                                                               \")\n",
    "print(\" E = f(n/d) multiplied by f-1(d). Noting d < n is key to recursion.            \")\n",
    "print(\"                                                                               \")\n",
    "\n",
    "test_n = 35\n",
    "F      = [ant.N(n) for n in range(test_n + 1)]      # F holds values of f in the recursive space from 1 to test_n\n",
    "                                                    #   Note the list indices match n and we ignore F[0]\n",
    "G      = [0, 1/F[1]]                                # G[] will be the inverse of F built recursively, also with 'ignore 0' indexing\n",
    "    \n",
    "print('The list length of G should be two elements: ' + str(len(G)))\n",
    "print()\n",
    "\n",
    "# here the recursion formula calculates values of G.\n",
    "#   each value is a sum (so there is an inner loop)\n",
    "#   and the recursion proceeds from 1 to test_n (outer loop)\n",
    "\n",
    "G = [0]\n",
    "G.append(1/F[1])\n",
    "leading_coefficient = -1/F[1]\n",
    "for n in range(2, test_n + 1):\n",
    "    S = 0\n",
    "    for d in ant.divisors_less_n(n): S += ant.N(n//d)*G[d]\n",
    "    G.append(leading_coefficient * S)\n",
    "    \n",
    "print(\"Confirmation of Dirichlet inverse:\")\n",
    "print()\n",
    "print(\"The function f(n) is in list F: \" + str(F[1:]))\n",
    "print(\"The function f-1(n) is in list G: \" + str(G[1:]))\n",
    "print()\n",
    "print(\"Confirm their product is I: \")\n",
    "print(ant.ScanningListDirichlet(F, G)[1:])\n",
    "print(ant.ScanningListDirichlet(G, F)[1:])\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69bcb476-794f-42b5-bba8-cc161bc30749",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First belabor: Show the unit function U and the Mobius function are Dirichlet inverses\n",
    "mu, u = [0], [0]\n",
    "for n in range(1, test_n + 1):\n",
    "    mu.append(ant.Mobius(n))\n",
    "    u.append(ant.U(n))\n",
    "    \n",
    "print(ant.ScanningListDirichlet(mu, u)[1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55045cf4-602a-4510-be9e-e2644f9809e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Second belabor: Show that the Mobius Inversion Formula returns f by substitution:\n",
    "#   Suppose some arithmetic function f is equal to the sum over the divisors of n of some function g(d).\n",
    "#   Then g(n) is f * mu. The plan is to calculate g and show it generates f.\n",
    "    \n",
    "GMob = [0]\n",
    "for n in range(1, test_n + 1):\n",
    "    GMob.append(ant.Dirichlet(ant.N, ant.Mobius, n))\n",
    "\n",
    "FRegen = [0]\n",
    "for n in range(1, test_n + 1):\n",
    "    FRegen.append(ant.ListDirichlet(GMob, u, n))\n",
    "\n",
    "print(F[1:])\n",
    "print(GMob[1:])\n",
    "print()\n",
    "print(\"        (by the way that looks an awful lot like the totient...)\")\n",
    "print()\n",
    "print(\"Confirm F is recovered as N(n): \")\n",
    "print(FRegen[1:])\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87e1aaa1-2769-4f20-ac86-66a6acdd12b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Third belabor: Demonstrate that the product form of the Mobius Inversion Formula works as advertised.\n",
    "#   First condition: f(n) > 0 for all n: Define f(n) = n as above.\n",
    "#   Second condition: a(n) is real and a(1) is not zero: Set a(n) = n\n",
    "#   Third condition: a-inv is the Dirichlet inverse of a(n): a-inv(n) = n times mu(n) (not a Dirichlet product).\n",
    "#\n",
    "#   Result to show: Define g as the product over divisors d: f(d) raised to the a(n/d); if and only if \n",
    "#       f is the product over divisors d of g(d) raised to the a-inv(n/d). \n",
    "\n",
    "def a(n): return ant.N(n)\n",
    "\n",
    "def a_inverse(n): return ant.N(n)*ant.Mobius(n)\n",
    "\n",
    "print(ant.Dirichlet(a, a_inverse, 3))      # Confirms result is I(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40460215-089b-4773-93f1-dfb99c407481",
   "metadata": {},
   "source": [
    "## 2.25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "c7cc28d9-7419-454d-bb65-45000c824513",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n     tot      mu      tot*mu      tm*m      I\n",
      "1       1       1       1       1       1\n",
      "2       1       -1       -1       0       0\n",
      "3       2       -1       -2       0       0\n",
      "4       2       0       0       1       0\n",
      "5       4       -1       -4       0       0\n",
      "6       2       1       2       0       0\n",
      "7       6       -1       -6       0       0\n",
      "8       4       0       0       2       0\n",
      "9       6       0       0       2       0\n",
      "10       4       1       4       0       0\n",
      "11       10       -1       -10       0       0\n",
      "12       4       0       0       0       0\n",
      "13       12       -1       -12       0       0\n",
      "14       6       1       6       0       0\n",
      "15       8       1       8       0       0\n",
      "16       8       0       0       4       0\n",
      "17       16       -1       -16       0       0\n",
      "18       6       0       0       0       0\n",
      "19       18       -1       -18       0       0\n",
      "20       8       0       0       0       0\n",
      "21       12       1       12       0       0\n",
      "22       10       1       10       0       0\n",
      "23       22       -1       -22       0       0\n",
      "24       8       0       0       0       0\n",
      "25       20       0       0       4       0\n",
      "26       12       1       12       0       0\n",
      "27       18       0       0       6       0\n",
      "28       12       0       0       0       0\n",
      "29       28       -1       -28       0       0\n",
      "30       8       -1       -8       0       0\n",
      "31       30       -1       -30       0       0\n",
      "32       16       0       0       8       0\n",
      "33       20       1       20       0       0\n",
      "34       16       1       16       0       0\n",
      "35       24       1       24       0       0\n",
      "36       12       0       0       2       0\n",
      "37       36       -1       -36       0       0\n",
      "38       18       1       18       0       0\n",
      "39       24       1       24       0       0\n",
      "40       16       0       0       0       0\n",
      "41       40       -1       -40       0       0\n",
      "42       12       -1       -12       0       0\n"
     ]
    }
   ],
   "source": [
    "print(\"n     tot      mu      tot*mu      tm*m      I\")\n",
    "for n in range(1, 43):\n",
    "    t = ant.Totient(n)\n",
    "    m = ant.Mobius(n)\n",
    "    d = ant.Dirichlet(ant.TotientMobiusProduct, ant.Totient, n)\n",
    "    I = 1 if n == 1 else 0\n",
    "    print(n, '     ', t, '     ', m, '     ', t*m, '     ', d, '     ', I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87084822-3919-4140-a0da-febef50e2b59",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}