Using a significantly enlarged 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 flat spectrum radio quasars, 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 modelled spectra, BL Lacs contribute similar to 20 per cent of the total extragalactic gamma-ray background (EGB) at E > 100 MeV, similar to 100 per cent of the EGB at E > 50 GeV, and the unresolved BL Lacs contribute similar to 20 per cent of the isotropic diffuse gamma-ray background at E > 100 MeV. A prediction of the TeV EGB spectra are given, which may be tested by the future detectors.