  Using a significantly enlarged \textit{Fermi}-LAT BL Lac objects (BL Lacs)
sample, we construct the gamma-ray luminosity function (GLF) of BL Lacs, by the
joint use of the space density distribution and source counts distribution. We
use three well-studied forms of the GLF, i.e., the forms of pure density
evolution (PDE), pure luminosity evolution (PLE) and luminosity-dependent
density evolution (LDDE). The Markov Chain Monte Carlo (MCMC) technique is used
to constrain model parameters. Our results suggest that LDDE model can give the
best description for the BL Lac GLF. And the model shows that the BL Lacs with
a harder GeV spectrum and a less luminosity evolve as strongly as FSRQs, and
the evolution decreases as increasing luminosity. We also model the average
photon spectra of BL Lacs with a double power laws model. Using this modeled
spectra, BL Lacs contribute $ \sim 20\%$ of the total extragalactic gamma-ray
background (EGB) at $E>100$ MeV, $\sim 100\%$ of the EGB at $E>50$ GeV, and the
unresolved BL Lacs contribute $\sim 20 \%$ of the isotropic diffuse gamma-ray
background (IGRB) at $E>100$ MeV. A prediction of the TeV EGB spectra are
given, which may be tested by the future detectors.
