curvedsky.rand_alm crashes with single-precision sims if lmax less than the ps lmax

[zatkins2]

• pixell version: 0.28.0
• numpy version: 1.26.4
• Python version: 3.10.16
• Operating System: rhel 8.10

Description

Drawing sims with curvedsky.rand_alm crashes under the following conditions:

• User supplies lmax that is smaller than the lmax that would correspond to the shape of ps
• The dtype is np.complex64

What I Did

Here's an example script that captures the bug.

from pixell import enmap, curvedsky, utils, enplot
import numpy as np

# make a geometry to put maps onto without aliasing
shape, wcs = enmap.fullsky_geometry(res=8*utils.arcmin, variant='fejer1')
print(shape, wcs)

# make an arbitrary power spectrum with an lmax of 1000
ps = np.exp(-np.arange(1001)/100 + 1)

# double prec sim
ps_dtype = np.float64
alm_dtype = np.complex128
map_dtype = np.float64

# if lmax of my realization is >= 1000, pads ps with 0 (fine)
omap = curvedsky.alm2map(curvedsky.rand_alm(ps.astype(ps_dtype), lmax=1300, seed=123, dtype=alm_dtype), enmap.empty(shape, wcs, map_dtype))

# if lmax of my realization is < 1000, still OK
omap = curvedsky.alm2map(curvedsky.rand_alm(ps.astype(ps_dtype), lmax=500, seed=123, dtype=alm_dtype), enmap.empty(shape, wcs, map_dtype))

# single prec sim
ps_dtype = np.float32
alm_dtype = np.complex64
map_dtype = np.float32

# if lmax of my realization is >= 1000, pads ps with 0 (fine)
omap = curvedsky.alm2map(curvedsky.rand_alm(ps.astype(ps_dtype), lmax=1300, seed=123, dtype=alm_dtype), enmap.empty(shape, wcs, map_dtype))

# if lmax of my realization is < 1000, now this crashes
omap = curvedsky.alm2map(curvedsky.rand_alm(ps.astype(ps_dtype), lmax=500, seed=123, dtype=alm_dtype), enmap.empty(shape, wcs, map_dtype))

# issue is specific to the alm_dtype being single prec
ps_dtype = np.float64
alm_dtype = np.complex64
map_dtype = np.float64

# if lmax of my realization is >= 1000, pads ps with 0 (fine)
omap = curvedsky.alm2map(curvedsky.rand_alm(ps.astype(ps_dtype), lmax=1300, seed=123, dtype=alm_dtype), enmap.empty(shape, wcs, map_dtype))

# if lmax of my realization is < 1000, this also crashes
omap = curvedsky.alm2map(curvedsky.rand_alm(ps.astype(ps_dtype), lmax=500, seed=123, dtype=alm_dtype), enmap.empty(shape, wcs, map_dtype))

The error in either of the crashing cases is

free(): corrupted unsorted chunks
Aborted (core dumped)

If you print the output of curvedsky.rand_alm rather than project it into a map, you get some weird looking elements in the first crashing case, before the same error message occurs:

# double prec, if lmax of my realization is >= 1000, pads ps with 0 (fine)
[-1.78990227+0.j  0.46422573+0.j -3.96110789+0.j ...  0.        +0.j
  0.        +0.j  0.        +0.j]

# double prec, if lmax of my realization is < 1000, still OK
[-1.78990227+0.j          0.46422573+0.j         -3.96110789+0.j
 ... -0.06780981-0.04358139j  0.05597249-0.04199121j
  0.13843811+0.0326628j ]

# single prec, if lmax of my realization is >= 1000, pads ps with 0 (fine)
[-1.7899024 +0.j  0.46422574+0.j -3.9611077 +0.j ...  0.        +0.j
  0.        +0.j  0.        +0.j]

# single prec, if lmax of my realization is < 1000, now this crashes
[-1.7899024e+00+0.0e+00j  4.6422574e-01+0.0e+00j -3.9611077e+00+0.0e+00j
 ... -3.5032462e-44-2.2e-44j  2.5223372e-44-1.8e-44j
  5.6051939e-45+1.4e-45j]