In Newtonian gravity, the energy is proportional to the sum over pairs of masses of $-m_1{m}_2/r_{1,2}$. Or in the continous case, $-\int \int \frac{d{m}_{1}d{m}_2}{|pos\left(m_1\right)-pos\left(m_2\right){|}_2}$. This does not actually strike me as very simple in whatever language makes Maxwell's equations and the Schrödinger equation simple.

Partially-formed idea 1:

Here's an outline of something that seems like it might be a simpler theory of gravity. Given a density function $\rho$ over physical space (possibly including Kronecker deltas for point masses), we need to come up with a gravitational potential $U$. First, convolve $\rho$ with some (radially symmetric) function $f:R^3\to R$. Write this $\rho \ast f$. Let $F\{\rho \ast f\}$ be the Fourier transform. I think taking the Fourier transform of something in position space gives something in momentum space. Then, you can take a "measurement" of the momentum squared by evaluating ${\int }_{R^3}F\{\rho \ast f\}(\omega {)}^2{\omega }^{2}d\omega$. And I think the momentum squared is energy-like. If $\rho$ is the sum of two delta functions that are a distance $d$ apart, surely there is some $f$, such that for large $d$, this looks like $c_1-c_2/d$.

It is possible that this could be made to resemble normal Newtonian gravity in the case of a few point masses, or a single mass (like the sun) that is much more massive than other nearby masses. But it would maybe also give different-looking answers when there is lots of mass evenly spread out of large distances, on the scale of, say, rotating galaxies.

Even-less-formed idea 2:

If you are experiencing gravity, you don't actually know exactly where the masses are which are causing this gravity. Doing some heavy handwaving with the gist of the Heisenberg Uncertainty Principle, to the extent you know the precise position of something, then from your perspective, it has lots of momentum. Maybe as you become more able to judge the precise position of nearby gravitating masses, your relative momentum increases to compensate. (This would look like speeding up as you approach a mass).

If you're around other masses, you don't detect the gradient of gravitational potential, because you fall down it, and free-fall feels like floating. But, if you have volume, you can detect the Hessian of gravitational potential. In the extreme case, this manifests as spaghettification near a black hole. (Perhaps a cost of having volume is that you don't know exactly where you are, so relative positions have some inherent uncertainty.) But in any case, your uncertainty about the relative location of the distribution of mass around you may be dominated by the fact that many different mass distributions could lead to the same Hessian of gravitational potential.

I'm not sure how to deal with your uncertainty about the total mass in the universe, but I have the intuition that treating it as known wouldn't cause anything to break.

So then you have the unknown potential function $U\left(\rho \right)$, and let $\tilde{\rho }\left\{x\right\}$ be the alternative density function when you are at position $x\in R^3$, and you have some tiny mass $m$, and the rest of the mass is distributed according to $\rho$. One thing to keep track of then is $H_{x}U\left(\tilde{\rho }\right\{x\left\}\right){|}_{x=0}/m$. (That's a Hessian, not a Hamiltonian). Denote this matrix $A_U\left\{\rho \right\}$, so $A_U:\Delta R^3\to R^{3\times 3}$, and it's totally determined by the function $U$. Now, we want to consider $P_M=\{\rho \in \Delta R^3:A_U\{\rho \}=M\}$ for a measured Hessian $M$; it's the set of viable mass distributions that would have given rise to $M$. Next, we need a probability distribution over these mass distributions. Maybe there's a sensible uniform distribution over this set. Or maybe it's better to consider the distribution that puts weight on a mass distribution $\rho$ proportional to $\exp (-{\int }_{R^3}F\{\rho \left\}\right(\omega {)}^2{\omega }^{2}d\omega )$, using similar intuitions as above, and a vague recollection from QM about exponentiation being right here. In any case, the distribution over mass distributions can be flattened into a single mass distribution ${\rho }_U^M$. And then, let's say, $U\left(\rho \right)={\int }_{R^3}F\left\{{\rho }_U^{H_{x}U\left(\tilde{\rho }\right\{x\left\}\right){|}_{x=0}/m}\right\}(\omega {)}^2{\omega }^{2}d\omega$. Now, we just solve for $U$...

None of that includes a possible trade-off between measurement error of $M$ and measurement error of your own position, so maybe that would have to be added.

I hope it goes without saying this is intensely hypothetical.

Genre of ideas:

Those ideas are probably very uninformed, but there may be some simple functions $U\left(\rho \right)$ out there that have the right behavior in the limit where $\rho$ is two point masses far apart. And maybe they would explain galaxy rotation curves too.