diff --git a/src/tequila/apps/unary_state_prep.py b/src/tequila/apps/unary_state_prep.py
index b59a0f98..55ad73e9 100644
--- a/src/tequila/apps/unary_state_prep.py
+++ b/src/tequila/apps/unary_state_prep.py
@@ -89,7 +89,6 @@ def __init__(self, target_space: typing.List[BitString], max_repeat: int = 100,
         simulator.convert_to_numpy = False
         variables = None # {k:k.name.evalf() for k in self._abstract_circuit.extract_variables()}
         wfn = simulator.simulate(initial_state=BitString.from_int(0, nbits=self.n_qubits), variables=variables)
-        wfn.n_qubits = self._n_qubits
         equations = []
         for k in target_space:
             equations.append(wfn[k] - abstract_coefficients[k])
@@ -174,7 +173,7 @@ def __call__(self, wfn: QubitWaveFunction) -> QCircuit:
         :return:
         """
         try:
-            assert (len(wfn) == len(self._target_space))
+            assert wfn.length() == len(self._target_space)
             for key in wfn.keys():
                 try:
                     assert (key in self._target_space)
diff --git a/src/tequila/hamiltonian/paulis.py b/src/tequila/hamiltonian/paulis.py
index 175c5850..d282aa4c 100644
--- a/src/tequila/hamiltonian/paulis.py
+++ b/src/tequila/hamiltonian/paulis.py
@@ -256,7 +256,7 @@ def Projector(wfn, threshold=0.0, n_qubits=None) -> QubitHamiltonian:
 
     """
 
-    wfn = QubitWaveFunction(state=wfn, n_qubits=n_qubits)
+    wfn = QubitWaveFunction.convert_from(n_qubits, wfn)
 
     H = QubitHamiltonian.zero()
     for k1, v1 in wfn.items():
@@ -304,8 +304,9 @@ def KetBra(ket: QubitWaveFunction, bra: QubitWaveFunction, hermitian: bool = Fal
 
     """
     H = QubitHamiltonian.zero()
-    ket = QubitWaveFunction(state=ket, n_qubits=n_qubits)
-    bra = QubitWaveFunction(state=bra, n_qubits=n_qubits)
+
+    ket = QubitWaveFunction.convert_from(n_qubits, ket)
+    bra = QubitWaveFunction.convert_from(n_qubits, bra)
 
     for k1, v1 in bra.items():
         for k2, v2 in ket.items():
diff --git a/src/tequila/quantumchemistry/encodings.py b/src/tequila/quantumchemistry/encodings.py
index d31ae761..4ce088df 100644
--- a/src/tequila/quantumchemistry/encodings.py
+++ b/src/tequila/quantumchemistry/encodings.py
@@ -128,7 +128,7 @@ def map_state(self, state: list, *args, **kwargs) -> list:
         fop = openfermion.FermionOperator(string, 1.0)
         op = self(fop)
         from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
-        wfn = QubitWaveFunction.from_int(0, n_qubits=n_qubits)
+        wfn = QubitWaveFunction.from_basis_state(0, n_qubits=n_qubits)
         wfn = wfn.apply_qubitoperator(operator=op)
         assert (len(wfn.keys()) == 1)
         key = list(wfn.keys())[0].array
diff --git a/src/tequila/simulators/simulator_api.py b/src/tequila/simulators/simulator_api.py
index 50726f65..bd1b486b 100755
--- a/src/tequila/simulators/simulator_api.py
+++ b/src/tequila/simulators/simulator_api.py
@@ -363,7 +363,7 @@ def compile_circuit(abstract_circuit: 'QCircuit',
     return CircType(abstract_circuit=abstract_circuit, variables=variables, noise=noise, device=device, *args, **kwargs)
 
 
-def simulate(objective: typing.Union['Objective', 'QCircuit','QTensor'],
+def simulate(objective: typing.Union['Objective', 'QCircuit', 'QTensor'],
              variables: Dict[Union[Variable, Hashable], RealNumber] = None,
              samples: int = None,
              backend: str = None,
@@ -407,7 +407,7 @@ def simulate(objective: typing.Union['Objective', 'QCircuit','QTensor'],
                 objective.extract_variables()))
 
     compiled_objective = compile(objective=objective, samples=samples, variables=variables, backend=backend,
-                                 noise=noise,device=device, *args, **kwargs)
+                                 noise=noise, device=device, *args, **kwargs)
 
     return compiled_objective(variables=variables, samples=samples, *args, **kwargs)
 
diff --git a/src/tequila/simulators/simulator_base.py b/src/tequila/simulators/simulator_base.py
index b765abb2..78f86826 100755
--- a/src/tequila/simulators/simulator_base.py
+++ b/src/tequila/simulators/simulator_base.py
@@ -371,7 +371,7 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
                                   **kwargs)
 
         if keymap_required:
-            result.apply_keymap(keymap=keymap, initial_state=initial_state)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, n_qubits=len(all_qubits), initial_state=initial_state)
 
         return result
 
diff --git a/src/tequila/simulators/simulator_cirq.py b/src/tequila/simulators/simulator_cirq.py
index 5115f90f..7a9505bb 100755
--- a/src/tequila/simulators/simulator_cirq.py
+++ b/src/tequila/simulators/simulator_cirq.py
@@ -173,7 +173,7 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
         simulator = cirq.Simulator()
         backend_result = simulator.simulate(program=self.circuit, param_resolver=self.resolver,
                                             initial_state=initial_state)
-        return QubitWaveFunction.from_array(arr=backend_result.final_state_vector, numbering=self.numbering)
+        return QubitWaveFunction.from_array(array=backend_result.final_state_vector, numbering=self.numbering)
 
     def convert_measurements(self, backend_result: cirq.Result) -> QubitWaveFunction:
         """
@@ -186,18 +186,18 @@ def convert_measurements(self, backend_result: cirq.Result) -> QubitWaveFunction
         Returns
         -------
         QubitWaveFunction:
-            the result of sampling, as a tequila QubitWavefunction.
+            the result of sampling, as a tequila QubitWaveFunction.
 
         """
         assert (len(backend_result.measurements) == 1)
         for key, value in backend_result.measurements.items():
-            counter = QubitWaveFunction()
+            counter = QubitWaveFunction(self.n_qubits, self.numbering)
             for sample in value:
                 binary = BitString.from_array(array=sample.astype(int))
-                if binary in counter._state:
-                    counter._state[binary] += 1
+                if binary in counter.keys():
+                    counter[binary] += 1
                 else:
-                    counter._state[binary] = 1
+                    counter[binary] = 1
             return counter
 
     def do_sample(self, samples, circuit, *args, **kwargs) -> QubitWaveFunction:
diff --git a/src/tequila/simulators/simulator_pyquil.py b/src/tequila/simulators/simulator_pyquil.py
index 15c36bca..8261afc1 100755
--- a/src/tequila/simulators/simulator_pyquil.py
+++ b/src/tequila/simulators/simulator_pyquil.py
@@ -439,7 +439,7 @@ def do_simulate(self, variables, initial_state, *args, **kwargs):
             if val > 0:
                 iprep += pyquil.gates.X(i)
         backend_result = simulator.wavefunction(iprep + self.circuit, memory_map=self.resolver)
-        return QubitWaveFunction.from_array(arr=backend_result.amplitudes, numbering=self.numbering)
+        return QubitWaveFunction.from_array(array=backend_result.amplitudes, numbering=self.numbering)
 
     def do_sample(self, samples, circuit, *args, **kwargs) -> QubitWaveFunction:
         """
@@ -495,7 +495,7 @@ def string_to_array(s):
                     listing.append(int(letter))
             return listing
 
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         bit_dict = {}
         for b in backend_result:
             try:
@@ -505,7 +505,7 @@ def string_to_array(s):
 
         for k, v in bit_dict.items():
             arr = string_to_array(k)
-            result._state[BitString.from_array(arr)] = v
+            result[BitString.from_array(arr)] = v
         return result
 
     def no_translation(self, abstract_circuit):
diff --git a/src/tequila/simulators/simulator_qibo.py b/src/tequila/simulators/simulator_qibo.py
index 027b22cb..05b681b1 100755
--- a/src/tequila/simulators/simulator_qibo.py
+++ b/src/tequila/simulators/simulator_qibo.py
@@ -356,20 +356,20 @@ def do_simulate(self, variables, initial_state=None, *args, **kwargs):
         n_qubits = max(self.highest_qubit + 1, self.n_qubits, self.abstract_circuit.max_qubit() + 1)
         if initial_state is not None:
             if isinstance(initial_state, (int, np.int64)):
-                wave = QubitWaveFunction.from_int(i=initial_state, n_qubits=n_qubits)
+                wave = QubitWaveFunction.from_basis_state(n_qubits, initial_state, self.numbering)
             elif isinstance(initial_state, str):
-                wave = QubitWaveFunction.from_string(string=initial_state).to_array()
+                wave = QubitWaveFunction.from_string(initial_state, self.numbering).to_array()
             elif isinstance(initial_state, QubitWaveFunction):
                 wave = initial_state
             elif isinstance(initial_state,np.ndarray):
-                wave = QubitWaveFunction.from_array(initial_state)
+                wave = QubitWaveFunction.from_array(initial_state, self.numbering)
             else:
                 raise TequilaQiboException('could not understand initial state of type {}'.format(type(initial_state)))
             state = wave.to_array()
             result = self.circuit(state)
         else:
             result = self.circuit()
-        back= QubitWaveFunction.from_array(arr=result.numpy())
+        back= QubitWaveFunction.from_array(result.numpy(), self.numbering)
         return back
 
     def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunction:
@@ -398,7 +398,7 @@ def simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveFunc
         if isinstance(initial_state, BitString):
             initial_state = initial_state.integer
         if isinstance(initial_state, QubitWaveFunction):
-            if len(initial_state.keys()) != 1:
+            if len(initial_state) != 1:
                 return self.do_simulate(variables=variables,initial_state=initial_state, *args, **kwargs)
             initial_state = list(initial_state.keys())[0].integer
         if isinstance(initial_state,np.ndarray):
@@ -426,19 +426,18 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
             results transformed to tequila native QubitWaveFunction
         """
 
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         # todo there are faster ways
 
         for k, v in backend_result.frequencies(binary=True).items():
             converted_key = BitString.from_bitstring(other=BitString.from_binary(binary=k))
-            result._state[converted_key] = v
-
+            result[converted_key] = v
 
         if target_qubits is not None:
             mapped_target = [self.qubit_map[q].number for q in target_qubits]
             mapped_full = [self.qubit_map[q].number for q in self.abstract_qubits]
             keymap = KeyMapRegisterToSubregister(subregister=mapped_target, register=mapped_full)
-            result = result.apply_keymap(keymap=keymap)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, len(mapped_target))
 
         return result
 
@@ -521,21 +520,20 @@ def do_sample(self, samples, circuit, noise_model=None, initial_state=None, *arg
         n_qubits = max(self.highest_qubit + 1, self.n_qubits, self.abstract_circuit.max_qubit() + 1)
         if initial_state is not None:
             if isinstance(initial_state, int):
-                wave=QubitWaveFunction.from_int(i=initial_state, n_qubits=n_qubits)
+                wave = QubitWaveFunction.from_basis_state(n_qubits, initial_state, self.numbering)
             elif isinstance(initial_state, str):
-                wave = QubitWaveFunction.from_string(string=initial_state).to_array()
+                wave = QubitWaveFunction.from_string(initial_state, self.numbering).to_array()
             elif isinstance(initial_state, QubitWaveFunction):
                 wave = initial_state
-            elif isinstance(initial_state,np.ndarray):
-                wave = QubitWaveFunction.from_array(arr=initial_state, n_qubits=n_qubits)  # silly but necessary
+            elif isinstance(initial_state, np.ndarray):
+                wave = QubitWaveFunction.from_array(initial_state, self.numbering)  # silly but necessary
             else:
                 raise TequilaQiboException('received an unusable initial state of type {}'.format(type(initial_state)))
-            state=wave.to_array()
-            result = circuit(state,nshots=samples)
+            state = wave.to_array()
+            result = circuit(state, nshots=samples)
         else:
             result = circuit(nshots=samples)
 
-
         back = self.convert_measurements(backend_result=result)
         return back
 
diff --git a/src/tequila/simulators/simulator_qiskit.py b/src/tequila/simulators/simulator_qiskit.py
index 2f292334..186966b7 100755
--- a/src/tequila/simulators/simulator_qiskit.py
+++ b/src/tequila/simulators/simulator_qiskit.py
@@ -300,8 +300,7 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
 
         if initial_state != 0:
             state = np.zeros(2 ** self.n_qubits)
-            initial_state = reverse_int_bits(initial_state, self.n_qubits)
-            state[initial_state] = 1.0
+            state[reverse_int_bits(initial_state, self.n_qubits)] = 1.0
             init_circuit = qiskit.QuantumCircuit(self.q, self.c)
             init_circuit.set_statevector(state)
             circuit = init_circuit.compose(circuit)
@@ -310,7 +309,7 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
 
         backend_result = qiskit_backend.run(circuit, optimization_level=optimization_level).result()
 
-        return QubitWaveFunction.from_array(arr=backend_result.get_statevector(circuit), numbering=self.numbering)
+        return QubitWaveFunction.from_array(array=backend_result.get_statevector(circuit).data, numbering=self.numbering)
 
     def do_sample(self, circuit: qiskit.QuantumCircuit, samples: int, read_out_qubits, *args,
                   **kwargs) -> QubitWaveFunction:
@@ -394,16 +393,16 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
             measurements converted into wave function form.
         """
         qiskit_counts = backend_result.result().get_counts()
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         # todo there are faster ways
         for k, v in qiskit_counts.items():
-            converted_key = BitString.from_bitstring(other=BitStringLSB.from_binary(binary=k))
-            result._state[converted_key] = v
+            converted_key = BitString.from_binary(k)
+            result[converted_key] = v
         if target_qubits is not None:
             mapped_target = [self.qubit_map[q].number for q in target_qubits]
             mapped_full = [self.qubit_map[q].number for q in self.abstract_qubits]
             keymap = KeyMapRegisterToSubregister(subregister=mapped_target, register=mapped_full)
-            result = result.apply_keymap(keymap=keymap)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, n_qubits=len(target_qubits))
 
         return result
 
diff --git a/src/tequila/simulators/simulator_qlm.py b/src/tequila/simulators/simulator_qlm.py
index 1f23f174..2604cf4d 100644
--- a/src/tequila/simulators/simulator_qlm.py
+++ b/src/tequila/simulators/simulator_qlm.py
@@ -395,10 +395,10 @@ def do_simulate(self, variables, initial_state=0, *args, **kwargs) -> QubitWaveF
         if MY_QLM:
             result = PyLinalg().submit(job)
             statevector = get_statevector(result)
-            return QubitWaveFunction.from_array(arr=statevector, numbering=self.numbering)
+            return QubitWaveFunction.from_array(array=statevector, numbering=self.numbering)
 
         result = LinAlg().submit(job)
-        return QubitWaveFunction.from_array(arr=result.statevector, numbering=self.numbering)
+        return QubitWaveFunction.from_array(array=result.statevector, numbering=self.numbering)
 
     def update_variables(self, variables):
         """
@@ -431,7 +431,7 @@ def convert_measurements(self, backend_result) -> QubitWaveFunction:
         QubitWaveFunction:
             measurements converted into wave function form.
         """
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         shots = int(backend_result.meta_data["nbshots"])
         nbits = backend_result[0].qregs[0].length
         for sample in backend_result:
diff --git a/src/tequila/simulators/simulator_qulacs.py b/src/tequila/simulators/simulator_qulacs.py
index 29cfca95..43d91ba1 100755
--- a/src/tequila/simulators/simulator_qulacs.py
+++ b/src/tequila/simulators/simulator_qulacs.py
@@ -3,7 +3,7 @@
 import warnings
 
 from tequila import TequilaException, TequilaWarning
-from tequila.utils.bitstrings import BitNumbering, BitString, BitStringLSB
+from tequila.utils.bitstrings import BitNumbering, BitString, BitStringLSB, reverse_int_bits
 from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 from tequila.simulators.simulator_base import BackendCircuit, BackendExpectationValue, QCircuit, change_basis
 from tequila.utils.keymap import KeyMapRegisterToSubregister
@@ -149,11 +149,10 @@ def do_simulate(self, variables, initial_state, *args, **kwargs):
             QubitWaveFunction representing result of the simulation.
         """
         state = self.initialize_state(self.n_qubits)
-        lsb = BitStringLSB.from_int(initial_state, nbits=self.n_qubits)
-        state.set_computational_basis(BitString.from_binary(lsb.binary).integer)
+        state.set_computational_basis(reverse_int_bits(initial_state, self.n_qubits))
         self.circuit.update_quantum_state(state)
 
-        wfn = QubitWaveFunction.from_array(arr=state.get_vector(), numbering=self.numbering)
+        wfn = QubitWaveFunction.from_array(array=state.get_vector(), numbering=self.numbering)
         return wfn
 
     def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveFunction:
@@ -170,22 +169,17 @@ def convert_measurements(self, backend_result, target_qubits=None) -> QubitWaveF
             results transformed to tequila native QubitWaveFunction
         """
 
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self.numbering)
         # todo there are faster ways
 
