Source

shor-algo-ocaml / quantum.ml

open Log;;
open Printf;;
(* Quelques fonctions sur les complexes {{{1 *)
let complex_of_float f = {
  Complex.re = f;
  Complex.im = 0.;
}
and complex_of_int i = { 
  Complex.re = float_of_int i;
  Complex.im = 0.;
}
let ( *! ) = Complex.mul
(* }}} *)
let pi = acos (-1.);;

module ComplexMatrix = MatrixFactory.Make( 
  struct
    type t = Complex.t
    let zero = Complex.zero
    let one = Complex.one
    let prod = Complex.mul
    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
  end
);;
open ComplexMatrix;;


(* Register Class {{{1 *)
exception Quant_Bad_Access of string;;
class register n = object
  val size = n
  val state = new vector ~rows:(1 lsl n) () (* Equivalent à 2^n *)
  method size () = size
  method nbStates () = 1 lsl n
  method normalize () = state#normalize ()
  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";
    for i = 1 to state#rows () do
      printf "État %i : %f.\n" (i-1) (Complex.norm(state#row i));
    done
  (* Met les n premiers états dans un état de superposition uniforme *)
  method setUniformSuperposition n =
    if n > state#rows () then raise (Quant_Bad_Access "setUniformSuperposition")
    else begin
      let prob = complex_of_float (sqrt(1. /. (float_of_int n))) in
      for i = 1 to n do
        state#rowset i prob
      done
    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) (
          Complex.exp (
            (complex_of_float(2. *. pi /. (float_of_int q) *. float_of_int(a * c)  ) )
            *! Complex.i
          )
        )
      done
    done;
    for c = 1 to n do
      state#rowset c (state#row c *! cvals.(c) *! d)
    done
  (* }}} *)
  (* measureState {{{2 *)
  method measureState () =
    (* On procède ainsi, imaginons l'état suivant :
     *             00              01     10   11
     * |          0.5          |  0.25  | 0 | O.25 |
     * |  ------------ alea ------->|0.63  
     * La probabilité de chaque état est bien respecté   *)
    let measured      = ref false
    and stateMeasured = ref (-1)
    and alea          = float_of_int (Random.int max_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 begin state#rowset i Complex.zero end;
        bottom := !bottom +. norm 
      end else state#rowset i Complex.zero
    done;
    !stateMeasured
  (* }}} *)
end
(* }}} *)