diff --git a/src/core/circuit.h b/src/core/circuit.h
index 9fe91147..a199c8b2 100644
--- a/src/core/circuit.h
+++ b/src/core/circuit.h
@@ -33,10 +33,10 @@ namespace qx
    {
       private:
 
-         uint64_t            n_qubit;
+         size_t              n_qubit;
          std::vector<gate *> gates;
          std::string         name;
-         uint64_t            iteration;
+         size_t              iteration;
          double              time;
 
       public:
@@ -44,7 +44,7 @@ namespace qx
          /**
           * \brief circuit constructor
           */
-         circuit(uint64_t n_qubit, std::string name = "", uint32_t iteration=1) : n_qubit(n_qubit), name(name), iteration(iteration)
+         circuit(size_t n_qubit, std::string name = "", size_t iteration=1) : n_qubit(n_qubit), name(name), iteration(iteration)
          {
          }
 
@@ -91,7 +91,7 @@ namespace qx
          /**
           * \brief return gate <i>
           */
-         qx::gate * operator [] (uint64_t i)
+         qx::gate * operator [] (size_t i)
          {
             assert(i < gates.size());
             return gates[i];
@@ -100,7 +100,7 @@ namespace qx
          /**
           * \brief set iterations number
           */
-         void set_iterations(uint64_t n)
+         void set_iterations(size_t n)
          {
             iteration = n;
          }
@@ -108,7 +108,7 @@ namespace qx
          /**
           * \brief set iterations number
           */
-         uint64_t get_iterations()
+         size_t get_iterations()
          {
             return iteration;
          }
@@ -116,30 +116,32 @@ namespace qx
          /**
           * \brief return gate <i>
           */
-         qx::gate * get(uint64_t i)
+         qx::gate * get(size_t i)
          {
             assert(i < gates.size());
             return gates[i];
          }
 
-         void execute(qu_register& reg, bool verbose=false, bool only_binary=false)
+         void execute(qu_register& reg, bool verbose=false, bool silent=false, bool only_binary=false)
          {
-            uint64_t it = iteration;
+            size_t it = iteration;
 
 #ifdef XPU_TIMER
-            println("[+] executing circuit '" << name << "' (" << it << " iter) ...");
             xpu::timer tmr;
-            tmr.start();
+            if (!silent)
+            {
+               println("[+] executing circuit '" << name << "' (" << it << " iter) ...");
+               tmr.start();
+            }
 #endif
-
             while (it--)
             {
                if (!verbose) 
-                  for (uint64_t i=0; i<gates.size(); ++i)
+                  for (size_t i=0; i<gates.size(); ++i)
                      gates[i]->apply(reg);
                else
                {
-                  for (uint64_t i=0; i<gates.size(); ++i)
+                  for (size_t i=0; i<gates.size(); ++i)
                   {
                      println("[-] executing gate " << i << "...");
                      gates[i]->dump();
@@ -149,15 +151,18 @@ namespace qx
                } 
             }
 #ifdef XPU_TIMER
-            tmr.stop();
-            println("[+] circuit execution time: " << tmr.elapsed() << " sec.");
+            if (!silent)
+            {
+               tmr.stop();
+               println("[+] circuit execution time: " << tmr.elapsed() << " sec.");
+            }
 #endif // XPU_TIMER
          }
 
          /**
           * \return gates count
           */
-         uint64_t size()
+         size_t size()
          {
             return gates.size();
          }
@@ -166,7 +171,7 @@ namespace qx
           * insert a gate at the specified 
           * position
           */
-         void insert(uint64_t pos, qx::gate * g)
+         void insert(size_t pos, qx::gate * g)
          {
             gates.insert(gates.begin()+pos,g);
          }
@@ -183,19 +188,19 @@ namespace qx
          void dump()
          {
             println("[+++] circuit '" << name << "' :");
-            for (uint64_t i=0; i<gates.size(); ++i)
+            for (size_t i=0; i<gates.size(); ++i)
             {
                print("  |--" << i << "--");
                gates[i]->dump();
             }
          }
 
-         uint64_t get_qubit_count()
+         size_t get_qubit_count()
          {
             return n_qubit;
          }
 
-         void set_qubit_count(uint64_t n)
+         void set_qubit_count(size_t n)
          {
             n_qubit = n;
             // check gate validity (target/ctrl qubits < n_qubit)
diff --git a/src/core/depolarizing_channel.h b/src/core/depolarizing_channel.h
index ea02aa20..4e805222 100644
--- a/src/core/depolarizing_channel.h
+++ b/src/core/depolarizing_channel.h
@@ -7,7 +7,7 @@
 namespace qx
 {
 
-   typedef std::vector<uint64_t>   histogram_t;
+   typedef std::vector<size_t>   histogram_t;
 
 
    #define __third__ (0.33333333333333333333333333333f)
@@ -24,427 +24,422 @@ namespace qx
       public:
 
         /**
-	 * ctor
-	 */
-        depolarizing_channel(qx::circuit * c, uint64_t nq, double pe) : nrg(0,__third__), c(c), nq(nq),
-							                pe(pe), 
-	                                                                xp(__third__), 
-							                yp(__third__), 
-							                zp(__third__),
-									overall_error_probability(0),
-									total_errors(0),
-									error_histogram(nq+1,0),
-									error_recording(false)
-	{
-
-	   for (uint64_t i=1; i<(nq+1); ++i)
-	   {
-	      double pp = error_probability(nq,i,pe);
-	      overall_error_probability += pp;
-	   }
-
-	   x_errors = 0;
-	   z_errors = 0;
-	   y_errors = 0;
-
-	   srand48(xpu::timer().current());
-	   qx::circuit * noisy_c = new qx::circuit(nq,c->id() + "(noisy)");
-	}
+	      * ctor
+	      */
+        depolarizing_channel(qx::circuit * c, size_t nq, double pe) : nrg(0,__third__), c(c), nq(nq),
+                                                                      pe(pe), 
+                                                                      xp(__third__), 
+                                                                      yp(__third__), 
+                                                                      zp(__third__),
+                                                                      overall_error_probability(0),
+                                                                      total_errors(0),
+                                                                      error_histogram(nq+1,0),
+                                                                      error_recording(false)
+         {
+            for (size_t i=1; i<(nq+1); ++i)
+            {
+               double pp = error_probability(nq,i,pe);
+               overall_error_probability += pp;
+            }
+
+            x_errors = 0;
+            z_errors = 0;
+            y_errors = 0;
+
+            srand48(xpu::timer().current());
+            qx::circuit * noisy_c = new qx::circuit(nq,c->id() + "(noisy)");
+         }
         
 
-	/**
-	 * ctor
-	 */
-	depolarizing_channel(qx::circuit * c, uint64_t nq, double pe, double xp, double yp, double zp) : nrg(0,__third__), c(c), nq(nq), 
-	                                                                                                 pe(pe), 
-											                 xp(xp), 
-											                 yp(yp), 
-											                 zp(zp),
-													 overall_error_probability(0),
-									                                 total_errors(0),
-													 error_histogram(nq+1,0),
-													 error_recording(false)
-	{
-	   for (uint64_t i=1; i<(nq+1); ++i)
-	   {
-	      double pp = error_probability(nq,i,pe);
-	      overall_error_probability += pp;
-	   }
-
-	   x_errors = 0;
-	   z_errors = 0;
-	   y_errors = 0;
-
-	   srand48(xpu::timer().current());
-	   qx::circuit * noisy_c = new qx::circuit(nq,c->id() + "(noisy)");
-	}
-
-	/**
-	 * check whether a qubit is idle
-	 */
-	bool is_idle(uint32_t q, std::vector<uint32_t> used)
-	{
-	   for (int i=0; i<used.size(); ++i)
-	      if (q == used[i])
-		 return false;
-	   return true;
-	}
-
-	/**
-	 * return idle qubits
-	 */
-	std::vector<uint32_t> idle_qubits(size_t n, std::vector<uint32_t> & used)
-	{
-	   std::vector<uint32_t> r;
-	   for (uint32_t i=0; i<n; i++)
-	   {
-	      if (is_idle(i,used))
-		 r.push_back(i);
-	   }
-	   return r;
-	}
-
-	bool is_measurement(qx::gate * g, uint32_t q)
-	{
-	   // println("is_measurement on " << q);
-	   if (g->type() == __measure_reg_gate__)
-	   {
-	      //println(" => register measurement !");
-	      return true;
-	   }
-	   if (g->type() == __measure_gate__)
-	   {
-	      // println(" => qubit measurement on " << g->qubits()[0]);
-	      if (g->qubits()[0] == q)
-		 return true;
-	   }
-	   if (g->type() == __parallel_gate__)
-	   {
-	      // println(" => parallel gates ...");
-	      std::vector<qx::gate *> gates = ((qx::parallel_gates*)g)->get_gates();
-	      for (uint32_t i=0; i<gates.size(); ++i)
-	      {
-		 if (is_measurement(gates[i],q))
-		    return true;
-	      }
-	   }
-	   return false;
-	}
-
-	/**
-	 * \brief enable error recording
-	 */
-	void enable_error_recording()
-	{
-	   error_recording = true;
-	}
-
-
-	/**
-	 * \brief disable error recording
-	 */
-	void disable_error_recording()
-	{
-	   error_recording = false;
-	}
-
-	const char error_type[3] = { 'x', 'z', 'y' };
-
-	/**
-	 * \brief report injected errors 
-	 */
-	void report_errors()
-	{
-	   if (errors.empty())
-	      return;
-	   println("[x] injected errors report : ");
-	   for (size_t i=0; i<errors.size(); ++i)
-	   {
-	      println("     " << error_location[i] << " : "<< error_type[errors[i].first] << " on qubit " << errors[i].second);
-	   }
-	   println("[+] stats : " << x_errors << " x | " << z_errors << " z | " << y_errors << " y ");
-	}
-
-	#define __verbose__ if (verbose)
-
-
-	/**
-	 * \brief inject errors
-	 * \return noisy circuit 
-	 */
-	qx::circuit *  inject(bool verbose=false)
-	{
-	   x_errors = 0;
-	   z_errors = 0;
-	   y_errors = 0;
-	  
-
-	   __verbose__ println("    [e] depolarizing_channel : injecting errors in circuit '" << c->id() << "'...");
-	   uint64_t steps = c->size();
-	   qx::circuit * noisy_c = new qx::circuit(nq,c->id() + "(noisy)");
-	   
-	   errors.clear();
-	   error_location.clear();
-	   total_errors = 0;
-
-	   __verbose__ println("    [+] circuit steps : " << steps);
-	   for (uint64_t p=0; p<steps; ++p)
-	   {
-	      qx::gate_type_t gt = c->get(p)->type();
-	      if ((gt==qx::__display__) || (gt==qx::__display_binary__))
-	      {
-		 noisy_c->add(c->get(p));
-		 continue;
-	      }
-	      // std::vector<uint32_t> used = c->get(p)->qubits();
-	      // std::vector<uint32_t> idle = idle_qubits(nq,used);
-	      // uint32_t idle_nq=idle.size();
-	      
-	      double x = uniform_rand();
-	      
-	      if (x<overall_error_probability)
-	      {
-		 uint32_t affected_qubits = 1;
-		 
-		 for (uint32_t i=2; i<(nq+1); ++i)
-		 {
-		    if (x>error_probability(nq,i,pe))
-		       break;
-		    affected_qubits++;
-		 }
-		 // affected_qubits = (affected_qubits > idle_nq ? idle_nq : affected_qubits);
-		 total_errors += affected_qubits;
-		 error_histogram[affected_qubits]++;
-
-		 __verbose__ println("    [>>>] error injection step " << p << " : number of affected qubits: " << affected_qubits);
-		 if (error_recording)
-		 {
-		    for (size_t e=0; e<affected_qubits; ++e)
-		       error_location.push_back(p);
-		 }
-
-		 if (affected_qubits==1)
-		 {
-		    // uint64_t q = idle[(rand()%idle_nq)];
-		    uint64_t q = rand()%nq;
-		    
-		    if (is_measurement(c->get(p),q))
-		    {
-		       __verbose__  print("      |--> error on qubit " << q << " (potential measurement error) ");
-		       noisy_c->add(measurement_error(q,verbose));
-		       noisy_c->add(c->get(p));
-		    }
-		    else
-		    {
-		       __verbose__  print("      |--> error on qubit  " << q);
-		       noisy_c->add(c->get(p));
-		       noisy_c->add(single_qubit_error(q,verbose));
-		    }
-		    
-
-		 }
-		 else
-		 {
-		    qx::parallel_gates * g = new qx::parallel_gates();
-		    bool measurement = false;
-
-		    std::vector<bool> v(nq,0);
-		    for (uint64_t i=0; i<affected_qubits; ++i)
-		    { 
-		       // uint64_t q = idle[(rand()%idle_nq)];
-		       uint64_t q = rand()%nq;
-		       while (v[q])
-			  q = rand()%nq;
-			  // q = idle[(rand()%idle_nq)];
-		       v[q] = 1;
-
-		       if (is_measurement(c->get(p),q))
-			  measurement = true;
-		       __verbose__  print("      |--> error on qubit  " << q << (measurement ? "(potential measurement error)" : ""));
-		       g->add(single_qubit_error(q,verbose));
-		    }
-		    
-		    if (measurement)
-		    {
-		       //__verbose__  println("      |--> (!) potential measurement error.");
-		       noisy_c->add(g);
-		       noisy_c->add(c->get(p));
-		    }
-		    else
-		    {
-		       noisy_c->add(c->get(p));
-		       noisy_c->add(g);
-		    }
-		 }
-	      }
-	      else
-	      {
-		 noisy_c->add(c->get(p));
-	      }
-	   }
-	   __verbose__ println("    [+] total injected errors in circuit '" <<  c->id() << "': " << total_errors);
-
-	   return noisy_c;
-	}
-
-
-	/**
-	 * \brief total errors
-	 * \return total errors
-	 */
-	 uint32_t get_total_errors()
-	 {
-	    return total_errors;
-	 }
-
-	 /**
-	  * overall probability of error
-	  */
-	  double get_overall_error_probability()
-	  {
-	     return overall_error_probability;
-	  }
-
-
-	/**
-	 * dump
-	 */
-	void dump()
-	{
-	   println("   [+] depolarizing channel :");
-	   println("   [-] single qubit error probability : " << pe);
-	   for (uint64_t i=1; i<(nq+1); ++i)
-	   {
-	      double pp = error_probability(nq,i,pe);
-	      println("   [i] simultaneous error(s) probability of " << i << " qubits out of " << nq << " : " << pp);
-	   }
-	   println("   [-] overall probability of errors: " << overall_error_probability);
-	   println("   [-] probability of (x) errors: " << xp);
-	   println("   [-] probability of (z) errors: " << yp);
-	   println("   [-] probability of (y) errors: " << zp);
-	}
+        /**
+         * ctor
+         */
+        depolarizing_channel(qx::circuit * c, size_t nq, double pe, double xp, double yp, double zp) : nrg(0,__third__), c(c), nq(nq), 
+	                                                                                               pe(pe), 
+                                                                                                  xp(xp), 
+                                                                                                  yp(yp), 
+                                                                                                  zp(zp),
+                                                                                                  overall_error_probability(0),
+                                                                                                  total_errors(0),
+                                                                                                  error_histogram(nq+1,0),
+                                                                                                  error_recording(false)
+        {
+           for (size_t i=1; i<(nq+1); ++i)
+           {
+              double pp = error_probability(nq,i,pe);
+              overall_error_probability += pp;
+           }
+
+           x_errors = 0;
+           z_errors = 0;
+           y_errors = 0;
+
+           srand48(xpu::timer().current());
+           qx::circuit * noisy_c = new qx::circuit(nq,c->id() + "(noisy)");
+        }
+
+        /**
+         * check whether a qubit is idle
+         */
+        bool is_idle(size_t q, std::vector<size_t> used)
+        {
+           for (int i=0; i<used.size(); ++i)
+              if (q == used[i])
+                 return false;
+           return true;
+        }
+
+        /**
+         * return idle qubits
+         */
+        std::vector<size_t> idle_qubits(size_t n, std::vector<size_t> & used)
+        {
+           std::vector<size_t> r;
+           for (size_t i=0; i<n; i++)
+           {
+              if (is_idle(i,used))
+                 r.push_back(i);
+           }
+           return r;
+        }
+
+        bool is_measurement(qx::gate * g, size_t q)
+        {
+           if (g->type() == __measure_reg_gate__)
+           {
+              return true;
+           }
+           if (g->type() == __measure_gate__)
+           {
+              if (g->qubits()[0] == q)
+                 return true;
+           }
+           if (g->type() == __parallel_gate__)
+           {
+              std::vector<qx::gate *> gates = ((qx::parallel_gates*)g)->get_gates();
+              for (size_t i=0; i<gates.size(); ++i)
+              {
+                 if (is_measurement(gates[i],q))
+                    return true;
+              }
+           }
+           return false;
+        }
+
+        /**
+         * \brief enable error recording
+         */
+        void enable_error_recording()
+        {
+           error_recording = true;
+        }
+
+
+        /**
+         * \brief disable error recording
+         */
+        void disable_error_recording()
+        {
+           error_recording = false;
+        }
+
+        const char error_type[3] = { 'x', 'z', 'y' };
+
+        /**
+         * \brief report injected errors 
+         */
+        void report_errors()
+        {
+           if (errors.empty())
+              return;
+           println("[x] injected errors report : ");
+           for (size_t i=0; i<errors.size(); ++i)
+           {
+              println("     " << error_location[i] << " : "<< error_type[errors[i].first] << " on qubit " << errors[i].second);
+           }
+           println("[+] stats : " << x_errors << " x | " << z_errors << " z | " << y_errors << " y ");
+        }
+
+#define __verbose__ if (verbose)
+
+
+        /**
+         * \brief inject errors
+         * \return noisy circuit 
+         */
+        qx::circuit *  inject(bool verbose=false)
+        {
+           x_errors = 0;
+           z_errors = 0;
+           y_errors = 0;
+
+
+           __verbose__ println("    [e] depolarizing_channel : injecting errors in circuit '" << c->id() << "'...");
+           size_t steps = c->size();
+           qx::circuit * noisy_c = new qx::circuit(nq,c->id() + "(noisy)");
+
+           errors.clear();
+           error_location.clear();
+           total_errors = 0;
+
+           __verbose__ println("    [+] circuit steps : " << steps);
+           for (size_t p=0; p<steps; ++p)
+           {
+              qx::gate_type_t gt = c->get(p)->type();
+              if ((gt==qx::__display__) || (gt==qx::__display_binary__))
+              {
+                 noisy_c->add(c->get(p));
+                 continue;
+              }
+              // std::vector<size_t> used = c->get(p)->qubits();
+              // std::vector<size_t> idle = idle_qubits(nq,used);
+              // size_t idle_nq=idle.size();
+
+              double x = uniform_rand();
+
+              if (x<overall_error_probability)
+              {
+                 size_t affected_qubits = 1;
+
+                 for (size_t i=2; i<(nq+1); ++i)
+                 {
+                    if (x>error_probability(nq,i,pe))
+                       break;
+                    affected_qubits++;
+                 }
+                 // affected_qubits = (affected_qubits > idle_nq ? idle_nq : affected_qubits);
+                 total_errors += affected_qubits;
+                 error_histogram[affected_qubits]++;
+
+                 __verbose__ println("    [>>>] error injection step " << p << " : number of affected qubits: " << affected_qubits);
+                 if (error_recording)
+                 {
+                    for (size_t e=0; e<affected_qubits; ++e)
+                       error_location.push_back(p);
+                 }
+
+                 if (affected_qubits==1)
+                 {
+                    // size_t q = idle[(rand()%idle_nq)];
+                    size_t q = rand()%nq;
+
+                    if (is_measurement(c->get(p),q))
+                    {
+                       __verbose__  print("      |--> error on qubit " << q << " (potential measurement error) ");
+                       noisy_c->add(measurement_error(q,verbose));
+                       noisy_c->add(c->get(p));
+                    }
+                    else
+                    {
+                       __verbose__  print("      |--> error on qubit  " << q);
+                       noisy_c->add(c->get(p));
+                       noisy_c->add(single_qubit_error(q,verbose));
+                    }
+
+
+                 }
+                 else
+                 {
+                    qx::parallel_gates * g = new qx::parallel_gates();
+                    bool measurement = false;
+
+                    std::vector<bool> v(nq,0);
+                    for (size_t i=0; i<affected_qubits; ++i)
+                    { 
+                       // size_t q = idle[(rand()%idle_nq)];
+                       size_t q = rand()%nq;
+                       while (v[q])
+                          q = rand()%nq;
+                       // q = idle[(rand()%idle_nq)];
+                       v[q] = 1;
+
+                       if (is_measurement(c->get(p),q))
+                          measurement = true;
+                       __verbose__  print("      |--> error on qubit  " << q << (measurement ? "(potential measurement error)" : ""));
+                       g->add(single_qubit_error(q,verbose));
+                    }
+
+                    if (measurement)
+                    {
+                       //__verbose__  println("      |--> (!) potential measurement error.");
+                       noisy_c->add(g);
+                       noisy_c->add(c->get(p));
+                    }
+                    else
+                    {
+                       noisy_c->add(c->get(p));
+                       noisy_c->add(g);
+                    }
+                 }
+              }
+              else
+              {
+                 noisy_c->add(c->get(p));
+              }
+           }
+           __verbose__ println("    [+] total injected errors in circuit '" <<  c->id() << "': " << total_errors);
+
+           return noisy_c;
+        }
+
+
+        /**
+         * \brief total errors
+         * \return total errors
+         */
+        size_t get_total_errors()
+        {
+           return total_errors;
+        }
+
+        /**
+         * overall probability of error
+         */
+        double get_overall_error_probability()
+        {
+           return overall_error_probability;
+        }
+
+
+        /**
+         * dump
+         */
+        void dump()
+        {
+           println("   [+] depolarizing channel :");
+           println("   [-] single qubit error probability : " << pe);
+           for (size_t i=1; i<(nq+1); ++i)
+           {
+              double pp = error_probability(nq,i,pe);
+              println("   [i] simultaneous error(s) probability of " << i << " qubits out of " << nq << " : " << pp);
+           }
+           println("   [-] overall probability of errors: " << overall_error_probability);
+           println("   [-] probability of (x) errors: " << xp);
+           println("   [-] probability of (z) errors: " << yp);
+           println("   [-] probability of (y) errors: " << zp);
+        }
 
       private:
 
 
-	/**
-	 * uniform random number generator
-	 */
-	double uniform_rand()
-	{
-	   return urg.next();
-	}
-
-	/**
-	 * normal random number generator
-	 */
-	double normal_rand()
-	{
-	   double r = 1-std::abs(nrg.next()); 
-	   r = (r < 0 ? 0 : (r > __limit__ ? __limit__ : r));
-	   return r; 
-	}
-
-       /**
-        * single qubit error 
-	*/
-	qx::gate * single_qubit_error(uint32_t q, bool verbose=false)
-	{
-	   double p = uniform_rand();
-	   if (p<xp)
-	   {
-	      __verbose__ println(" (x error) ");
-	      if (error_recording)
-		 errors.push_back(error_t(__x_error__,q));
-	      x_errors++;
-	      return new qx::pauli_x(q);
-	   }
-	   else if (p<(zp+xp))
-	   {
-	      __verbose__ println(" (z error) ");
-	      if (error_recording)
-		 errors.push_back(error_t(__z_error__,q));
-	      z_errors++;
-	      return new qx::pauli_z(q);
-	   }
-	   else
-	   {
-	      __verbose__ println(" (y error) ");
-	      if (error_recording)
-		 errors.push_back(error_t(__y_error__,q));
-	      y_errors++;
-	      return new qx::pauli_y(q);
-	   }
-	}
-       
-       /**
-        * measurement error 
-	*/
-	qx::gate * measurement_error(uint32_t q, bool verbose=false)
-	{
-	   __verbose__ println(" (measurement error) ");
-	   if (error_recording)
-	      errors.push_back(error_t(__x_error__,q));
-	   return new qx::pauli_x(q);
-	}
-
-
-	/**
-	 * factorial
-	 */
-	uint64_t factorial(uint64_t n)
-	{
-	   return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
-	}
-
-	/**
-	 * combinations
-	 */
-	uint64_t combinations(uint64_t n, uint64_t k)
-	{
-	   return (factorial(n)/(factorial(k)*factorial(n-k)));
-	}
-
-	/**
-	 * nq : total number of qubits
-	 * ne : number of simultaneously affected qubits
-	 * p  : probability of single error
-	 */
-	double error_probability(uint64_t nq, uint64_t ne, double p)
-	{
-	   double c = (double)combinations(nq,ne);
-	   double res = c*std::pow(p,ne)*std::pow(1-p,nq-ne);
-	   return res;
-	}
+        /**
+         * uniform random number generator
+         */
+        double uniform_rand()
+        {
+           return urg.next();
+        }
+
+        /**
+         * normal random number generator
+         */
+        double normal_rand()
+        {
+           double r = 1-std::abs(nrg.next()); 
+           r = (r < 0 ? 0 : (r > __limit__ ? __limit__ : r));
+           return r; 
+        }
+
+        /**
+         * single qubit error 
+         */
+        qx::gate * single_qubit_error(size_t q, bool verbose=false)
+        {
+           double p = uniform_rand();
+           if (p<xp)
+           {
+              __verbose__ println(" (x error) ");
+              if (error_recording)
+                 errors.push_back(error_t(__x_error__,q));
+              x_errors++;
+              return new qx::pauli_x(q);
+           }
+           else if (p<(zp+xp))
+           {
+              __verbose__ println(" (z error) ");
+              if (error_recording)
+                 errors.push_back(error_t(__z_error__,q));
+              z_errors++;
+              return new qx::pauli_z(q);
+           }
+           else
+           {
+              __verbose__ println(" (y error) ");
+              if (error_recording)
+                 errors.push_back(error_t(__y_error__,q));
+              y_errors++;
+              return new qx::pauli_y(q);
+           }
+        }
+
+        /**
+         * measurement error 
+         */
+        qx::gate * measurement_error(size_t q, bool verbose=false)
+        {
+           __verbose__ println(" (measurement error) ");
+           if (error_recording)
+              errors.push_back(error_t(__x_error__,q));
+           return new qx::pauli_x(q);
+        }
+
+
+        /**
+         * factorial
+         */
+        size_t factorial(size_t n)
+        {
+           return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
+        }
+
+        /**
+         * combinations
+         */
+        size_t combinations(size_t n, size_t k)
+        {
+           return (factorial(n)/(factorial(k)*factorial(n-k)));
+        }
+
+        /**
+         * nq : total number of qubits
+         * ne : number of simultaneously affected qubits
+         * p  : probability of single error
+         */
+        double error_probability(size_t nq, size_t ne, double p)
+        {
+           double c = (double)combinations(nq,ne);
+           double res = c*std::pow(p,ne)*std::pow(1-p,nq-ne);
+           return res;
+        }
 
       private:
 
-	/**
-	 * parameters
-	 */
-	 
-	 qx::normal_random_number_generator  nrg;
-	 qx::uniform_random_number_generator urg;
-	 
-	 qx::circuit *        c;
-	 qx::circuit *        noisy_c;
-	 uint64_t             nq;
-	 
-	 double               pe;
-	 double               xp;
-	 double               yp;
-	 double               zp;
-	 
-	 double               overall_error_probability;
-	 uint64_t             total_errors;
-	 histogram_t          error_histogram;
-
-	 bool                 error_recording;
-	 std::vector<error_t> errors;
-	 std::vector<size_t>  error_location;
-	 size_t               x_errors;
-	 size_t               z_errors;
-	 size_t               y_errors;
+        /**
+         * parameters
+         */
+
+        qx::normal_random_number_generator  nrg;
+        qx::uniform_random_number_generator urg;
+
+        qx::circuit *        c;
+        qx::circuit *        noisy_c;
+        size_t               nq;
+
+        double               pe;
+        double               xp;
+        double               yp;
+        double               zp;
+
+        double               overall_error_probability;
+        size_t               total_errors;
+        histogram_t          error_histogram;
+
+        bool                 error_recording;
+        std::vector<error_t> errors;
+        std::vector<size_t>  error_location;
+        size_t               x_errors;
+        size_t               z_errors;
+        size_t               y_errors;
 
    };
 }
diff --git a/src/core/error_model.h b/src/core/error_model.h
index f8ccfcc3..3f113ce7 100644
--- a/src/core/error_model.h
+++ b/src/core/error_model.h
@@ -16,6 +16,19 @@ namespace qx
       __unknown_error_model__
    } error_model_t;
 
+
+   qx::circuit * noisy_dep_ch(qx::circuit * c, double p, size_t& total_errors)
+   {
+      if (c)
+      {
+         qx::depolarizing_channel dep_ch(c, c->get_qubit_count(), p);
+         qx::circuit * noisy_c = dep_ch.inject(false); //true);
+         total_errors += dep_ch.get_total_errors();
+         return noisy_c;
+      }
+      return NULL;
+   }
+
 };
 
 
diff --git a/src/core/register.cc b/src/core/register.cc
index ecb849a1..383a2fd3 100644
--- a/src/core/register.cc
+++ b/src/core/register.cc
@@ -232,47 +232,54 @@ int64_t qx::qu_register::measure()
  */
 void qx::qu_register::dump(bool only_binary=false)
 {
-   if (!only_binary)
-   {
-      println("--------------[quantum state]-------------- ");
-      std::cout << std::fixed;
-      for (int i=0; i<data.size(); ++i)
-      {
-         if ((std::abs(data[i].re) > __amp_epsilon__) ||
-             (std::abs(data[i].im) > __amp_epsilon__))
-	 {
-	    print("   " << std::showpos << std::setw(14) << data[i] << " |"); to_binary(i,n_qubits); println("> +");
-	 }
-      }
-   }
-   if (measurement_averaging_enabled)
-   {
-      println("------------------------------------------- ");
-      print("[>>] measurement averaging (ground state) :");
-      print(" ");
-      for (int i=measurement_averaging.size()-1; i>=0; --i)
-      {
-	 double gs = measurement_averaging[i].ground_states;
-	 double es = measurement_averaging[i].exited_states;
-	 // println("(" << gs << "," << es << ")");
-	 double av = ((es+gs) != 0. ? (gs/(es+gs)) : 0.);
-	 print(" | " << av);  
-      }
-      println(" |");
-   }
-   println("------------------------------------------- ");
-   print("[>>] measurement prediction:");
-   print(" ");
-   for (int i=measurement_prediction.size()-1; i>=0; --i)
-      print(" | " << __format_bin(measurement_prediction[i]));  
-   println(" |");
-   println("------------------------------------------- ");
-   print("[>>] measurement register  :");
-   print(" ");
-   for (int i=measurement_register.size()-1; i>=0; --i)
-      print(" | " << (measurement_register[i] ? '1' : '0'));  
-   println(" |");
-   println("------------------------------------------- ");
+    if (!only_binary)
+    {
+        println("--------------[quantum state]-------------- ");
+        std::streamsize stream_size = std::cout.precision();
+        // std::cout.precision(std::numeric_limits<double>::digits10);
+        std::cout.precision(7); //std::numeric_limits<double>::digits10);
+        std::cout << std::fixed;
+        for (int i=0; i<data.size(); ++i)
+        {
+            if ((std::abs(data[i].re) > __amp_epsilon__) ||
+                (std::abs(data[i].im) > __amp_epsilon__))
+            {
+                print("  [p = " << std::showpos << data[i].norm() << "]");
+                print("  " << std::showpos << data[i] << " |");
+                to_binary(i,n_qubits);
+                println("> +");
+            }
+        }
+        std::cout.precision(stream_size);
+    }
+    if (measurement_averaging_enabled)
+    {
+        std::setprecision(9);
+        println("------------------------------------------- ");
+        print("[>>] measurement averaging (ground state) :");
+        print(" ");
+        for (int i=measurement_averaging.size()-1; i>=0; --i)
+        {
+            double gs = measurement_averaging[i].ground_states;
+            double es = measurement_averaging[i].exited_states;
+            double av = ((es+gs) != 0. ? (gs/(es+gs)) : 0.);
+            print(" | " << std::setw(9) << av);  
+        }
+        println(" |");
+    }
+    println("------------------------------------------- ");
+    print("[>>] measurement prediction               :");
+    print(" ");
+    for (int i=measurement_prediction.size()-1; i>=0; --i)
+        print(" | " <<  std::setw(9) << __format_bin(measurement_prediction[i]));  
+    println(" |");
+    println("------------------------------------------- ");
+    print("[>>] measurement register                 :");
+    print(" ");
+    for (int i=measurement_register.size()-1; i>=0; --i)
+        print(" | " <<  std::setw(9) << (measurement_register[i] ? '1' : '0'));  
+    println(" |");
+    println("------------------------------------------- ");
 }
 
 /**
@@ -409,4 +416,3 @@ double fidelity(qu_register& s1, qu_register& s2)
    
    return f;
 }
-
diff --git a/src/libqasm_interface.h b/src/libqasm_interface.h
index 0a5a6113..5f27bee5 100644
--- a/src/libqasm_interface.h
+++ b/src/libqasm_interface.h
@@ -4,7 +4,7 @@
 #include <core/circuit.h>
 #include <qasm_ast.hpp>
 
-#define ITER_FOR_IN(e, l) for (auto e = l.begin(); e != l.end(); e++)
+#define __for_in(e, l) for (auto e = l.begin(); e != l.end(); e++)
 
 #define __ret_gate_1(__g) \
 {\
@@ -148,8 +148,10 @@ qx::gate *gateLookup(compiler::Operation &operation)
    }
 
    ///////// common sq gates //////
+   if (type == "i")
+      __ret_gate_1(qx::identity)
    if (type == "x")
-      __ret_gate_1(qx::pauli_x)
+	  __ret_gate_1(qx::pauli_x)
    if (type == "y")
       __ret_gate_1(qx::pauli_y)
    if (type == "z")
@@ -435,7 +437,7 @@ qx::gate *gateLookup(compiler::Operation &operation)
    return NULL;
 }
 
-qx::circuit *qxCircuitFromCQasm(uint64_t qubits_count, compiler::SubCircuit &subcircuit)
+qx::circuit * load_cqasm_code(uint64_t qubits_count, compiler::SubCircuit &subcircuit)
 {
   uint64_t iterations = subcircuit.numberIterations();
   std::string name = subcircuit.nameSubCircuit();
@@ -446,11 +448,11 @@ qx::circuit *qxCircuitFromCQasm(uint64_t qubits_count, compiler::SubCircuit &sub
   const std::vector<compiler::OperationsCluster*>& clusters
     = subcircuit.getOperationsCluster();
 
-  ITER_FOR_IN(p_cluster, clusters)
+  __for_in(p_cluster, clusters)
   {
     const std::vector<compiler::Operation*> operations
       = (*p_cluster)->getOperations();
-    ITER_FOR_IN(p_operation, operations)
+    __for_in(p_operation, operations)
     {
        qx::gate * g;
        try
diff --git a/src/qx_representation.h b/src/qx_representation.h
index e82675fa..e54676cb 100644
--- a/src/qx_representation.h
+++ b/src/qx_representation.h
@@ -8,56 +8,56 @@
 namespace qx
 {
 
-class QxRepresentation
-{
-  public:
-    QxRepresentation(uint32_t count) : qubits_count(count) {}
-    ~QxRepresentation()
-    {
-      for (auto circuit = circuits_.begin(); circuit != circuits_.end(); ++circuit)
-        delete(*circuit);
-    }
-
-
-    uint32_t getQubitsCount()
-    {
-      return qubits_count;
-    }
-
-
-    error_model_t getErrorModel(void)
-    {
-      return error_model;
-    }
-
-    void setErrorModel(error_model_t model)
-    {
-      error_model = model;
-    }
-
-
-    void setErrorProbability(double probability)
-    {
-      error_probability = probability;
-    }
-
-    double getErrorProbability(void)
-    {
-      return error_probability;
-    }
-
-
-    std::vector<qx::circuit*> &circuits(void)
-    {
-      return circuits_;
-    }
-
-  private:
-    uint32_t      qubits_count;
-    error_model_t error_model =  __unknown_error_model__;
-    double        error_probability;
-    std::vector<qx::circuit*> circuits_;
-};
+   class QxRepresentation
+   {
+      public:
+         QxRepresentation(size_t count) : qubits_count(count) {}
+         ~QxRepresentation()
+         {
+            for (auto circuit = circuits_.begin(); circuit != circuits_.end(); ++circuit)
+               delete(*circuit);
+         }
+
+
+         size_t getQubitsCount()
+         {
+            return qubits_count;
+         }
+
+
+         error_model_t getErrorModel(void)
+         {
+            return error_model;
+         }
+
+         void setErrorModel(error_model_t model)
+         {
+            error_model = model;
+         }
+
+
+         void setErrorProbability(double probability)
+         {
+            error_probability = probability;
+         }
+
+         double getErrorProbability(void)
+         {
+            return error_probability;
+         }
+
+
+         std::vector<qx::circuit*> &circuits(void)
+         {
+            return circuits_;
+         }
+
+      private:
+         size_t      qubits_count;
+         error_model_t error_model =  __unknown_error_model__;
+         double        error_probability;
+         std::vector<qx::circuit*> circuits_;
+   };
 
 } // namepsace qx
 
diff --git a/src/simulator.cc b/src/simulator.cc
index c04e043b..550eae4a 100644
--- a/src/simulator.cc
+++ b/src/simulator.cc
@@ -36,126 +36,178 @@ void print_banner() {
  */
 int main(int argc, char **argv)
 {
-  std::string file_path;
-  size_t ncpu = 0;
-  print_banner();
+   std::string file_path;
+   size_t ncpu = 0;
+   size_t navg = 0;
+   print_banner();
 
-  if (!(argc == 2 || argc == 3))
-  {
+   if (!(argc == 2 || argc == 3 || argc == 4))
+   {
       println("error : you must specify a circuit file !");
-      println("usage: \n   " << argv[0] << " file.qc");
+      println("usage: \n   " << argv[0] << " file.qc [iterations] [num_cpu]");
       return -1;
-  }
-    
-  // parse arguments and initialise xpu cores
-  file_path = argv[1];
-  if (argc == 3) ncpu = (atoi(argv[2]));
-  if (ncpu && ncpu < 128) xpu::init(ncpu);
-  else xpu::init();
-
-  // parse file and create abstract syntax tree
-  println("[+] loading circuit from '" << file_path << "' ...");
-  FILE *qasm_file = fopen(file_path.c_str(), "r");
-  if (!qasm_file)
-  {
-    std::cerr << "[x] error: could not open " << file_path << std::endl;
-    xpu::clean();
-    return -1;
-  }
-
-  // construct libqasm parser and safely parse input file
-  compiler::QasmSemanticChecker *parser;
-  compiler::QasmRepresentation ast;
-  try
-  {
-    parser = new compiler::QasmSemanticChecker(qasm_file);
-    ast = parser->getQasmRepresentation();
-  }
-  catch (std::exception &e)
-  {
-    std::cerr << "error while parsing file " << file_path << ": " << std::endl;
-    std::cerr << e.what() << std::endl;
-    xpu::clean();
-    return -1;
-  }
-
-  
-  // Create qx circuits and register
-  std::vector<qx::circuit*> circuits;
-  uint32_t qubits = ast.numQubits();
-  println("[+] creating quantum register of " << qubits << " qubits... ");
-  qx::qu_register *reg = NULL;
-
-  try {
-	   reg = new qx::qu_register(qubits);
-  } catch(std::bad_alloc& exception) {
-	   std::cerr << "[x] not enough memory, aborting" << std::endl;
-	   xpu::clean();
-	   return -1;
-  } catch(std::exception& exception) {
-	   std::cerr << "[x] unexpected exception (" << exception.what() << "), aborting" << std::endl;
-	   xpu::clean();
-	   return -1;
-  }
-
-  // Convert libqasm ast to qx internal representation
-  qx::QxRepresentation qxr = qx::QxRepresentation(qubits);
-  std::vector<compiler::SubCircuit> subcircuits
-    = ast.getSubCircuits().getAllSubCircuits();
-  ITER_FOR_IN(subcircuit, subcircuits)
-  {
-    try
-    {
-      circuits.push_back(qxCircuitFromCQasm(qubits, *subcircuit));
-    }
-    catch (std::string type)
-    {
-      std::cerr << "[x] Encountered unsuported gate: " << type << std::endl;
+   }
+
+   // parse arguments and initialise xpu cores
+   file_path = argv[1];
+   if (argc > 2) navg = (atoi(argv[2]));
+   if (argc > 3) ncpu = (atoi(argv[3]));
+   if (ncpu && ncpu < 128) xpu::init(ncpu);
+   else xpu::init();
+
+   // parse file and create abstract syntax tree
+   println("[+] loading circuit from '" << file_path << "' ...");
+   FILE * qasm_file = fopen(file_path.c_str(), "r");
+   if (!qasm_file)
+   {
+      std::cerr << "[x] error: could not open " << file_path << std::endl;
       xpu::clean();
       return -1;
-    }
-  }
-
-  // Execute qx code
-  uint32_t total_errors = 0;
-
-  // check whether an error model is specified
-  if (qxr.getErrorModel() == qx::__depolarizing_channel__)
-  {
-     println("[+] generating noisy circuits...");
-     double error_probability = qxr.getErrorProbability(); 
-     std::vector<qx::circuit*> perfect_circuits = qxr.circuits();
-     for (uint32_t i=0; i<perfect_circuits.size(); i++)
-     {
-        println("[>] processing circuit '" << perfect_circuits[i]->id() << "'...");
-  uint32_t iterations = perfect_circuits[i]->get_iterations();
-  if (iterations > 1)
-  {
-     for (uint32_t it=0; it<iterations; ++it)
-     {
-        qx::depolarizing_channel dep_ch(perfect_circuits[i], qubits, error_probability);
-        qx::circuit * noisy_c = dep_ch.inject(true);
-        println("[+] noisy circuit : " << std::hex << noisy_c);
-        circuits.push_back(noisy_c);
-        total_errors += dep_ch.get_total_errors();
-     }
-  }
-  else
-  {
-     qx::depolarizing_channel dep_ch(perfect_circuits[i], qubits, error_probability);
-     qx::circuit * noisy_c = dep_ch.inject(true);
-     circuits.push_back(noisy_c);
-     total_errors += dep_ch.get_total_errors();
-  }
-     }
-     println("[+] total errors injected in all circuits : " << total_errors);
-
-  }
- 
-  for (uint32_t i=0; i<circuits.size(); i++)
-     circuits[i]->execute(*reg);
- 
-  xpu::clean();
-
-  return 0;
+   }
+
+   // construct libqasm parser and safely parse input file
+   compiler::QasmSemanticChecker * parser;
+   compiler::QasmRepresentation ast;
+   try
+   {
+      parser = new compiler::QasmSemanticChecker(qasm_file);
+      ast = parser->getQasmRepresentation();
+   }
+   catch (std::exception &e)
+   {
+      std::cerr << "error while parsing file " << file_path << ": " << std::endl;
+      std::cerr << e.what() << std::endl;
+      xpu::clean();
+      return -1;
+   }
+
+   // quantum state and circuits
+   size_t                     qubits = ast.numQubits();
+   qx::qu_register *          reg = NULL;
+   std::vector<qx::circuit*>  circuits;
+   std::vector<qx::circuit*>  noisy_circuits;
+   std::vector<qx::circuit *> perfect_circuits; 
+   // error model parameters
+   size_t                     total_errors      = 0;
+   double                     error_probability = 0;
+   qx::error_model_t          error_model       = qx::__unknown_error_model__;
+
+   // create the quantum state
+   println("[+] creating quantum register of " << qubits << " qubits... ");
+   try {
+      reg = new qx::qu_register(qubits);
+   } catch(std::bad_alloc& exception) {
+      std::cerr << "[x] not enough memory, aborting" << std::endl;
+      xpu::clean();
+      return -1;
+   } catch(std::exception& exception) {
+      std::cerr << "[x] unexpected exception (" << exception.what() << "), aborting" << std::endl;
+      xpu::clean();
+      return -1;
+   }
+
+   // convert libqasm ast to qx internal representation
+   // qx::QxRepresentation qxr = qx::QxRepresentation(qubits);
+   std::vector<compiler::SubCircuit> subcircuits = ast.getSubCircuits().getAllSubCircuits();
+   __for_in(subcircuit, subcircuits)
+   {
+      try
+      {
+         // qxr.circuits().push_back(load_cqasm_code(qubits, *subcircuit));
+         perfect_circuits.push_back(load_cqasm_code(qubits, * subcircuit));
+      }
+      catch (std::string type)
+      {
+         std::cerr << "[x] encountered unsuported gate: " << type << std::endl;
+         xpu::clean();
+         return -1;
+      }
+   }
+
+   println("[i] loaded " << perfect_circuits.size() << " circuits.");
+
+   // check whether an error model is specified
+   if (ast.getErrorModelType() == "depolarizing_channel")
+   {
+      error_probability = ast.getErrorModelProbability();
+      error_model       = qx::__depolarizing_channel__;
+   }
+
+   // measurement averaging
+   if (navg)
+   {
+      if (error_model == qx::__depolarizing_channel__)
+      {
+         qx::measure m;
+         for (size_t s=0; s<navg; ++s)
+         {
+            reg->reset();
+            for (size_t i=0; i<perfect_circuits.size(); i++)
+            {
+               if (perfect_circuits[i]->size() == 0)
+                  continue;
+               size_t iterations = perfect_circuits[i]->get_iterations();
+               if (iterations > 1)
+               {
+                  for (size_t it=0; it<iterations; ++it)
+                  {
+                     qx::noisy_dep_ch(perfect_circuits[i],error_probability,total_errors)->execute(*reg,false,true);
+                  }
+               } else
+                     qx::noisy_dep_ch(perfect_circuits[i],error_probability,total_errors)->execute(*reg,false,true);
+            }
+            m.apply(*reg);
+         }
+      }
+      else
+      {
+         qx::measure m;
+         for (size_t s=0; s<navg; ++s)
+         {
+            reg->reset();
+            for (size_t i=0; i<perfect_circuits.size(); i++)
+               perfect_circuits[i]->execute(*reg,false,true);
+            m.apply(*reg);
+         }
+      }
+      
+      println("[+] average measurement after " << navg << " shots:");
+      reg->dump(true);
+   }
+   else
+   {
+      // if (qxr.getErrorModel() == qx::__depolarizing_channel__)
+      if (error_model == qx::__depolarizing_channel__)
+      {
+         // println("[+] generating noisy circuits (p=" << qxr.getErrorProbability() << ")...");
+         for (size_t i=0; i<perfect_circuits.size(); i++)
+         {
+            if (perfect_circuits[i]->size() == 0)
+               continue;
+            // println("[>] processing circuit '" << perfect_circuits[i]->id() << "'...");
+            size_t iterations = perfect_circuits[i]->get_iterations();
+            if (iterations > 1)
+            {
+               for (size_t it=0; it<iterations; ++it)
+                  circuits.push_back(qx::noisy_dep_ch(perfect_circuits[i],error_probability,total_errors));
+            }
+            else
+            {
+               circuits.push_back(qx::noisy_dep_ch(perfect_circuits[i],error_probability,total_errors));
+            }
+         }
+         // println("[+] total errors injected in all circuits : " << total_errors);
+      }
+      else 
+         circuits = perfect_circuits; // qxr.circuits();
+
+
+      for (size_t i=0; i<circuits.size(); i++)
+         circuits[i]->execute(*reg);
+   }
+
+   // exit(0);
+   xpu::clean();
+
+   return 0;
 }