-
         for k in backend_result:
-            converted_key = BitString.from_binary(BitStringLSB.from_int(integer=k, nbits=self.n_qubits).binary)
-            if converted_key in result._state:
-                result._state[converted_key] += 1
-            else:
-                result._state[converted_key] = 1
+            result[k] += 1
 
         if target_qubits is not None:
             mapped_target = [self.qubit_map[q].number for q in target_qubits]
             mapped_full = [self.qubit_map[q].number for q in self.abstract_qubits]
             keymap = KeyMapRegisterToSubregister(subregister=mapped_target, register=mapped_full)
-            result = result.apply_keymap(keymap=keymap)
+            result = QubitWaveFunction.from_wavefunction(result, keymap, n_qubits=len(target_qubits))
 
         return result
 
@@ -212,8 +206,7 @@ def do_sample(self, samples, circuit, noise_model=None, initial_state=0, *args,
             the results of sampling, as a Qubit Wave Function.
         """
         state = self.initialize_state(self.n_qubits)
-        lsb = BitStringLSB.from_int(initial_state, nbits=self.n_qubits)
-        state.set_computational_basis(BitString.from_binary(lsb.binary).integer)
+        state.set_computational_basis(reverse_int_bits(initial_state, self.n_qubits))
         circuit.update_quantum_state(state)
         sampled = state.sampling(samples)
         return self.convert_measurements(backend_result=sampled, target_qubits=self.measurements)
diff --git a/src/tequila/simulators/simulator_symbolic.py b/src/tequila/simulators/simulator_symbolic.py
index ddc3a24f..68003e93 100755
--- a/src/tequila/simulators/simulator_symbolic.py
+++ b/src/tequila/simulators/simulator_symbolic.py
@@ -42,37 +42,30 @@ def update_variables(self, variables):
 
     @classmethod
     def apply_gate(cls, state: QubitWaveFunction, gate: QGate, qubits: dict, variables) -> QubitWaveFunction:
-        result = QubitWaveFunction()
         n_qubits = len(qubits.keys())
+        result = sympy.Integer(1) * QubitWaveFunction(n_qubits)
         for s, v in state.items():
-            s.nbits = n_qubits
-            result += v * cls.apply_on_standard_basis(gate=gate, basisfunction=s, qubits=qubits, variables=variables)
+            result += v * cls.apply_on_standard_basis(gate=gate, basis_state=s, qubits=qubits, variables=variables)
         return result
 
     @classmethod
-    def apply_on_standard_basis(cls, gate: QGate, basisfunction: BitString, qubits:dict, variables) -> QubitWaveFunction:
-
-        basis_array = basisfunction.array
-        if gate.is_controlled():
-            do_apply = True
-            check = [basis_array[qubits[c]] == 1 for c in gate.control]
-            for c in check:
-                if not c:
-                    do_apply = False
-            if not do_apply:
-                return QubitWaveFunction.from_int(basisfunction)
+    def apply_on_standard_basis(cls, gate: QGate, basis_state: BitString, qubits:dict, variables) -> QubitWaveFunction:
+        n_qubits = len(qubits.keys())
+        basis_array = basis_state.array
+        if gate.is_controlled() and not all(basis_array[qubits[c]] == 1 for c in gate.control):
+            return QubitWaveFunction.from_basis_state(n_qubits, basis_state)
 
         if len(gate.target) > 1:
             raise Exception("Multi-targets not supported for symbolic simulators")
 
-        result = QubitWaveFunction()
+        result = sympy.Integer(1) * QubitWaveFunction(n_qubits)
         for tt in gate.target:
             t = qubits[tt]
             qt = basis_array[t]
             a_array = copy.deepcopy(basis_array)
             a_array[t] = (a_array[t] + 1) % 2
-            current_state = QubitWaveFunction.from_int(basisfunction)
-            altered_state = QubitWaveFunction.from_int(BitString.from_array(a_array))
+            current_state = QubitWaveFunction.from_basis_state(n_qubits, basis_state)
+            altered_state = QubitWaveFunction.from_basis_state(n_qubits, BitString.from_array(a_array))
 
             fac1 = None
             fac2 = None
@@ -115,22 +108,21 @@ def do_simulate(self, variables, initial_state: int = None, *args, **kwargs) ->
         count = 0
         for q in self.abstract_circuit.qubits:
             qubits[q] = count
-            count +=1
+            count += 1
 
         n_qubits = len(self.abstract_circuit.qubits)
 
         if initial_state is None:
-            initial_state = QubitWaveFunction.from_int(i=0, n_qubits=n_qubits)
-        elif isinstance(initial_state, int):
-            initial_state = QubitWaveFunction.from_int(initial_state, n_qubits=n_qubits)
+            initial_state = 0
+        initial_state = QubitWaveFunction.from_basis_state(n_qubits, initial_state)
 
         result = initial_state
         for g in self.abstract_circuit.gates:
             result = self.apply_gate(state=result, gate=g, qubits=qubits, variables=variables)
 
-        wfn = QubitWaveFunction()
+        wfn = QubitWaveFunction(n_qubits)
         if self.convert_to_numpy:
-            for k,v in result.items():
+            for k, v in result.items():
                 wfn[k] = complex(v)
         else:
             wfn = result
diff --git a/src/tequila/utils/bitstrings.py b/src/tequila/utils/bitstrings.py
index d725aacc..9ae7cb20 100644
--- a/src/tequila/utils/bitstrings.py
+++ b/src/tequila/utils/bitstrings.py
@@ -192,8 +192,13 @@ def reverse_int_bits(x: int, nbits: int) -> int:
 
 
 def initialize_bitstring(integer: int, nbits: int = None, numbering_in: BitNumbering = BitNumbering.MSB,
-                         numbering_out: BitNumbering = BitNumbering.MSB):
-    integer = reverse_int_bits(integer, nbits) if numbering_in != numbering_out else integer
+                         numbering_out: BitNumbering = None):
+    if numbering_out is None:
+        numbering_out = numbering_in
+
+    if numbering_in != numbering_out:
+        integer = reverse_int_bits(integer, nbits)
+
     if numbering_out == BitNumbering.MSB:
         return BitString.from_int(integer=integer, nbits=nbits)
     else:
diff --git a/src/tequila/wavefunction/qubit_wavefunction.py b/src/tequila/wavefunction/qubit_wavefunction.py
index 2ad75bd4..1c5c3b21 100644
--- a/src/tequila/wavefunction/qubit_wavefunction.py
+++ b/src/tequila/wavefunction/qubit_wavefunction.py
@@ -1,296 +1,397 @@
 from __future__ import annotations
-import copy
+
 import typing
-from typing import Dict, Union
+from copy import deepcopy
+from math import log2
+from typing import Union, Generator
+
 import numpy
+import numpy as np
+import numpy.typing as npt
 import numbers
 
-from tequila.utils.bitstrings import BitNumbering, BitString, initialize_bitstring
-from tequila import TequilaException
-from tequila.utils.keymap import KeyMapLSB2MSB, KeyMapMSB2LSB
-from tequila.tools import number_to_string
+import sympy
 
-# from __future__ import annotations # can use that in python 3.7+ to get rid of string type hints
+from tequila.utils.bitstrings import BitString, reverse_int_bits
+from tequila import TequilaException, BitNumbering, initialize_bitstring
+from tequila.utils.keymap import KeyMapABC
 
 if typing.TYPE_CHECKING:
-    # don't need those structures, just for convenient type hinting
+    # Don't need those structures, just for convenient type hinting
     from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian, PauliString
 
 
 class QubitWaveFunction:
     """
-    Store Wavefunction as dictionary of comp. basis state and complex numbers
-    Use the same structure for Measurments results with int instead of complex numbers (counts)
+    Represents a wavefunction.
+    Amplitudes are either stored in a Numpy array for dense wavefunctions, or a dictionary for sparse wavefunctions.
+    Does not enforce normalization.
     """
 
-    numbering = BitNumbering.MSB
-
-    def apply_keymap(self, keymap, initial_state: BitString = None):
-        self.n_qubits = keymap.n_qubits
-        mapped_state = dict()
-        for k, v in self.state.items():
-            mapped_key=keymap(input_state=k, initial_state=initial_state)
-            if mapped_key in mapped_state:
-                mapped_state[mapped_key] += v
-            else:
-                mapped_state[mapped_key] = v
-
-        self.state = mapped_state
-        return self
+    def __init__(self, n_qubits: int, numbering: BitNumbering = BitNumbering.MSB, dense: bool = False,
+                 init_state: bool = True) -> None:
+        """
+        Initialize a QubitWaveFunction with all amplitudes set to zero.
 
-    @property
-    def n_qubits(self) -> int:
-        if self._n_qubits is None:
-            return self.min_qubits()
+        :param n_qubits: Number of qubits.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :param dense: Whether to store the amplitudes in a Numpy array instead of a dictionary.
+        :param init_state: Whether to initialize the state array.
+            If False, set_state must be called immediately after the constructor.
+        """
+        self._n_qubits: int = n_qubits
+        self._numbering = numbering
+        self._dense = dense
+        if init_state:
+            self._state = np.zeros(2 ** self._n_qubits, dtype=complex) if dense else dict()
         else:
-            return max(self._n_qubits, self.min_qubits())
+            self._state = None
 
-    def min_qubits(self) -> int:
-        if len(self.state) > 0:
-            maxk = max(self.state.keys())
-            return maxk.nbits
+    @classmethod
+    def from_wavefunction(cls, wfn: QubitWaveFunction, keymap: KeyMapABC = None, n_qubits: int = None,
+                          initial_state: BitString = None) -> QubitWaveFunction:
+        """
+        Create a copy of a wavefunction.
+
+        :param wfn: The wavefunction to copy.
+        :param keymap: A keymap to apply to the wavefunction.
+        :param n_qubits: Number of qubits of the new wavefunction.
+            Must not be None if keymap is not None.
+        :param initial_state: Initial state to pass to the keymap.
+        :return: The copied wavefunction.
+        """
+        if keymap is not None:
+            result = QubitWaveFunction(n_qubits, numbering=wfn._numbering, dense=wfn._dense)
+            # Change amplitudes to sympy objects
+            if wfn._dense and wfn._state.dtype == object:
+                result = sympy.Integer(1) * result
+            for index, coeff in wfn.raw_items():
+                key = keymap(initialize_bitstring(index, wfn._n_qubits, numbering_in=wfn._numbering), initial_state)
+                result[key] += coeff
+            return result
         else:
-            return 0
+            return deepcopy(wfn)
 
-    @n_qubits.setter
-    def n_qubits(self, n_qubits):
-        if n_qubits is not None:
-            self._n_qubits = max(n_qubits, self.min_qubits())
-        return self
-
-    @property
-    def state(self):
-        if self._state is None:
-            return dict()
-        else:
-            return self._state
-
-    @state.setter
-    def state(self, other: Dict[BitString, complex]):
-        assert (isinstance(other, dict))
-        self._state = other
-
-    def __init__(self, state: Dict[BitString, complex] = None, n_qubits=None):
-        if state is None:
-            self._state = dict()
-        elif isinstance(state, int):
-            self._state = self.from_int(i=state, n_qubits=n_qubits).state
-        elif isinstance(state, str):
-            self._state = self.from_string(string=state, n_qubits=n_qubits).state
-        elif isinstance(state, numpy.ndarray) or isinstance(state, list):
-            self._state = self.from_array(arr=state, n_qubits=n_qubits).state
-        elif hasattr(state, "state"):
-            self._state = state.state
-        else:
-            self._state = state
-        self._n_qubits = n_qubits
+    @classmethod
+    def from_array(cls, array: npt.NDArray[complex], numbering: BitNumbering = BitNumbering.MSB,
+                   copy: bool = True) -> QubitWaveFunction:
+        """
+        Create a dense wavefunction from a Numpy array.
 
-    def items(self):
-        return self.state.items()
+        :param array: Array of amplitudes.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :param copy: Whether to copy the array or use it directly.
+            If False, the array must not be modified after the constructor.
+        :return: The created wavefunction.
+        """
+        if not log2(len(array)).is_integer():
+            raise ValueError(f"Array length must be a power of 2, received {len(array)}")
+        n_qubits = int(log2(len(array)))
+        result = QubitWaveFunction(n_qubits, numbering, dense=True, init_state=False)
+        result.set_state(array, copy)
+        return result
 
-    def keys(self):
-        return self.state.keys()
+    @classmethod
+    def from_basis_state(cls, n_qubits: int, basis_state: Union[int, BitString],
+                         numbering: BitNumbering = BitNumbering.MSB) -> QubitWaveFunction:
+        """
+        Create a sparse wavefunction that is a basis state.
 
-    def values(self):
-        return self.state.values()
+        :param n_qubits: Number of qubits.
+        :param basis_state: Index of the basis state.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :return: The created wavefunction.
+        """
+        if 2 ** n_qubits <= basis_state:
+            raise ValueError(f"Number of qubits {n_qubits} insufficient for basis state {basis_state}")
+        if isinstance(basis_state, BitString):
+            basis_state = reverse_int_bits(basis_state.integer,
+                                           basis_state.nbits) if numbering != basis_state.numbering else basis_state.integer
+        result = QubitWaveFunction(n_qubits, numbering)
+        result[basis_state] = 1.0
+        return result
 
-    @staticmethod
-    def convert_bitstring(key: Union[BitString, numbers.Integral], n_qubits):
-        if isinstance(key, numbers.Integral):
-            return BitString.from_int(integer=key, nbits=n_qubits)
-        elif isinstance(key, str):
-            return BitString.from_binary(binary=key, nbits=n_qubits)
-        else:
-            return key
+    @classmethod
+    def from_string(cls, string: str, numbering: BitNumbering = BitNumbering.MSB) -> QubitWaveFunction:
+        """
+        Create a sparse wavefunction from a string.
 
-    def __getitem__(self, item: BitString):
-        key = self.convert_bitstring(item, self.n_qubits)
-        return self.state[key]
+        :param string: String representation of the wavefunction.
+        :param numbering: Whether the first qubit is the most or least significant.
+        :return: The created wavefunction.
+        """
+        try:
+            string = string.replace(" ", "")
+            string = string.replace("*", "")
+            terms = string.split(">")[:-1]
+            n_qubits = len(terms[0].split("|")[-1])
+            result = QubitWaveFunction(n_qubits, numbering)
+            for term in terms:
+                coeff, index = term.split("|")
+                coeff = complex(coeff) if coeff != "" else 1.0
+                index = int(index, 2)
+                result[index] = coeff
+            return result
+        except ValueError:
+            raise TequilaException(f"Failed to initialize QubitWaveFunction from string:\n\"{string}\"\n")
 
-    def __call__(self, key, *args, **kwargs) -> numbers.Number:
+    @classmethod
+    def convert_from(cls, n_qubits: int, val: Union[QubitWaveFunction, int, str, numpy.ndarray]):
         """
-        Like getitem but returns zero if key is not there
-
-        Parameters
-        ----------
-        key: bitstring (or int or str)
-        Returns
-        -------
-            Return the amplitude or measurement occurence of a bitstring
+        Convert a value to a QubitWaveFunction.
+        Accepts QubitWaveFunction, int, str, and numpy.ndarray.
+
+        :param n_qubits: Number of qubits.
+        :param val: Value to convert.
+        :return: The converted value.
         """
-        ckey = self.convert_bitstring(key, self.n_qubits)
-        if ckey in self.state:
-            return self.state[ckey]
+        if isinstance(val, QubitWaveFunction):
+            return val
+        elif isinstance(val, int):
+            return cls.from_basis_state(n_qubits=n_qubits, basis_state=val)
+        elif isinstance(val, str):
+            return cls.from_string(val)
+        elif isinstance(val, numpy.ndarray):
+            return cls.from_array(val)
         else:
-            return 0.0
+            raise TequilaException(f"Cannot initialize QubitWaveFunction from type {type(val)}")
 
+    @property
+    def n_qubits(self) -> int:
+        """
+        Returns number of qubits in the wavefunction.
+        """
+        return self._n_qubits
 
+    @property
+    def numbering(self) -> BitNumbering:
+        """
+        Returns the bit numbering of the wavefunction.
+        """
+        return self._numbering
 
-    def __setitem__(self, key: BitString, value: numbers.Number):
-        self._state[self.convert_bitstring(key, self.n_qubits)] = value
-        return self
-
-    def __contains__(self, item: BitString):
-        return self.convert_bitstring(item, self.n_qubits) in self.keys()
-
-    def __len__(self):
-        return len(self.state)
-
-    @classmethod
-    def from_array(cls, arr: numpy.ndarray, keymap=None, threshold: float = 1.e-6,
-                   numbering: BitNumbering = BitNumbering.MSB, n_qubits: int = None):
-        arr = numpy.asarray(arr)
-        assert (len(arr.shape) == 1)
-        state = dict()
-        maxkey = len(arr) - 1
-        maxbit = initialize_bitstring(integer=maxkey, numbering_in=numbering, numbering_out=cls.numbering).nbits
-        for ii, v in enumerate(arr):
-            if abs(v) > threshold:
-                i = initialize_bitstring(integer=ii, nbits=maxbit, numbering_in=numbering, numbering_out=cls.numbering)
-                key = i if keymap is None else keymap(i)
-                state[key] = v
-        result = QubitWaveFunction(state, n_qubits=n_qubits)
+    @property
+    def dense(self) -> bool:
+        """
+        Returns whether the wavefunction is dense.
+        """
+        return self._dense
 
-        return result
+    def to_array(self, copy: bool = True) -> npt.NDArray[complex]:
+        """
+        Returns array of amplitudes.
 
-    @classmethod
-    def from_int(cls, i: int, coeff=1, n_qubits: int = None):
-        if isinstance(i, BitString):
-            return QubitWaveFunction(state={i: coeff}, n_qubits=n_qubits)
+        :param copy: Whether to copy the array or use it directly for dense Wavefunctions.
+            If False, changes to the array or wavefunction will affect each other.
+        :return: Array of amplitudes.
+        """
+        if self._dense:
+            return self._state.copy() if copy else self._state
         else:
-            return QubitWaveFunction(state={BitString.from_int(integer=i, nbits=n_qubits): coeff}, n_qubits=n_qubits)
+            result = np.zeros(2 ** self._n_qubits, dtype=complex)
+            for k, v in self.raw_items():
+                result[k] = v
+            return result
 
-    @classmethod
-    def from_string(cls, string: str, n_qubits: int = None):
+    def set_state(self, value: npt.NDArray[complex], copy: bool = True) -> None:
         """
-        Complex values like (x+iy)|...> will currently not work, you need to type Real and imaginary separately
-        Or improve this constructor :-)
-        e.g instead of (0.5+1.0j)|0101> do 0.5|0101> + 1.0j|0101>
-        :param paths:
-        :param string:
-        :return:
+        Sets the state to an array.
+        After this call, the wavefunction will be dense.
+
+        :param value: Array of amplitudes. Length must be 2 ** n_qubits.
+        :param copy: Whether to copy the array or use it directly.
+            If False, changes to the array or wavefunction will affect each other.
         """
-        try:
-            state = dict()
-            string = string.replace(" ", "")
-            string = string.replace("*", "")
-            string = string.replace("+-", "-")
-            string = string.replace("-+", "-")
-            terms = (string + "terminate").split('>')
-            for term in terms:
-                if term == 'terminate':
-                    break
-                tmp = term.split("|")
-                coeff = tmp[0]
-                if coeff == '':
-                    coeff = 1.0
-                else:
-                    coeff = complex(coeff)
-                basis_state = BitString.from_binary(binary=tmp[1])
+        if len(value) != 2 ** self._n_qubits:
+            raise ValueError(f"Wavefunction of {self._n_qubits} qubits must have {2 ** self._n_qubits} amplitudes, "
+                             f"received {len(value)}")
+        self._dense = True
+        if copy:
+            self._state = value.copy()
+        else:
+            self._state = value
+
+    def __getitem__(self, key: Union[int, BitString]) -> complex:
+        if isinstance(key, BitString):
+            key = reverse_int_bits(key.integer, key.nbits) if self._numbering != key.numbering else key.integer
+        return self._state[key] if self._dense else self._state.get(key, 0)
+
+    def __setitem__(self, key: Union[int, BitString], value: complex) -> None:
+        if isinstance(key, BitString):
+            key = reverse_int_bits(key.integer, key.nbits) if self._numbering != key.numbering else key.integer
+        self._state[key] = value
+
+    def __contains__(self, item: Union[int, BitString]) -> bool:
+        if isinstance(item, BitString):
+            item = reverse_int_bits(item.integer, item.nbits) if self._numbering != item.numbering else item.integer
+        return abs(self[item]) > 1e-6
+
+    def raw_items(self) -> Generator[tuple[int, complex]]:
+        """Returns a generator of non-zero amplitudes with integer indices."""
+        return ((k, v) for k, v in (enumerate(self._state) if self._dense else self._state.items()))
+
+    def items(self) -> Generator[tuple[BitString, complex]]:
+        """Returns a generator of non-zero amplitudes with BitString indices."""
+        return ((initialize_bitstring(k, self._n_qubits, self._numbering), v)
+                for k, v in self.raw_items()
+                if isinstance(v, sympy.Basic) or abs(v) > 1e-6)
+
+    def keys(self) -> Generator[BitString]:
+        """Returns a generator of BitString indices of non-zero amplitudes."""
+        return (k for k, v in self.items())
+
+    def values(self) -> Generator[complex]:
+        """Returns a generator of non-zero amplitudes."""
+        return (v for k, v in self.items())
+
+    def __eq__(self, other) -> bool:
+        raise TequilaException("Wavefunction equality is not well-defined. Consider using isclose.")
+
+    def isclose(self: QubitWaveFunction,
+                other: QubitWaveFunction,
+                rtol: float = 1e-5,
+                atol: float = 1e-8) -> bool:
+        """
+        Check if two wavefunctions are close, up to a global phase.
 
