{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# M1 Cryptographie - TD 3 et 4 : Arithmétique modulaire, nombres premiers\n",
    "\n",
    "## PGCD récursif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def gcd(a,b):\n",
    "    if b==0: return a\n",
    "    return gcd(b,a%b)\n",
    "\n",
    "gcd(12,18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notons que si $b>a$, le premier appel récursif sera \n",
    "<pre>\n",
    "gcd(a,b%a)\n",
    "</pre>\n",
    "Il n'y a donc pas besoin de distinguer ce cas.\n",
    "L'algorithme étendu peut aussi se coder de manière récursive :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, -3, 2)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def xgcd(a,b):\n",
    "    if b==0:\n",
    "        return (a,1,0)\n",
    "    else:\n",
    "        d,u,v = xgcd(b,a % b)\n",
    "        return (d,v,u-(a//b)*v)\n",
    "xgcd(17,26)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 2, -3)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xgcd(26,17)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explication\n",
    "Si $a=bq+r$ et $d=a\\wedge b = b\\wedge r$ vérifie $d = u'b+v'r$,\n",
    "alors \n",
    "$$\n",
    "d = u'b + v'(a-bq) = v'a +(u'-qv')b.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inversemod(a,n):\n",
    "    d,u,v = xgcd(a,n)\n",
    "    if d==1: return u%n\n",
    "    else: return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "170141183460469231731687303715884105727\n",
      "100909560761089108909934662113775107751\n"
     ]
    }
   ],
   "source": [
    "p = 2**127-1\n",
    "print (p)\n",
    "print (inversemod(666,p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "None\n"
     ]
    }
   ],
   "source": [
    "print (inversemod(8,26))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Crible d'Eratosthène\n",
    "\n",
    "Pour trouver les nombres premiers entre 2 et n, on part d'une liste `pp` contenant seulement 2, et d'une liste `mm`\n",
    "contenant tous les nombres entre 2 et n qui ne sont pas multiple de 2. `mm[0]` est donc 3, qui est forcément premier\n",
    "puisque pas multiple des nombres premiers qui lui sont inférieurs. On fait donc passer 3 à la fin de `pp` et on élimine\n",
    "de `mm` tous les multiples de 3. On recommence, 5 passe dans `pp` et les multiples de 5 sont éliminés. A chaque\n",
    "étape, `mm[0]` est un nombre premier $p$, on l'ajoute à `pp` et on élimine ses multiples de `mm`, jusqu'à ce\n",
    "que `mm` soit vide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def primes(n):\n",
    "    if n<2: return []\n",
    "    pp = [2]; mm = [m for m in range(3,n) if m%2]\n",
    "    while mm:\n",
    "        p = mm[0]\n",
    "        pp.append(p)\n",
    "        mm = [m for m in mm[1:] if m%p]\n",
    "    return pp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]\n"
     ]
    }
   ],
   "source": [
    "print (primes(1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9592"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ll = primes(100000)\n",
    "len(ll)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponentiation rapide\n",
    "\n",
    "A chaque étape, `m%2` est un bit de l'exposant $m$. Le suivant est obtenu en remplaçant `m` par `m//2 = m>>1`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "9>>1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def powermod(a,m,n):\n",
    "    res = 1; x = a\n",
    "    while m:\n",
    "        if m & 1: res = (res*x) % n\n",
    "        x = (x*x) % n\n",
    "        m >>=1           # m //= 2 \n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "powermod(3,2**107-2,2**107-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test par division\n",
    "\n",
    "Un entier $<30$ qui n'est pas premier est forcément divisible par 2, 3, ou 5. On repart de 7.\n",
    "\n",
    "Partant de cette observation, on pourrait améliorer l'algorithme pour faire moins de tests ([wheel factorization](https://en.wikipedia.org/wiki/Wheel_factorization)).\n",
    "\n",
    "*Exercice* : le faire ...\n",
    "\n",
    "*Solution* : fichier `ent3.py` disponible la semaine prochaine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trial_division(n):\n",
    "    if n%2 == 0: return 2\n",
    "    p = 3\n",
    "    while p*p <= n:\n",
    "        if n%p == 0: return p\n",
    "        p +=2\n",
    "    return n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "641"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trial_division(2**32+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Les 200000 premiers nombres premiers ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "ll = [2,3]\n",
    "p = 5\n",
    "i = 2\n",
    "while i<200000:\n",
    "    if trial_division(p) == p: \n",
    "        ll.append(p)\n",
    "        i+=1\n",
    "    p+=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2749919, 2749921, 2749991, 2750021, 2750029, 2750053, 2750071, 2750123, 2750131, 2750159]\n"
     ]
    }
   ],
   "source": [
    "print (ll[-10:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test de Fermat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import randrange\n",
    "\n",
    "\n",
    "def is_pseudo_prime(a,n):\n",
    "    return powermod(a,n-1,n)==1\n",
    "\n",
    "\n",
    "def fermat(n,tests=4):\n",
    "    for i in range(tests):\n",
    "        a = randrange(2,n-1)\n",
    "        if not is_pseudo_prime(a,n): return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "mm = [2,3]\n",
    "p = 5\n",
    "i = 2\n",
    "while i<200000:\n",
    "    if fermat(p): \n",
    "        mm.append(p)\n",
    "        i+=1\n",
    "    p+=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = set(ll); M=set(mm)\n",
    "failed = M.difference(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{2113921, 1193221, 294409, 1909001, 950797, 1152271, 552721, 62745, 399001, 114589, 29341, 2465, 876961, 736291, 488881, 1569457, 252601, 1461241, 1615681, 46657, 2433601, 314821, 334153, 512461, 670033, 340561, 1857241, 15841, 35425, 75361, 162401, 63973, 525673, 2628073, 852841, 2508013, 91001}\n"
     ]
    }
   ],
   "source": [
    "print (failed)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "37"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(failed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test de Miller-Rabin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def test_base(a,n):\n",
    "    m = n-1; k = 0\n",
    "    while not m%2:\n",
    "        k+=1; m//=2\n",
    "    b = powermod(a,m,n)\n",
    "    if b==1: return True\n",
    "    for i in range(k):\n",
    "        x = b\n",
    "        b = pow(b,2,n)\n",
    "        if b==1:\n",
    "            if x != n-1: return False\n",
    "            else: return True\n",
    "    return False\n",
    "\n",
    "def miller_rabin(n, tests=4):\n",
    "    if n in [2,3]: return True\n",
    "    if not n%2: return False\n",
    "    for i in range(tests):\n",
    "        a = randrange(2,n-1)\n",
    "        if not test_base(a,n): return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn = [2,3]\n",
    "p = 5\n",
    "i = 2\n",
    "while i<200000:\n",
    "    if miller_rabin(p): \n",
    "        nn.append(p)\n",
    "        i+=1\n",
    "    p+=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set()\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "N = set(nn)\n",
    "failed = L.difference(N)\n",
    "print (failed)\n",
    "print (len(failed))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restes chinois"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "\n",
    "def solve_chinese(yy,mm):\n",
    "    assert len(yy) == len(mm), \"Arguments must have same length.\"\n",
    "    k = len(yy); m = reduce(lambda x,y:x*y, mm); x=0\n",
    "    for i in range(k):\n",
    "        u = reduce(lambda x,y:x*y, [mm[j] for j in range(k) if j!=i])\n",
    "        v = inversemod(u,mm[i])\n",
    "        x = (x + yy[i]*v*u) % m\n",
    "    return x\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "785"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve_chinese([3,4,5],[17,11,6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test de Solovay-Strassen\n",
    "\n",
    "Lorsque $n$\n",
    "(impair) n'est pas premier, $\\left(\\frac{m}{n}\\right)$ est appelé [symbole de Jacobi](https://fr.wikipedia.org/wiki/Symbole_de_Jacobi). Son calcul est basé sur la [loi de réciprocité quadratique](https://fr.wikipedia.org/wiki/Loi_de_r%C3%A9ciprocit%C3%A9_quadratique) de Gauss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jacobi(m,n):\n",
    "    assert n%2, \"n must be odd\"\n",
    "    m = m%n\n",
    "    if m == 0: return 0\n",
    "    if m == 1: return 1\n",
    "    \n",
    "    if m%2 == 0:\n",
    "        m //=2\n",
    "        if n%8 in [1,7]: return jacobi(m,n)\n",
    "        else: return -jacobi(m,n)\n",
    "    elif m%4 == 3 and n%4 == 3: return -jacobi(n,m)  \n",
    "    else: return jacobi(n,m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobi(1001,9907)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_sol_stra(a,n):\n",
    "    j = jacobi(a,n)\n",
    "    if j==0: return False\n",
    "    else: return (pow(a,(n-1)//2,n)-jacobi(a,n)) % n  == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 2**(2**5)+1\n",
    "test_sol_stra(3,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solovay_strassen(n, tests=4):\n",
    "    for i in range(tests):\n",
    "        a = randrange(2,n-1)\n",
    "        if not test_sol_stra(a,n): return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "aa = [2,3]\n",
    "p = 5\n",
    "i = 2\n",
    "while i<200000:\n",
    "    if solovay_strassen(p): \n",
    "        aa.append(p)\n",
    "        i+=1\n",
    "    p+=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set()\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "A = set(aa)\n",
    "failed = L.difference(A)\n",
    "print (failed)\n",
    "print (len(failed))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}