We study the phase behavior of a mixture of colloidal hard rods with a length-to-diameter ratio of L/sigma(c)=5 and nonadsorbing ideal polymer. We map our binary mixture onto an effective one-component system by integrating out the degrees of freedom of the polymer coils. We derive a formal expression for the exact effective Hamiltonian of the colloidal rods, i.e., it includes all effective many-body interactions and it is related to the exact free volume available for the polymer. We determine numerically on a grid the free volume available for the ideal polymer coils "on the fly" for each colloidal rod configuration during our Monte Carlo simulations. This allows us to go beyond first-order perturbation theory, which employs the pure hard-rod system as reference state. We perform free energy calculations for the isotropic, nematic, smectic, and crystal phase using thermodynamic integration and common tangent constructions are used at fixed polymer fugacities to map out the phase diagram. The phase behavior is determined for size ratios q=sigma(p)/sigma(c)=0.15, 0.5, and 1, where sigma(p) is the diameter of the polymer coils. The phase diagrams based on the full effective Hamiltonian are compared with those obtained from first-order perturbation theory, from simulations using the effective pair potential approximation to the effective Hamiltonian, and with those based on an empiric effective depletion potential for the rods. We find that the many-body character of the effective interactions stabilizes the nematic and smectic phases for large q, while the effective pair potential description overestimates the attractive interactions and favors, hence, a broad isotropic-crystal coexistence.