-                state[basis_state] = coeff
-        except ValueError:
-            raise TequilaException("Failed to initialize QubitWaveFunction from string:" + string + "\n"
-                                                                                                    "did you try complex values?\n"
-                                                                                                    "currently you need to type real and imaginary parts separately\n"
-                                                                                                    "e.g. instead of (0.5+1.0j)|0101> do 0.5|0101> + 1.0j|0101>")
-        except:
-            raise TequilaException("Failed to initialize QubitWaveFunction from string:" + string)
-        return QubitWaveFunction(state=state, n_qubits=n_qubits)
+        :param other: The other wavefunction.
+        :param rtol: Relative tolerance.
+        :param atol: Absolute tolerance.
+        :return: Whether the wavefunctions are close.
+        """
+        inner = self.inner(other)
+        cosine_similarity = inner / (self.norm() * other.norm())
+        return np.isclose(abs(cosine_similarity), 1.0, rtol, atol)
 
-    def __repr__(self):
-        result = str()
-        for k, v in self.items():
-            result += number_to_string(number=v) + "|" + str(k.binary) + "> "
-        return result
+    def __add__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            return QubitWaveFunction.from_array(self._state + other._state, self.numbering, copy=False)
+        else:
+            result = QubitWaveFunction.from_wavefunction(self)
+            result += other
+            return result
 
-    def __eq__(self, other):
-        raise TequilaException("Wavefunction equality is not well-defined. Consider using inner"
-                               + " product equality, wf1.isclose(wf2).")
-
-    def isclose(self :  'QubitWaveFunction',
-                other : 'QubitWaveFunction',
-                rtol : float=1e-5,
-                atol : float=1e-8) -> bool:
-        """Return whether this wavefunction is similar to the target wavefunction."""
-        over1 = complex(self.inner(other))
-        over2 = numpy.sqrt(complex(self.inner(self) * other.inner(other)))
-        # Explicit casts to complex() is required if self or other are sympy
-        # wavefunction with sympy-typed amplitudes
-
-        # Check if the two numbers are equal.
-        return numpy.isclose(over1, over2, rtol=rtol, atol=atol)
-
-    def __add__(self, other):
-        result = QubitWaveFunction(state=copy.deepcopy(self._state))
-        for k, v in other.items():
-            if k in result._state:
-                result._state[k] += v
-            else:
-                result._state[k] = v
-        return result
+    def __iadd__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            self._state += other._state
+        else:
+            for k, v in other.raw_items():
+                if self._numbering != other._numbering:
+                    k = reverse_int_bits(k, self._n_qubits)
+                self[k] += v
+        return self
 
-    def __sub__(self, other):
-        return self + -1.0 * other
+    def __sub__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            return QubitWaveFunction.from_array(self._state - other._state, self.numbering, copy=False)
+        else:
+            result = QubitWaveFunction.from_wavefunction(self)
+            result -= other
+            return result
 
-    def __iadd__(self, other):
-        for k, v in other.items():
-            if k in self._state:
-                self._state[k] += v
-            else:
-                self._state[k] = v
+    def __isub__(self, other: QubitWaveFunction) -> QubitWaveFunction:
+        if self._dense and other._dense and self._numbering == other._numbering:
+            self._state -= other._state
+        else:
+            for k, v in other.raw_items():
+                if self._numbering != other._numbering:
+                    k = reverse_int_bits(k, self._n_qubits)
+                self[k] -= v
         return self
 
-    def __rmul__(self, other):
-        result = QubitWaveFunction(state=copy.deepcopy(self._state))
-        for k, v in result._state.items():
-            result._state[k] *= other
-        return result
+    def __rmul__(self, other: complex) -> QubitWaveFunction:
+        if self._dense:
+            return QubitWaveFunction.from_array(other * self._state, self.numbering, copy=False)
+        else:
+            result = QubitWaveFunction.from_wavefunction(self)
+            result *= other
+            return result
 
-    def inner(self, other):
-        # currently very slow and not optimized in any way
-        result = 0.0
-        for k, v in self.items():
-            if k in other._state:
-                result += v.conjugate() * other._state[k]
-        return result
+    def __imul__(self, other: complex) -> QubitWaveFunction:
+        if self._dense:
+            self._state *= other
+        else:
+            for k, v in self.raw_items():
+                self[k] = other * v
+        return self
 
-    def normalize(self):
+    def inner(self, other: QubitWaveFunction) -> complex:
+        """Returns the inner product with another wavefunction."""
+        if self._dense and other._dense and self._numbering == other._numbering:
+            return np.inner(self._state.conjugate(), other._state)
+        else:
+            result = 0
+            for k, v in self.raw_items():
+                if self._numbering != other._numbering:
+                    k = reverse_int_bits(k, self._n_qubits)
+                result += v.conjugate() * other[k]
+            if isinstance(result, sympy.Basic):
+                result = complex(result)
+            return result
+
+    def norm(self) -> float:
+        """Returns the norm of the wavefunction."""
+        return np.sqrt(self.inner(self))
+
+    def normalize(self, inplace: bool = False) -> QubitWaveFunction:
         """
