Commits

Gabriel Pichot committed bf71366

Graphe de l'état après transformation de Fourier

Comments (0)

Files changed (6)

 # }}}
 
 DIRS =
-SOURCES = log.ml matrixFactory.ml quantum.ml arithmetik.ml main.ml
+SOURCES = log.ml matrixFactory.ml quantum.ml print.ml arithmetik.ml main.ml
 EXEC = shor
 
 
 # The list of Caml libraries needed by the program
 
 # LIBS=$(WITHGRAPHICS) $(WITHUNIX) $(WITHSTR) $(WITHNUMS) $(WITHTHREADS) $(WITHDBM)
-
+LIBS=$(WITHGRAPHICS)
 # Should be set to -custom if you use any of the libraries above
 # or if any C code have to be linked with your program
 # (irrelevant for ocamlopt)
 open Log;;
 open Printf;;
 open Quantum;;
+open Print;;
 
 (* primeTest {{{1
  * Détermine si un nombre est premier ou non (exact).
   | 0 -> 0
   | n -> 1 + (nb_bits (n lsr 1))
 
-let order p n = 
+let order ~print ~texPrint p n = 
   log "\n---- Recherche de l'ordre de %i dans %i ----\n" p n; 
   let l = getQ n in
   let q = pow 2 l in
   (* ... mais on veut des probabilités donc on normalise *)
   reg2#normalize ();
   (* On mesure le second registre *)
-  let value = 2 in (*reg2#measureState () in *)
+  let value = 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 oublie pas de normaliser *)
   reg1#normalize ();
-  reg1#dump ();
   (* Maintenant on applique la transformée de Fourier *)
   log "Transformée de Fourier.\n";
-  let a = reg1#qft () in
-  Array.iteri (fun i j -> Printf.printf "Etat %i vaut %s.\n" i (c_to_string j)) a;
-  printf "vrai dft";
-  reg1#dft q;
-  reg1#dump ();
+  reg1#fft (); (* ou reg1#dft q() mais moins performant :) *)
   log "Fin de la transformée de Fourier.\n";
   reg1#normalize ();
+  (* Si on doit afficher le résultat de la transformée de Fourier *)
+  if print then showProb (reg1#state ());
+  if texPrint then printAsTex (reg1#state ()) p n;
   (* On mesure c sur le premier registre *)
   let c = reg1#measureState () in
   log "On trouve pour c : %i.\n" c;

latex/tipe_dossier_ens.tex

 % Just Fuck Listings Package which do not support UTF-8 O_o
 %\usepackage{listings}
 \usepackage{lgrind}
+\usepackage{tikz}
+\usetikzlibrary{calc}
 \usepackage{amsmath}
 \usepackage{amssymb}
 \usepackage{amsfonts}
 
 \subsection{Approximation de réels par des fractions continues}
 
+\input{../output/shor_11_15.tex}
+\input{../output/shor_8_15.tex}
+\input{../output/shor_7_99.tex}
+
 \appendix
 \section{Implémentation}
 \subsection{Fichier main.ml}
   end
 in
 printf "Le nombre qui va nous servir (et qui ne divise pas n) est : %i.\n" p;
+(* On peut aussi afficher les graphes d'après de transformée de Fourier avec
+ * l'option --print mise en dernière *)
+let printState = nb_args >= 4 && argv.(3) = "--print" in
+let texPrint   = nb_args >= 4 && argv.(3) = "--tex" in
 (* }}} *)
 
 (* Vérifications {{{1 *) 
 (* 1) N n'est pas pair *)
 if n mod 2 = 0 then begin
-  print "Erreur : Le nombre à factoriser doit être impair !";
+  print "Erreur : Le nombre à factoriser doit être impair !\n";
   exit 0
 end;
 (* 2) Si le nombre est premier *)
 if primeTest n then begin
-  print "Erreur : Le nombre ne doit pas être premier !";
+  print "Erreur : Le nombre ne doit pas être premier !\n";
   exit 0
 end;
 (* @TODO
  * 3) Teste si le nombre est la puissance d'un nombre premier (ne fait rien 
  * pour le moment en vérité) *)
 if primePowTest n then begin
-  print "Erreur : le nombre ne doit être la puisance d'un nombre premier !";
+  print "Erreur : le nombre ne doit être la puisance d'un nombre premier !\n";
   exit 0
 end;
 (* Teste si 1 <= p <= n - 1 *)
 if 1 > p || p >= n then begin
-  printf "Erreur : l'entier p doit être compris entre 1 et le nombre à factoriser strictement";
+  printf "Erreur : l'entier p doit être compris entre 1 et le nombre à factoriser strictement.\n";
   exit 0
 end;
 (* Teste si n et p sont bien premiers entre eux *)
 
 (* Ce qui suit peut paraître dérisoire mais c'est pourtant le coeur de notre
  * algorithme et tout notre sujet :) *)
-let r     = order p n in
+let r     = order ~print:printState ~texPrint p n in
 let ppow  = pow p (r / 2) in
 
   (* Tests sur l'ordre {{{1 *)
 (* Si r est impair *)
 if r mod 2 = 1 then begin
-  printf "Erreur : l'ordre est impair (%i), réessayez avec un autre nombre !" r;
+  printf "Erreur : l'ordre est impair (%i), réessayez avec un autre nombre !\n" r;
   exit 0
 end;
 (* Si r est congru a -1 modulo n *)
 if ppow mod n = n - 1 then begin
-  print "Erreur : x^(r/2) est congru à -1 modulo n !, un autre entier ?";
+  print "Erreur : x^(r/2) est congru à -1 modulo n !, un autre entier ?\n";
   exit 0
 end;
 (* }}} *)
+open Log;;
+open Graphics;;
+open Quantum;;
+
+let showProb u = 
+  let n = u#rows () in
+  let base = 800 / n in
+  let swidth = base * n
+  and sheight = 300 in
+  open_graph (Printf.sprintf " %ix%i" swidth sheight);
+  for i = 0 to n - 1 do
+    let height = int_of_float (Complex.norm2 (u#row (i + 1)) *. (float_of_int sheight)) in
+    fill_rect (i * base) 0 base height
+  done;
+  let _ = read_key () in
+  close_graph ()
+;;
+
+let printAsTex u n p =
+  (* On récupère la liste des données qui peuvent être vues *)
+  let rec obtainList = function
+    | 0 -> []
+    | k -> let y = Complex.norm2 (u#row k) in
+      if y > 0.0001 then (k,y) :: (obtainList (k - 1))
+      else obtainList (k - 1)
+  and max = List.fold_left (fun a b -> if snd a > snd b then a else b) (0,0.) in
+  let data = obtainList (u#rows ()) in
+  let max = snd (max data) in
+  let m = u#rows () in
+  let s = List.fold_left begin fun s (i,y) -> 
+    let x = float_of_int i *. 16. /. (float_of_int m) +. 0.01
+    and y = y /. max *. 2. in
+    s ^ (Printf.sprintf "(%f,%f)\n" x y) 
+  end "" data in
+  let file = open_out (Printf.sprintf "output/shor_%i_%i.tex" n p) in
+  let out = ref (Printf.sprintf "\\begin{tikzpicture}
+    \\draw (0,0) -- (0,3);
+    \\draw (0,0) -- (16,0);\n" ) in
+  for i = 0 to 8 do
+    out := !out ^ (Printf.sprintf "  \\draw(%i,0)node[below]{%i};\n" (i * 2) (i * m / 8))
+  done;
+  out := !out ^ "  \\foreach \\x in {0,2,...,16} \\draw (\\x,0) -- (\\x,0.2);
+    \\draw[very thick] plot[ycomb] coordinates {\n";
+  output_string file (!out ^ s ^ "};\n\\end{tikzpicture}")
+;;
+   
+
   val size = n
   val state = new vector ~rows:(1 lsl n) () (* Equivalent à 2^n *)
   method size () = size
+  method state () = state
   method nbStates () = 1 lsl n
   method norm () = state#norm ()
   method normalize () = state#normalize ()
     done;
     self#setState dftvals
   (* }}} *)
-  (* Transformée de Fourier rapide (car q est une puissance de 2 !!! {{{2 *)
-  method qft () =
-    let rec qft_aux ?(step=1) ?(start=0) () =
+  (* Transformée de Fourier rapide (car q est une puissance de 2 !!!) {{{2 *)
+  method fft () =
+    let rec fft_aux ?(step=1) ?(start=0) () =
       let n = state#rows () / step in
       if n = 1 then [| state#row (start + 1) |]
       else begin
-        let even = qft_aux ~step:(step * 2) ~start ()
-        and odd  = qft_aux ~step:(step * 2) ~start:(start + step) () in
+        let even = fft_aux ~step:(step * 2) ~start ()
+        and odd  = fft_aux ~step:(step * 2) ~start:(start + step) () in
         let c  = ref Complex.one
         and w  = w n
         and u' = Array.make n Complex.zero in
         done;
         u'
       end
-    in qft_aux ()
+    in self#setState (fft_aux ())
     (* }}} *)
   (* measureState {{{2 *)
   method measureState () =