Algorithm used to compute one row of the covariance , using the HEALPix tools for a fixed ℓ′ and varying ℓ. Square, diamond, or circle boxes are arrays representing maps, power spectra, or spherical harmonic coefficients, respectively. Operations are symbolized by arrows and described alongside. The indices of the arrays are indicated on the side of the corresponding boxes. For example, near the top of the diagram, the HEALPIX.map2alm operation applied on the product Y_l’m′ W produces the array ₀I_LMl′m′, with indices L, M. At the bottom of the diagram, the operations before the final summation produce the array with indices ℓ,m for a fixed ℓ’,m’ pair. This part of the algorithm scales as O(ℓ’³) because this is the scaling of HEALPix operations (which we applied three times) when the resolution of the map is comparable to the maximum multipole index considered. The last operation, summing over the indices m, m′, requires repeating the precedent steps for all m′ ∊ [−ℓ′, ℓ′], thus repeating them 2ℓ′ + 1 times. Therefore, the final complexity for producing a single row with fixed multipole index ℓ′ is O(ℓ′⁴). Computing the covariance matrix for all ℓ′ then increases the computational time to O(ℓ′⁵).