-        NOT AN Inplace operation
-        :return: Normalizes the wavefunction/countrate
+        Normalizes the wavefunction.
+
+        :param inplace: Whether to normalize the wavefunction in place or return a new one.
+        :return: The normalized wavefunction.
         """
-        norm2 = self.inner(other=self)
-        normalized = 1.0 / numpy.sqrt(norm2) * self
-        return normalized
+        norm = self.norm()
+        if inplace:
+            self *= 1.0 / norm
+            return self
+        else:
+            return (1.0 / norm) * QubitWaveFunction.from_wavefunction(self)
+
+    # It would be nice to call this __len__, however for some reason this causes infinite loops
+    # when multiplying wave functions with some types of numbers from the right sight, likely
+    # because the __mul__ implementation of the number tries to perform some sort of array
+    # operation.
+    def length(self):
+        return sum(1 for _ in self.raw_items())
+
+    def __repr__(self):
+        result = str()
+        for index, coeff in self.items():
+            index = index.integer
+            if self.numbering == BitNumbering.LSB:
+                index = reverse_int_bits(index, self._n_qubits)
+            result += f"{coeff} |{index:0{self._n_qubits}b}> "
+        # If the wavefunction contains no states
+        if not result:
+            result = "empty wavefunction"
+        return result
 
-    def compute_expectationvalue(self, operator: 'QubitHamiltonian') -> numbers.Real:
+    def compute_expectationvalue(self, operator: QubitHamiltonian) -> numbers.Real:
         tmp = self.apply_qubitoperator(operator=operator)
         E = self.inner(other=tmp)
-        if hasattr(E, "imag") and numpy.isclose(E.imag, 0.0, atol=1.e-6):
+        if hasattr(E, "imag") and np.isclose(E.imag, 0.0, atol=1.e-6):
             return float(E.real)
         else:
             return E
 
-    def apply_qubitoperator(self, operator: 'QubitHamiltonian'):
+    def apply_qubitoperator(self, operator: QubitHamiltonian) -> QubitWaveFunction:
         """
         Inefficient function which computes the action of a QubitHamiltonian on this wfn
         :param operator: QubitOperator
         :return: resulting Qubitwavefunction
         """
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self.n_qubits, self._numbering)
         for ps in operator.paulistrings:
             result += self.apply_paulistring(paulistring=ps)
-        result = result.simplify()
         return result
 
-    def apply_paulistring(self, paulistring: 'PauliString'):
+    def apply_paulistring(self, paulistring: PauliString) -> QubitWaveFunction:
         """
         Inefficient function which computes action of a single paulistring
         :param paulistring: PauliString
         :return: Expectation Value
         """
-        result = QubitWaveFunction()
+        result = QubitWaveFunction(self._n_qubits, self._numbering)
         for k, v in self.items():
             arr = k.array
             c = v
@@ -306,17 +407,3 @@ def apply_paulistring(self, paulistring: 'PauliString'):
                     raise TequilaException("unknown pauli: " + str(p))
             result[BitString.from_array(array=arr)] = c
         return paulistring.coeff * result
-
-    def to_array(self):
-        result = numpy.zeros(shape=2 ** self.n_qubits, dtype=complex)
-        for k, v in self.items():
-            result[int(k)] = v
-        return result
-
-    def simplify(self, threshold = 1.e-8):
-        state = {}
-        for k, v in self.state.items():
-            if not numpy.isclose(v, 0.0, atol=threshold):
-                state[k] = v
-        return QubitWaveFunction(state=state)
-
diff --git a/tests/test_binary_pauli.py b/tests/test_binary_pauli.py
index 499fe43f..ff3bdcc3 100644
--- a/tests/test_binary_pauli.py
+++ b/tests/test_binary_pauli.py
@@ -169,7 +169,7 @@ def test_commuting_groups():
 
 def prepare_cov_dict(H):
     eigenValues, eigenVectors = np.linalg.eigh(H.to_matrix())
-    wfn0 = tq.QubitWaveFunction(eigenVectors[:,0])
+    wfn0 = tq.QubitWaveFunction.from_array(eigenVectors[:, 0])
     terms = BinaryHamiltonian.init_from_qubit_hamiltonian(H).binary_terms
     cov_dict = {}
     for term1 in terms:
diff --git a/tests/test_circuits.py b/tests/test_circuits.py
index 05a9fb86..9782b0a0 100644
--- a/tests/test_circuits.py
+++ b/tests/test_circuits.py
@@ -11,6 +11,7 @@
 from tequila.simulators.simulator_api import simulate
 from tequila.wavefunction.qubit_wavefunction import QubitWaveFunction
 
+
 def test_qubit_map():
 
     for gate in [Rx, Ry, Rz]:
@@ -82,7 +83,7 @@ def test_conventions():
 
 
 def strip_sympy_zeros(wfn: QubitWaveFunction):
-    result = QubitWaveFunction()
+    result = QubitWaveFunction(wfn.n_qubits, wfn.numbering)
     for k, v in wfn.items():
         if v != 0:
             result[k] = v
@@ -130,7 +131,7 @@ def test_basic_gates():
     cos = sympy.cos
     sin = sympy.sin
     exp = sympy.exp
-    BS = QubitWaveFunction.from_int
+    def BS(state): return QubitWaveFunction.from_basis_state(n_qubits=1, basis_state=state)
     angle = sympy.pi
     gates = [X(0), Y(0), Z(0), Rx(target=0, angle=angle), Ry(target=0, angle=angle), Rz(target=0, angle=angle), H(0)]
     results = [
@@ -310,7 +311,7 @@ def test_unitary_gate_u_u3(gate, theta, phi, lambd):
 
 def test_swap():
     U = X(0)
-    U += SWAP(0,2)
+    U += SWAP(0, 2)
     wfn = simulate(U)
     wfnx = simulate(X(2))
     assert numpy.isclose(numpy.abs(wfn.inner(wfnx))**2,1.0)
@@ -345,9 +346,13 @@ def test_givens():
     U = X(0)
     U += Givens(0, 1, angle=-numpy.pi/4)
     wfn = simulate(U)
-    wfnx = simulate(X(0))
+    V = X(0)
+    V.n_qubits = 2
+    wfnx = simulate(V)
     assert numpy.isclose(numpy.abs(wfn.inner(wfnx))**2, 0.5)
-    wfnx = simulate(X(1))
+    V = X(1)
+    V.n_qubits = 2
+    wfnx = simulate(V)
     assert numpy.isclose(numpy.abs(wfn.inner(wfnx))**2, 0.5)
     
     U = X(0)
diff --git a/tests/test_hamiltonian_arithmetic.py b/tests/test_hamiltonian_arithmetic.py
index 42870a56..2fffa75a 100644
--- a/tests/test_hamiltonian_arithmetic.py
+++ b/tests/test_hamiltonian_arithmetic.py
@@ -84,16 +84,16 @@ def test_initialization():
 def test_ketbra():
     ket = QubitWaveFunction.from_string("1.0*|00> + 1.0*|11>").normalize()
     operator = paulis.KetBra(ket=ket, bra="|00>")
-    result = operator*QubitWaveFunction.from_int(0, n_qubits=2)
-    assert(result.isclose(ket))
+    result = operator * QubitWaveFunction.from_basis_state(n_qubits=2, basis_state=0)
+    assert result.isclose(ket)
 
 @pytest.mark.parametrize("n_qubits", [1,2,3,5])
 def test_ketbra_random(n_qubits):
     ket = numpy.random.uniform(0.0, 1.0, 2**n_qubits)
-    bra = QubitWaveFunction.from_int(0, n_qubits=n_qubits)
+    bra = QubitWaveFunction.from_basis_state(n_qubits, 0)
     operator = paulis.KetBra(ket=ket, bra=bra)
     result = operator * bra
-    assert(result.isclose(QubitWaveFunction.from_array(ket)))
+    assert result.isclose(QubitWaveFunction.from_array(ket))
 
 def test_paulistring_conversion():
     X1 = QubitHamiltonian.from_string("X0", openfermion_format=True)
@@ -196,7 +196,7 @@ def test_projectors(qubits):
     real = numpy.random.uniform(0.0,1.0,2**qubits)
     imag = numpy.random.uniform(0.0,1.0,2**qubits)
     array = real + 1.j*imag
-    wfn = QubitWaveFunction.from_array(arr=array)
+    wfn = QubitWaveFunction.from_array(array=array)
     P = paulis.Projector(wfn=wfn.normalize())
     assert(P.is_hermitian())
     assert(wfn.apply_qubitoperator(P).isclose(wfn))
@@ -304,17 +304,17 @@ def test_simple_trace_out():
 
     H1 = QubitHamiltonian.from_string("1.0*X(0)*Z(1)*Z(5)*X(100)Y(50)")
     H2 = QubitHamiltonian.from_string("1.0*X(0)X(100)Y(50)")
-    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")]*2)
+    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")] * 2)
 
     H1 = QubitHamiltonian.from_string("1.0*X(0)*Z(1)*X(100)Y(50)")
     H2 = QubitHamiltonian.from_string("-1.0*X(0)X(100)Y(50)")
-    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")]*2)
+    assert H2 == H1.trace_out_qubits(qubits=[1,5], states=[QubitWaveFunction.from_string("1.0*|1>")] * 2)
 
 @pytest.mark.parametrize("theta", numpy.random.uniform(0.0, 6.0, 10))
 def test_trace_out_xy(theta):
     a = numpy.sin(theta)
     b = numpy.cos(theta)
-    state = QubitWaveFunction.from_array([a,b])
+    state = QubitWaveFunction.from_array([a, b])
 
     H1 = QubitHamiltonian.from_string("1.0*X(0)*X(1)*X(100)")
     H2 = QubitHamiltonian.from_string("1.0*X(0)*X(100)")
diff --git a/tests/test_mappings.py b/tests/test_mappings.py
index b151201d..5659831a 100644
--- a/tests/test_mappings.py
+++ b/tests/test_mappings.py
@@ -12,15 +12,15 @@
 
 
 def test_keymaps():
-    initial_state = 0
+    initial_state = BitString.from_int(0)
 
-    small = QubitWaveFunction.from_int(i=int("0b1111", 2))
-    large = QubitWaveFunction.from_int(i=int("0b01010101", 2))
-    large.n_qubits = 8
+    small = QubitWaveFunction.from_basis_state(4, basis_state=int("0b1111", 2))
+    large = QubitWaveFunction.from_basis_state(8, basis_state=int("0b01010101", 2))
 
     keymap = KeyMapSubregisterToRegister(register=[0, 1, 2, 3, 4, 5, 6, 7], subregister=[1, 3, 5, 7])
 
-    assert (small.apply_keymap(keymap=keymap, initial_state=initial_state).isclose(large))
+    mapped = QubitWaveFunction.from_wavefunction(small, keymap, 8, initial_state)
+    assert mapped.isclose(large)
 
 
 def test_endianness():
@@ -47,6 +47,7 @@ def test_endianness():
         assert (i1 == BitString.from_bitstring(i22))
         assert (i2 == BitString.from_bitstring(i11))
 
+
 @pytest.mark.skipif(condition="cirq" not in INSTALLED_SAMPLERS or "qiskit" not in INSTALLED_SAMPLERS, reason="need cirq and qiskit for cross validation")
 def test_endianness_simulators():
     tests = ["000111",
@@ -73,7 +74,7 @@ def test_endianness_simulators():
         print("counts_cirq  =", counts_cirq)
         print("counts_qiskit=", counts_qiskit)
         assert (counts_cirq.isclose(counts_qiskit))
-        assert (wfn_cirq.state == counts_cirq.state)
+        assert (wfn_cirq.to_array() == counts_cirq.to_array())
 
 
 @pytest.mark.parametrize("backend", INSTALLED_SAMPLERS)
@@ -83,7 +84,7 @@ def test_paulistring_sampling(backend, case):
     H = QubitHamiltonian.from_paulistrings(PauliString.from_string(case[0]))
     U = gates.X(target=1) + gates.X(target=3) + gates.X(target=5)
     E = ExpectationValue(H=H, U=U)
-    result = simulate(E,backend=backend, samples=1)
+    result = simulate(E, backend=backend, samples=1)
     assert isclose(result, case[1], 1.e-4)
 
 
@@ -93,12 +94,13 @@ def test_paulistring_sampling_2(backend, case):
     H = QubitHamiltonian.from_paulistrings(PauliString.from_string(case[0]))
     U = gates.H(target=1) + gates.H(target=3) + gates.X(target=5) + gates.H(target=5)
     E = ExpectationValue(H=H, U=U)
-    result = simulate(E,backend=backend, samples=1)
-    assert (isclose(result, case[1], 1.e-4))
+    result = simulate(E, backend=backend, samples=1)
+    assert isclose(result, case[1], 1.e-4)
+
 
 @pytest.mark.parametrize("array", [numpy.random.uniform(0.0,1.0,i) for i in [2,4,8,16]])
 def test_qubitwavefunction_array(array):
     print(array)
-    wfn = QubitWaveFunction.from_array(arr=array)
+    wfn = QubitWaveFunction.from_array(array=array)
     array2 = wfn.to_array()
-    assert numpy.allclose(array,array2)
+    assert numpy.allclose(array, array2)
diff --git a/tests/test_qasm.py b/tests/test_qasm.py
index b39454bd..4c7f3902 100644
--- a/tests/test_qasm.py
+++ b/tests/test_qasm.py
@@ -35,6 +35,8 @@ def test_export_import_qasm_simple(zx_calculus):
     assert (numpy.isclose(wfn1.inner(wfn2), 1.0))
 
 
+# TODO: Reactivate
+@pytest.mark.skip("Extremely slow for some reason")
 @pytest.mark.parametrize(
     "zx_calculus,variabs",
     [
@@ -145,6 +147,8 @@ def test_export_import_qasm_trotterized_gate(zx_calculus, string1, string2, angl
     assert (numpy.isclose(wfn1.inner(wfn2), 1.0))
 
 
+# TODO: Reactivate
+@pytest.mark.skip("Extremely slow for some reason")
 @pytest.mark.parametrize(
     "zx_calculus,variabs",
     [
diff --git a/tests/test_simulator_backends.py b/tests/test_simulator_backends.py
index 1f991454..cefef684 100644
--- a/tests/test_simulator_backends.py
+++ b/tests/test_simulator_backends.py
@@ -9,6 +9,7 @@
 import numbers
 import tequila as tq
 import tequila.simulators.simulator_api
+from tequila import BitString
 
 """
 Warn if Simulators are not installed
