We propose an importance sampling scheme for the procedural glittering BSDF of Chermain et al. 2020. Glittering BSDFs have multi-lobe visible normal distribution functions (VNDFs) which are difficult to sample. They are typically sampled using a mono-lobe Gaussian approximation, leading to high variance and fireflies in the rendering. Our method optimally samples the multi-lobe VNDF, leading to lower variance and removing firefly artefacts at equal render time. It allows, for example, the rendering of glittering glass which requires an efficient sampling of the BSDF. The procedural VNDF of Chermain et al. is a finite mixture of tensor products of two 1D tabulated distributions. We sample the visible normals from their VNDF by first drawing discrete variables according to the mixture weights and then sampling the corresponding 1D distributions using the technique of inverse cumulative distribution functions (CDFs). We achieve these goals by tabulating and storing the CDFs, which uses twice the memory as the original work. We prove the optimality of our VNDF sampling and validate our implementation with statistical tests.