diff --git a/src/qsylvan_simulator.c b/src/qsylvan_simulator.c
index 3fc7e94..0b1bee0 100644
--- a/src/qsylvan_simulator.c
+++ b/src/qsylvan_simulator.c
@@ -359,17 +359,18 @@ check_ctrls_before_targ(BDDVAR *c1, BDDVAR *c2, BDDVAR *c3, BDDVAR t)
 TASK_IMPL_3(QMDD, qmdd_gate, QMDD, qmdd, gate_id_t, gate, BDDVAR, target)
 {
     qmdd_do_before_gate(&qmdd);
-    return qmdd_gate_rec(qmdd, gate, target);
+    aadd_refs_push(qmdd);
+    QMDD res = qmdd_gate_rec(qmdd, gate, target);
+    aadd_refs_pop(1);
+    return res;
 }
 
 /* Wrapper for applying controlled gates with 1, 2, or 3 control qubits. */
 QMDD _qmdd_cgate(QMDD state, gate_id_t gate, BDDVAR c1, BDDVAR c2, BDDVAR c3, BDDVAR t, BDDVAR n)
 {
-    qmdd_do_before_gate(&state);
-
     if (check_ctrls_before_targ(&c1, &c2, &c3, t)) {
         BDDVAR cs[4] = {c1, c2, c3, AADD_INVALID_VAR}; // last pos is to mark end
-        return qmdd_cgate_rec(state, gate, cs, t);
+        return RUN(qmdd_cgate, state, gate, cs, t);//qmdd_cgate_rec(state, gate, cs, t);
     }
     else {
         assert(n != 0 && "ERROR: when ctrls > targ, nqubits must be passed to cgate() function.");
@@ -377,6 +378,14 @@ QMDD _qmdd_cgate(QMDD state, gate_id_t gate, BDDVAR c1, BDDVAR c2, BDDVAR c3, BD
         return aadd_matvec_mult(gate_matrix, state, n);
     }
 }
+TASK_IMPL_4(QMDD, qmdd_cgate, QMDD, state, gate_id_t, gate, BDDVAR*, cs, BDDVAR, t)
+{
+    qmdd_do_before_gate(&state);
+    aadd_refs_push(state);
+    QMDD res = qmdd_cgate_rec(state, gate, cs, t);
+    aadd_refs_pop(1);
+    return res;
+}
 
 /* Wrapper for applying a controlled gate where the controls are a range. */
 TASK_IMPL_5(QMDD, qmdd_cgate_range, QMDD, qmdd, gate_id_t, gate, BDDVAR, c_first, BDDVAR, c_last, BDDVAR, t)
diff --git a/src/qsylvan_simulator.h b/src/qsylvan_simulator.h
index dc460dc..2e6ca63 100644
--- a/src/qsylvan_simulator.h
+++ b/src/qsylvan_simulator.h
@@ -272,6 +272,7 @@ TASK_DECL_3(QMDD, qmdd_gate, QMDD, gate_id_t, BDDVAR);
 #define qmdd_cgate2(state,gate,c1,c2,t,...) (_qmdd_cgate(state,gate,c1,c2,AADD_INVALID_VAR,t,(0, ##__VA_ARGS__)))
 #define qmdd_cgate3(state,gate,c1,c2,c3,t,...) (_qmdd_cgate(state,gate,c1,c2,c3,t,(0, ##__VA_ARGS__)))
 QMDD _qmdd_cgate(QMDD state, gate_id_t gate, BDDVAR c1, BDDVAR c2, BDDVAR c3, BDDVAR t, BDDVAR n);
+TASK_DECL_4(QMDD, qmdd_cgate, QMDD, gate_id_t, BDDVAR*, BDDVAR);
 
 /* Applies given controlled gate to |q>. */
 #define qmdd_cgate_range(qmdd,gate,c_first,c_last,t) (RUN(qmdd_cgate_range,qmdd,gate,c_first,c_last,t))
diff --git a/src/sylvan_aadd.c b/src/sylvan_aadd.c
index 1f23374..0be6349 100644
--- a/src/sylvan_aadd.c
+++ b/src/sylvan_aadd.c
@@ -649,14 +649,20 @@ TASK_IMPL_2(AADD, aadd_plus, AADD, a, AADD, b)
 TASK_IMPL_3(AADD, aadd_matvec_mult, AADD, mat, AADD, vec, BDDVAR, nvars)
 {
     aadd_do_before_mult(&mat, &vec);
-    return CALL(aadd_matvec_mult_rec, mat, vec, nvars, 0);
+    aadd_refs_push(mat); aadd_refs_push(vec);
+    AADD res = CALL(aadd_matvec_mult_rec, mat, vec, nvars, 0);
+    aadd_refs_pop(2);
+    return res;
 }
 
 /* Wrapper for matrix vector multiplication. */
 TASK_IMPL_3(AADD, aadd_matmat_mult, AADD, a, AADD, b, BDDVAR, nvars)
 {
     aadd_do_before_mult(&a, &b);
-    return CALL(aadd_matmat_mult_rec, a, b, nvars, 0);
+    aadd_refs_push(a); aadd_refs_push(b);
+    AADD res = CALL(aadd_matmat_mult_rec, a, b, nvars, 0);
+    aadd_refs_pop(2);
+    return res;
 }
 
 TASK_IMPL_4(AADD, aadd_matvec_mult_rec, AADD, mat, AADD, vec, BDDVAR, nvars, BDDVAR, nextvar)