{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "from random import *\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on calcule a^n mod N grâce à l'exponentiation rapide \n",
    "def Base(n,b): #Ecriture de n en base b\n",
    "    L=[]\n",
    "    while n>0:\n",
    "        L=L+[n%b]\n",
    "        n=n//b\n",
    "    return L  \n",
    "\n",
    "def ExpRapide(a,n,N): #Calcul de a^n modulo N par exponentiation rapide. \n",
    "    L=Base(n,2)\n",
    "    P=[a%N]\n",
    "    i=0\n",
    "    while i<len(L)-1: #Calcul des carrés itérés de a modulo N\n",
    "        b=(P[i]*P[i])%N\n",
    "        P=P+[b]\n",
    "        i=i+1\n",
    "    p=1\n",
    "    for i in range(len(L)): #Calcul le produit des carrés itérés de a associés aux coefficients non nuls de L\n",
    "        if L[i]!=0:\n",
    "            p=(p*P[i])%N\n",
    "    return p    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie True a^{n-1} = 1 mod n pour tout a premier à n (cf condition (1) du TP \n",
    "# et renvoie False sinon\n",
    "def TestFermat(n):\n",
    "    for a in range(1,n):\n",
    "        if math.gcd(a,n)==1 and ExpRapide(a,n-1,n) != 1:\n",
    "            return False\n",
    "    return True\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 True\n",
      "3 True\n",
      "4 False\n",
      "5 True\n",
      "6 False\n",
      "7 True\n",
      "8 False\n",
      "9 False\n",
      "10 False\n",
      "11 True\n",
      "12 False\n",
      "13 True\n",
      "14 False\n",
      "15 False\n",
      "16 False\n",
      "17 True\n",
      "18 False\n"
     ]
    }
   ],
   "source": [
    "for n in range(2,19):\n",
    "    print(n,TestFermat(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie la liste des entiers n <= x tels que TestFermat(n) = True\n",
    "def PremierFermat(x):\n",
    "    L = []  # liste des entiers à renvoyer\n",
    "    for n in range(2,x+1): # on parcourt les entiers entre 2 et x\n",
    "        if TestFermat(n):  # si TestFermat(n) = True on ajoute n à notre liste\n",
    "            L = L+[n]\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]"
      ]
     },
     "execution_count": 239,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PremierFermat(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on supprime l'élément e (avec répétitions éventuelles) de la liste L\n",
    "# sous-algo utilisé dans crible\n",
    "def supprElemListe(L,e):\n",
    "    N = L.count(e)  # on compte le nombre de fois où e apparaît dans la liste L\n",
    "    for i in range(N): # on supprime e le nombre de fois où il apparaît\n",
    "        L.remove(e)\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie la liste des nombres premiers <= n en utilisant le crible d'Eratosthène\n",
    "def crible(n):\n",
    "    d=list(range(n+1)) # d est la liste des entiers entre 0 et n\n",
    "    if n > 0: # 1 n'est pas premier\n",
    "        d[1]=0\n",
    "    r=2 #on initialise à p=2\n",
    "    while r<math.sqrt(n)+1: # on peut arrêter le crible dès qu'on est plus grand que racine de n (1ère optimisation)\n",
    "        if d[r]!=0: # on vérifie que le d[r] n'a pas déjà été supprimé (2ème optimisation)\n",
    "            for i in range(r*r,n+1,r): #on \"supprime\" les multiples de r (en commençant à r^2) en mettant leur valeur à 0\n",
    "                d[i]=0\n",
    "        r=r+1   \n",
    "    return supprElemListe(d,0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "435\n",
      "430\n"
     ]
    }
   ],
   "source": [
    "# la liste PremierFermat(N) contient tous les nombres premiers, donc contient crible(N)\n",
    "# On remarque que pour N = 3000 il y a 5 entiers qui sont dans PremierFermat(N) mais \n",
    "# pas dans crible(N). Ces entiers ne sont donc pas des nombres premiers!\n",
    "\n",
    "N = 3000\n",
    "print(len(PremierFermat(N)))\n",
    "print(len(crible(N)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie la liste des k premiers contre-exemples (qui ne sont pas Fermatpremiers mais pas premiers)\n",
    "# horrible temps de calcul, plusieurs minutes pour calculer pour k=2\n",
    "# cet algorithme n'est pas recommandé, il y a beaucoup de répétition de calculs\n",
    "def contreExemples(k):\n",
    "    n = 2\n",
    "    L = []  # liste des contre-exemples\n",
    "    j= 0   # compte le nombre de contre-exemples trouvés\n",
    "    for j in range(k): \n",
    "        # on compare simplement la taille des listes\n",
    "        # dès qu'il y a un saut de taille, c'est qu'on vient de trouver un contre-exemple\n",
    "        while len(PremierFermat(n)) == len(crible(n))+j: \n",
    "            n =  n + 1  \n",
    "        L = L+[n]    \n",
    "    return L\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[561]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contreExemples(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "# On teste si un nombre est premier\n",
    "# sous-algo utilisé dans contreExemples2\n",
    "def isPrime(n):\n",
    "    for i in range(2,int(n**0.5)+1):\n",
    "        if n%i==0:\n",
    "            return False       \n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [],
   "source": [
    "# On code une deuxième version de PremierFermat qui renvoie un couple de booléens (b1,b2)\n",
    "# b1 = True si n est premier, False sinon, et b2 = True si est premier de Fermat, False sinon\n",
    "# En faisant cela on optimise un peu l'algo contreExemples2 ci-dessous\n",
    "\n",
    "\n",
    "def TestFermat2(n):\n",
    "    Prime = True  # change en false si n n'est pas premier\n",
    "    for a in range(1,n):\n",
    "        if math.gcd(a,n) != 1:\n",
    "            Prime = False\n",
    "        elif ExpRapide(a,n-1,n) != 1: # si n n'est pas premier de Fermat, il n'est pas non plus premier\n",
    "            return (False,False)\n",
    "    return (Prime,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie la liste des k premiers contre-exemples (qui ne sont pas premiers) en un temps \n",
    "# \"raisonnable\"\n",
    "def contreExemples2(k):\n",
    "    L = []  # liste des contre-exemples à renvoyer\n",
    "    n = 2\n",
    "    while len(L) < k:  # tant qu'on n'a pas atteint le bon nombre de contre-exemples on continue\n",
    "        if TestFermat2(n) == (False,True): # si n n'est pas premier mais Fermat-premier on l'ajoute\n",
    "            L = L+[n]\n",
    "        n = n+1\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[561, 1105, 1729, 2465, 2821]"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contreExemples2(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on tire un nombre a au hasard et on teste la condition (1) du TP pour ce nombre\n",
    "def TestFermatP(n):\n",
    "    a = randint(1,n-1) # entier aléatoire entre 1 et n-1\n",
    "    if math.gcd(a,n) !=1 or ExpRapide(a,n-1,n) != 1:\n",
    "        return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 248,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TestFermatP(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.5014634132385254\n",
      "1.3357200622558594\n"
     ]
    }
   ],
   "source": [
    "n = 10**2022\n",
    "\n",
    "# on mesure le temps de calcul pour l'algo TestFermat\n",
    "start=time.time()\n",
    "\n",
    "TestFermat(n)\n",
    "\n",
    "end=time.time()\n",
    "duree=end-start\n",
    "print(duree)\n",
    "\n",
    "# on mesure maintenant le temps de calcul pour TestFermatP\n",
    "start=time.time()\n",
    "\n",
    "TestFermatP(n)\n",
    "\n",
    "end=time.time()\n",
    "duree=end-start\n",
    "print(duree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [],
   "source": [
    "# algo intermédiaire utilisé dans TestMillerRabin\n",
    "#  Entrée : entier n impair.\n",
    "# On renvoie h et N tels que n = 1 + 2^h*m avec m impair\n",
    "def decompoMiller(n):\n",
    "    h = 0\n",
    "    m = n-1\n",
    "    while m%2 == 0:\n",
    "        h = h+1\n",
    "        m = m//2\n",
    "    return (h,m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {},
   "outputs": [],
   "source": [
    "# algo intermédiaire utilisé dans TestMillerRabin\n",
    "# On construit la liste des a^{2^k*m} mod n avec k entre 0 et h (où m et h sont donnés par l'algo précédent)\n",
    "\n",
    "def CreationListe(n,a):\n",
    "    (h,m) = decompoMiller(n)  # rappel : n = 1+2^h*m\n",
    "    L = [ExpRapide(a,m,n)]    # on initialise en calculant a^m mod n\n",
    "    for i in range(h):        # le terme L[i+1] est simplement L[i]^2\n",
    "        L = L+[(L[i]**2)%n]\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie True si la liste ne contient que des 1, False sinon\n",
    "def testListe1(L):\n",
    "     return L.count(1) == len(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [],
   "source": [
    "# algo récursif pour tester si la liste est de la forme [*,...,*,n-1,1,...,1]\n",
    "def testListe2(L,n):\n",
    "    # si on a la bonne forme, L termine au minimum par n-1,1, donc sa longueur est >= 2\n",
    "    if len(L) < 2:     \n",
    "        return False\n",
    "    # Si L commence par n-1 et ensuite ne contient que des 1, c'est bon\n",
    "    if L[0] == n-1 and L.count(1) == len(L) - 1:\n",
    "        return True\n",
    "    # Sinon, on supprime le premier terme et on applique l'algo à la liste tronquée\n",
    "    G = L[1:]  # subtilité ici : on fait une copie de L pour éviter que notre algo ne modifie L!\n",
    "    return testListe2(G, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remarquons que si a^{2^k*m} = n-1 mod n, alors toutes les valeurs suivantes valent 1\n",
    "# Donc pour tester si la liste CreationListe(n,a) a la forme 2, il suffit simplement\n",
    "# de vérifier que n-1 appartient à la liste ! Cela donne le programme simplifié suivant :\n",
    "\n",
    "# algo récursif pour tester si la liste est de la forme [*,...,*,n-1,1,...,1]\n",
    "def testListe2bis(L,n):\n",
    "    G = L[:-1]  # On fait une copie de L en supprimant le dernier élément pour éviter que notre algo ne modifie L\n",
    "    return n-1 in G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prend en entrée n. Renvoie True si n vérifie la condition (2) du TP, False sinon\n",
    "\n",
    "def TestMillerRabin(n):\n",
    "    if n%2 == 0 and n > 2: # Si n est pair > 2, on renvoie directement False\n",
    "        return False\n",
    "    a = randint(1,n-1)\n",
    "    L = CreationListe(n,a)\n",
    "    if math.gcd(a,n) != 1:\n",
    "        return False\n",
    "    return (math.gcd(a,n) == 1) and (testListe1(L) or testListe2(L,n))\n",
    "        \n",
    "    # si pgcd(a,n) != 1, renvoie False.\n",
    "    # Si pgcd(a,n) = 1, renvoie True dès que L est du type 1 ou du type 1, renvoie False sinon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "print(TestMillerRabin(11))  # TestMillerRabin(n) est toujours vrai pour n premier\n",
    "print(TestMillerRabin(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 263,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testListe2([4,1,1,5,4,1,1,1],5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8279\n"
     ]
    }
   ],
   "source": [
    "n = 100129*200257\n",
    "\n",
    "N = 10000 # nombre de répétitions\n",
    "nbFalse = 0 # on compte les fois où on obtient False\n",
    "\n",
    "for k in range(N):\n",
    "    if not TestMillerRabin(n):\n",
    "        nbFalse = nbFalse + 1\n",
    "\n",
    "print(nbFalse/N)  # On fait la moyenne. On constate que cette moyenne est >= 0.75"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie True si n passe le test de Miller-Rabin pour k répétitions indépendantes\n",
    "# on se sert de cet algo dans CribleMR\n",
    "def TestMillerRabinR(n,k):\n",
    "    j = 0\n",
    "    while TestMillerRabin(n) and j < k:\n",
    "        j = j+1\n",
    "    return j == k # on renvoie True si on a passé les k tests\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie la liste des entiers <= N qui ont passé k fois le test de Miller-Rabin avec succès\n",
    "def CribleMR(N,k):\n",
    "    L = []\n",
    "    for n in range(2,N+1):\n",
    "        if TestMillerRabinR(n,k):\n",
    "            L = L+[n]\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Comparaison entre crible(1000) et CribleMR(1000,k)\n",
      "On affiche k et la différence des longueurs\n",
      "1 9\n",
      "2 1\n",
      "3 0\n",
      "4 0\n",
      "5 0\n",
      "6 0\n",
      "7 0\n",
      "8 0\n",
      "9 0\n",
      "10 0\n"
     ]
    }
   ],
   "source": [
    "# CribleMR(N,k) contient la liste Crible(N) des premiers <= N\n",
    "# Pour comparer ces listes, il suffit donc de comparer leur nombre d'éléments\n",
    "\n",
    "print(\"Comparaison entre crible(1000) et CribleMR(1000,k)\")\n",
    "print(\"On affiche k et la différence des longueurs\")\n",
    "\n",
    "N = 1000\n",
    "nbPrime = len(crible(N))\n",
    "\n",
    "for k in range(1,11):\n",
    "    print(k, len(CribleMR(N,k)) - nbPrime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {},
   "outputs": [],
   "source": [
    "# renvoie le plus petit p >= m passant avec succès 10 fois le test de Miller-Rabin\n",
    "def PremierSuivant(m):\n",
    "    p = m\n",
    "    while not TestMillerRabinR(p,10):\n",
    "        p = p+1\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 269,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1009"
      ]
     },
     "execution_count": 269,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PremierSuivant(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on suppose ici n impair\n",
    "# on renvoie la (r,L), où L est la liste des a ne vérifiant pas la condition (2) du TP\n",
    "# et r = longueur de L (donc = nombre de témoins)\n",
    "\n",
    "def TemoinsMR(n):\n",
    "    G = []   # liste des a ne vérifiant pas (2)\n",
    "    for a in range(1,n):\n",
    "        L = CreationListe(n,a)\n",
    "        if not (testListe1(L) or testListe2bis(L,n)):\n",
    "            G = G+[a]\n",
    "    return len(G),G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on suppose ici n impair\n",
    "# on renvoie la (r,L), où L est la liste des a ne vérifiant pas la condition (2) du TP\n",
    "# et r = longueur de L (donc = nombre de témoins)\n",
    "# c'est le même algo que précédemment mais un peu plus optimisé sans faire appel aux sous-algo\n",
    "\n",
    "def TemoinsMR2(n):\n",
    "    (h,m) = decompoMiller(n)\n",
    "    L = []   # liste des a ne vérifiant pas (2)\n",
    "    for a in range(1,n):\n",
    "        M = ExpRapide(a,m,n)\n",
    "        if M != 1:  # on teste si la première condition de (2) est fausse\n",
    "            k = 0\n",
    "            while M != n-1 and k < h:  # on teste si la 2ème condition de (2) est fausse\n",
    "                M = M**2%n\n",
    "                k = k+1\n",
    "            if k == h:    # si k = h c'est que (2) est fausse pour tout k < h\n",
    "                L = L+[a]\n",
    "    return len(L),L\n",
    "                \n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(12, [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])\n",
      "(12, [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])\n"
     ]
    }
   ],
   "source": [
    "n = 15\n",
    "print(TemoinsMR(n))\n",
    "print(TemoinsMR2(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Liste des témoins pour n=25\n",
      "(20, [2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 19, 20, 21, 22, 23])\n",
      " \n",
      "Les non-témoins sont 7, 18 et 24\n",
      " \n",
      "On affiche à la main la liste des a^{2^km} mod n\n",
      "1 [1, 1, 1, 1]\n",
      "2 [8, 14, 21, 16]\n",
      "3 [2, 4, 16, 6]\n",
      "4 [14, 21, 16, 6]\n",
      "5 [0, 0, 0, 0]\n",
      "6 [16, 6, 11, 21]\n",
      "7 [18, 24, 1, 1]\n",
      "8 [12, 19, 11, 21]\n",
      "9 [4, 16, 6, 11]\n",
      "10 [0, 0, 0, 0]\n",
      "11 [6, 11, 21, 16]\n",
      "12 [3, 9, 6, 11]\n",
      "13 [22, 9, 6, 11]\n",
      "14 [19, 11, 21, 16]\n",
      "15 [0, 0, 0, 0]\n",
      "16 [21, 16, 6, 11]\n",
      "17 [13, 19, 11, 21]\n",
      "18 [7, 24, 1, 1]\n",
      "19 [9, 6, 11, 21]\n",
      "20 [0, 0, 0, 0]\n",
      "21 [11, 21, 16, 6]\n",
      "22 [23, 4, 16, 6]\n",
      "23 [17, 14, 21, 16]\n",
      "24 [24, 1, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "n = 25\n",
    "print(\"Liste des témoins pour n=25\")\n",
    "print(TemoinsMR(n))\n",
    "\n",
    "print(\" \")\n",
    "print(\"Les non-témoins sont 7, 18 et 24\")\n",
    "\n",
    "print(\" \")\n",
    "print(\"On affiche à la main la liste des a^{2^km} mod n\")\n",
    "for a in range(1,n):\n",
    "    print(a,CreationListe(n,a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}