We present an in-depth and systematic variability study of a sample of 20 powerful blazars, including 12 BL Lacs and 8 flat-spectrum radio quasars, applying various analysis tools such as flux distribution, symmetry analysis, and time-series analysis on the decade-long Fermi/LAT observations. The results show that blazars with steeper γ-ray spectral indexes are found to be more variable, and the γ-ray flux distribution closely resembles a log-normal probability distribution function. The statistical variability properties of the sources as studied by power spectral density analysis are consistent with flicker noise (P(ν) ∝ 1/ν)-an indication of long-memory processes at work. Statistical analysis of the distribution of flux rise and decay rates in the light curves of the sources, aimed at distinguishing between particle acceleration and energy-dissipation timescales, counterintuitively suggests that both kinds of rates follow a similar distribution and the derived mean variability timescales are on the order of a few weeks. The corresponding emission region size is used to constrain the location of γ-ray production sites in the sources to be a few parsecs. Additionally, using Lomb-Scargle periodogram and weighted wavelet z-transform methods and extensive Monte Carlo simulations, we detected year-timescale quasi-periodic oscillations in the sources S5 0716+714, Mrk 421, ON +325, PKS 1424-418, and PKS 2155-304. The detection significance was computed taking proper account of the red noise and other artifacts inherent in the observations. We explain the results in light of current blazar models with relativistic shocks propagating down the jet viewed close to the line of sight.