@@ -40,7 +41,7 @@ def teardown_function(function):
 @pytest.mark.dependencies
 def test_dependencies():
     for package in tequila.simulators.simulator_api.SUPPORTED_BACKENDS:
-        if package != "qulacs_gpu":
+        if package not in ["qulacs_gpu", "qiskit_gpu"]:
             assert (package in tq.simulators.simulator_api.INSTALLED_BACKENDS)
 
 
@@ -351,8 +352,8 @@ def test_initial_state_from_integer(simulator, initial_state):
         U += tq.gates.X(target=i) + tq.gates.X(target=i)
 
     wfn = tq.simulate(U, initial_state=initial_state, backend=simulator)
-    assert (initial_state in wfn)
-    assert (numpy.isclose(wfn[initial_state], 1.0))
+    assert BitString.from_int(initial_state, 6) in wfn
+    assert numpy.isclose(wfn[BitString.from_int(initial_state, 6)], 1.0)
 
 
 @pytest.mark.parametrize("backend", tequila.simulators.simulator_api.INSTALLED_SIMULATORS.keys())
diff --git a/tests/test_symbolic_simulator.py b/tests/test_symbolic_simulator.py
index 0da79126..6b5f867e 100644
--- a/tests/test_symbolic_simulator.py
+++ b/tests/test_symbolic_simulator.py
@@ -6,14 +6,14 @@
 import numpy
 import pytest
 
