Repository: qBraid/pyqasm

Issue:
[BUG] `pyqasm.unroll()` handling of include statements
Right now, `pyqasm.unroll()` seems to be adding an `include "stdgates.inc"` statement to all programs, even if the original qasm program didn't have it. This is problematic when using together with a library that doesn't support "include" statements, e.g. Amazon Braket:

```python
from braket.circuits import Circuit
from braket.circuits.serialization import IRType
from braket.ir.openqasm import Program

from pyqasm import dumps, loads

bk_bell = Circuit().h(0).cnot(0, 1).measure(0).measure(1)

qasm_str = bk_bell.to_ir(IRType.OPENQASM).source

program = Program(source=qasm_str)

# Reconstruct the original Circuit from QASM string - this works as expected
circuit = Circuit.from_ir(source=program.source, inputs=program.inputs)

module = loads(qasm_str)

module.unroll()

qasm_unrolled = dumps(module)

program_unrolled = Program(source=qasm_unrolled)

# Attempt to reconstruct a Circuit from the unrolled QASM string - this raises a FileNotFoundError
# Error: FileNotFoundError: [Errno 2] No such file or directory: 'stdgates.inc'
circuit_unrolled = Circuit.from_ir(source=program_unrolled.source, inputs=program_unrolled.inputs)
```

The `include "stdgates.inc"` statement should only be added to the unrolled qasm if it was present in the original, input program.

Another related issue, let's say I unroll a qasm string with 2 custom include statements:

```python
from pyqasm import dumps, loads

custom_qasm = """
OPENQASM 3.0;
include "customgates.inc";
include "othergates.inc";
bit[2] b;
qubit[2] q;
h q[0];
cnot q[0], q[1];
b[0] = measure q[0];
b[1] = measure q[1];
"""

module = loads(custom_qasm)

module.unroll()

qasm_unrolled = dumps(module)
```
```
OPENQASM 3.0;
include "stdgates.inc";
bit[2] b;
qubit[2] q;
h q[0];
cx q[0], q[1];
b[0] = measure q[0];
b[1] = measure q[1];
```

Using `pyqasm.unroll()` all existing include statements are overwritten by `include "stdgates.inc";`, which shouldn't happen.

Discussion:

Solution patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c63219d..fb46b5e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,107 +15,17 @@ Types of changes:
 ## Unreleased
 
 ### Added
