# Computing a linear combination of all reduced density matrices of a given pure state

**URL:** <https://discourse.julialang.org/t/computing-a-linear-combination-of-all-reduced-density-matrices-of-a-given-pure-state/72070>\
**Category:** Quantum\
**Created:** [November 25, 2021, 4:07pm UTC](https://discourse.julialang.org/t/computing-a-linear-combination-of-all-reduced-density-matrices-of-a-given-pure-state/72070 "2021-11-25T16:07:41Z")\
**Posts on this page:** 1\
**Showing post:** 8

<div class="post-metadata">

**Author:** ![rmedinar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rmedinar/32/25263_2.png) [@rmedinar](https://discourse.julialang.org/u/rmedinar)\
**Post date:** [November 29, 2021, 9:47pm UTC](https://discourse.julialang.org/t/computing-a-linear-combination-of-all-reduced-density-matrices-of-a-given-pure-state/72070/8 "2021-11-29T21:47:17Z")

</div>

Hi once again. Besides the function from QuantumOptics I manage with the help of a coworker to implement the following solution to the above issue. It uses the "parallel bit deposit and extract” instructions, which are part of BMI2. See [this](https://discourse.julialang.org/t/bit-manipulation-instruction-set/5368). This is the working code: it receives the reduced density matrix (`2^nA x 2^nA`) as well as another matrix ( `2^N x 2^N` ) which is then modified:

```julia
function get_mask(rA::Vector{Int}, N::Int)
    mask = UInt64(0)
    for a in rA
        mask = mask | (1 << (N-a))
    end 
    return mask
end

pdep(x::UInt64, y::UInt64) = ccall("llvm.x86.bmi.pdep.64", llvmcall, UInt64, (UInt64, UInt64), x, y)

function embed_rdm(rho:: Matrix{ComplexF64}, N::Int, rA::Vector{Int}, rhoToModify::Matrix{ComplexF64})
    nA = length(rA)
    nB = N - nA;

    maskA = get_mask(rA, N);
    maskB = get_mask(setdiff(collect(1:N), rA), N)

    for iB in 0:(2^nB - 1)
        sB = pdep(UInt64(iB), maskB)
        for iA1 in 0:(2^nA - 1)
            sA1 = pdep(UInt64(iA1), maskA)
            for iA2 in 0:(2^nA - 1)
                sA2 = pdep(UInt64(iA2), maskA)
                rhoToModify[(sB | sA1) + 1, (sB | sA2) + 1] += rho[iA1 + 1, iA2 + 1]/(2^nB)
            end
        end
    end
end

```

I hope is useful to someone. Thanks for all the help 🙂

---

_[View the full topic](https://discourse.julialang.org/t/computing-a-linear-combination-of-all-reduced-density-matrices-of-a-given-pure-state/72070)._