{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The largest prime (so far)\n",
    "\n",
    "A new record for the largest prime has been found lately. Let explore this number and know more about it and how to deal with it using some Python, Numpy and finally using we will take a loog at Cython and GMP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from datetime import datetime\n",
    "\n",
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This number of was found in 2016 using the Great Internet Mersenne Prime Search (GIMPS). This prime is one of the Mersenne Prime whice are of the defined by the following formula:\n",
    "\n",
    "$$2^p-1$$\n",
    "\n",
    "Where $p$ is a prime number. This is called a psudo-prime meaning not all Mersenne Prime are primes.\n",
    "\n",
    "$2^2-1=3$\n",
    "\n",
    "$2^3-1=7$\n",
    "\n",
    "$2^5-1=31$\n",
    "\n",
    "$2^7-1=127$   ...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2**2-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2**3-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "31"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2**5-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "127"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2**7-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The smallest none-prime (composite) number is $2^{11}-1=2047$ which is a composite number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2047"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2**11-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The largest prime number is that was found lately is:\n",
    "$$2^{74,207,281}-1$$\n",
    "\n",
    "The reason why I didn't put the actual value is that the number if very large (over 22 million digits)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p = 74207281\n",
    "the_number = (2 ** p) - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lucas–Lehmer test\n",
    "\n",
    "A very easy way to test a Mersenne Prime is the Lucas-Lehmer test which uses the Lucas series.\n",
    "\n",
    "$4, 14, 37634 ... $\n",
    "\n",
    "The series starts with 4 then uses the last value in this formula:\n",
    "\n",
    "$S_i=S_{i-1}^2-2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "S = 4\n",
    "print(S)\n",
    "print(len(str(S)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\n",
      "2 digits\n"
     ]
    }
   ],
   "source": [
    "S = S ** 2 - 2\n",
    "print(S)\n",
    "print(len(str(S)), \"digits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "194\n",
      "3 digits\n"
     ]
    }
   ],
   "source": [
    "S = S ** 2 - 2\n",
    "print(S)\n",
    "print(len(str(S)), \"digits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "37634\n",
      "5 digits\n"
     ]
    }
   ],
   "source": [
    "S = S ** 2 - 2\n",
    "print(S)\n",
    "print(len(str(S)), \"digits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1416317954\n",
      "10 digits\n"
     ]
    }
   ],
   "source": [
    "S = S ** 2 - 2\n",
    "print(S)\n",
    "print(len(str(S)), \"digits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2005956546822746114\n",
      "19 digits\n"
     ]
    }
   ],
   "source": [
    "S = S ** 2 - 2\n",
    "print(S)\n",
    "print(len(str(S)), \"digits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22338618 digits\n"
     ]
    }
   ],
   "source": [
    "# The largest prime (so far)\n",
    "the_number = 2 ** 74207281 - 1\n",
    "print(int(math.log10(the_number))+1, \"digits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test is if the Mersenne Prime number of $n^2-1$ is a divided by $S_{n-2}$ you should get a remaining of 0. The problem is this series goes very large very fast. So there is another way to do this test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0:00:00.000423 \n",
      "1 0:00:00.000017 \n",
      "2 0:00:00.000012 \n",
      "3 0:00:00.000014 \n",
      "4 0:00:00.000014 \n",
      "5 0:00:00.000013 \n",
      "6 0:00:00.000014 \n",
      "7 0:00:00.000014 \n",
      "8 0:00:00.000015 \n",
      "9 0:00:00.000014 \n",
      "10 0:00:00.000036 \n",
      "11 0:00:00.000051 \n",
      "12 0:00:00.000127 \n",
      "13 0:00:00.000358 \n",
      "14 0:00:00.001040 \n",
      "15 0:00:00.002735 \n",
      "16 0:00:00.004968 \n",
      "17 0:00:00.014782 \n",
      "18 0:00:00.044420 \n",
      "19 0:00:00.132890 \n",
      "20 0:00:00.398497 \n",
      "21 0:00:01.198042 \n",
      "22 0:00:03.592832 \n",
      "23 0:00:10.796848 \n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-20-4c550569e183>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mtime_stamp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdatetime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m     \u001b[0mS\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mS\u001b[0m \u001b[1;33m**\u001b[0m \u001b[1;36m2\u001b[0m \u001b[1;33m-\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m%\u001b[0m \u001b[0mthe_number\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      8\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      9\u001b[0m         \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdatetime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mtime_stamp\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m\"\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "p = 74207281\n",
    "the_number = (2 ** p) - 1\n",
    "\n",
    "S = 4\n",
    "time_stamp = datetime.now()\n",
    "for i in range(p-2):\n",
    "    S = (S ** 2 - 2) % the_number\n",
    "    if i % 1 == 0:\n",
    "        print(i, datetime.now() - time_stamp,\"\")\n",
    "        time_stamp = datetime.now()\n",
    "if S == 0:\n",
    "    print(\"PRIME\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "10\n",
      "20\n",
      "30\n",
      "40\n",
      "50\n",
      "PRIME\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.4/dist-packages/IPython/utils/path.py:264: UserWarning: get_ipython_cache_dir has moved to the IPython.paths module\n",
      "  warn(\"get_ipython_cache_dir has moved to the IPython.paths module\")\n"
     ]
    }
   ],
   "source": [
    "%%cython\n",
    "cdef unsigned long p = 61\n",
    "cdef unsigned long P = (2 ** p) - 1\n",
    "\n",
    "S = 4\n",
    "for i in range(p-2):\n",
    "    S = S ** 2 - 2\n",
    "    S = S % P\n",
    "    if i % 10 == 0:\n",
    "        print(i)\n",
    "if S == 0:\n",
    "    print(\"PRIME\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1000\n",
      "2000\n",
      "3000\n",
      "4000\n",
      "5000\n",
      "6000\n",
      "7000\n",
      "8000\n",
      "9000\n",
      "10000\n",
      "11000\n",
      "12000\n",
      "13000\n",
      "14000\n",
      "15000\n",
      "16000\n",
      "17000"
     ]
    }
   ],
   "source": [
    "%%cython --link-args=-lgmp\n",
    "\n",
    "cdef extern from \"gmp.h\":\n",
    "    ctypedef struct mpz_t:\n",
    "        pass\n",
    "    \n",
    "    ctypedef unsigned long mp_bitcnt_t\n",
    "    \n",
    "    cdef void mpz_init(mpz_t)\n",
    "    cdef void mpz_init_set_ui(mpz_t, unsigned int)\n",
    "    \n",
    "    cdef void mpz_add(mpz_t, mpz_t, mpz_t)\n",
    "    cdef void mpz_add_ui(mpz_t, const mpz_t, unsigned long int)\n",
    "    \n",
    "    cdef void mpz_sub (mpz_t, const mpz_t, const mpz_t)\n",
    "    cdef void mpz_sub_ui (mpz_t, const mpz_t, unsigned long int)\n",
    "    cdef void mpz_ui_sub (mpz_t, unsigned long int, const mpz_t)\n",
    "    \n",
    "    cdef void mpz_mul (mpz_t, const mpz_t, const mpz_t)\n",
    "    cdef void mpz_mul_si (mpz_t, const mpz_t, long int)\n",
    "    cdef void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int)\n",
    "    \n",
    "    cdef void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t)\n",
    "    \n",
    "    cdef void mpz_mod (mpz_t, const mpz_t, const mpz_t)\n",
    "    \n",
    "    cdef unsigned long int mpz_get_ui(const mpz_t)\n",
    "\n",
    "#cdef unsigned long p = 61\n",
    "cdef mp_bitcnt_t p = 74207281\n",
    "cdef mpz_t t # = 1\n",
    "cdef mpz_t a # = 1\n",
    "cdef mpz_t P # = (2 ** p) - 1\n",
    "cdef mpz_t S # = 4\n",
    "\n",
    "mpz_init(t)\n",
    "mpz_init_set_ui(t, 1)\n",
    "\n",
    "mpz_init(a)\n",
    "mpz_init_set_ui(a, 2)\n",
    "\n",
    "mpz_init(P)\n",
    "mpz_mul_2exp(P,t,p)\n",
    "mpz_sub_ui(P,P,1)\n",
    "\n",
    "mpz_init(S)\n",
    "mpz_init_set_ui(S, 4)\n",
    "\n",
    "for i in range(p-2):\n",
    "    #S = S ** 2 - 2\n",
    "    mpz_mul(S,S,S)\n",
    "    mpz_sub_ui(S,S,2)\n",
    "    \n",
    "    #S = S % P\n",
    "    mpz_mod(S,S,P)\n",
    "    \n",
    "    if i % 1000 == 0:\n",
    "        print(i)\n",
    "if mpz_get_ui(S) == 0:\n",
    "    print(\"PRIME\")\n",
    "else:\n",
    "    print(\"COMPOSITE\")\n",
    "\n",
    "#print(mpz_get_ui(P))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}