{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ed824e72",
   "metadata": {},
   "source": [
    "# Séance 8 : factorisation d'entiers -- crible quadratique\n",
    "\n",
    "\n",
    "## Exercice\n",
    "\n",
    "Dans cet exercice, on propose d'implémenter (de façon naïve et non-optimisée) l'algorithme de crible quadratique qui trouve un diviseur propre d'un entier $N$.\n",
    "\n",
    "\n",
    "**Question 1.** Écrire une fonction `calcule_base_facteurs(N, B)` qui calcule une base de facteurs plus petits que `B` dans le but de factoriser un entier `N`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70ca831c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcule_base_facteurs(N, B):\n",
    "    L = [ 2 ]\n",
    "    p = 3\n",
    "    while p <= B:\n",
    "        if jacobi_symbol(N, p) == 1:\n",
    "            L.append(p)\n",
    "        p = next_prime(p)\n",
    "    return L\n",
    "    \n",
    "print(calcule_base_facteurs(N=369713, B=21))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6df0ba3",
   "metadata": {},
   "source": [
    "**Question 2.** Écrire une fonction de `quad_sieve(FB, N, A)` qui prend en entrée une base de facteurs `FB`, un entier `N` à factoriser et une borne `A`, et qui effectue l'étape d'effritement (crible)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a836e22f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def quad_sieve(FB, N, A): \n",
    "    S = ceil(sqrt(N))\n",
    "    Q = lambda x: x**2 + 2*S*x + S**2 - N\n",
    "    T = [ Q(x) for x in range(A) ]\n",
    "    for p in FB:\n",
    "        a = mod(N, p)\n",
    "        if a.is_square():\n",
    "            roots = a.sqrt(all=True)\n",
    "            roots = [Integer(r-S) for r in roots ]\n",
    "            roots = [r for r in roots if r < A ]\n",
    "            e = 1\n",
    "            while roots != []:\n",
    "                for x in roots:\n",
    "                    while x < A:\n",
    "                        T[x] = T[x]//p\n",
    "                        x = x + p**e\n",
    "                e += 1\n",
    "                a = mod(N, p**e)\n",
    "                roots = a.sqrt(all=True)\n",
    "                roots = [Integer(r-S) for r in roots ]\n",
    "                roots = [r for r in roots if r < A ]\n",
    "    return [ [x,Q(x)] for x in range(A) if T[x] == 1 ]\n",
    "    \n",
    "N = 369713\n",
    "FB = calcule_base_facteurs(N=369713, B=21)\n",
    "A = ceil(sqrt(N))//3\n",
    "L = quad_sieve(FB, N, A)\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d205573c",
   "metadata": {},
   "source": [
    "**Question 3.** Écrire une fonction de `build_matrix(FB, sieve)` qui prend en entrée une base de facteurs `FB` et une liste d'entiers effrités `sieve`, et qui construit la matrice associée."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9627d5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_matrix(FB, res_sieve):\n",
    "    n = len(FB)\n",
    "    m = len(res_sieve)\n",
    "    M = zero_matrix(ZZ, m, n)\n",
    "    for i in range(m):\n",
    "        r = res_sieve[i][1]\n",
    "        for j in range(n):\n",
    "            p = FB[j]\n",
    "            while (r != 0) and ((r % p) == 0):\n",
    "                M[i,j] = M[i,j] + 1\n",
    "                r = r//p\n",
    "    return M\n",
    "    \n",
    "M = build_matrix(FB, L)\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b51baed1",
   "metadata": {},
   "source": [
    "**Question 4.** Assembler les fonctions pour obtenir une fonction qui, étant donné un entier `N` en entrée, retourne un diviseur de `N` grâce à l'algorithme du crible quadratrique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e895ded",
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadratic_sieve(N, verbose=False):\n",
    "\n",
    "    # Initialisation des constantes\n",
    "    sqN = ceil(sqrt(N))\n",
    "    Q = lambda x: x**2 + 2*sqN*x + sqN**2 - N\n",
    "    A = sqN//3 # A = sqN\n",
    "    B = floor(2**(0.5*sqrt(log(N, 2) * log(log(N, 2), 2))))\n",
    "    B_large_enough = False\n",
    "    \n",
    "    # Collecte de relations (en avoir + que de facteurs)\n",
    "    #  et solution du système \n",
    "    while not(B_large_enough):\n",
    "        FB = calcule_base_facteurs(N, B)\n",
    "        s = len(FB)\n",
    "        sieve = quad_sieve(FB, N, A)\n",
    "        M = build_matrix(FB, sieve)\n",
    "        M2 = matrix(GF(2), M)\n",
    "        B_large_enough = (M2.kernel().dimension() > 0)\n",
    "        if not(B_large_enough) and verbose:\n",
    "            print(\"on augmente B\")\n",
    "            B *= 2\n",
    "        \n",
    "    # Reconstruction d'un diviseur\n",
    "    factor_found = False\n",
    "    while not(factor_found):\n",
    "        u = M2.left_kernel().random_element()\n",
    "        a = 1\n",
    "        b = 1 \n",
    "        vec = vector(ZZ, [0 for j in range(s)])\n",
    "        for i in range(len(u)):\n",
    "            if u[i] == 1:\n",
    "                a = (a * (sieve[i][0] + sqN)) % N\n",
    "                vec += M[i]\n",
    "        vec = vec/2\n",
    "        for j in range(s):\n",
    "            b = (b * FB[j]**vec[j]) % N\n",
    "        p, q = gcd(a-b, N), gcd(a+b, N)\n",
    "        factor_found = (p not in  [1, N]) or (q not in  [1, N])\n",
    "        if not(factor_found) and verbose:\n",
    "            print(\"facteur non trouvé, on réessaie\")\n",
    "    return p, q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c757445",
   "metadata": {},
   "source": [
    "**Question 5.** Testez votre fonction avec $N =$ par exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44c2c904",
   "metadata": {},
   "outputs": [],
   "source": [
    "# N1 = 73217\n",
    "N2 = 369713\n",
    "# N3 = 7029871\n",
    "\n",
    "# print(quadratic_sieve(N1))\n",
    "print(quadratic_sieve(N2))\n",
    "# print(quadratic_sieve(N3))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Sagemath",
   "language": "python",
   "name": "sagemath"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}