diff --git a/src/algorithms/solvers.jl b/src/algorithms/solvers.jl
index b9f717c..866c14f 100644
--- a/src/algorithms/solvers.jl
+++ b/src/algorithms/solvers.jl
@@ -187,7 +187,8 @@ function _core_msolve(
         la_option::Int=2,                     # linear algebra option
         get_param::Int=1,                     # get rational parametrization
         info_level::Int=0,                    # info level for print outs
-        precision::Int=32                     # precision of the solution set
+        precision::Int=32,                    # precision of the solution set
+        worker_pool::Union{Nothing,AbstractWorkerPool}=nothing, # worker pool for parallel computation
         )
 
     F = I.gens
@@ -207,9 +208,14 @@ function _core_msolve(
     # convert Singular ideal to flattened arrays of ints
     lens, cfs, exps, nr_gens = _convert_to_msolve(F)
 
-    jl_len, jl_vnames, jl_cf_lf, jl_cf, jl_sols_num, jl_sols_den =
-        _core_msolve_array(lens, cfs, exps, variable_names, field_char;
-            initial_hts, nr_thrds, max_nr_pairs, la_option, get_param, info_level, precision)
+    run_core_array = () -> _core_msolve_array(lens, cfs, exps, variable_names, field_char;
+        initial_hts, nr_thrds, max_nr_pairs, la_option, get_param, info_level, precision)
+
+    jl_len, jl_vnames, jl_cf_lf, jl_cf, jl_sols_num, jl_sols_den = if isnothing(worker_pool)
+        run_core_array()
+    else
+        remotecall_fetch(run_core_array, worker_pool)
+    end
 
     # convert to julia array, also give memory management to julia
     jl_ld = Int32(length(jl_len))
@@ -490,7 +496,8 @@ function real_solutions(
         la_option::Int=2,                     # linear algebra option
         info_level::Int=0,                    # info level for print outs
         precision::Int=32,                    # precision of the solution set
-        interval::Bool=false                  # return real solutions as intervals
+        interval::Bool=false,                  # return real solutions as intervals
+        worker_pool::Union{Nothing,AbstractWorkerPool}=nothing, # worker pool for parallel computation
         )
 
     isdefined(I, :real_sols) ||
@@ -500,7 +507,8 @@ function real_solutions(
                  max_nr_pairs = max_nr_pairs,
                  la_option = la_option,
                  info_level = info_level,
-                 precision = precision)
+                 precision = precision,
+                 worker_pool = worker_pool)
 
     if interval
         return I.inter_sols
diff --git a/test/algorithms/solvers.jl b/test/algorithms/solvers.jl
index a6fe62a..43e0bf8 100644
--- a/test/algorithms/solvers.jl
+++ b/test/algorithms/solvers.jl
@@ -1,3 +1,5 @@
+using Distributed
+
 @testset "Algorithms -> Solvers" begin
     R, (x1,x2,x3,x4) = polynomial_ring(QQ,["x1","x2","x3","x4"], internal_ordering=:degrevlex)
     I = Ideal([x1 + 2*x2 + 2*x3 + 2*x4 - 1,
@@ -122,3 +124,40 @@ end
     @test (res_len, res_vnames, res_cf_lf, res_cf, res_sols_num, res_sols_den) == AlgebraicSolving._core_msolve_array(
         lens, cfs, exps, variable_names, field_char; get_param=2)
 end
+
+@testset "Algorithms -> Solvers with workers" begin
+    R, (x, y) = polynomial_ring(QQ, ["x", "y"])
+
+    F1 = [x - 1, y + 2, R(0)]
+    H1 = Vector{QQFieldElem}[
+        [1, -2]
+    ]
+
+    F2 = [x^2 - 1, y]
+    H2 = Vector{QQFieldElem}[
+        [-1, 0],
+        [1, 0]
+    ]
+
+    F3 = [x^2 - 1, y^2]
+    H3 = H2
+
+    nb_tests = 42
+    F = [F1, F2, F3]
+    G = Vector{Vector{Vector{QQFieldElem}}}(undef, nb_tests)
+    H = [H1, H2, H3]
+
+    nb_workers = 2
+    worker_ids = addprocs(nb_workers; exeflags="--project=$(Base.active_project())")
+    @everywhere worker_ids using AlgebraicSolving
+    worker_pool = WorkerPool(worker_ids)
+
+    Threads.@threads for i in 1:nb_tests
+        G[i] = real_solutions(Ideal(F[i%3+1]), worker_pool=worker_pool)
+    end
+    for i in 1:nb_tests
+        @test issetequal(G[i], H[i%3+1])
+    end
+
+    rmprocs(worker_ids)
+end