A, by now, well-known signal-to-noise problem plagues Monte Carlo calculations of quantum-information-theoretic observables in systems of interacting fermions, particularly the Renyi entanglement entropies Sn, even in many cases where the infamous sign problem does not appear. Several methods have been put forward to circumvent this affliction including ensemble-switching techniques using auxiliary partition-function ratios. This dissertation presents an algorithm that modifies the recently proposed free-fermion decomposition in an essential way: we incorporate the entanglement-sensitive correlations directly into the probability measure in a natural way. Implementing this algorithm, we demonstrate that it is compatible with the hybrid Monte Carlo algorithm, the workhorse of the lattice quantum chromodynamics community and an essential tool for studying gauge theories that contain dynamical fermions. By studying a simple one-dimensional Hubbard model, we demonstrate that our method does not exhibit the same debilitating numerical difficulties that naive attempts to study entanglement often encounter. Following that, we illustrate some key probabilistic insights, using intuition derived from the previous method and its successes to construct a simpler, better behaved, and more elegant algorithm. Using this method, in combination with new identities which allow us to avoid seemingly necessary numerical difficulties, the inversion of the restricted one-body density matrices, we compute high order Renyi entropies and perform a thorough comparison to this new algorithm's predecessor using the Hubbard model mentioned before. Finally, we characterize non-perturbatively the Renyi entropies of degree n=2,3,4, and 5 of three-dimensional, strongly coupled many-fermion systems in the scale-invariant regime of short interaction range and large scattering length, i.e. in the unitary limit using the algorithms detailed herein. We also detail an exact, few-body projective method which we use to characterize the entanglement properties of the two-body sector across a broad range of attractive couplings. In the many-body case, we determine universal scaling properties of this system, and for the two-body case, we compute the entanglement spectrum exactly, successfully characterizing a vast range of entanglement behavior across the BCS-BEC crossover.