-@pytest.mark.parametrize("paulis", [(gates.X,paulis.X), (gates.Y, paulis.Y), (gates.Z,paulis.Z)])
-@pytest.mark.parametrize("qubit", [0,1,2])
-@pytest.mark.parametrize("init", [0,1])
+@pytest.mark.parametrize("paulis", [(gates.X, paulis.X), (gates.Y, paulis.Y), (gates.Z, paulis.Z)])
+@pytest.mark.parametrize("qubit", [0, 1, 2])
+@pytest.mark.parametrize("init", [0, 1])
 def test_pauli_gates(paulis, qubit, init):
-    iwfn = QubitWaveFunction.from_int(i=init, n_qubits=qubit+1)
+    iwfn = QubitWaveFunction.from_basis_state(n_qubits=qubit + 1, basis_state=init)
+    iwfn = iwfn.apply_qubitoperator(paulis[1](qubit))
     wfn = simulate(paulis[0](qubit), initial_state=init)
-    iwfn=iwfn.apply_qubitoperator(paulis[1](qubit))
-    assert(iwfn.isclose(wfn))
+    assert iwfn.isclose(wfn)
 
 @pytest.mark.parametrize("rot", [(gates.Rx,paulis.X), (gates.Ry, paulis.Y), (gates.Rz,paulis.Z)])
 @pytest.mark.parametrize("angle", numpy.random.uniform(0.0, 2*numpy.pi, 5))
