Many of the existing stochastic models of gene expression contain the first-order decay reaction term that may describe active protein degradation or dilution. If the model variable is interpreted as the molecule number, and not concentration, the decay term may also approximate the loss of protein molecules due to cell division as a continuous degradation process. The seminal model of that kind leads to gamma distributions of protein levels, whose parameters are defined by the mean frequency of protein bursts and mean burst size. However, such models (whether interpreted in terms of molecule numbers or concentrations) do not correctly account for the noise due to protein partitioning between daughter cells. We present an exactly solvable stochastic model of gene expression in cells dividing at random times, where we assume description in terms of molecule numbers with a constant mean protein burst size. The model is based on a population balance equation supplemented with protein production in random bursts. If protein molecules are partitioned equally between daughter cells, we obtain at steady state the analytical expressions for probability distributions similar in shape to gamma distributions, yet with quite different values of mean burst size and mean burst frequency than would result from fitting of the classical continuous-decay model to these distributions. For random partitioning of protein molecules between daughter cells, we obtain the moment equations for the protein number distribution and thus the analytical formulas for the squared coefficient of variation.