Commits

Gabriel Pichot committed fc60c12

Création du dépôt (make ne marchera pas)

Comments (0)

Files changed (8)

+
+#	{{{ Documentation
+# Advanced usage:
+# ---------------
+
+# If you want to fix the name of the executable, set the variable
+# EXEC, for instance, if you want to obtain a program named my_prog:
+# EXEC = my_prog
+
+# If you need special libraries provided with the Caml system,
+# (graphics, arbitrary precision numbers, regular expression on strings, ...),
+# you must set the variable LIBS to the proper set of libraries. For
+# instance, to use the graphics library set LIBS to $(WITHGRAPHICS):
+# LIBS=$(WITHGRAPHICS)
+
+# You may use any of the following predefined variable
+# WITHGRAPHICS : provides the graphics library
+# WITHUNIX : provides the Unix interface library
+# WITHSTR : provides the regular expression string manipulation library
+# WITHNUMS : provides the arbitrary precision arithmetic package
+# WITHTHREADS : provides the byte-code threads library
+# WITHDBM : provides the Data Base Manager library
+#
+#
+# }}}
+
+DIRS =
+SOURCES = log.ml arithmetik.ml matrixFactory.ml misc.ml quantum.ml main.ml
+EXEC = shor
+
+
+# {{{ The Caml compilers.
+CAMLC = ocamlc
+CAMLOPT = ocamlopt
+CAMLDEP = ocamldep
+CAMLLEX = ocamllex
+CAMLYACC = ocamlyacc
+# }}}
+
+# The list of Caml libraries needed by the program
+
+# LIBS=$(WITHGRAPHICS) $(WITHUNIX) $(WITHSTR) $(WITHNUMS) $(WITHTHREADS) $(WITHDBM)
+
+# 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)
+
+CUSTOM=-custom
+
+# {{{ WITH Settings
+# Default setting of the WITH* variables. Should be changed if your
+# local libraries are not found by the compiler.
+WITHGRAPHICS =graphics.cma -cclib -lgraphics -cclib -L/usr/X11R6/lib -cclib -lX11
+WITHUNIX =unix.cma -cclib -lunix
+WITHSTR =str.cma
+WITHNUMS =nums.cma -cclib -lnums
+WITHTHREADS =threads.cma -cclib -lthreads
+WITHDBM =dbm.cma -cclib -lmldbm -cclib -lndbm
+WITHXMLLIGHT =-I +xml-light xml-light.cma
+WITHLABLGTK =-I +lablgtk2 lablgtk.cma gtkInit.cmo
+WITHGLADE =-I +lablgtk2 lablglade.cma
+# }}}
+
+################ This part should be generic
+################ Nothing to set up or fix here
+
+# {{{
+
+all: depend $(EXEC)
+
+opt : $(EXEC).opt
+
+#ocamlc -custom other options graphics.cma other files -cclib -lgraphics -cclib -lX11
+#ocamlc -thread -custom other options threads.cma other files -cclib -lthreads
+#ocamlc -custom other options str.cma other files -cclib -lstr
+#ocamlc -custom other options nums.cma other files -cclib -lnums
+#ocamlc -custom other options unix.cma other files -cclib -lunix
+#ocamlc -custom other options dbm.cma other files -cclib -lmldbm -cclib -lndbm
+
+SOURCES1 = $(SOURCES:.mly=.ml)
+SOURCES2 = $(SOURCES1:.mll=.ml)
+OBJS = $(SOURCES2:.ml=.cmo)
+OPTOBJS = $(SOURCES2:.ml=.cmx)
+
+$(EXEC): $(OBJS)
+	$(CAMLC) $(CUSTOM) -o $(EXEC) $(LIBS) $(OBJS)
+
+$(EXEC).opt: $(OPTOBJS)
+	$(CAMLOPT) -o $(EXEC) $(LIBS:.cma=.cmxa) $(OPTOBJS)
+
+.SUFFIXES:
+.SUFFIXES: .ml .mli .cmo .cmi .cmx .mll .mly
+
+.ml.cmo:
+	$(CAMLC) $(LIBS) -annot -c $<
+
+.mli.cmi:
+	$(CAMLC) -c $<
+
+.ml.cmx:
+	$(CAMLOPT) -c $<
+
+.mll.cmo:
+	$(CAMLLEX) $<
+	$(CAMLC) -c $*.ml
+
+.mll.cmx:
+	$(CAMLLEX) $<
+	$(CAMLOPT) -c $*.ml
+
+.mly.cmo:
+	$(CAMLYACC) $<
+	$(CAMLC) -c $*.mli
+	$(CAMLC) -c $*.ml
+
+.mly.cmx:
+	$(CAMLYACC) $<
+	$(CAMLOPT) -c $*.mli
+	$(CAMLOPT) -c $*.ml
+
+.mly.cmi:
+	$(CAMLYACC) $<
+	$(CAMLC) -c $*.mli
+
+.mll.ml:
+	$(CAMLLEX) $<
+
+.mly.ml:
+	$(CAMLYACC) $<
+
+clean:
+	rm -f *.cm[iox] ui/*.cm[iox] *.annot .*~ .*~ #*#
+	rm -f $(EXEC)
+	rm -f $(EXEC).opt
+
+depend: $(SOURCES2)
+	$(CAMLDEP) *.mli *.ml > .depend
+
+include .depend
+
+# }}}
+# Une implémentation de l'algorithme de Shor
+
+Vous trouverez ci-joint une implémentation de l'algorithme de Shor
+que j'ai écrite dans le cadre de mes TIPE en OCaml.
+
+## Liste des fichiers
+
+ - main.ml : fichier principal fait des vérifications, et autres 
+interactions avec l'utilisateur.
+ - misc.ml : contient diverses fonctions peu utilisées
+ - arithmetik.ml : toutes les fonctions arithmétiques dont la 
+recherche de l'ordre.
+ - quantum.ml : implémentation de registre quantique.
+ - matrixFactory.ml : module de matrice générique.
+ - log.ml : gestion des entrées et sorties utilisateurs, spécifie
+un mode de debogagge ou non (très simplet).
+
+## Makefile
+
+Pour générer l'exécutable il suffit de faire un :
+    
+    make
+
+suivit d'un :
+
+    ./shor
+
+pour exécuter l'algorithme, noter toutefois qu'il est aussi possible
+de rentrer directement les deux nombres choisis directement en entrée :
+
+    ./shor 15 8
+
+par exemple pour factoriser 15 quand le nombre dont on chercher l'ordre 
+est 8.
+(*
+ * Ce fichier contient sans doute la partie la plus intéressante :
+ * recherche de l'ordre
+ * puis diverses fonctions arithmétiques : pgcd, test si premier...
+ *)
+
+open Log;;
+
+(* primeTest {{{1
+ * Détermine si un nombre est premier ou non (exact).
+ * Renvoit true s'il l'est *)
+let primeTest n = 
+  let racine      = int_of_float (sqrt (float_of_int n)) in
+  let isNotPrime  = ref (n mod 2 = 0 && n <> 2) in
+  let j           = ref 1 in
+  while not !isNotPrime && !j <= racine / 2 do
+    if n mod (2 * !j + 1) = 0 then 
+      isNotPrime := true;
+    incr j
+  done;
+  not !isNotPrime
+(* }}} *)
+(* primePowTest {{{1 @TODO
+ * Teste si l'entier n passé en argument n'est pas la puissance d'un nombre
+ * premier *)
+let primePowTest n = false
+(* }}} *)
+(* Pgcd {{{1 *)
+let rec pgcd a b = match b with 
+  | 0 -> a
+  | _ -> pgcd b (a mod b)
+(* }}} *)
+(* isANonTrivialDivisor {{{1 *)
+let isANonTrivialDivisor p n = p <> 1 && p <> n && p mod n = 0
+(* }}} *)
+(* pow {{{1 *)
+let rec pow x n = match n with
+  | 0 -> 1
+  | _ -> let r = pow x (n / 2) in
+    if n mod 2 = 0 then r * r else r * r * x
+(* }}} *)
+
+(* Recherche de l'ordre de p dans Z/nZ {{{1 *)
+(* On créé plusieurs petites fonctions auxiliaires *)
+(* getQ : récupère l'entier tel que n^2 <= q < 2n^2 *)
+let getQ n = 
+  let n2 = n * n and q = ref 1 in
+  let p = ref 2 in
+  while !p < n2 do
+    incr q;
+    p := !p * 2
+  done;
+  !q
+
+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;
+  l
+open Printf;;
+
+
+
+(** Vars **)
+(* Open the logfile *)
+let logfile = stdout (* open_out "shor.log" *)
+(* Enabled debug mod *)
+let mod_debug = ref false
+
+(** Methods **)
+(* Active, unactive debug mode *)
+let set_debug b = mod_debug := b
+(* Log the fmt into the file *)
+(* let log fmt = fprintf logfile fmt *)
+let log fmt = printf fmt
+(* Log the debug fmt if necessary *)
+let debug fmt = if !mod_debug then fprintf logfile fmt else ifprintf logfile fmt
+(* Show in the console (do not forget to use the ! to flush the data) *)
+let print = printf "%s%!" 
+(*
+ * Ce fichier contient le squelette de l'application, vérifications etc.
+ * @TODO :
+ *  - faire la fonction primePowTest
+ *  -? ne pas obliger l'utilisateur à rentrer l'entier p. 
+ *)
+open Log;;        (* Log, debug *)
+open Quantum;;  
+open Sys;;        (* Arguments *)
+open Misc;;       (* Scanf *)
+open Printf;;     (* Printf *)
+open Arithmetik;; (* primeTest *)
+
+(* Initialisation (debug, args, random...) {{{1 *)
+print "Simulation de l'algorithme de Shor\n";
+Random.init 42;
+set_debug true;
+if !mod_debug then print "Mode DEBUG : on.\n";
+
+
+(* Le nombre que l'on souhaite factoriser *)
+let nb_args = Array.length argv in
+let n = 
+  if nb_args >= 2 then int_of_string argv.(1)
+  else begin
+    print "Veuillez entrer un nombre à factoriser : ";
+    read_int ();
+  end
+in
+printf "Le nombre à factoriser est %i.\n" n;
+(* Un nombre qui ne divise pas n *)
+let p = 
+  if nb_args >= 3 then int_of_string argv.(2)
+  else begin
+    print "Veuillez entrer un nombre ne divisant pas n : "; 
+    read_int ();
+  end
+in
+printf "Le nombre qui va nous servir (et qui ne divise pas n) est : %i.\n" p;
+(* }}} *)
+
+(* Vérifications {{{1 *) 
+(* 1) N n'est pas pair *)
+if n mod 2 = 0 then begin
+  print "Erreur : Le nombre à factoriser doit être impair !";
+  exit 0
+end;
+(* 2) Si le nombre est premier *)
+if primeTest n then begin
+  print "Erreur : Le nombre ne doit pas être premier !";
+  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 !";
+  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";
+  exit 0
+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;
+  exit 0
+end;
+(* }}} *)
+
+(* 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 ppow  = pow p (r / 2) in
+
+(* 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 !";
+  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 ?";
+  exit 0
+end;
+(* }}} *)
+
+(* It's done !! {{{1 *)
+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
+else if isANonTrivialDivisor r2 n then
+  printf "\nSuccès : %i divise %i." r2 n
+else
+  print "Impossible de factoriser ce nombre, désolé.";
+(* }}} *)
+
+
+
+let rec iterator ?(i=0) n f acc = match i with
+  | _ when i > n || i = n -> acc
+  | _ -> iterator ~i:(i+1) n f (f i acc)
+let rec viterator ?(i=0) n f = match i with
+  | _ when i > n || i = n -> ()
+  | _ -> viterator ~i:(i+1) n f
+
+(* Ev Module {{{1 *)
+module type K = sig
+  type t
+  val zero : t
+  val one : t
+  val prod : t -> t -> t
+  val div : t -> t -> t
+  val sqrt : t -> t 
+  val add : t -> t -> t
+  val norm2 : t -> float
+  val to_string : t -> string
+end
+(* }}} *)
+
+(* Matrix Module {{{1 *)
+module Make = functor (Ev : K) ->  struct
+  (* Exceptions {{{2 *)
+  exception Dims_Not_Matching of string
+  (* }}} *)
+  (* Matrix Class {{{2 *)
+  class matrix ?data ?(dims=(0,0)) () = object
+    val rows = match data with | Some x -> Array.length x     | None -> fst dims
+    val cols = match data with | Some x -> Array.length x.(0) | None -> snd dims
+    val data = match data with 
+      | Some x -> x 
+      | None -> Array.make_matrix (fst dims) (snd dims) Ev.zero
+    method rows () = rows
+    method cols () = cols
+    method get row col = data.(row-1).(col-1)
+    method set row col value = data.(row-1).(col-1) <- value
+    method to_string () =
+      Array.fold_left begin fun s row ->
+        s ^ (Array.fold_left (fun t coef -> t ^ " " ^ (Ev.to_string coef)) "(" row) ^ " )\n"
+      end "" data
+  end
+  (* }}} *)
+  (* Square Matrix Class {{{2 *)
+  class square_matrix ?data n = object(self)
+    inherit matrix ~dims:(n,n) ?data () as super
+    method trace () = iterator n (fun i sum -> Ev.add (self#get i i) sum) Ev.zero
+  end
+  (* }}} *)
+  (* Vector Class {{{2 *)
+  class vector ?data n = object(self)
+    inherit matrix ~dims:(n,1) ?data () as super
+    method row n = super#get n 0
+    method rowset n = super#set n 0
+    method norm () = sqrt (iterator n begin fun i sum -> 
+     (Ev.norm2 (self#row i)) +. sum
+    end 0.)
+    method normalize () = let norm = self#norm () in
+      viterator n (fun i -> self#rowset i (Ev.div (self#row i) norm))
+  end
+  (* }}} *)
+  (* Operators Module {{{2 *)
+  module Op = struct
+    let build_matrix_array f rows cols = 
+      Array.init rows begin fun i ->
+        Array.init cols (fun j -> f i j)
+      end
+
+    let add a b = if a#cols <> b#cols && a#rows <> b#rows then
+        raise (Dims_Not_Matching "Op.add")
+      else begin
+        new matrix ~data:begin
+          build_matrix_array (fun i j -> Ev.add (a#get i j) (b#get i j)) a#rows a#cols
+        end
+      end
+
+    let mul k a = new matrix ~data:begin
+      build_matrix_array (fun i j -> Ev.prod k (a#get i j)) a#rows a#cols
+    end
+    
+    let transpose a = new matrix ~data:begin
+      build_matrix_array (fun i j -> a#get j i) a#cols a#rows
+    end
+    
+  let prod a b = 
+    if a#cols <> b#rows && a#rows <> a#cols then
+      raise (Dims_Not_Matching "Op.prod")
+    else begin
+      new matrix ~data:begin
+        build_matrix_array begin fun i j -> 
+          iterator a#cols (fun k s -> Ev.add s (Ev.add (a#get i k) (b#get k j))) Ev.zero 
+        end a#rows b#cols
+      end
+    end
+  end
+  (* }}} *)
+end
+(* }}} *)
+
+
+(* Input reading *)
+let read_int () = Scanf.scanf " %i" (fun x -> x)
+
+(* Module Matrix : initialisation {{{1 *)
+module ComplexMatrix = MatrixFactory.Make (
+  struct
+    type t = Complex.t
+    let zero = Complex.zero
+    let one = Complex.one
+    let prod = Complex.mul
+    let div = Complex.div
+    let sqrt = Complex.sqrt
+    let add = Complex.add
+    let to_string = fun (c:Complex.t) ->
+      Printf.sprintf "%.2f + %.2fI " c.Complex.re c.Complex.im
+  end
+);;
+(* }}} *)
+(* Quelques fonctions sur les complexes {{{1 *)
+let complex_of_float f = {
+  Complex.re = f;
+  Complex.im = 0.;
+}
+let ( *: )  = Complex.mul
+let ( /: )  = Complex.div
+let ( +: )  = Complex.add
+(* }}} *)
+(*
+(* Replace function {{{1 *)
+let str_replace search mixed = 
+  Str.global_replace (Str.regexp_string search) mixed
+in 
+(* Translator function *)
+let exec str l = 
+  let rec _exec str l i = 
+    match l with
+    | [] -> str
+    | t :: q -> 
+      let search = "%" ^ (string_of_int i) in
+      _exec (str_replace search t str) q (i+1)
+  in _exec str l 0
+and str = " %0 %1 %2 "
+in
+
+log (exec str ["test"; string_of_int 5; string_of_float 5.]);;
+(* }}} *)
+*)
+
+
+
+module ComplexMatrix = MatrixFactory.Make( 
+  struct
+    type t = Complex.t
+    let zero = Complex.zero
+    let one = Complex.one
+    let prod = Complex.mul
+    let div = Complex.div
+    let sqrt = Complex.sqrt
+    let add = Complex.add
+    let norm2 = Comlex.norm2
+    let to_string = fun (c:Complex.t) ->
+      Printf.sprintf "%.2f + %.2fI " c.Complex.re c.Complex.im
+  end
+);;
+open ComplexMatrix
+
+
+(* Register Class {{{1 *)
+exception State_Does_Not_Exist of string;;
+class register n = object
+  val size = n
+  val state = new vector (1 lsl n) (* Equivalent à 2^n *)
+  (* Non destructif *)
+  method getStateProbability s = 
+    if s >= state#rows () then raise State_Does_Not_Exist "getStateProbability"
+    else state#row s
+  method size () = size
+  (* Destructif *)
+  method measureState () =
+    if s >= state#rows () then raise State_Does_Not_Exist "measureState"
+    else begin
+      (* FIXME trop naif ? *)
+      let measured      = ref false
+      and stateMeasured = ref -1
+      and alea          = float_of_int (Random.int ()) /. float_of_int max_int 
+      and bottom        = ref 0. 
+      and top           = ref 0. in  
+      for i = 1 to state#rows do
+        if not measured then begin
+          let norm = (Complex.norm2 (state#row i)) in
+          top := !top +. norm;
+          if !bottom < alea && alea < !top then begin
+            measured := true;
+            stateMeasured := i
+            state#rowset i Complex.one
+          end else state#rowset i Complex.zero
+          bottom := !bottom +. norm 
+        end else state#rowset i Complex.zero
+      done;
+      !stateMeasured
+    end
+end
+(* }}} *)
+(*
+(* Circuit Class {{{1 *)
+class circuit n = object
+  (* Nombre de qubits du circuit *)
+  val nb_qubits = n
+  (* Portes du circuit, les circuits étant linéaires *)
+  val mutable gates = []
+  (* Ajout d'une porte *)
+  method add_gate (gate:gate) = gates = gate :: gates
+end
+(* }}} *)
+ *)