A computational scheme for the evaluation of the static first (beta) and second (gamma) hyperpolarizability tensors of systems periodic in 1D (polymers), 2D (slabs), 3D (crystals), and, as a limiting case, 0D (molecules) has been implemented, within the coupled perturbed Hartree-Fock framework (CPHF), in the CRYSTAL code, which uses a Gaussian type basis set. This generalizes to 2D and 3D the work by Bishop et al. (J. Chem. Phys. 114, 7633 (2001)). CPHF is applied for beta and gamma (the polarizability tensor alpha is also reported for completeness) of LiF in different aggregation states: finite and infinite chains, slabs, and cubic crystal. Correctness of the computational scheme and its numerical efficiency are documented by the trend of beta and gamma for increasing dimensionality: for a finite linear chain containing N LiF units, the hyperpolarizability tends to the infinite chain value at large N, N parallel chains give the slab value when N is sufficiently large, and N superimposed slabs tend to the bulk value. High numerical accuracy can be achieved at relatively low cost, with a dependence on the computational parameters similar to that observed for field-free self-consistent field (SCF) calculations.