diff --git a/algorithms/quantum_linear_solvers/hhl/hhl.ipynb b/algorithms/quantum_linear_solvers/hhl/hhl.ipynb
index 18c633182..89a1c7b0a 100644
--- a/algorithms/quantum_linear_solvers/hhl/hhl.ipynb
+++ b/algorithms/quantum_linear_solvers/hhl/hhl.ipynb
@@ -589,22 +589,9 @@
     "qmod_hhl_exact = hhl_model(main, backend_preferences)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 10,
-   "id": "31",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from classiq import write_qmod\n",
-    "\n",
-    "# Save qmod file\n",
-    "write_qmod(qmod_hhl_exact, \"hhl\", decimal_precision=20, symbolic_only=False)"
-   ]
-  },
   {
    "cell_type": "markdown",
-   "id": "32",
+   "id": "31",
    "metadata": {},
    "source": [
     "#### Synthesizing the Model (exact)"
@@ -613,7 +600,7 @@
   {
    "cell_type": "code",
    "execution_count": 11,
-   "id": "33",
+   "id": "32",
    "metadata": {
     "colab": {
      "base_uri": "https://localhost:8080/"
@@ -646,7 +633,7 @@
     }
    },
    "cell_type": "markdown",
-   "id": "34",
+   "id": "33",
    "metadata": {},
    "source": [
     "![image.png](attachment:70e1642f-7684-4c35-82e7-81f3fecb9476.png)\n",
@@ -657,7 +644,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "35",
+   "id": "34",
    "metadata": {},
    "source": [
     "#### Statevector Simulation (exact)\n",
@@ -668,7 +655,7 @@
   {
    "cell_type": "code",
    "execution_count": 12,
-   "id": "36",
+   "id": "35",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -679,7 +666,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "37",
+   "id": "36",
    "metadata": {},
    "source": [
     "#### Results (exact)"
@@ -688,7 +675,7 @@
   {
    "cell_type": "code",
    "execution_count": 13,
-   "id": "38",
+   "id": "37",
    "metadata": {
     "colab": {
      "base_uri": "https://localhost:8080/",
@@ -726,7 +713,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "39",
+   "id": "38",
    "metadata": {},
    "source": [
     "### Example: Approximated Hamiltonian Evolution (Suzuki-Trotter)\n",
@@ -736,7 +723,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "40",
+   "id": "39",
    "metadata": {},
    "source": [
     "For the QPE we are going to use Classiq's `suzuki-trotter` function of order one for Hamiltonian evolution $e^{-i H t}$ [[3](#Trotter)]. This function is an approximated one, where its `repetitions` parameter controls its error. For a QPE algorithm with estimator size $m$ a series of controlled-unitaries $U^{2^k}$ with $0 \\leq k \\leq n-1$ are applied, for each of which we would like to pass a different `repetitions` parameter depth (to keep a roughly same error, the repetitions for approximating $U=e^{2\\pi i 2^k A }$ is expected to be $\\sim 2^k$ times the repetitions of $U=e^{2\\pi i A }$). In Classiq this can be done by working with a [`qpe_flexible`](https://github.com/Classiq/classiq-library/blob/main/functions/qmod_library_reference/classiq_open_library/qpe/qpe.ipynb), and passing a \"rule\" for how to exponentiate each step within the QPE in `repetitions` parameter."
@@ -745,7 +732,7 @@
   {
    "cell_type": "code",
    "execution_count": 14,
-   "id": "41",
+   "id": "40",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -771,7 +758,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "42",
+   "id": "41",
    "metadata": {},
    "source": [
     "The parameters `R0` and `REPS_SCALING_FACTOR` dictate the number of repititions in the Suzuki-Trotter approximation. The specific number of repititions depend on the Hamiltonian, and therefore, were chosen by trial and error. For other examples, one would need to use different values for these parameters, please compare with specific example in [Flexible QPE tutorial](https://github.com/Classiq/classiq-library/blob/main/tutorials/advanced_tutorials/high_level_modeling_flexible_qpe/high_level_modeling_flexible_qpe.ipynb). The relevant literature that discusses the errors of product formulas is available in Ref. [[5](#trotter-error)]."
@@ -780,7 +767,7 @@
   {
    "cell_type": "code",
    "execution_count": 15,
-   "id": "43",
+   "id": "42",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -826,7 +813,7 @@
   {
    "cell_type": "code",
    "execution_count": 17,
-   "id": "44",
+   "id": "43",
    "metadata": {
     "colab": {
      "base_uri": "https://localhost:8080/",
@@ -878,7 +865,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "45",
+   "id": "44",
    "metadata": {},
    "source": [
     "We explored the HHL algorithm for solving linear systems using exact and approximated Hamiltonian simulations. The exact method, with a smaller circuit depth, is computationally less intensive but lacks scalability. In contrast, the approximated method, with a greater circuit depth, offers flexibility and can handle larger, more complex systems. This trade-off underscores the importance of selecting the appropriate method based on the problem's size and complexity."
@@ -886,7 +873,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "46",
+   "id": "45",
    "metadata": {},
    "source": [
     "## Technical Notes\n",
@@ -899,7 +886,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "47",
+   "id": "46",
    "metadata": {},
    "source": [
     "### Generalizations\n",
@@ -962,7 +949,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "48",
+   "id": "47",
    "metadata": {},
    "source": [
     "## References\n",
@@ -1007,7 +994,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.7"
+   "version": "3.11.9"
   }
  },
  "nbformat": 4,
diff --git a/algorithms/quantum_linear_solvers/hhl/hhl.qmod b/algorithms/quantum_linear_solvers/hhl/hhl.qmod
index 543749bb2..2d3e6f435 100644
--- a/algorithms/quantum_linear_solvers/hhl/hhl.qmod
+++ b/algorithms/quantum_linear_solvers/hhl/hhl.qmod
@@ -1,164 +1,113 @@
-qfunc load_b_expanded___0(amplitudes: real[], output memory: qbit[2]) {
-  prepare_amplitudes(amplitudes, 0.0, memory);
-}
+// Harrow–Hassidim–Lloyd (HHL) Algorithm for Solving Linear Systems
 
-qfunc apply_to_all_expanded___0(target: qbit[4]) {
-  repeat (index: 4) {
-    H(target[index]);
-  }
-}
+// This example implements the HHL quantum linear-systems algorithm to solve
+//   A |x⟩ = |b⟩
+// where A is a Hermitian matrix and |b⟩ is the right-hand-side vector prepared
+// as a quantum state.
 
-qfunc hamiltonian_evolution_with_power_0_lambda___0_0_expanded___0(pw: int, target: qbit[2]) {
-  power (pw) {
-    unitary([
-      [
-        ((-0.09406240950199857) + 0.8149069223122054j),
-        (0.03521871946675126 - 0.029763534641642615j),
-        ((-0.018800717000078293) - 0.16142879795007106j),
-        (0.43769245930764733 + 0.32705554908759304j)
-      ],
-      [
-        (0.03521871946675127 - 0.029763534641642633j),
-        ((-0.15347248298890326) - 0.17275282472948233j),
-        (0.23117644455908531 + 0.8872069971297389j),
-        (0.23971825754883572 + 0.21548267921288933j)
-      ],
-      [
-        ((-0.01880071700007826) - 0.16142879795007103j),
-        (0.23117644455908523 + 0.8872069971297386j),
-        ((-0.1219131720516462) + 0.13200138126428373j),
-        (0.29584069101495575 + 0.11488938733473114j)
-      ],
-      [
-        (0.43769245930764744 + 0.32705554908759327j),
-        (0.2397182575488357 + 0.21548267921288933j),
-        (0.29584069101495586 + 0.1148893873347311j),
-        ((-0.6563827949579104) + 0.25690988991104674j)
-      ]
-    ], target);
-  }
-}
 
-qfunc unitary_with_power_0_lambda___0_0_expanded___0(k: int, memory_captured__hhl__1: qbit[2]) {
-  hamiltonian_evolution_with_power_0_lambda___0_0_expanded___0(k, memory_captured__hhl__1);
-}
+// In this example, the matrix A is given as a sum of Pauli strings:
+//   A = 0.03*X(0)*X(1) + 0.05*Z(0)*Z(1) + 0.02*Y(0)*Y(1) + 0.08*I
+// This Pauli decomposition corresponds to the matrix
+//   A = [[0.13, 0.00, 0.00, 0.01],
+//        [0.00, 0.03, 0.05, 0.00],
+//        [0.00, 0.05, 0.03, 0.00],
+//        [0.01, 0.00, 0.00, 0.13]]
 
-qfunc qft_no_swap_expanded___0(qbv: qbit[4]) {
-  repeat (i: 4) {
-    H(qbv[i]);
-    repeat (j: (4 - i) - 1) {
-      CPHASE(pi / (2 ** (j + 1)), qbv[(i + j) + 1], qbv[i]);
-    }
-  }
-}
+hamiltonian: SparsePauliOp = SparsePauliOp {
+    terms=[
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=1, index=0},
+          IndexedPauli {pauli=1, index=1}
+        ],
+        coefficient=0.03
+      },
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=3, index=1},
+          IndexedPauli {pauli=3, index=0}
+        ],
+        coefficient=0.05
+      },
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=2, index=0},
+          IndexedPauli {pauli=2, index=1}
+        ],
+        coefficient=0.02
+      },
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=0, index=0}
+        ],
+        coefficient=0.08
+      }
+    ],
+    num_qubits=2
+    };
 
-qfunc qft_expanded___0(target: qbit[4]) {
-  repeat (index: 2.0) {
-    SWAP(target[index], target[3 - index]);
-  }
-  qft_no_swap_expanded___0(target);
-}
+// For the right-hand-side vector |b⟩, we use a normalized vector of size 4.
+rhs_vector: real[] = [sqrt(1/6), sqrt(1/6), sqrt(1/6), sqrt(1/2)];
 
-qfunc qpe_flexible_expanded___0(phase: qbit[4], memory_captured__hhl__1: qbit[2]) {
-  apply_to_all_expanded___0(phase);
-  repeat (index: 4) {
-    control (phase[index]) {
-      unitary_with_power_0_lambda___0_0_expanded___0(2 ** index, memory_captured__hhl__1);
-    }
-  }
-  invert {
-    qft_expanded___0(phase);
-  }
+// In our example the matrix A is a sum of commuting Pauli strings.
+// We define how to take powers of the Hamiltonian evolution exp(2π A), which is applied
+// within the Quantum Phase Estimation. We simply multiply the evolution time by the power value,
+// keeping the number of repetitions at 1, since the terms commute.
+qfunc powered_suzuki_trotter_commuting_terms(k: int, hamiltonian: SparsePauliOp, qba: qbit[]) {
+  suzuki_trotter(hamiltonian, -2*pi* k, 1, 1, qba);
 }
 
-qfunc assign_amplitude_table_expanded___0(const index: qbit[4], indicator: qbit) {
-  RY(0.49069646327199606, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.159660988575539, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.215103068840025, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.114786680603947, indicator);
-  skip_control {
-    CX(index[2], indicator);
-  }
-  RY(0.0734201414091853, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.222322935969005, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.181317761599329, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.22482930086683403, indicator);
-  skip_control {
-    CX(index[3], indicator);
-  }
-  RY(0.192520246080629, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.184296564729869, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.22312206827512002, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.0676093870745949, indicator);
-  skip_control {
-    CX(index[2], indicator);
-  }
-  RY(0.0978641117835104, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.216731023385388, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.168241915420623, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.309069995704199, indicator);
-  skip_control {
-    CX(index[3], indicator);
-  }
+// We define a quantum function that assigns eigenvalue inversion to the amplitudes,
+// conditioned on an additional indicator qubit being in the |1⟩ state:
+// |value⟩|0⟩ -> (C / value)|value⟩|1⟩ + sqrt(1 − (C / value)^2)|value⟩|0⟩.
+// The variable `value` is a signed qnum of size 4 with 4 fractional bits,
+// and the inversion amplitudes are normalized using the smallest possible value,
+// C = 2^(−value.size).
+qfunc assign_inversion_size4(const value: qbit[4], indicator: qbit) {
+  RY(-0.0157, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(-0.0157, indicator);
+  skip_control { CX(value[1], indicator); }
+  RY(-0.3379, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(0.3066, indicator);
+  skip_control { CX(value[2], indicator); }
+  RY(-0.1047, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(-0.1047, indicator);
+  skip_control { CX(value[1], indicator); }
+  RY(-0.3181, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(0.4649, indicator);
+  skip_control { CX(value[3], indicator); }
+  RY(-0.0475, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(-0.0475, indicator);
+  skip_control { CX(value[1], indicator); }
+  RY(-0.3407, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(0.2457, indicator);
+  skip_control { CX(value[2], indicator); }
+  RY(-0.0939, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(-0.0939, indicator);
+  skip_control { CX(value[1], indicator); }
+  RY(-0.3122, indicator);
+  skip_control { CX(value[0], indicator); }
+  RY(0.8154, indicator);
+  skip_control { CX(value[3], indicator); }
 }
 
-qfunc hhl_expanded___0(rhs_vector: real[], output memory: qbit[2], output estimator: qnum<4, False, 4>, output indicator: qbit) {
-  allocate(4, False, 4, estimator);
-  load_b_expanded___0([
-    0.18257418583505536,
-    0.3651483716701107,
-    0.7302967433402214,
-    0.5477225575051661
-  ], memory);
-  allocate(1, indicator);
+qfunc main(output solution: qnum<2, False, 0>, output phase_var: qnum<4, True, 4>, output indicator: qbit) {
+  allocate(phase_var);
+  allocate(indicator);
+  prepare_amplitudes(rhs_vector, 0.0, solution);
   within {
-    qpe_flexible_expanded___0(estimator, memory);
+    qpe_flexible(lambda(k) {
+      powered_suzuki_trotter_commuting_terms(k, hamiltonian, solution);
+    }, phase_var);
   } apply {
-    assign_amplitude_table_expanded___0(estimator, indicator);
+    assign_inversion_size4(phase_var, indicator);
   }
 }
-
-qfunc main(output res: qnum<2, False, 0>, output estimator_var: qnum<4, False, 4>, output indicator: qbit) {
-  hhl_expanded___0([
-    0.18257418583505536,
-    0.3651483716701107,
-    0.7302967433402214,
-    0.5477225575051661
-  ], res, estimator_var, indicator);
-}