numberphile.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A prime number generator\n",
    "__MATH 420__<br>\n",
    "_Spring 2021_<br>\n",
    "\n",
    "See https://www.youtube.com/watch?v=_gCKX6VMvmU\n",
    "\n",
    "Let's write a Julia function that implements the prime number generator that is described in the Numberphile video. It's a good candidate for a memoizing function. Here it is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Memoize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mygenprime (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@memoize function mygenprime(n::Integer)\n",
    "    if n == 1  \n",
    "        2.920050977316\n",
    "    else\n",
    "        f = mygenprime(n-1)\n",
    "        floor(f)*(f - floor(f) + 1)\n",
    "    end\n",
    "end            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use this function to build an array of the first 15 primes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15-element Vector{Tuple{Int64, Int64}}:\n",
       " (1, 2)\n",
       " (2, 3)\n",
       " (3, 5)\n",
       " (4, 7)\n",
       " (5, 11)\n",
       " (6, 13)\n",
       " (7, 17)\n",
       " (8, 19)\n",
       " (9, 23)\n",
       " (10, 29)\n",
       " (11, 31)\n",
       " (12, 37)\n",
       " (13, 40)\n",
       " (14, 42)\n",
       " (15, 43)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[(k, floor(Int, mygenprime(k))) for k = 1 : 15]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oops! Forty isn't a prime. Neither is forty-two. What's the story? It's crazy to think that _all_ primes can be encoded into the number $2.920050977316$. Almost surely, to generate the primes past $37$, we need a more digits of the constant $2.920050977316$.  \n",
    "\n",
    "Fortunately, the video tells us how to compute this constant. Not surprisingly, the constant is expressed as an _infinite series_. The best we can do is sum its first few terms.\n",
    "\n",
    "Here is a function `Q` that sums the first `n` terms of the infinite series for this constant. So that we don't have to worry about rounding errors, we'll sum this using exact rational numbers.  The optional package `Primes` has a function that generates the primes for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Primes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Q (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function Q(n::Integer)\n",
    "    prod = BigInt(1)\n",
    "    s = BigInt(0)\n",
    "    for k = 1 : n\n",
    "        pk = prime(k)\n",
    "        s += (pk-1)//prod\n",
    "        prod *= pk\n",
    "    end\n",
    "    s\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Summing the first $64$ terms, we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2982076387136897711868778837562533128923101076886460857776351024737697552705029412066629838751794837984142947810782961//1021241207877052728809092149097626662644537663247243355111401451312772219284716228359734476103267543069998638899157773"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q(64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Converting to a decimal, we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.92005097731613471209256291711201946800272789932142671977268253310773377212777"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "BigFloat(Q(64))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Summing the first $128$ terms is no problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5227554455733777871260305151323808032492772521481472963485838824360748775617077180526329978266811309524859073074134253177854209904470343077182696221368015025622622213062490474563836827606091357779777088425776324911260290461579448439440774443274564890853890530370099364635758182172299843563765//1790227121493100199266413788084742281796363638940326181979063572262335989939146297006758045608216781053380619178432229671020652188884002829412501970104427630569751221728626885262586272945339984183616494370041526619082760890373018744558545629914917132431516486188023601797767156032558566683497"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q(128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.92005097731613471209256291711201946800272789932142671977268253310773377212777"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "BigFloat(Q(128))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference $Q(128) - Q(64)$ is small--and I mean tiny."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.60409249577861863914907299171977958705313392093787080488841313278988811892698e-123"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "BigFloat(Q(128)-Q(64))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the value of $Q(128)$ to generate primes. Let's rewrite our function and make the value of its first term an argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mygenprime2 (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@memoize function mygenprime2(n::Integer, f0::Rational)\n",
    "    if n == 1  \n",
    "       f0\n",
    "    else\n",
    "       f = mygenprime2(n-1,f0)\n",
    "       floor(f)*(f - floor(f) + 1)\n",
    "    end\n",
    "end            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = Q(128);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the 127th prime is correct, but the 128th is wrong:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "128-element Vector{Tuple{Int64, Int64, Int64}}:\n",
       " (1, 2, 2)\n",
       " (2, 3, 3)\n",
       " (3, 5, 5)\n",
       " (4, 7, 7)\n",
       " (5, 11, 11)\n",
       " (6, 13, 13)\n",
       " (7, 17, 17)\n",
       " (8, 19, 19)\n",
       " (9, 23, 23)\n",
       " (10, 29, 29)\n",
       " ⋮\n",
       " (120, 659, 659)\n",
       " (121, 661, 661)\n",
       " (122, 673, 673)\n",
       " (123, 677, 677)\n",
       " (124, 683, 683)\n",
       " (125, 691, 691)\n",
       " (126, 701, 701)\n",
       " (127, 709, 709)\n",
       " (128, 719, 718)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[(k, prime(k), floor(Int, mygenprime2(k,f0))) for k = 1 : 128]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a guess, in general using $Q(n)$ as the first term, we can correctly find the first $n-1$ primes, but not the nth prime; example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1669, 1669)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "(prime(263), floor(Int, mygenprime2(263, Q(264))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1693, 1692)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "(prime(264), floor(Int, mygenprime2(264, Q(264))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Question__: Could we use the same to scheme to generate terms of other sequences? Say we have an integer-valued sequence $p_1, p_2, p_3, \\dots $.  We calculate the constant the same way, but instead of the primes, we use the terms of some other sequence; thus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QQ (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function QQ(F::Function, n::Integer)\n",
    "    prod = BigInt(1)\n",
    "    s = BigInt(0)\n",
    "    for k = 1 : n\n",
    "        pk = F(k)\n",
    "        s += (pk-1)//prod\n",
    "        prod *= pk\n",
    "    end\n",
    "    s\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try everyone's favorite sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fib (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@memoize function fib(n::Integer)\n",
    "    if n == 1 || n == 2\n",
    "        1\n",
    "    else\n",
    "        fib(n-2) + fib(n-1)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ha! We get the terms of the Fibonacci sequence, but the sequence is _shifted_ by two terms. Actually, since the constant is\n",
    "$$\n",
    "   (p_1 -1) + \\frac{p_2 -1}{p_1}  + \\cdots,\n",
    "$$\n",
    "the first two terms of the series both vanish, so the values of $p_1$ and $p_2$ don't get incorporated into the constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = QQ(fib, 128);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24-element Vector{Tuple{Int64, Int64, BigInt}}:\n",
       " (1, 1, 2)\n",
       " (2, 1, 3)\n",
       " (3, 2, 5)\n",
       " (4, 3, 8)\n",
       " (5, 5, 13)\n",
       " (6, 8, 21)\n",
       " (7, 13, 34)\n",
       " (8, 21, 55)\n",
       " (9, 34, 89)\n",
       " (10, 55, 144)\n",
       " ⋮\n",
       " (16, 987, 2584)\n",
       " (17, 1597, 4181)\n",
       " (18, 2584, 6765)\n",
       " (19, 4181, 10946)\n",
       " (20, 6765, 17711)\n",
       " (21, 10946, 28657)\n",
       " (22, 17711, 46368)\n",
       " (23, 28657, 75025)\n",
       " (24, 46368, 121393)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[(k, fib(k), floor(BigInt, mygenprime2(k,f0))) for k = 1 : 24]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the conjecture that maybe we want to exclude $1$ from the range of the sequence, let's try $k \\mapsto k+1$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "70911686272574852599660141113652666333760118218640746369131464891915518607035687428261//26086951518479474527844035636340764967694390660885113527714205258062233600000000000000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xxx = QQ(x -> x+1, 64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK--this seems to work OK:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24-element Vector{Tuple{Int64, BigInt}}:\n",
       " (2, 2)\n",
       " (3, 3)\n",
       " (4, 4)\n",
       " (5, 5)\n",
       " (6, 6)\n",
       " (7, 7)\n",
       " (8, 8)\n",
       " (9, 9)\n",
       " (10, 10)\n",
       " (11, 11)\n",
       " ⋮\n",
       " (17, 17)\n",
       " (18, 18)\n",
       " (19, 19)\n",
       " (20, 20)\n",
       " (21, 21)\n",
       " (22, 22)\n",
       " (23, 23)\n",
       " (24, 24)\n",
       " (25, 25)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[(k+1, floor(BigInt, mygenprime2(k, xxx))) for k = 1 : 24]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about a sequence of squares?  We'll try $2^2, 3^3, \\dots$. First, we compute the constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "xxx = QQ(x -> (1+x)^2, 128);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This fails from the start--maybe there is an upper bound on the growth of the sequence that allows the method to work. And maybe there is some condition on the difference between consecutive terms as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16-element Vector{Tuple{Int64, Int64, BigInt}}:\n",
       " (1, 4, 5)\n",
       " (2, 9, 7)\n",
       " (3, 16, 9)\n",
       " (4, 25, 10)\n",
       " (5, 36, 11)\n",
       " (6, 49, 18)\n",
       " (7, 64, 31)\n",
       " (8, 81, 54)\n",
       " (9, 100, 93)\n",
       " (10, 121, 141)\n",
       " (11, 144, 147)\n",
       " (12, 169, 277)\n",
       " (13, 196, 356)\n",
       " (14, 225, 484)\n",
       " (15, 256, 758)\n",
       " (16, 289, 1512)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[(k, (1+k)^2, floor(BigInt, mygenprime2(k, xxx))) for k = 1 : 16]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Searching `oeis` (https://oeis.org/) for the sequence $5,7,9,10,11,18,31, \\dots $ doesn't find a match.  You might like to try sequences with other growth rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.0",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}