At a quantum critical point, the low-energy physics of a quantum spin chain is described by conformal field theory (CFT). Given the Hamiltonian of a critical spin chain, we propose a variational method to build an approximate lattice representation ϕ_{α} of the corresponding primary CFT operators ϕ_{α}^{CFT}. We then show how to numerically compute the operator product expansion coefficients C_{αβγ}^{CFT} governing the fusion of two primary fields. In this way, we complete the implementation of Cardy's program, outlined in the 1980s, which aims to extract the universality class of a phase transition, as encoded in the so-called conformal data of the underlying CFT, starting from a microscopic description. Our approach, demonstrated here for the critical quantum Ising model, only requires a generic (i.e., in general, nonintegrable) critical lattice Hamiltonian as its input.