diff --git a/python/README.md b/python/README.md
index b986585..af8a8de 100644
--- a/python/README.md
+++ b/python/README.md
@@ -74,6 +74,7 @@ print("Status:", m.Status)
 print("Objective:", m.ObjVal)
 print("Primal solution:", m.X)
 print("Dual solution:", m.Pi)
+print("Reduced cost:", m.RC)
 ```
 
 ## Modeling
@@ -190,6 +191,7 @@ After calling `m.optimize()`, the solver stores results in a set of read-only at
 | `RelGap` | float | Relative primal-dual gap. |
 | `X` | numpy.ndarray | Primal solution vector \(x\). May be `None` if no feasible solution was found. |
 | `Pi` | numpy.ndarray | Dual solution vector (Lagrange multipliers). |
+| `RC` | numpy.ndarray | Reduced Cost vector. |
 | `IterCount` | int | Number of iterations performed. |
 | `Runtime` | float | Total wall-clock runtime in seconds. |
 | `RescalingTime` | float | Time spent on preprocessing and rescaling (seconds). |
@@ -214,6 +216,7 @@ print("Iterations:", m.IterCount, " Runtime (s):", m.Runtime)
 # Access solutions
 print("Primal solution:", m.X)
 print("Dual solution:", m.Pi)
+print("Reduced cost:", m.RC)
 
 # Check residuals
 print("Primal residual:", m.RelPrimalResidual)
diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py
index a7630af..5302353 100644
--- a/python/cupdlpx/model.py
+++ b/python/cupdlpx/model.py
@@ -128,6 +128,7 @@ def __init__(
         # initialize solution attributes
         self._x: Optional[np.ndarray] = None # primal solution
         self._y: Optional[np.ndarray] = None # dual solution
+        self._rc: Optional[np.ndarray] = None # reduced costs
         self._objval: Optional[float] = None # objective value
         self._dualobj: Optional[float] = None # dual objective value
         self._gap: Optional[float] = None # primal-dual gap
@@ -376,6 +377,7 @@ def optimize(self):
         # solutions
         self._x = np.asarray(info.get("X")) if info.get("X") is not None else None
         self._y = np.asarray(info.get("Pi")) if info.get("Pi") is not None else None
+        self._rc = np.asarray(info.get("RC")) if info.get("RC") is not None else None
         # objectives & gaps
         primal_obj_eff = info.get("PrimalObj")
         dual_obj_eff   = info.get("DualObj")
@@ -404,7 +406,7 @@ def _clear_solution_cache(self) -> None:
         """
         Clear cached solution attributes.
         """
-        self._x = self._y = None
+        self._x = self._y = self._rc = None
         self._objval = self._dualobj = None
         self._gap = self._rel_gap = None
         self._status = None
@@ -423,6 +425,10 @@ def X(self) -> Optional[np.ndarray]:
     @property
     def Pi(self) -> Optional[np.ndarray]:
         return self._y
+    
+    @property
+    def RC(self) -> Optional[np.ndarray]:
+        return self._rc
 
     @property
     def ObjVal(self) -> Optional[float]:
diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp
index e0d4072..c8b4376 100644
--- a/python_bindings/_core_bindings.cpp
+++ b/python_bindings/_core_bindings.cpp
@@ -518,16 +518,19 @@ static py::dict solve_once(
     const int m_out = res->num_constraints;
     py::array_t<double> x({n_out});
     py::array_t<double> y({m_out});
+    py::array_t<double> rc({n_out});
     {
-        auto xb = x.request(), yb = y.request();
+        auto xb = x.request(), yb = y.request(), rcb = rc.request();
         std::memcpy(xb.ptr, res->primal_solution, sizeof(double) * n_out);
         std::memcpy(yb.ptr, res->dual_solution, sizeof(double) * m_out);
+        std::memcpy(rcb.ptr, res->reduced_costs, sizeof(double) * n_out);
     }
     // build info dict
     py::dict info;
     // solution
     info["X"] = x;
     info["Pi"] = y;
+    info["RC"] = rc;
     // objectives and gaps
     info["PrimalObj"] = res->primal_objective_value;
     info["DualObj"] = res->dual_objective_value;
diff --git a/src/presolve.c b/src/presolve.c
index 102ed0d..f5933dc 100644
--- a/src/presolve.c
+++ b/src/presolve.c
@@ -201,6 +201,18 @@ void pslp_postsolve(const cupdlpx_presolve_info_t *info,
     memcpy(result->primal_solution, info->presolver->sol->x, original_prob->num_variables * sizeof(double));
     memcpy(result->dual_solution, info->presolver->sol->y, original_prob->num_constraints * sizeof(double));
     memcpy(result->reduced_cost, info->presolver->sol->z, original_prob->num_variables * sizeof(double));
+    // TODO: this can be removed if PSLP implements it.
+    for (int i=0; i <original_prob->num_variables; i++)
+    {
+        if (!isfinite(original_prob->variable_lower_bound[i]))
+        {
+            result->reduced_cost[i] = fmin(result->reduced_cost[i], 0.0);
+        }
+        if (!isfinite(original_prob->variable_upper_bound[i]))
+        {
+            result->reduced_cost[i] = fmax(result->reduced_cost[i], 0.0);
+        }
+    }
     result->presolve_time = info->presolve_time;
     if (info->presolver->reduced_prob->n == 0)
     {
diff --git a/src/solver.cu b/src/solver.cu
index 062ec0f..40de675 100644
--- a/src/solver.cu
+++ b/src/solver.cu
@@ -56,7 +56,6 @@ __global__ void compute_delta_solution_kernel(
     const double *__restrict__ initial_dual,
     const double *__restrict__ pdhg_dual, double *__restrict__ delta_dual,
     int n_vars, int n_cons);
-static void rescale_solution(pdhg_solver_state_t *state);
 static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state, const lp_problem_t *original_problem);
 static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params);
 static void initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, const pdhg_parameters_t *params);
@@ -290,12 +289,24 @@ __global__ void compute_and_rescale_reduced_cost_kernel(
     const double *__restrict__ variable_rescaling,
     const double objective_vector_rescaling,
     const double constraint_bound_rescaling,
+    const double *__restrict__ variable_lower_bound,
+    const double *__restrict__ variable_upper_bound,
     int n_vars)
 {
     int i = blockIdx.x * blockDim.x + threadIdx.x;
     if (i < n_vars)
     {
-        reduced_cost[i] = (objective[i] - dual_product[i]) * variable_rescaling[i] / objective_vector_rescaling;
+        double rc = (objective[i] - dual_product[i]) * variable_rescaling[i] / objective_vector_rescaling;
+
+        if (!isfinite(variable_lower_bound[i]))
+        {
+            rc = fmin(rc, 0.0);
+        }
+        if (!isfinite(variable_upper_bound[i]))
+        {
+            rc = fmax(rc, 0.0);
+        }
+        reduced_cost[i] = rc;
     }
 }
 
@@ -902,15 +913,6 @@ static void compute_next_dual_solution(pdhg_solver_state_t *state,const int k_of
     }
 }
 
-static void rescale_solution(pdhg_solver_state_t *state)
-{
-    rescale_solution_kernel<<<state->num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>(
-        state->pdhg_primal_solution, state->pdhg_dual_solution,
-        state->variable_rescaling, state->constraint_rescaling,
-        state->objective_vector_rescaling, state->constraint_bound_rescaling,
-        state->num_variables, state->num_constraints);
-}
-
 static void perform_restart(pdhg_solver_state_t *state,
                             const pdhg_parameters_t *params)
 {
@@ -1202,9 +1204,15 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state, co
         state->variable_rescaling,
         state->objective_vector_rescaling,
         state->constraint_bound_rescaling,
+        state->variable_lower_bound,
+        state->variable_upper_bound,
         state->num_variables);
 
-    rescale_solution(state);
+    rescale_solution_kernel<<<state->num_blocks_primal_dual, THREADS_PER_BLOCK, 0, state->stream>>>(
+        state->pdhg_primal_solution, state->pdhg_dual_solution,
+        state->variable_rescaling, state->constraint_rescaling,
+        state->objective_vector_rescaling, state->constraint_bound_rescaling,
+        state->num_variables, state->num_constraints);
 
     results->primal_solution =
         (double *)safe_malloc(state->num_variables * sizeof(double));