We study the gamma -ray variability of 13 blazars observed with the Fermi/Large Area Telescope (LAT). These blazars have the most complete light curves collected during the first four years of the Fermi sky survey. We model them with the Ornstein-Uhlenbeck (OU) process or a mixture of the OU processes. The OU process has power spectral density (PSD) proportional to 1/f super( alpha ) with alpha changing at a characteristic timescale, tau sub(0), from 0(tau >> tau sub(0)) to 2 (tau Lt tau sub(0)). The PSD of the mixed OU process has two characteristic timescales and an additional intermediate region with 0 < alpha < 2. We show that the OU model provides a good description of the Fermi/LAT light curves of three blazars in our sample. For the first time, we constrain a characteristic gamma -ray timescale of variability in two BL Lac sources, 3C 66A and PKS 2155-304 (tau sub(0) Asymptotically = to 25 days and tau sub(0) Asymptotically = to 43 days, respectively, in the observer's frame), which are longer than the soft X-ray timescales detected in blazars and Seyfert galaxies. We find that the mixed OU process approximates the light curves of the remaining 10 blazars better than the OU process. We derive limits on their long and short characteristic timescales, and infer that their Fermi/LAT PSD resemble power-law functions. We constrain the PSD slopes for all but one source in the sample. We find hints for sub-hour Fermi/LAT variability in four flat spectrum radio quasars. We discuss the implications of our results for theoretical models of blazar variability.