We give a detailed account of the calculation of nonrenormalizable terms in the free fermionic formulation of the heterotic string in four dimensions. In this formulation such terms can be calculated exactly in σ-model and remain uncorrected to all orders in string loops, therefore considerably enhancing the predictive power of models built in this context. We obtain general results for the types of nonvanishing N-point function contributions to the superpotential, and present explicit computations of sample quartic and quintic terms, and corrections to a typical D-term, in a specific model. We also discuss the connection between our results and results obtained using the orbifold approach. Finally, we elucidate the Ashkin-Teller-model method of computation of arbitrary Ising model correlators needed in the evaluation of nonrenormalizable terms in realistic models.