Gabriel Pichot avatar Gabriel Pichot committed 065286f

iSauvegarde car marche mais on va passer par la QFT

Comments (0)

Files changed (4)

 # }}}
 
 DIRS =
-SOURCES = log.ml matrixFactory.ml testMatrix.ml quantum.ml arithmetik.ml main.ml
+SOURCES = log.ml matrixFactory.ml quantum.ml arithmetik.ml main.ml
 EXEC = shor
 
 
   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 *)
+(* Les fractions continues fournissent une bonne approximation *)
 let approx x q = 
   let rec reduite (p0,p1) (q0,q1) r precis =
     if abs_float (p1 /. q1 -. x) < precis then (p1, q1)
   | n -> 1 + (nb_bits (n lsr 1))
 
 let order p n = 
-  log "order %i in %i\n" 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
-  log "On trouve 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 (nb_bits n) in
-  printf "Création de deux registres de tailles respectives %i,  %i.\n" 
+  log "Création de deux registres de tailles respectives %i et %i.\n" 
     (reg1#size()) (reg2#size());
   (* On met le premier registre dans un état de superposition uniforme *)
   reg1#setUniformSuperposition q;
-  (* On calcule x ^ a mod n pour tous les a *)
+  (* On calcule x ^ a mod n pour tous les a, on doit cependant
+   * penser à les conserver (ce que ne peut pas faire un calculateur 
+   * quantique) pour les "superposer" avec le premier registre ensuite. *)
   let expModTemp = Array.make n Complex.zero in
   let sauvExpMod = Array.make q 0 in
   for a = 0 to q - 1 do
     expModTemp.(state) <- expModTemp.(state) +! Complex.one
   done;
   reg2#setState expModTemp;
-  (* ... mais on veut des probabilités ! *)
+  (* ... mais on veut des probabilités donc on normalise *)
   reg2#normalize ();
   (* On mesure le second registre *)
   let value = 2 in (*reg2#measureState () in *)
   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 *)
+  (* 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 ();
+  log "Fin de la transformée de Fourier.\n";
   reg1#normalize ();
   (* On mesure c sur le premier registre *)
-  (*reg1#dump ();*)
   let c = reg1#measureState () in
-  printf "On trouve pour c : %i.\n" c;
+  log "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;
+  log "---- Un ordre possible est donc %i. ----\n" r;
   r
   
 
 end;
 (* Teste si n et p sont bien premiers entre eux *)
 if pgcd n p <> 1 then begin
-  printf "Pseudo-Succès :  %i | %i (sans passer l'algorithme de Shor).\n" p n;
+  printf "Pseudo-Succès :  %i | %i (sans passer l'algorithme de Shor).\n" (pgcd p n) n;
   exit 0
 end;
 (* }}} *)
 }
 let ( *! ) = Complex.mul
 let ( +! ) = Complex.add
+let ( -! ) = Complex.sub
 (* }}} *)
 let pi = acos (-1.);;
 let c_to_string = fun (c:Complex.t) ->
       Printf.sprintf "%.9F + %.9FI " c.Complex.re c.Complex.im
+let w n = Complex.exp (
+  (complex_of_float(2. *. pi /. (float_of_int n)  ) ) *! Complex.i
+)
 
 module ComplexMatrix = MatrixFactory.Make( 
   struct
     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) () =
+      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 c  = ref Complex.one
+        and w  = w n
+        and u' = Array.make n Complex.zero in
+        for k = 0 to n / 2 - 1 do 
+          u'.(k) <- even.(k) +! !c *! odd.(k);
+          u'.(k + n / 2) <- even.(k) -! !c *! odd.(k);
+          c := !c *! w
+        done;
+        u'
+      end
+    in qft_aux ()
+    (* }}} *)
   (* measureState {{{2 *)
   method measureState () =
     (* On procède ainsi, imaginons l'état suivant :
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.