  We present a comprehensive study for the gamma-ray luminosity function (GLF)
of blazars and their contribution to the extragalactic diffuse gamma-ray
background (EGRB). Radio and gamma-ray luminosity correlation is introduced to
take into account the radio detectability which is important for the blazar
identification. Previous studies considered only pure luminosity evolution
(PLE) or pure density evolution, but we introduce the luminosity dependent
density evolution (LDDE) model, which is favored from the evolution of X-ray
luminosity function (XLF) of AGNs. The model parameters are constrained by
likelihood analyses about the observed redshift and gamma-ray flux
distributions of the EGRET blazars. We find that the LDDE model gives a better
fit to the observed distributions than the PLE model, indicating that the LDDE
model is also appropriate for gamma-ray blazars, and that the jet activity is
universally correlated with the accretion history of AGNs. The normalization
between the GLF and XLF is consistent with the unified picture of AGNs, when
the beaming and a reasonable duty cycle of jet activity are taken into account.
We then find that only 25--50% of the EGRB can be explained by unresolved
blazars with the best-fit LDDE parameters. Unresolved blazars can account for
all the EGRB only with a steeper index of the faint-end slope of the GLF, which
is marginally consistent with the EGRET data but inconsistent with that of the
XLF. Therefore unresolved AGNs cannot be the dominant source of the EGRB,
unless there is a new population of gamma-ray emitting AGNs that evolves
differently from the XLF of AGNs. Predictions for the GLAST mission are made,
and we find that the best-fit LDDE model predicts about 3000 blazars in the
entire sky, which is considerably fewer than a previous estimate.
