I think "buckets" is a histogram-approach driven constraint. Histogram is also known as empirical-PDF, or ePDF.
If (yes, big word) you have sufficient samples (whatever that means) then you could contrive you distribution in empirical CDF instead. In that space, you can interpolate using appropriate (whatever that means) basis functions and using uniform random number generators create a bottomless pile of numbers drawn from a distribution that is similar enough, based on sample sufficiency, to the actual generating distribution. In the CDF domain, using slopes, you are not limited to bins, or finite number of digits. This also applies to the mode.
A more clever person would look at the slopes of the interpolating function, because in the limit it is the density. One could look at Gaussian Mixture models as equivalent to linear interpolation in (I think it is called) the probit domain (not the same as ECDF) and leverage the existing work showing their properties as interpolating functions for continuous probabilities, and so the normsinv as an interpolant in eCDF domain. This would allow handling of multi-modal distributions. Mode is about returning the locations of peaks.
I prefer to LSQ fit in eCDF domain, because repeated values are removed from the dataset, so it can be faster, and because the sum acts to smooth out centered noise. There are "textbook" ways to look at samples, error, and number of components, so that one selects the correct number of components. One could also look at a Kernel based approach, if they wanted. It might allow a 2-pass approach for the estimate.
