We present an accurate and computationally efficient first-principles methodology to calculate natural optical activity. Our approach is based on the long-wave density-functional perturbation theory and includes self-consistent field terms naturally in the formalism, which are found to be of crucial importance. The final result is expressed exclusively in terms of response functions to uniform field perturbations and avoids troublesome summations over empty states. Our strategy is validated by computing the natural optical activity tensor in representative chiral crystals (trigonal Se, α-HgS, and α-SiO_{2}) and molecules (C_{4}H_{4}O_{2}), finding excellent agreement with experiment and previous theoretical calculations.