Commits

Gabriel Pichot committed 197ad92

It Works !!! (15 8, or 27 3)

Comments (0)

Files changed (4)

  *)
 
 open Log;;
+open Printf;;
 open Quantum;;
 
 (* primeTest {{{1
   | _ -> pgcd b (a mod b)
 (* }}} *)
 (* isANonTrivialDivisor {{{1 *)
-let isANonTrivialDivisor p n = p <> 1 && p <> n && p mod n = 0
+let isANonTrivialDivisor p n = p <> 1 && p <> n && n mod p = 0
 (* }}} *)
 (* pow {{{1 *)
 let rec pow x n = match n with
       in expMod_aux (xi * xi mod n) (ai lsr 1) r
   in expMod_aux x a 1
 (* }}} *)
+(* Approximation du réél a par une fraction {{{ *) 
+(* Les fractions continues fournissent une bonne approximation, on veut
+ * toutefois un dénominateur inférieur à n *)
+let approx x q = 
+  let rec reduite (p0,p1) (q0,q1) r precis =
+    if abs_float (p1 /. q1 -. x) < precis then (p1, q1)
+    else begin
+      let a = floor r in
+      let d = r -. a in
+      reduite (p1, a *. p1 +. p0) (q1, a *. q1 +. q0) (1. /. d) precis
+    end
+  in let (p,q) = reduite (1., floor x) (0., 1.) (1. /. (x -. floor x)) (0.5 /. q)
+  in (int_of_float p, int_of_float q)
+(* }}} *)
 
 (* Recherche de l'ordre de p dans Z/nZ {{{1 *)
 (* On créé plusieurs petites fonctions auxiliaires *)
     p := !p * 2
   done;
   !q
+(* Donne le nombre de bits nécessaire pour représenter l'entier *)
+let rec nb_bits = function 
+  | 0 -> 0
+  | n -> 1 + (nb_bits (n lsr 1))
 
 let order p n = 
   log "order %i in %i\n" p n; 
   let l = getQ n in
   let q = pow 2 l in
-  log "q = %i = 2 ^ %i\n" q l;
+  log "On trouve q = %i = 2 ^ %i\n" q l;
   let reg1 = new register l
-  and reg2 = new register l in
+  and reg2 = new register (nb_bits n) in
+  printf "Création de deux registres de tailles respectives %i,  %i.\n" 
+    (reg1#size()) (reg2#size());
   (* On met le premier registre dans un état de superposition uniforme *)
   reg1#setUniformSuperposition q;
-  (* On calcule x ^ i mod n dans le second... *)
-  for i = 0 to reg2#nbStates () - 1 do
-    reg2#setStateProbability i (complex_of_int(expModulaire p i n))
+  (* On calcule x ^ a mod n pour tous les a *)
+  let expModTemp = Array.make n Complex.zero in
+  let sauvExpMod = Array.make q 0 in
+  for a = 0 to q - 1 do
+    let state = expModulaire p a n in
+    sauvExpMod.(a) <- state;
+    expModTemp.(state) <- expModTemp.(state) +! Complex.one
   done;
+  reg2#setState expModTemp;
   (* ... mais on veut des probabilités ! *)
   reg2#normalize ();
-
+  (* On mesure le second registre *)
+  let value = 2 in (*reg2#measureState () in *)
+  (* On doit répercuter cette mesure sur le premier registre *)
+  for a = 0 to q - 1 do
+    reg1#setStateProbability a (if sauvExpMod.(a) = value then Complex.one else Complex.zero)
+  done;
+  (* Et on oublis pas de normaliser *)
+  reg1#normalize ();
   (* Maintenant on applique la transformée de Fourier *)
   reg1#dft q;
-  let r = reg1#measureState () in
+  reg1#normalize ();
+  (* On mesure c sur le premier registre *)
+  (*reg1#dump ();*)
+  let c = reg1#measureState () in
+  printf "On trouve pour c : %i.\n" c;
+  (* On approche le réel grâce aux fractions continues *)
+  let s =  (float_of_int c) /. (float_of_int q) in
+  let (d,r) = approx s (float_of_int q) in
+  printf "On trouve pour approximation de %f ~ d / r = %i / %i.\n" s d r;
   r
   
 

latex/tipe_dossier_ens.tex

 \usepackage[utf8x]{inputenc}
 \usepackage[T1]{fontenc}
 \usepackage[french]{babel}
-
+\usepackage{braket}
 % Just Fuck Listings Package which do not support UTF-8 O_o
 %\usepackage{listings}
 \usepackage{lgrind}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{amsfonts}
+\usepackage{textcomp}
+\usepackage{stmaryrd}
+\usepackage{mathabx}
 \usepackage{xcolor}
 \usepackage{lscape}
 \usepackage[top=2cm,bottom=2cm,left=1.5cm,right=1.5cm]{geometry}
 
 %opening
-\title{}
-\author{}
+\title{L'algorithme de \textsc{Shor}}
+\author{Gabriel Pichot}
 
-\begin{document}
+\newcommand{\abs}[1]{\lvert #1 \rvert}
+\newcommand{\C}{\ensuremath{\mathbb{C}}}
+\newcommand{\intEntier}[1]{\ensuremath{\llbracket #1 \rrbracket}}
+\renewcommand{\O}{\ensuremath{\mathcal{O}}}
+\newcommand{\code}[1]{#1}
 
+\begin{document}
 \maketitle
 
 \begin{abstract}
-
+Le présent document présente les éléments principaux intervenant dans
+l'algorithme de \textsc{Shor} et donc dans la recherche de l'ordre d'un élément.
+Une implémentation en langage OCaml est proposée.
 \end{abstract}
 
-\section{}
+\newpage
+Trouver des algorithmes pour factoriser des nombres quelconques est en soi une
+problématique. On sait que tout nombre est se décompose en produit de facteur(s)
+premier(s). Le soucis est bien entendu de trouver ses facteurs. Pour des petits
+nombres on y arrive facilement, pour de plus grands par des propriétés
+remarquables (carré parfait, multiple de deux, congruences...). Cependant pour
+des grands nombre cela vient plus difficile un outil informatique
+
+\section{Tirer parti d'un calculateur quantique}
+\subsection{Classique et quantique}
+Dans un ordinateur classique, l'information minimale est d'ordinaire stockée
+sous la forme d'un bit qui peut prendre soit la valeur $0$ soit la valeur $1$.
+Le nombre de bits requis pour stocker une dépendant alors de la manière de
+représenter cette information. 
+De la même façon, dans un ordinateur quantique on considère des qubits (pour
+quantum bit). Ceux-ci, a contrario des précédents sont dans une superposition
+des deux états possibles, $0$ et $1$. Ainsi l'information codée n'est plus
+stockée de façon binaire soit $0$ soit $1$ mais au contraire comme composition
+de ces deux composantes. Ainsi un qubit peut aisément être représenté par un
+vecteur de $\C^2$. Appelant $\ket 0$ et $\ket 1$ les deux éléments formant une
+base de cet espace on obtient que tout qubit peut s'écrire formellement :
+\[ a \cdot \ket 0 + b \cdot \ket 1 \text{ avec } a,b\in\C \text{ et } \abs a
+^2 + \abs b ^2 = 1 \]
+Cette dernière condition traduit une condition supplémentaire : en vérité, on
+ne peut pas accéder à l'information stockée dans un qubit. Ainsi il nous est
+tout simplement impossible de mesurer les deux composantes, car lors de la
+mesure le qubit se \og fige \fg{} dans un état, soit $0$ soit $1$. Cet état
+dépend tout de même de $a$ et de $b$, dont les modules au carré représentent la
+probabilité de tomber dans cet état.
+Par extension, on parle de registre pour un ensemble de qubits, cette fois-ci
+l'état du registre est donc un vecteur, dans le cas d'un registre de $n$
+qubits, de $\C^{2n}$, ce qui induit donc la représentation formelle :
+\[ \sum_{i=0}^{n-1} a_i \ket {U_i} \text{ avec } \forall i \in
+\intEntier{0,n-1} , a_i \in\C \text{ et } \sum_{i=0}^{n-1} \abs{a_i} ^2 =
+1 \]
+
+Si on apprends à l'école à diviser des nombres on s'arrange vite 
+\subsection{Intrication quantique}
+\section{À la recherche de l'ordre}
+\subsection{Exponentiation modulaire}
+Dans le second registre, on doit calculer $x ^ a \mod n$, pour cela il existe
+plusieurs méthodes certaines plus astucieuses que d'autres. Une première
+méthode plutôt naïve consisterait à calculer les résidus modulo $n$ successifs
+des puissances de $x$ successives (on élimine directement l'idée de calculer
+$x^a$ puis son résidu, car on peut facilement dépasser la limite des entiers
+définie pour le langage (\code{max\_int} en OCaml)). Ce qui donne un
+algorithme d'un ordre de complexité en $\O(a)$. Cependant une solution beaucoup
+plus efficace est de considérer un algorithme d’exponentiation rapide.
+Si on écrit $a = \sum_{i=0}^{n} a_i 2^i$, c'est à dire $a$ dans son écriture
+binaire $a_n a_{n-1} \dots a_1 a_0$. On en déduit donc 
+\[ x^a = x ^ {\sum_{i=0}^{n} a_i 2^i} = \prod_{i=0}^{n} x ^ {a_i 2^i}
+ =  \prod_{i=0}^{n} \left( x ^ {2^i}\right) ^{a_i} \]
+L'opération de congruence possédant des propriétés remarquables par rapport au
+produit (et donc au puissance), on remarque qu'il suffit de considérer
+uniquement des puissances de $x$ de la forme $x ^ {2^i}$. Encore mieux, lorsque
+$a_i = 0$ le résidu est simplement $1$. Par suite il vient l'algorithme suivant
+:
+\begin{verbatim}
+let expModulaire x a n =
+  let residu  = ref 1
+  and a_bin   = ref a
+  and x_puis  = ref x in
+  while !a_bin > 0 do
+    if a_bin land 1 = 1 then begin(* Si le bit de poids faible est 1 *)
+      residu := !residu * x_puis % n
+    end;
+    x_puis := !x_puis * x_puis % n; (* Contient les puissance de 2 de x *)
+    a_bin  := a_bin lsr 1; (* On décale vers la droite (supprime le bit de
+poids faible *)
+  done;
+  !residu
+\end{verbatim}
+Cette algorithme présente l'avantage d'être très rapide ainsi pour des entiers
+codés sur $32$ bits, on fait dans le pire des cas $32$ tours de boucle. En
+fait, on montre (même principe que pour l'exponentiation rapide) que la
+fonction ci-dessus a une complexité en $\O(log_2(a))$. 
+
+\subsection{Approximation de réels par des fractions continues}
 
+\appendix
+\section{Implémentation}
+\subsection{Fichier main.ml}
 \lgrindfile{main.tex}
+\subsection{Module Arithmetik}
 \lgrindfile{arithmetik.tex}
+\subsection{Module Quantum.ml}
 \lgrindfile{quantum.tex}
+\subsection{Module MatrixFactory}
 \lgrindfile{matrixFactory.tex}
 
 \end{document}
 
 (* Initialisation (debug, args, random...) {{{1 *)
 print "Simulation de l'algorithme de Shor\n";
-Random.init 42;
+Random.self_init ();
 set_debug true;
 if !mod_debug then print "Mode DEBUG : on.\n";
 
 let r     = order p n in
 let ppow  = pow p (r / 2) in
 
-(* Tests sur l'ordre {{{1 *)
+  (* Tests sur l'ordre {{{1 *)
 (* Si r est impair *)
 if r mod 2 = 1 then begin
-  print "Erreur : l'ordre est impair, réessayez avec un autre nombre !";
+  printf "Erreur : l'ordre est impair (%i), réessayez avec un autre nombre !" r;
   exit 0
 end;
 (* Si r est congru a -1 modulo n *)
 let r1 = pgcd (ppow - 1) n
 and r2 = pgcd (ppow + 1) n in
 if isANonTrivialDivisor r1 n then
-  printf "\nSuccès : %i divise %i." r1 n
+  printf "\nSuccès : %i = %i x %i.\n" n r1 (n / r1)
 else if isANonTrivialDivisor r2 n then
-  printf "\nSuccès : %i divise %i." r2 n
+  printf "\nSuccès : %i = %i x %i.\n" n r2 (n / r2)
 else
-  print "Impossible de factoriser ce nombre, désolé.";
+  printf "Impossible de factoriser ce nombre, désolé (r = %i).\n" r;
 (* }}} *)
 
 
   Complex.im = 0.;
 }
 let ( *! ) = Complex.mul
+let ( +! ) = Complex.add
 (* }}} *)
 let pi = acos (-1.);;
+let c_to_string = fun (c:Complex.t) ->
+      Printf.sprintf "%.9F + %.9FI " c.Complex.re c.Complex.im
 
 module ComplexMatrix = MatrixFactory.Make( 
   struct
     let mul a b = Complex.mul a (complex_of_float b) 
     let add = Complex.add
     let norm2 = Complex.norm2
-    let to_string = fun (c:Complex.t) ->
-      Printf.sprintf "%.2f + %.2fI " c.Complex.re c.Complex.im
+    let to_string = c_to_string
   end
 );;
 open ComplexMatrix;;
 
 (* Register Class {{{1 *)
 exception Quant_Bad_Access of string;;
-class register n = object
+class register n = object(self)
   val size = n
   val state = new vector ~rows:(1 lsl n) () (* Equivalent à 2^n *)
   method size () = size
   method nbStates () = 1 lsl n
+  method norm () = state#norm ()
   method normalize () = state#normalize ()
+  method setState s = 
+    if Array.length s > state#rows () then raise (Quant_Bad_Access "setState")
+    else begin
+      for i = 0 to Array.length s - 1 do
+        state#rowset (i + 1) s.(i)
+      done;
+    end
   method setStateProbability s v = state#rowset (s+1) v
   method getStateProbability s = 
     if s > state#rows () then raise (Quant_Bad_Access "getStateProbability")
     else state#row (s+1)
   method dump () =
-    print "Le registre est dans l'état :\n";
+    printf "Le registre est dans l'état (norme %f):\n" (self#norm () );
     for i = 1 to state#rows () do
-      printf "État %i : %f.\n" (i-1) (Complex.norm(state#row i));
+      printf "État %i : %s.\n" (i-1) (c_to_string (state#row i));
     done
   (* Met les n premiers états dans un état de superposition uniforme *)
   method setUniformSuperposition n =
     end
   (* Transformation de Fourier discrète sur les q premiers états {{{2 *)
   method dft q = 
-    let d     = complex_of_float(float_of_int q ** (-0.5)) in
-    let cvals = Array.make q Complex.zero in 
-    for a = 1 to n do
-      for c = 1 to n do
-        cvals.(c) <- Complex.add cvals.(c) (
+    let d       = complex_of_float(float_of_int q ** (-0.5)) in
+    let dftvals = Array.make q Complex.zero in 
+    for a = 0 to q - 1 do
+      for c = 0 to q - 1 do
+        dftvals.(c) <- dftvals.(c) +! (
           Complex.exp (
-            (complex_of_float(2. *. pi /. (float_of_int q) *. float_of_int(a * c)  ) )
-            *! Complex.i
-          )
+            (complex_of_float(2. *. pi /. (float_of_int q) *. float_of_int(a * c)  ) ) *! Complex.i
+          ) *! d *! (state#row (a + 1))
         )
       done
     done;
-    for c = 1 to n do
-      state#rowset c (state#row c *! cvals.(c) *! d)
-    done
+    self#setState dftvals
   (* }}} *)
   (* measureState {{{2 *)
   method measureState () =
         top := !top +. norm;
         if !bottom < alea && alea < !top then begin
           measured := true;
-          stateMeasured := i;
+          stateMeasured := i - 1;
           state#rowset i Complex.one
         end else begin state#rowset i Complex.zero end;
         bottom := !bottom +. norm