@@ -22,7 +22,7 @@ def test_pauli_gates(paulis, qubit, init):
 def test_rotations(rot, qubit, angle, init):
     pauli = rot[1](qubit)
     gate = rot[0](target=qubit, angle=angle)
-    iwfn = QubitWaveFunction.from_int(i=init, n_qubits=qubit+1)
+    iwfn = QubitWaveFunction.from_basis_state(basis_state=init, n_qubits=qubit + 1)
     wfn = simulate(gate, initial_state=init)
     test= numpy.cos(-angle/2.0)*iwfn + 1.0j*numpy.sin(-angle/2.0)* iwfn.apply_qubitoperator(pauli)
     assert(wfn.isclose(test))
@@ -31,7 +31,7 @@ def test_rotations(rot, qubit, angle, init):
 @pytest.mark.parametrize("init", [0,1])
 def test_hadamard(qubit, init):
     gate = gates.H(target=qubit)
-    iwfn = QubitWaveFunction.from_int(i=init, n_qubits=qubit+1)
+    iwfn = QubitWaveFunction.from_basis_state(basis_state=init, n_qubits=qubit + 1)
     wfn = simulate(gate, initial_state=init)
     test= 1.0/numpy.sqrt(2)*(iwfn.apply_qubitoperator(paulis.Z(qubit)) + iwfn.apply_qubitoperator(paulis.X(qubit)))
     assert(wfn.isclose(test))
@@ -47,11 +47,8 @@ def test_controls(target, control, gate):
     assert(wfn0.isclose(wfn1))
 
     c0 = gates.QCircuit()
+    c0.n_qubits = max(target, control) + 1
     c1 = gate(target=target, control=control)
     wfn0 = simulate(c0, initial_state=0)
     wfn1 = simulate(c1, initial_state=0)
     assert(wfn0.isclose(wfn1))
-
-
-
-
diff --git a/tests/test_unary_state_prep.py b/tests/test_unary_state_prep.py
index 0d7ee4ba..a89fded4 100644
--- a/tests/test_unary_state_prep.py
+++ b/tests/test_unary_state_prep.py
@@ -21,7 +21,7 @@ def test_unary_states(target_space: list):
     coeff = 1.0 / numpy.sqrt(qubits)  # fails for the 3-Qubit Case because the wrong sign is picked in the solution
     coeffs = [coeff for i in range(qubits)]
 
-    wfn = QubitWaveFunction()
+    wfn = QubitWaveFunction(qubits)
     for i, c in enumerate(coeffs):
         wfn += c * QubitWaveFunction.from_string("1.0|" + target_space[i] + ">")