-- Added a `dumps` and `formatted_qasm` method to the `QasmModule` class to allow for the conversion of a `QasmModule` object to a string representation of the QASM code ([#71](https://github.com/qBraid/pyqasm/pull/71))
-- Added the `populate_idle_qubits` method to the `QasmModule` class to populate idle qubits with an `id` gate ([#72](https://github.com/qBraid/pyqasm/pull/72))
-- Added gate definitions for "c3sqrtx", "u1", "rxx", "cu3", "csx", "rccx" , "ch" , "cry", "cp", "cu", "cu1", "rzz" in `maps.py` ([#74](https://github.com/qBraid/pyqasm/pull/74))
-- Added support for skipping the unrolling for externally linked gates. The `QasmModule.unroll()` method now accepts an `external_gates` parameter which is a list of gate names that should not be unrolled ([#59](https://github.com/qBraid/pyqasm/pull/59)). Usage - 
-
-```python
-In [30]: import pyqasm
-
-In [31]: qasm_str = """OPENQASM 3.0;
-    ...:     include "stdgates.inc";
-    ...:     gate custom q1, q2, q3{
-    ...:         x q1;
-    ...:         y q2;
-    ...:         z q3;
-    ...:     }
-    ...:
-    ...:     qubit[4] q;
-    ...:     custom q[0], q[1], q[2];
-    ...:     cx q[1], q[2];"""
-
-In [32]: module = pyqasm.loads(qasm_str)
-
-In [33]: module.unroll(external_gates= ["custom"])
-
-In [34]: pyqasm.dumps(module).splitlines()
-Out[34]:
-['OPENQASM 3.0;',
- 'include "stdgates.inc";',
- 'qubit[4] q;',
- 'custom q[0], q[1], q[2];',
- 'cx q[1], q[2];']
-```
-- **Major Change**: Added the `load`, `loads`, `dump`, and `dumps` functions to the `pyqasm` module to allow for the loading and dumping of QASM code ([#76](https://github.com/qBraid/pyqasm/pull/76)). Usage - 
-
-```python
-In [18]: import pyqasm
-
-In [19]: qasm_str = """OPENQASM 3.0;
-    ...:     include "stdgates.inc";
-    ...:     qreg q1[2];
-    ...:     qubit[2] q2;"""
-
-In [20]: module = pyqasm.loads(qasm_str)
-
-In [21]: print(pyqasm.dumps(module))
-OPENQASM 3.0;
-include "stdgates.inc";
-qubit[2] q1;
-qubit[2] q2;
-
-
-In [22]: file_path = "test.qasm"
-
-In [23]: pyqasm.dump(module, file_path)
-
-In [24]: module = pyqasm.load(file_path)
-
-In [25]: print(pyqasm.dumps(module))
-OPENQASM 3.0;
-include "stdgates.inc";
-qubit[2] q1;
-qubit[2] q2;
-```
-- Added definitions for various gates in `maps.py` and tests for qasm formatting functions of the qbraid-sdk ([#82](https://github.com/qBraid/pyqasm/pull/82), [#84](https://github.com/qBraid/pyqasm/pull/84))
-- Added `pyqasm.accelerate` module to hold `.pyx` files with Cython-based optimizations for computationally intensive functions  ([#83](https://github.com/qBraid/pyqasm/pull/83))
-- Added `has_barriers` method for checking if a `QasmModule` object contains barriers ([#85](https://github.com/qBraid/pyqasm/pull/85))
-- Added `pyqasm.cli` module with `typer` integration to enable using pyqasm as a command-line tool ([#87](https://github.com/qBraid/pyqasm/pull/87))
-
-```bash
-$ pip install 'pyqasm[cli]'
-$ pyqasm --help
-Usage: pyqasm [OPTIONS] COMMAND [ARGS]...
-$ pyqasm --version
-pyqasm/0.1.0a1
-$ pyqasm validate tests/cli/resources
-tests/cli/resources/invalid1.qasm: error: Index 2 out of range for register of size 1 in qubit [validation]
-Found errors in 1 file (checked 3 source files)
-$ pyqasm validate tests/cli/resources --skip tests/cli/resources/invalid1.qasm
-Success: no issues found in 2 source files
-```
+- Added support for `gphase`, `toffoli`, `not`, `c3sx` and `c4x` gates ([#86](https://github.com/qBraid/pyqasm/pull/86))
 
 ### Improved / Modified
-- Changed the `__init__` method for the `QasmModule` class to only accept an `openqasm3.ast.Program` object as input ([#71](https://github.com/qBraid/pyqasm/pull/71))
-- Changed `DepthNode`, `QubitDepthNode`, `ClbitDepthNode`, and `Variable` to dataclasses. `__repr__` method is therefore handled automatically and you don't need all of the redundant private / public attribute and setters ([#79](https://github.com/qBraid/pyqasm/pull/79))
-- Simplified `map_qasm_op_to_callable` redundant `KeyError` handling with loop ([#79](https://github.com/qBraid/pyqasm/pull/79))
-- The `load` function has been renamed to `loads` and `load` is now used to load a QASM file. `QasmModule.dumps()` has been replaced with `__str__` method ([#76](https://github.com/qBraid/pyqasm/pull/76))
-- Experimental Cython integration:  ([#83](https://github.com/qBraid/pyqasm/pull/83))
-    - Migrated `pyqasm.linalg._kronecker_factor` to `pyqasm.linalg_cy` with ~60% speedup
-    - Migrated `pyqasm.linalg._so4_to_so2()` to to `pyqasm.linalg_cy` with ~5% speedup
-- Changed source code directory from `./pyqasm` to `./src/pyqasm` to prevents conflicts between the local source directory and the installed package in site-packages, ensuring Python's module resolution prioritizes the correct version. Required for local testing with new Cython build step ([#83](https://github.com/qBraid/pyqasm/pull/83))
-- Updated the build process for `pyqasm` due to Cython integration. Wheels are now built for each platform and uploaded to PyPI along with the source distributions ([#88](https://github.com/qBraid/pyqasm/pull/88))
+ - Refactored the initialization of `QasmModule` to remove default include statements. Only user supplied include statements are now added to the generated QASM code ([#86](https://github.com/qBraid/pyqasm/pull/86))
 
 ### Deprecated
 
 ### Removed
-- Removed the `from_program` method from the `QasmModule` class ([#71](https://github.com/qBraid/pyqasm/pull/71))
-- `QasmModule.formatted_qasm()` method has been removed ([#76](https://github.com/qBraid/pyqasm/pull/76))
 
 ### Fixed
-- Updated docs custom CSS used for sphinx to make version stable/latest drop-down visible. Previously was set white so blended into background and wasn't visible. ([#78](https://github.com/qBraid/pyqasm/pull/78))
-- Fixed bug in `pyqasm.linalg.so_bidiagonalize()` in final dot product order ([#83](https://github.com/qBraid/pyqasm/pull/83))
+- Fixed bugs in implementations of `gpi2` and `prx` gates ([#86](https://github.com/qBraid/pyqasm/pull/86))
 
 ### Dependencies
 
@@ -123,6 +33,7 @@ Success: no issues found in 2 source files
 
 Archive of changelog entries from previous releases:
 
+- [v0.1.0-alpha](https://github.com/qBraid/pyqasm/releases/tag/v0.1.0-alpha)
 - [v0.0.3](https://github.com/qBraid/pyqasm/releases/tag/v0.0.3)
 - [v0.0.2](https://github.com/qBraid/pyqasm/releases/tag/v0.0.2)
 - [v0.0.1](https://github.com/qBraid/pyqasm/releases/tag/v0.0.1)
diff --git a/src/pyqasm/maps.py b/src/pyqasm/maps.py
index a10dd47..89752a7 100644
--- a/src/pyqasm/maps.py
+++ b/src/pyqasm/maps.py
@@ -32,6 +32,7 @@
     QuantumGate,
     QuantumGateDefinition,
     QuantumMeasurementStatement,
+    QuantumPhase,
     QuantumReset,
     QubitDeclaration,
     SubroutineDefinition,
@@ -159,6 +160,24 @@ def u2_inv_gate(phi, lam, qubits) -> list[QuantumGate]:
     return u3_inv_gate(CONSTANTS_MAP["pi"] / 2, phi, lam, qubits)
 
 
+def global_phase_gate(theta: float, qubit_list: list[IndexedIdentifier]) -> list[QuantumPhase]:
+    """
+    Builds a global phase gate with the given theta and qubit list.
+
+    Args:
+        theta (float): The phase angle.
+        qubit_list (list[IndexedIdentifier]): The list of qubits on which to apply the phase.
+
+    Returns:
+        list[QuantumPhase]: A QuantumPhase object representing the global phase gate.
+    """
+    return [
+        QuantumPhase(
+            argument=FloatLiteral(value=theta), qubits=qubit_list, modifiers=[]  # type: ignore
+        )
+    ]
+
+
 def sxdg_gate_op(qubit_id) -> list[QuantumGate]:
     """
     Implements the conjugate transpose of the Sqrt(X) gate as a decomposition of other gates.
@@ -593,7 +612,7 @@ def cu1_gate(
         q_0: ┤ U(0,0,theta/2) ├──■───────────────────────■────────────────────
              └────────────────┘┌─┴─┐┌─────────────────┐┌─┴─┐┌────────────────┐
         q_1: ──────────────────┤ X ├┤ U(0,0,-theta/2) ├┤ X ├┤ U(0,0,theta/2) ├
-                             └───┘└─────────────────┘└───┘└────────────────┘
+                               └───┘└─────────────────┘└───┘└────────────────┘
     """
     result: list[QuantumGate] = []
     result.extend(u3_gate(0, 0, theta / 2, qubit0))
@@ -645,13 +664,13 @@ def csx_gate(qubit0: IndexedIdentifier, qubit1: IndexedIdentifier) -> list[Quant
 
 def rxx_gate(
     theta: Union[int, float], qubit0: IndexedIdentifier, qubit1: IndexedIdentifier
-) -> list[QuantumGate]:
+) -> list[Union[QuantumGate, QuantumPhase]]:
     """
     Implements the RXX gate as a decomposition of other gates.
     """
-    # TODO: to implement global phase of -theta / 2
 
-    result: list[QuantumGate] = []
+    result: list[Union[QuantumGate, QuantumPhase]] = []
+    result.extend(global_phase_gate(-theta / 2, [qubit0, qubit1]))
     result.extend(one_qubit_gate_op("h", qubit0))
     result.extend(one_qubit_gate_op("h", qubit1))
     result.extend(two_qubit_gate_op("cx", qubit0, qubit1))
@@ -682,7 +701,7 @@ def rccx_gate(
 
 def rzz_gate(
     theta: Union[int, float], qubit0: IndexedIdentifier, qubit1: IndexedIdentifier
-) -> list[QuantumGate]:
+) -> list[Union[QuantumGate, QuantumPhase]]:
     """
     Implements the RZZ gate as a decomposition of other gates.
 
@@ -705,11 +724,9 @@ def rzz_gate(
         q_1: ┤ X ├┤ U3(0,0,theta) ├┤ X ├
              └───┘└───────────────┘└───┘
     """
-    result: list[QuantumGate] = []
+    result: list[Union[QuantumGate, QuantumPhase]] = []
 
-    # TODO: verify implementation of the global phase
-    # result.extend(phaseshift_gate(-theta/2, qubit0))
-    # result.extend(phaseshift_gate(-theta/2, qubit1))
+    result.extend(global_phase_gate(-theta / 2, [qubit0, qubit1]))
     result.extend(two_qubit_gate_op("cx", qubit0, qubit1))
     result.extend(u3_gate(0, 0, theta, qubit1))
     result.extend(two_qubit_gate_op("cx", qubit0, qubit1))
@@ -818,9 +835,12 @@ def gpi2_gate(phi, qubit_id) -> list[QuantumGate]:
     """
     Implements the gpi2 gate as a decomposition of other gates.
     """
+    # Reference:
+    # https://amazon-braket-sdk-python.readthedocs.io/en/latest/_apidoc/braket.circuits.circuit.html#braket.circuits.circuit.Circuit.gpi2
+    # https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.U3Gate#u3gate
     theta_0 = CONSTANTS_MAP["pi"] / 2
-    phi_0 = phi + 3 * CONSTANTS_MAP["pi"] / 2
-    lambda_0 = -phi_0 + CONSTANTS_MAP["pi"] / 2
+    phi_0 = phi - CONSTANTS_MAP["pi"] / 2
+    lambda_0 = CONSTANTS_MAP["pi"] / 2 - phi
     return u3_gate(theta_0, phi_0, lambda_0, qubit_id)
 
 
@@ -903,13 +923,106 @@ def ecr_gate(qubit0: IndexedIdentifier, qubit1: IndexedIdentifier) -> list[Quant
     return result
 
 
+def c3sx_gate(
+    qubit0: IndexedIdentifier,
+    qubit1: IndexedIdentifier,
+    qubit2: IndexedIdentifier,
+    qubit3: IndexedIdentifier,
+) -> list[QuantumGate]:
+    """
+    Implements the c3sx gate as a decomposition of other gates.
+
+    Uses the following qiskit decomposition -
+
+    In [15]: qc.draw()
+    Out[15]:
+
+    q_0: ──■───
+           │
+    q_1: ──■───
+           │
+    q_2: ──■───
+         ┌─┴──┐
+    q_3: ┤ Sx ├
+         └────┘
+
+    In [16]: qc.decompose().draw()
+    Out[16]:
+
+    q_0: ──────■──────────■────────────────────■────────────────────────────────────────■────────
+               │        ┌─┴─┐                ┌─┴─┐                                      │
+    q_1: ──────┼────────┤ X ├──────■─────────┤ X ├──────■──────────■────────────────────┼────────
+               │        └───┘      │         └───┘      │        ┌─┴─┐                ┌─┴─┐
+    q_2: ──────┼───────────────────┼────────────────────┼────────┤ X ├──────■─────────┤ X ├──────
+         ┌───┐ │U1(π/8) ┌───┐┌───┐ │U1(-π/8) ┌───┐┌───┐ │U1(π/8) ├───┤┌───┐ │U1(-π/8) ├───┤┌───┐
+    q_3: ┤ H ├─■────────┤ H ├┤ H ├─■─────────┤ H ├┤ H ├─■────────┤ H ├┤ H ├─■─────────┤ H ├┤ H ├─
+         └───┘          └───┘└───┘           └───┘└───┘          └───┘└───┘           └───┘└───┘
+    «
+    «q_0:─────────────────────────────────■──────────────────────
+    «                                     │
+    «q_1:────────────■────────────────────┼──────────────────────
+    «              ┌─┴─┐                ┌─┴─┐
+    «q_2:─■────────┤ X ├──────■─────────┤ X ├──────■─────────────
+    «     │U1(π/8) ├───┤┌───┐ │U1(-π/8) ├───┤┌───┐ │U1(π/8) ┌───┐
+    «q_3:─■────────┤ H ├┤ H ├─■─────────┤ H ├┤ H ├─■────────┤ H ├
+    «              └───┘└───┘           └───┘└───┘          └───┘
+    """
+
+    result: list[QuantumGate] = []
+    result.extend(one_qubit_gate_op("h", qubit3))
+    result.extend(cu1_gate(CONSTANTS_MAP["pi"] / 8, qubit0, qubit3))
+    result.extend(two_qubit_gate_op("cx", qubit0, qubit1))
+    # h(q[3]) * h (q[3]) = Identity
+    result.extend(cu1_gate(-CONSTANTS_MAP["pi"] / 8, qubit1, qubit3))
+    result.extend(two_qubit_gate_op("cx", qubit0, qubit1))
+    # h(q[3]) * h (q[3]) = Identity
+    result.extend(cu1_gate(CONSTANTS_MAP["pi"] / 8, qubit1, qubit3))
+    result.extend(two_qubit_gate_op("cx", qubit1, qubit2))
+    # h(q[3]) * h (q[3]) = Identity
+    result.extend(cu1_gate(-CONSTANTS_MAP["pi"] / 8, qubit2, qubit3))
+    result.extend(two_qubit_gate_op("cx", qubit0, qubit2))
+    # h(q[3]) * h (q[3]) = Identity
+    result.extend(cu1_gate(CONSTANTS_MAP["pi"] / 8, qubit2, qubit3))
+    result.extend(two_qubit_gate_op("cx", qubit1, qubit2))
+    # h(q[3]) * h (q[3]) = Identity
+    result.extend(cu1_gate(-CONSTANTS_MAP["pi"] / 8, qubit2, qubit3))
+    result.extend(two_qubit_gate_op("cx", qubit0, qubit2))
+    # h(q[3]) * h (q[3]) = Identity
+    result.extend(cu1_gate(CONSTANTS_MAP["pi"] / 8, qubit2, qubit3))
+    result.extend(one_qubit_gate_op("h", qubit3))
+
+    return result
+
+
+def c4x_gate(
+    qubit0: IndexedIdentifier,
+    qubit1: IndexedIdentifier,
+    qubit2: IndexedIdentifier,
+    qubit3: IndexedIdentifier,
+) -> list[QuantumGate]:
+    """
+    Implements the c4x gate
+    """
+    return [
+        QuantumGate(
+            modifiers=[],
+            name=Identifier(name="c4x"),
+            arguments=[],
+            qubits=[qubit0, qubit1, qubit2, qubit3],
+        )
+    ]
+
+
 def prx_gate(theta, phi, qubit_id) -> list[QuantumGate]:
     """
     Implements the PRX gate as a decomposition of other gates.
     """
+    # Reference:
+    # https://amazon-braket-sdk-python.readthedocs.io/en/latest/_apidoc/braket.circuits.circuit.html#braket.circuits.circuit.Circuit.prx
+    # https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.U3Gate#u3gate
     theta_0 = theta
-    phi_0 = CONSTANTS_MAP["pi"] / 2 - phi
-    lambda_0 = -phi_0
+    phi_0 = phi - CONSTANTS_MAP["pi"] / 2
+    lambda_0 = CONSTANTS_MAP["pi"] / 2 - phi
     return u3_gate(theta_0, phi_0, lambda_0, qubit_id)
 
 
@@ -954,6 +1067,7 @@ def two_qubit_gate_op(
     "id": lambda qubit_id: one_qubit_gate_op("id", qubit_id),
     "h": lambda qubit_id: one_qubit_gate_op("h", qubit_id),
     "x": lambda qubit_id: one_qubit_gate_op("x", qubit_id),
+    "not": lambda qubit_id: one_qubit_gate_op("x", qubit_id),
     "y": lambda qubit_id: one_qubit_gate_op("y", qubit_id),
     "z": lambda qubit_id: one_qubit_gate_op("z", qubit_id),
     "s": lambda qubit_id: one_qubit_gate_op("s", qubit_id),
@@ -999,12 +1113,12 @@ def two_qubit_gate_op(
     "ch": ch_gate,
     "xx": rxx_gate,
     "rxx": rxx_gate,
-    "xy": xy_gate,
-    "xx_plus_yy": xx_plus_yy_gate,
     "yy": ryy_gate,
     "ryy": ryy_gate,
-    "zz": zz_gate,
+    "zz": rzz_gate,
     "rzz": rzz_gate,
+    "xy": xy_gate,
+    "xx_plus_yy": xx_plus_yy_gate,
     "pswap": pswap_gate,
     "iswap": iswap_gate,
     "cp": cphaseshift_gate,
@@ -1028,11 +1142,18 @@ def two_qubit_gate_op(
 
 THREE_QUBIT_OP_MAP = {
     "ccx": ccx_gate_op,
+    "toffoli": ccx_gate_op,
     "ccnot": ccx_gate_op,
     "cswap": cswap_gate,
     "rccx": rccx_gate,
 }
 
+FOUR_QUBIT_OP_MAP = {"c3sx": c3sx_gate, "c3sqrtx": c3sx_gate}
+
+FIVE_QUBIT_OP_MAP = {
+    "c4x": c4x_gate,
+}
+
 
 def map_qasm_op_to_callable(op_name: str) -> tuple[Callable, int]:
     """
@@ -1052,6 +1173,8 @@ def map_qasm_op_to_callable(op_name: str) -> tuple[Callable, int]:
         (ONE_QUBIT_ROTATION_MAP, 1),
         (TWO_QUBIT_OP_MAP, 2),
         (THREE_QUBIT_OP_MAP, 3),
+        (FOUR_QUBIT_OP_MAP, 4),
+        (FIVE_QUBIT_OP_MAP, 5),
     ]
 
     for op_map, qubit_count in op_maps:
diff --git a/src/pyqasm/modules/base.py b/src/pyqasm/modules/base.py
index a99a8ce..fc08b29 100644
--- a/src/pyqasm/modules/base.py
+++ b/src/pyqasm/modules/base.py
@@ -17,7 +17,7 @@
 from typing import Optional
 
 import openqasm3.ast as qasm3_ast
-from openqasm3.ast import Include, Program
+from openqasm3.ast import Program
 
 from pyqasm.analyzer import Qasm3Analyzer
 from pyqasm.elements import ClbitDepthNode, QubitDepthNode
@@ -48,7 +48,7 @@ def __init__(self, name: str, program: Program):
         self._has_measurements: Optional[bool] = None
         self._has_barriers: Optional[bool] = None
         self._validated_program = False
-        self._unrolled_ast = Program(statements=[Include("stdgates.inc")])
+        self._unrolled_ast = Program(statements=[])
 
     @property
     def name(self) -> str:
@@ -128,7 +128,7 @@ def has_measurements(self):
             # the presence of measurements
             stmts_to_check = (
                 self._unrolled_ast.statements
-                if len(self._unrolled_ast.statements) > 1
+                if len(self._unrolled_ast.statements) > 0
                 else self._statements
             )
             for stmt in stmts_to_check:
@@ -148,7 +148,7 @@ def remove_measurements(self, in_place: bool = True):
         """
         stmt_list = (
             self._statements
-            if len(self._unrolled_ast.statements) == 1
+            if len(self._unrolled_ast.statements) == 0
             else self._unrolled_ast.statements
         )
         stmts_without_meas = [
@@ -187,7 +187,7 @@ def has_barriers(self):
             # the presence of barriers
             stmts_to_check = (
                 self._unrolled_ast.statements
-                if len(self._unrolled_ast.statements) > 1
+                if len(self._unrolled_ast.statements) > 0
                 else self._statements
             )
             for stmt in stmts_to_check:
@@ -207,7 +207,7 @@ def remove_barriers(self, in_place: bool = True):
         """
         stmt_list = (
             self._statements
-            if len(self._unrolled_ast.statements) == 1
+            if len(self._unrolled_ast.statements) == 0
             else self._unrolled_ast.statements
         )
         stmts_without_barriers = [
@@ -393,6 +393,7 @@ def remove_idle_qubits(self, in_place: bool = True):
                     if isinstance(stmt, qasm3_ast.QubitDeclaration):
                         if stmt.qubit.name == reg_name:
                             qasm_module._unrolled_ast.statements.remove(stmt)
+                            qasm_module._statements.remove(stmt)
                             break
 
                 del qasm_module._qubit_registers[reg_name]
@@ -441,7 +442,7 @@ def reverse_qubit_order(self, in_place=True):
         #    the depth maps here
 
         # 2. replace each qubit index in the Quantum Operations with the new index
-        for operation in qasm_module.unrolled_ast.statements:
+        for operation in qasm_module._unrolled_ast.statements:
             if isinstance(operation, QUANTUM_STATEMENTS):
                 bit_list = Qasm3Analyzer.get_op_bit_list(operation)
                 for bit in bit_list:
@@ -457,7 +458,7 @@ def reverse_qubit_order(self, in_place=True):
                     bit.indices[0][0].value = -1 * new_reg_idx - 1
 
         # remove the -ve marker
-        for operation in qasm_module.unrolled_ast.statements:
+        for operation in qasm_module._unrolled_ast.statements:
             if isinstance(operation, QUANTUM_STATEMENTS):
                 bit_list = Qasm3Analyzer.get_op_bit_list(operation)
                 for bit in bit_list:
@@ -495,9 +496,7 @@ def unroll(self, **kwargs):
         except (ValidationError, UnrollError) as err:
             # reset the unrolled ast and qasm
             self.num_qubits, self.num_clbits = -1, -1
-            self._unrolled_ast = Program(
-                statements=[Include("stdgates.inc")], version=self.original_program.version
-            )
+            self._unrolled_ast = Program(statements=[], version=self.original_program.version)
             raise err
 
     def __str__(self) -> str:
@@ -507,7 +506,7 @@ def __str__(self) -> str:
             str: The string representation of the module
         """
 
-        if len(self.unrolled_ast.statements) > 1:
+        if len(self._unrolled_ast.statements) > 1:
             return self._qasm_ast_to_str(self.unrolled_ast)
         return self._qasm_ast_to_str(self.original_program)
 
diff --git a/src/pyqasm/modules/qasm2.py b/src/pyqasm/modules/qasm2.py
index b8644f2..f17e55e 100644
--- a/src/pyqasm/modules/qasm2.py
+++ b/src/pyqasm/modules/qasm2.py
@@ -37,7 +37,7 @@ class Qasm2Module(QasmModule):
 
     def __init__(self, name: str, program: Program):
         super().__init__(name, program)
-        self._unrolled_ast = Program(statements=[Include("qelib1.inc")], version="2.0")
+        self._unrolled_ast = Program(statements=[], version="2.0")
         self._whitelist_statements = {
             qasm3_ast.BranchingStatement,
             qasm3_ast.QubitDeclaration,
@@ -102,6 +102,6 @@ def accept(self, visitor):
         """
         self._filter_statements()
         unrolled_stmt_list = visitor.visit_basic_block(self._statements)
-        self.unrolled_ast.statements = [Include("qelib1.inc")]  # pylint: disable=W0201
-        self.unrolled_ast.statements.extend(unrolled_stmt_list)
-        # TODO: some finalizing method here probably
+        final_stmt_list = visitor.finalize(unrolled_stmt_list)
+
+        self.unrolled_ast.statements = final_stmt_list
diff --git a/src/pyqasm/modules/qasm3.py b/src/pyqasm/modules/qasm3.py
index a19182c..fb193c4 100644
--- a/src/pyqasm/modules/qasm3.py
+++ b/src/pyqasm/modules/qasm3.py
@@ -12,7 +12,7 @@
 Defines a module for handling OpenQASM 3.0 programs.
 """
 
-from openqasm3.ast import Include, Program
+from openqasm3.ast import Program
 from openqasm3.printer import dumps
 
 from pyqasm.modules.base import QasmModule
@@ -30,7 +30,7 @@ class Qasm3Module(QasmModule):
 
     def __init__(self, name: str, program: Program):
         super().__init__(name, program)
-        self._unrolled_ast = Program(statements=[Include("stdgates.inc")], version="3.0")
+        self._unrolled_ast = Program(statements=[], version="3.0")
 
     def _qasm_ast_to_str(self, qasm_ast):
         """Convert the qasm AST to a string"""
@@ -45,6 +45,6 @@ def accept(self, visitor):
             visitor (QasmVisitor): The visitor to accept
         """
         unrolled_stmt_list = visitor.visit_basic_block(self._statements)
-        self.unrolled_ast.statements = [Include("stdgates.inc")]  # pylint: disable=W0201
-        self.unrolled_ast.statements.extend(unrolled_stmt_list)
-        # TODO: some finalizing method here probably
+        final_stmt_list = visitor.finalize(unrolled_stmt_list)
+
+        self._unrolled_ast.statements = final_stmt_list
diff --git a/src/pyqasm/transformer.py b/src/pyqasm/transformer.py
index c3f3d70..7b5019d 100644
--- a/src/pyqasm/transformer.py
+++ b/src/pyqasm/transformer.py
@@ -12,6 +12,7 @@
 Module with transformation functions for QASM3 visitor
 
 """
+from copy import deepcopy
 from typing import Any, Optional, Union
 
 import numpy as np
@@ -30,6 +31,7 @@
 from openqasm3.ast import (
     QuantumBarrier,
     QuantumGate,
+    QuantumPhase,
     QuantumReset,
     RangeDefinition,
     UintType,
@@ -133,7 +135,7 @@ def get_qubits_from_range_definition(
 
     @staticmethod
     def transform_gate_qubits(
-        gate_op: QuantumGate, qubit_map: dict[str, IndexedIdentifier]
+        gate_op: Union[QuantumGate, QuantumPhase], qubit_map: dict[str, IndexedIdentifier]
     ) -> None:
         """Transform the qubits of a gate operation with a qubit map.
 
@@ -144,6 +146,9 @@ def transform_gate_qubits(
         Returns:
             None
         """
+        if isinstance(gate_op, QuantumPhase) and len(gate_op.qubits) == 0:
+            gate_op.qubits = deepcopy(list(qubit_map.values()))
+            return
         for i, qubit in enumerate(gate_op.qubits):
             if isinstance(qubit, IndexedIdentifier):
                 raise_qasm3_error(
@@ -195,7 +200,7 @@ def transform_expression(expression, variable_map: dict[str, Union[int, float, b
 
     @staticmethod
     def transform_gate_params(
-        gate_op: QuantumGate, param_map: dict[str, Union[int, float, bool]]
+        gate_op: Union[QuantumGate, QuantumPhase], param_map: dict[str, Union[int, float, bool]]
     ) -> None:
         """Transform the parameters of a gate operation with a parameter map.
 
@@ -211,10 +216,13 @@ def transform_gate_params(
 
         # param map is a "global dict for this gate" which contains the binding of the params
         # to the actual values used in the call
-        for i, actual_arg in enumerate(gate_op.arguments):
-            # recursively replace ALL instances of the parameter in the expression
-            # with the actual value
-            gate_op.arguments[i] = Qasm3Transformer.transform_expression(actual_arg, param_map)
+        if isinstance(gate_op, QuantumGate):
+            for i, actual_arg in enumerate(gate_op.arguments):
+                # recursively replace ALL instances of the parameter in the expression
+                # with the actual value
+                gate_op.arguments[i] = Qasm3Transformer.transform_expression(actual_arg, param_map)
+        else:
+            gate_op.argument = Qasm3Transformer.transform_expression(gate_op.argument, param_map)
 
     @staticmethod
     def get_branch_params(condition: Any) -> tuple[Optional[int], str, Optional[bool]]:
@@ -234,7 +242,6 @@ def get_branch_params(condition: Any) -> tuple[Optional[int], str, Optional[bool
                 span=condition.span,
             )
         if isinstance(condition, UnaryExpression):
-            print(condition)
             if condition.op != UnaryOperator["!"]:
                 raise_qasm3_error(
                     message="Only '!' supported in branching condition with classical register",
@@ -287,7 +294,7 @@ def get_branch_params(condition: Any) -> tuple[Optional[int], str, Optional[bool
     @classmethod
     def transform_function_qubits(
         cls,
-        q_op: Union[QuantumGate, QuantumBarrier, QuantumReset],
+        q_op: Union[QuantumGate, QuantumBarrier, QuantumReset, QuantumPhase],
         formal_qreg_sizes: dict[str, int],
         qubit_map: dict[tuple, tuple],
     ) -> list[IndexedIdentifier]:
diff --git a/src/pyqasm/visitor.py b/src/pyqasm/visitor.py
index 090cbdc..7267c77 100644
--- a/src/pyqasm/visitor.py
+++ b/src/pyqasm/visitor.py
@@ -58,6 +58,7 @@ def __init__(self, module, check_only: bool = False, external_gates: list[str] |
         self._module = module
         self._scope: deque = deque([{}])
         self._context: deque = deque([Context.GLOBAL])
+        self._included_files: set[str] = set()
         self._qubit_labels: dict[str, int] = {}
         self._clbit_labels: dict[str, int] = {}
         self._alias_qubit_labels: dict[tuple[str, int], tuple[str, int]] = {}
@@ -345,6 +346,15 @@ def _get_op_bits(
         if isinstance(operation, qasm3_ast.QuantumMeasurementStatement):
             assert operation.target is not None
             bit_list = [operation.measure.qubit] if qubits else [operation.target]
+        elif isinstance(operation, qasm3_ast.QuantumPhase) and operation.qubits is None:
+            for reg_name, reg_size in reg_size_map.items():
+                bit_list.append(
+                    qasm3_ast.IndexedIdentifier(
+                        qasm3_ast.Identifier(reg_name), [[qasm3_ast.IntegerLiteral(i)]]
+                    )
+                    for i in range(reg_size)
+                )
+            return bit_list
         else:
             bit_list = (
                 operation.qubits if isinstance(operation.qubits, list) else [operation.qubits]
@@ -646,7 +656,7 @@ def _broadcast_gate_operation(
             (All arguments of the callable should be qubits, i.e. all non-qubit arguments of the
             gate should already be evaluated, e.g. using functools.partial).
             all_targets (list[list[qasm3_ast.IndexedIdentifier]]):
-                The list of target of target qubits.
+                The list of target qubits.
                 The length of this list indicates the number of time the gate is invoked.
         Returns:
             List of all executed gates.
@@ -735,7 +745,7 @@ def _visit_basic_gate_operation(  # pylint: disable=too-many-locals
 
     def _visit_custom_gate_operation(
         self, operation: qasm3_ast.QuantumGate, inverse: bool = False
-    ) -> list[qasm3_ast.QuantumGate]:
+    ) -> list[Union[qasm3_ast.QuantumGate, qasm3_ast.QuantumPhase]]:
         """Visit a custom gate operation element recursively.
 
         Args:
@@ -782,11 +792,11 @@ def _visit_custom_gate_operation(
         self._push_context(Context.GATE)
         result = []
         for gate_op in gate_definition_ops:
-            if isinstance(gate_op, qasm3_ast.QuantumGate):
+            if isinstance(gate_op, (qasm3_ast.QuantumGate, qasm3_ast.QuantumPhase)):
+                gate_op_copy = copy.deepcopy(gate_op)
                 # necessary to avoid modifying the original gate definition
                 # in case the gate is reapplied
-                gate_op_copy = copy.deepcopy(gate_op)
-                if gate_op.name.name == gate_name:
+                if isinstance(gate_op, qasm3_ast.QuantumGate) and gate_op.name.name == gate_name:
                     raise_qasm3_error(
                         f"Recursive definitions not allowed for gate {gate_name}", span=gate_op.span
                     )
@@ -873,7 +883,45 @@ def gate_function(*qubits):
 
         return result
 
-    def _collapse_gate_modifiers(self, operation: qasm3_ast.QuantumGate) -> tuple:
+    def _visit_phase_operation(
+        self, operation: qasm3_ast.QuantumPhase, inverse: bool = False
+    ) -> list[qasm3_ast.QuantumPhase]:
+        """Visit a phase operation element.
+
+        Args:
+            operation (qasm3_ast.QuantumPhase): The phase operation to visit.
+            inverse (bool): Whether the operation is an inverse operation. Defaults to False.
+
+        Returns:
+            list[qasm3_ast.Statement]: The unrolled quantum phase operation.
+        """
+        logger.debug("Visiting phase operation '%s'", str(operation))
+
+        evaluated_arg = Qasm3ExprEvaluator.evaluate_expression(operation.argument)[0]
+        if inverse:
+            evaluated_arg = -1 * evaluated_arg
+        # remove the modifiers, as we have already applied the inverse
+        operation.modifiers = []
+
+        operation.argument = qasm3_ast.FloatLiteral(value=evaluated_arg)
+        # no qubit evaluation to be done here
+        # if args are provided in global scope, then we should raise error
+        if self._in_global_scope() and len(operation.qubits) != 0:
+            raise_qasm3_error(
+                f"Qubit arguments not allowed for phase operation {str(operation)} in global scope",
+                span=operation.span,
+            )
+
+        # if it were in function scope, then the args would have been evaluated and added to the
+        # qubit list
+        if self._check_only:
+            return []
+
+        return [operation]
+
+    def _collapse_gate_modifiers(
+        self, operation: Union[qasm3_ast.QuantumGate, qasm3_ast.QuantumPhase]
+    ) -> tuple:
         """Collapse the gate modifiers of a gate operation.
            Some analysis is required to get this result.
            The basic idea is that any power operation is multiplied and inversions are toggled.
@@ -908,8 +956,8 @@ def _collapse_gate_modifiers(self, operation: qasm3_ast.QuantumGate) -> tuple:
         return (power_value, inverse_value)
 
     def _visit_generic_gate_operation(
-        self, operation: qasm3_ast.QuantumGate
-    ) -> list[qasm3_ast.QuantumGate]:
+        self, operation: Union[qasm3_ast.QuantumGate, qasm3_ast.QuantumPhase]
+    ) -> list[Union[qasm3_ast.QuantumGate, qasm3_ast.QuantumPhase]]:
         """Visit a gate operation element.
 
         Args:
@@ -922,7 +970,11 @@ def _visit_generic_gate_operation(
         operation = copy.deepcopy(operation)
 
         # only needs to be done once for a gate operation
-        if not self._in_gate_scope() and len(self._function_qreg_size_map) > 0:
+        if (
+            len(operation.qubits) > 0
+            and not self._in_gate_scope()
+            and len(self._function_qreg_size_map) > 0
+        ):
             # we are in SOME function scope
             # transform qubits to use the global qreg identifiers
             operation.qubits = (
@@ -934,9 +986,11 @@ def _visit_generic_gate_operation(
             )
         # Applying the inverse first and then the power is same as
         # apply the power first and then inverting the result
-        result = []
+        result: list[Union[qasm3_ast.QuantumGate, qasm3_ast.QuantumPhase]] = []
         for _ in range(power_value):
-            if operation.name.name in self._external_gates:
+            if isinstance(operation, qasm3_ast.QuantumPhase):
+                result.extend(self._visit_phase_operation(operation, inverse_value))
+            elif operation.name.name in self._external_gates:
                 result.extend(self._visit_external_gate_operation(operation, inverse_value))
             elif operation.name.name in self._custom_gates:
                 result.extend(self._visit_custom_gate_operation(operation, inverse_value))
@@ -1695,6 +1749,24 @@ def _evaluate_case(statements):
             default_stmts = statement.default.statements
             return _evaluate_case(default_stmts)
 
+    def _visit_include(self, include: qasm3_ast.Include) -> list[qasm3_ast.Statement]:
+        """Visit an include statement element.
+
+        Args:
+            include (qasm3_ast.Include): The include statement to visit.
+
+        Returns:
+            None
+        """
+        filename = include.filename
+        if filename in self._included_files:
+            raise_qasm3_error(f"File '{filename}' already included", span=include.span)
+        self._included_files.add(filename)
+        if self._check_only:
+            return []
+
+        return [include]
+
     def visit_statement(self, statement: qasm3_ast.Statement) -> list[qasm3_ast.Statement]:
         """Visit a statement element.
 
@@ -1707,13 +1779,14 @@ def visit_statement(self, statement: qasm3_ast.Statement) -> list[qasm3_ast.Stat
         logger.debug("Visiting statement '%s'", str(statement))
         result = []
         visit_map = {
-            qasm3_ast.Include: lambda x: [],  # No operation
+            qasm3_ast.Include: self._visit_include,  # No operation
             qasm3_ast.QuantumMeasurementStatement: self._visit_measurement,
             qasm3_ast.QuantumReset: self._visit_reset,
             qasm3_ast.QuantumBarrier: self._visit_barrier,
             qasm3_ast.QubitDeclaration: self._visit_quantum_register,
             qasm3_ast.QuantumGateDefinition: self._visit_gate_definition,
             qasm3_ast.QuantumGate: self._visit_generic_gate_operation,
+            qasm3_ast.QuantumPhase: self._visit_generic_gate_operation,
             qasm3_ast.ClassicalDeclaration: self._visit_classical_declaration,
             qasm3_ast.ClassicalAssignment: self._visit_classical_assignment,
             qasm3_ast.ConstantDeclaration: self._visit_constant_declaration,
@@ -1754,3 +1827,24 @@ def visit_basic_block(self, stmt_list: list[qasm3_ast.Statement]) -> list[qasm3_
         for stmt in stmt_list:
             result.extend(self.visit_statement(stmt))
         return result
+
+    def finalize(self, unrolled_stmts):
+        """Finalize the unrolled statements.
+        Rules:
+        - Remove qubit args from phase operations if ALL qubits are used
+        To add more rules if needed
+
+        Args:
+            unrolled_stmts (list[qasm3_ast.Statement]): The list of unrolled statements.
+
+        Returns:
+            list[qasm3_ast.Statement]: The list of finalized statements.
+
+        """
+        # remove the gphase qubits if they use ALL qubits
+        for stmt in unrolled_stmts:
+            # Rule 1
+            if isinstance(stmt, qasm3_ast.QuantumPhase):
+                if len(stmt.qubits) == len(self._qubit_labels):
+                    stmt.qubits = []
+        return unrolled_stmts


Test patch:
diff --git a/tests/qasm2/test_declarations.py b/tests/qasm2/test_declarations.py
index 46ea5b6..2fe48cd 100644
--- a/tests/qasm2/test_declarations.py
+++ b/tests/qasm2/test_declarations.py
@@ -30,8 +30,8 @@ def test_qubit_declarations():
     
     """
 
-    expected_qasm = """OPENQASM 2.0;
-    include 'qelib1.inc';
+    expected_qasm = """
+    OPENQASM 2.0;
     qreg q3[3];
     qreg q[1];
     """
@@ -83,7 +83,8 @@ def test_qubit_clbit_declarations():
     creg c2[2];
     """
 
-    expected_qasm = """OPENQASM 2.0;
+    expected_qasm = """
+    OPENQASM 2.0;
     include 'qelib1.inc';
     qreg q1[1];
     qreg q2[2];
diff --git a/tests/qasm3/resources/gates.py b/tests/qasm3/resources/gates.py
index 21c8108..1d6bbed 100644
--- a/tests/qasm3/resources/gates.py
+++ b/tests/qasm3/resources/gates.py
@@ -19,6 +19,7 @@
 import pytest
 
 from pyqasm.maps import (
+    FOUR_QUBIT_OP_MAP,
     ONE_QUBIT_OP_MAP,
     ONE_QUBIT_ROTATION_MAP,
     THREE_QUBIT_OP_MAP,
@@ -30,6 +31,14 @@
 
 RESOURCES_DIR = os.path.join(os.path.dirname(__file__), "qasm")
 
+VALID_GATE_NAMES = set(
+    ONE_QUBIT_OP_MAP.keys()
+    | TWO_QUBIT_OP_MAP.keys()
+    | ONE_QUBIT_ROTATION_MAP.keys()
+    | THREE_QUBIT_OP_MAP.keys()
+    | FOUR_QUBIT_OP_MAP.keys()
+)
+
 
 def resources_file(filename: str) -> str:
     return os.path.join(RESOURCES_DIR, f"{filename}")
@@ -39,22 +48,10 @@ def _fixture_name(s: str) -> str:
     return f"Fixture_{s}"
 
 
-def _validate_gate_name(gate_name: str) -> str:
-    if gate_name in ONE_QUBIT_OP_MAP:
-        return True
-    if gate_name in TWO_QUBIT_OP_MAP:
-        return True
-    if gate_name in ONE_QUBIT_ROTATION_MAP:
-        return True
-    if gate_name in THREE_QUBIT_OP_MAP:
-        return True
-    return False
-
-
 def _generate_one_qubit_fixture(gate_name: str):
     @pytest.fixture()
     def test_fixture():
-        if not _validate_gate_name(gate_name):
+        if gate_name not in VALID_GATE_NAMES:
             raise ValueError(f"Unknown qasm3 gate {gate_name}")
         qasm3_string = f"""
         OPENQASM 3;
@@ -79,7 +76,7 @@ def test_fixture():
 def _generate_rotation_fixture(gate_name: str):
     @pytest.fixture()
     def test_fixture():
-        if not _validate_gate_name(gate_name):
+        if gate_name not in VALID_GATE_NAMES:
             raise ValueError(f"Unknown qasm3 gate {gate_name}")
         qasm3_string = f"""
         OPENQASM 3;
@@ -103,7 +100,7 @@ def test_fixture():
 def _generate_two_qubit_fixture(gate_name: str):
     @pytest.fixture()
     def test_fixture():
-        if not _validate_gate_name(gate_name):
+        if gate_name not in VALID_GATE_NAMES:
             raise ValueError(f"Unknown qasm3 gate {gate_name}")
 
         params = ""
@@ -135,7 +132,7 @@ def test_fixture():
 def _generate_three_qubit_fixture(gate_name: str):
     @pytest.fixture()
     def test_fixture():
-        if not _validate_gate_name(gate_name):
+        if gate_name not in VALID_GATE_NAMES:
             raise ValueError(f"Unknown qasm3 gate {gate_name}")
         qasm3_string = f"""
         OPENQASM 3;
@@ -156,6 +153,30 @@ def test_fixture():
     locals()[name] = _generate_three_qubit_fixture(gate)
 
 
+def _generate_four_qubit_fixture(gate_name: str):
+    @pytest.fixture()
+    def test_fixture():
+        if gate_name not in VALID_GATE_NAMES:
+            raise ValueError(f"Unknown qasm3 gate {gate_name}")
+        qasm3_string = f"""
+        OPENQASM 3;
+        include "stdgates.inc";
+
+        qubit[4] q;
+        {gate_name} q[0], q[1], q[2], q[3];
+        {gate_name} q;
+        """
+        return qasm3_string
+
+    return test_fixture
+
+
+# Generate four-qubit gate fixtures
+for gate in FOUR_QUBIT_OP_MAP:
+    name = _fixture_name(gate)
+    locals()[name] = _generate_four_qubit_fixture(gate)
+
+
 def _generate_custom_op_fixture(op_name: str):
     print(os.getcwd())
 
@@ -175,7 +196,7 @@ def test_fixture():
     locals()[name] = _generate_custom_op_fixture(test_name)
 
 single_op_tests = [_fixture_name(s) for s in ONE_QUBIT_OP_MAP]
-already_tested_single_op = ["id", "si", "ti", "v", "sx", "vi", "sxdg"]
+already_tested_single_op = ["id", "si", "ti", "v", "sx", "vi", "sxdg", "not"]
 for gate in already_tested_single_op:
     single_op_tests.remove(_fixture_name(gate))
 
@@ -212,10 +233,12 @@ def test_fixture():
     double_op_tests.remove(_fixture_name(gate))
 
 triple_op_tests = [_fixture_name(s) for s in THREE_QUBIT_OP_MAP]
-already_tested_triple_op = ["ccnot", "cswap", "rccx"]
+already_tested_triple_op = ["ccnot", "cswap", "rccx", "toffoli"]
 for gate in already_tested_triple_op:
     triple_op_tests.remove(_fixture_name(gate))
 
+four_op_tests = [_fixture_name(s) for s in FOUR_QUBIT_OP_MAP]
+
 custom_op_tests = [_fixture_name(s) for s in CUSTOM_OPS]
 
 # qasm_input, expected_error
@@ -285,6 +308,14 @@ def test_fixture():
 
 # qasm_input, expected_error
 CUSTOM_GATE_INCORRECT_TESTS = {
+    "incorrect_gphase_usage": (
+        """
+        OPENQASM 3.0;
+        qubit q;
+        gphase(pi) q;
+        """,
+        r"Qubit arguments not allowed for phase operation",
+    ),
     "undeclared_custom": (
         """
         OPENQASM 3;
diff --git a/tests/qasm3/test_entrypoint.py b/tests/qasm3/test_entrypoint.py
index c8050f4..327c91c 100644
--- a/tests/qasm3/test_entrypoint.py
+++ b/tests/qasm3/test_entrypoint.py
@@ -63,3 +63,9 @@ def test_incorrect_module_dump():
 def test_incorrect_qasm_version():
     with pytest.raises(ValidationError, match=r"Unsupported OpenQASM version"):
         loads("OPENQASM 4.0;\n qubit q;")
+
+
+def test_incorrect_module_unroll_raises_error():
+    with pytest.raises(ValidationError):
+        module = loads("OPENQASM 3.0;\n qubit q; h q[2]")
+        module.unroll()
diff --git a/tests/qasm3/test_format.py b/tests/qasm3/test_format.py
index 67e7807..dab09eb 100644
--- a/tests/qasm3/test_format.py
+++ b/tests/qasm3/test_format.py
@@ -203,31 +203,3 @@ def test_simplify_redundant_parentheses(qasm_input, expected_result):
     module.unroll(external_gates=["cry", "crx", "ry"])
     actual_result = dumps(module)
     check_unrolled_qasm(actual_result, expected_result)
-
-
-@pytest.mark.parametrize(
-    "qasm3_without, qasm3_with",
-    [
-        (
-            """
-        OPENQASM 3.0;
-        qubit[1] q;
-        h q[0];
-        ry(pi/4) q[0];
-                """,
-            """
-        OPENQASM 3.0;
-        include "stdgates.inc";
-        qubit[1] q;
-        h q[0];
-        ry(0.7853981633974483) q[0];
-        """,
-        ),
-    ],
-)
-def test_add_stdgates_include(qasm3_without, qasm3_with):
-    """Test adding stdgates include to OpenQASM 3.0 string"""
-    module = loads(qasm3_without)
-    module.unroll()
-    actual_qasm3 = dumps(module)
-    check_unrolled_qasm(actual_qasm3, qasm3_with)
diff --git a/tests/qasm3/test_gates.py b/tests/qasm3/test_gates.py
index 507a18d..048b1c1 100644
--- a/tests/qasm3/test_gates.py
+++ b/tests/qasm3/test_gates.py
@@ -14,13 +14,14 @@
 """
 import pytest
 
-from pyqasm.entrypoint import loads
+from pyqasm.entrypoint import dumps, loads
 from pyqasm.exceptions import ValidationError
 from tests.qasm3.resources.gates import (
     CUSTOM_GATE_INCORRECT_TESTS,
     SINGLE_QUBIT_GATE_INCORRECT_TESTS,
     custom_op_tests,
     double_op_tests,
+    four_op_tests,
     rotation_tests,
     single_op_tests,
     triple_op_tests,
@@ -28,10 +29,12 @@
 from tests.utils import (
     check_custom_qasm_gate_op,
     check_custom_qasm_gate_op_with_external_gates,
+    check_four_qubit_gate_op,
     check_single_qubit_gate_op,
     check_single_qubit_rotation_op,
     check_three_qubit_gate_op,
     check_two_qubit_gate_op,
+    check_unrolled_qasm,
 )
 
 
@@ -90,6 +93,20 @@ def test_three_qubit_qasm3_gates(circuit_name, request):
     check_three_qubit_gate_op(result.unrolled_ast, 2, qubit_list, gate_name)
 
 
+@pytest.mark.parametrize("circuit_name", four_op_tests)
+def test_four_qubit_qasm3_gates(circuit_name, request):
+    qubit_list = [[0, 1, 2, 3], [0, 1, 2, 3]]
+    gate_name = circuit_name.removeprefix("Fixture_")
+
+    qasm3_string = request.getfixturevalue(circuit_name)
+    result = loads(qasm3_string)
+    # we do not want to validate every gate inside it
+    result.unroll(external_gates=[gate_name])
+    assert result.num_qubits == 4
+    assert result.num_clbits == 0
+    check_four_qubit_gate_op(result.unrolled_ast, 2, qubit_list, gate_name)
+
+
 def test_gate_body_param_expression():
     qasm3_string = """
     OPENQASM 3;
@@ -204,6 +221,100 @@ def test_custom_ops(test_name, request):
     check_custom_qasm_gate_op(result.unrolled_ast, gate_type)
 
 
+def test_global_phase_gate():
+    qasm3_string = """OPENQASM 3.0;
+    qubit[2] q;
+    gphase(pi/4);
+    """
+
+    qasm3_expected = """
+    OPENQASM 3.0;
+    qubit[2] q;
+    gphase(0.7853981633974483);
+    """
+    module = loads(qasm3_string)
+    module.unroll()
+
+    assert module.num_qubits == 2
+    assert module.num_clbits == 0
+
+    check_unrolled_qasm(dumps(module), qasm3_expected)
+
+
+def test_global_phase_qubits_retained():
+    """Test that global phase gate is retained when applied on specific qubits"""
+    qasm3_string = """OPENQASM 3.0;
+    gate custom a,b,c { 
+       gphase(pi/8);
+       h a;
+    }
+    qubit[23] q2;
+    custom q2[0:3];
+    """
+
+    qasm3_expected = """
+    OPENQASM 3.0;
+    qubit[23] q2;
+    gphase(0.39269908169872414) q2[0], q2[1], q2[2];
+    h q2[0];
+    """
+    module = loads(qasm3_string)
+    module.unroll()
+
+    assert module.num_qubits == 23
+    assert module.num_clbits == 0
+
+    check_unrolled_qasm(dumps(module), qasm3_expected)
+
+
+def test_global_phase_qubits_simplified():
+    """Test that the global phase gate is simplified when applied on all qubits"""
+    qasm3_string = """OPENQASM 3.0;
+    qubit[3] q2;
+    gate custom a,b,c {
+        gphase(pi/8) a, b, c;
+    }
+    custom q2;
+    """
+
+    qasm3_expected = """
+    OPENQASM 3.0;
+    qubit[3] q2;
+    gphase(0.39269908169872414);
+    """
+    module = loads(qasm3_string)
+    module.unroll()
+
+    assert module.num_qubits == 3
+    assert module.num_clbits == 0
+
+    check_unrolled_qasm(dumps(module), qasm3_expected)
+
+
+def test_inverse_global_phase():
+    """Test that the inverse of global phase gate is simplified"""
+    qasm3_string = """OPENQASM 3.0;
+    qubit[3] q2;
+    gate custom a,b,c {
+        inv @ gphase(pi/8) a, b, c;
+    }
+    custom q2;
+    """
+
+    qasm3_expected = """
+    OPENQASM 3.0;
+    qubit[3] q2;
+    gphase(-0.39269908169872414);
+    """
+    module = loads(qasm3_string)
+    module.unroll()
+
+    assert module.num_qubits == 3
+    assert module.num_clbits == 0
+
+    check_unrolled_qasm(dumps(module), qasm3_expected)
+
+
 @pytest.mark.parametrize("test_name", custom_op_tests)
 def test_custom_ops_with_external_gates(test_name, request):
     qasm3_string = request.getfixturevalue(test_name)
diff --git a/tests/qasm3/test_measurement.py b/tests/qasm3/test_measurement.py
index 12e048a..2872d5a 100644
--- a/tests/qasm3/test_measurement.py
+++ b/tests/qasm3/test_measurement.py
@@ -40,8 +40,8 @@ def test_measure():
 
     """
 
-    expected_qasm = """OPENQASM 3.0;
-    include "stdgates.inc";
+    expected_qasm = """
+    OPENQASM 3.0;
     qubit[2] q1;
     qubit[5] q2;
     qubit[1] q3;
@@ -112,8 +112,8 @@ def test_remove_measurement():
 
     """
 
-    expected_qasm = """OPENQASM 3.0;
-    include 'stdgates.inc';
+    expected_qasm = """
+    OPENQASM 3.0;
     qubit[2] q1;
     qubit[5] q2;
     qubit[1] q3;
diff --git a/tests/test_include.py b/tests/test_include.py
new file mode 100644
index 0000000..aab0779
--- /dev/null
+++ b/tests/test_include.py
@@ -0,0 +1,65 @@
+# Copyright (C) 2024 qBraid
+#
+# This file is part of pyqasm
+#
+# Pyqasm is free software released under the GNU General Public License v3
+# or later. You can redistribute and/or modify it under the terms of the GPL v3.
+# See the LICENSE file in the project root or <https://www.gnu.org/licenses/gpl-3.0.html>.
+#
+# THERE IS NO WARRANTY for pyqasm, as per Section 15 of the GPL v3.
+
+"""
+Module containing unit tests for linalg.py functions.
+
+"""
+import pytest
+
+from pyqasm import ValidationError, dumps, loads
+from tests.utils import check_unrolled_qasm
+
+
+def test_no_include_added():
+    qasm_str = """
+    OPENQASM 3.0;
+    include "random.qasm";
+    """
+    expected_qasm_str = """
+    OPENQASM 3.0;
+    include "random.qasm";
+    """
+    module = loads(qasm_str)
+    module.unroll()
+    check_unrolled_qasm(dumps(module), expected_qasm_str)
+
+
+def test_includes_preserved():
+    qasm_str = """
+    OPENQASM 3.0;
+    include "stdgates.inc";
+    include "random.qasm";
+
+    qubit[2] q;
+    h q;
+    """
+    expected_qasm_str = """
+    OPENQASM 3.0;
+    include "stdgates.inc";
+    include "random.qasm";
+    qubit[2] q;
+    h q[0];
+    h q[1];
+    """
+    module = loads(qasm_str)
+    module.unroll()
+    check_unrolled_qasm(dumps(module), expected_qasm_str)
+
+
+def test_repeated_include_raises_error():
+    qasm_str = """
+    OPENQASM 3.0;
+    include "stdgates.inc";
+    include "stdgates.inc";
+    """
+    with pytest.raises(ValidationError):
+        module = loads(qasm_str)
+        module.validate()
diff --git a/tests/utils.py b/tests/utils.py
index e5a2b6f..94b396c 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -54,6 +54,25 @@ def check_single_qubit_gate_op(unrolled_ast, num_gates, qubit_list, gate_name):
     assert gate_count == num_gates
 
 
+def check_global_phase_op(unrolled_ast, num_gates, qubit_list, phase_list):
+    # qubit_list is list [list]
+    # qubit_list = [ [], [0,1,2], [3,4,5] ] means gphase applied on ALL , 3 , 3 qubits  with given
+    # phase values
+    qubits_id, phase_id = 0, 0
+    gate_count = 0
+    for stmt in unrolled_ast.statements:
+        if isinstance(stmt, qasm3_ast.QuantumPhase):
+            assert len(stmt.qubits) == len(qubit_list[qubits_id])
+            assert stmt.argument.value == phase_list[phase_id]
+            for i, _ in enumerate(stmt.qubits):
+                assert stmt.qubits[i].indices[0][0].value == qubit_list[qubits_id][i]
+            qubits_id += 1
+            phase_id += 1
+            gate_count += 1
+
+    assert gate_count == num_gates
+
+
 def _check_ch_gate_op(unrolled_ast, num_gates, qubits):
     check_single_qubit_gate_op(unrolled_ast, num_gates, [qubits[1]] * num_gates, "s")
     check_single_qubit_gate_op(unrolled_ast, 2 * num_gates, [qubits[1]] * 2 * num_gates, "h")
@@ -121,11 +140,15 @@ def _check_cry_gate_op(unrolled_ast, num_gates, qubits, theta):
 
 
 def _check_rxx_gate_op(unrolled_ast, num_gates, qubits, theta):
+    num_global_phase = num_gates
     num_h_gates = 4 * num_gates
     # because H is applied as H(q0) H(q1) H(q1) H(q0) in  1 rxx gate
     h_qubit_list = (qubits + qubits[::-1]) * num_gates
     num_cx_gates = 2 * num_gates
     num_rz_gates = num_gates
+    check_global_phase_op(
+        unrolled_ast, num_global_phase, [[]] * num_gates, [-theta / 2] * num_gates
+    )
     check_single_qubit_gate_op(unrolled_ast, num_h_gates, h_qubit_list, "h")
     check_two_qubit_gate_op(unrolled_ast, num_cx_gates, [qubits] * num_cx_gates, "cx")
     check_single_qubit_rotation_op(
@@ -154,9 +177,13 @@ def _check_ryy_gate_op(unrolled_ast, num_gates, qubits, theta):
 
 
 def _check_rzz_gate_op(unrolled_ast, num_gates, qubits, theta):
+    num_global_phase = num_gates
     num_cx_gates = 2 * num_gates
     num_u3_gates = num_gates
 
+    check_global_phase_op(
+        unrolled_ast, num_global_phase, [[]] * num_gates, [-theta / 2] * num_gates
+    )
     check_two_qubit_gate_op(unrolled_ast, num_cx_gates, [qubits] * num_cx_gates, "cx")
     check_u3_gate_op(
         unrolled_ast,
@@ -244,6 +271,20 @@ def check_three_qubit_gate_op(unrolled_ast, num_gates, qubit_list, gate_name):
     assert gate_count == num_gates
 
 
+def check_four_qubit_gate_op(unrolled_ast, num_gates, qubit_list, gate_name):
+    qubit_id, gate_count = 0, 0
+
+    for stmt in unrolled_ast.statements:
+        if isinstance(stmt, qasm3_ast.QuantumGate) and stmt.name.name == gate_name.lower():
+            assert len(stmt.qubits) == 4
+            for i in range(4):
+                assert stmt.qubits[i].indices[0][0].value == qubit_list[qubit_id][i]
+            qubit_id += 1
+            gate_count += 1
+
+    assert gate_count == num_gates
+
+
 def check_u3_gate_op(unrolled_ast, num_gates, qubit_list, param_list):
     op_count = 0
     q_id = 0