  In this paper and upcoming ones, we initiate a systematic study of Bethe
ansatz equations for integrable models by modern computational algebraic
geometry. We show that algebraic geometry provides a natural mathematical
language and powerful tools for understanding the structure of solution space
of Bethe ansatz equations. In particular, we find novel efficient methods to
count the number of solutions of Bethe ansatz equations based on Gr\"obner
basis and quotient ring. We also develop analytical approach based on companion
matrix to perform the sum of on-shell quantities over all physical solutions
without solving Bethe ansatz equations explicitly. To demonstrate the power of
our method, we revisit the completeness problem of Bethe ansatz of Heisenberg
spin chain, and calculate the sum rules of OPE coefficients in planar
$\mathcal{N}=4$ super-Yang-Mills theory.
