We formalize the procedure of functional development, in a general theoretical framework. Expansion in a functional basis set, and fitting via an error functional to a data set, casts functional development as a variational problem to obtain the functional basis-set and data-set limits. Overfitting is avoided by defining the optimum number of parameters. We implement our theory for an investigation of first- and second-order generalized gradient approximations (GGA) to the exchange-correlation and kinetic energy functionals, within an ab initio model. A variety of functional basis sets, including a general finite-element representation, is constructed to represent both one-dimensional and multidimensional GGA enhancement factors. An extensible data set consisting of 429 atomic and diatomic, neutral and cationic species, at stretched and equilibrium geometries, is constructed from Moller-Plesset level exchange-correlation energies, and Hartree-Fock kinetic energies. The range of chemically relevant density and gradient variables is examined. Exhaustive fitting investigations are carried out, to determine the accuracy of the GGA representation of the ab initio models. In the exchange-correlation case we demonstrate that we can reach the functional basis-set and data-set limit, which correspond to a root-mean-square (rms) error of similar to 10 mH (6.3 kcal/mol). Changing the functional basis set, higher-order density variables such as the kinetic energy density, multidimensional enhancement factors, and exact exchange yield no significant improvement, and our fits represent an effective solution of the GGA problem for exchange-correlation, at the Moller-Plesset level. In the kinetic energy case, accurate functionals with rms errors of similar to 80 mH (50 kcal/mol) are developed. These exhibit a beautifully simple kinetic energy enhancement factor, and are a step towards orbital-free calculations. (C) 2000 American Institute of Physics. [S0021-9606(00)30110-6].