Proving Mertens' theorem using the prime number theorem Mertens' Theorem states that
$$\sum_{p \leq x}\frac{1}{p} = \log \log x + M + O(1/\log x).$$
This is weaker than the prime number theorem; in fact according to the Wikipedia page, the prime number theorem is equivalent to
$$\sum_{p \leq x}\frac{1}{p} = \log \log x + M + o(1/\log x).$$
However no reference or proof is given, it just simply says that "Although this equivalence is not explicitly mentioned there, it can for instance be easily derived from the material in chapter I.3 of: G. Tenenbaum. Introduction to analytic and probabilistic number theory."
So what is a proof/reference for this fact? I naturally tried to prove it myself using partial summation, and it is easy to see where the $\log \log x$ comes from, however the constant $M$ and $o(1/\log x)$ are a bit more mysterious.
Standard proofs in prime number theory often proceed by introducing a logarithmic weight via the Von Mangoldt function, then proving an asymptotic and going back again. I would prefer to avoid such an approach as I have a different problem in mind coming from counting rational points on varieties where I can't introduce a logarithmic weight.
 A: That's not exactly OP's question, but I find it illuminating to look at the explicit formula for $\sum 1/p$. Allow me to work, for simplicity, with
$$\sum_{n \le x} \frac{\Lambda(n)}{n\log n}$$
which is closely related to the original sum because it equals
$$\sum_{p \le x}\frac{1}{p} + \sum_{p^k \le x, \, k \ge 2} \frac{1}{k p^k }.$$
(The contribution of higher powers can easily be absorbed in the constant $M$ and in the error term.)
There is an explicit formula for ${\sum_{n \le x}}' \Lambda(n)/n^s$ due to Landau (as usual, the $\prime$ in the sum indicates the last term is counted with weight $1/2$ if $x$ is a positive integer). Integrating it with respect to $s$, one obtains
$${\sum_{n \le x}}' \frac{\Lambda(n)}{n\log n} = \log \log x +\gamma -\sum_{\rho} \int_{0}^{\infty}\frac{x^{\rho-1-t}}{\rho-1-t}dt$$
where the sum is over all zeros of $\zeta$. So any given zero-free region will give a corresponding error term. The integral can be approximated as
$$\int_{0}^{\infty}\frac{x^{\rho-1-t}}{\rho-1-t}dt = \frac{x^{\rho-1}}{(\rho-1)\log x} \left(1+O\left( \frac{1}{\log x}\right)\right).$$
For instance, the Vinogradov--Korobov zero-free region implies
$$\sum_{p \le x} \frac{1}{p} = \log \log x + M + O(\exp(-c(\log x)^{3/5}(\log \log x)^{-1/5})).$$
(To be precise, one needs a truncated version of the above formula to deduce this.)
This explicit formula is known to experts, but a bit tricky to find in the literature. Two references, that do exactly as mentioned above (integrate Landau's formula) but with more detail:

*

*Youness Lamzouri, "A bias in Mertens' product formula", Int. J. Number Theory 12 (2016), no. 1, 97–109. In Proposition 2.1 you can find an explicit formula for $\sum_{p \le x} \log(1-1/p)$.

*Corollary 4.2, with $s=0$, in my arXiv preprint: https://arxiv.org/abs/2211.08973 .

A: The method described in GH from MO's answer requires three uses of partial summation (perhaps I'm mischaracterizing the use of Axer's theorem). But one can go directly from $\pi(x)$ to $\sum_{p\le x} \frac1p$ with a single application of partial summation, namely
$$
\sum_{p\le x} \frac1p = \int_{2^-}^x \frac1t \, d\pi(t) = \frac{\pi(x)}x + \int_2^x \frac{\pi(t)}{t^2} \, dt.
$$
For the present purposes, it's more convenient to start with
\begin{align*}
\sum_{p\le x} \frac1p - \log\log x &= \int_{2^-}^x \frac1t \, d\pi(t) - \int_e^x \frac1{t\log t} \, dt \\
&= \int_{2^-}^x \frac1t \, d\bigl( \pi(t) - \mathop{\rm li}(t) \bigr) - \log\log 2 \\
&= \frac{\pi(x) - \mathop{\rm li}(x)}x + \frac{\mathop{\rm li}(2)}2 - \log\log 2 + \int_2^x \frac{\pi(t)-\mathop{\rm li}(t)}{t^2} \,dt.
\end{align*}
One step that the OP might have missed is to write
$$
\int_2^x \frac{\pi(t)-\mathop{\rm li}(t)}{t^2} \,dt = \int_2^\infty \frac{\pi(t)-\mathop{\rm li}(t)}{t^2} \,dt - \int_x^\infty \frac{\pi(t)-\mathop{\rm li}(t)}{t^2} \,dt
$$
and remember to consider the convergent integral $\int_2^\infty \frac{\pi(t)-\mathop{\rm li}(t)}{t^2} \,dt$ as part of the constant $M$ rather than part of the error term.
A: Well, one can always say that the PNT is equivalent to
$$\sum_{p \leq x}\frac{1}{p} = \log \log x + M + o\left(\frac{1}{\log x}\right),\tag{$\ast$}$$
because both results are true (with better error terms). This is of course not what is meant by the Wikipedia page. Instead, the idea is that the equivalence PNT$\,\Leftrightarrow(\ast)$ can be established in a simpler way than either PNT or $(\ast)$. On the other hand, "simpler" is a subjective word, e.g. I usually find Tauberian arguments tricky.
At any rate, the first three exercises in Section 8.1.1 of Montgomery-Vaughan: Multiplicative number theory I address this question. For example, the PNT easily implies the relation $\psi(x)\sim x$ (logarithmic weights!), which then implies rather nontrivially (using Theorem 8.1 = Axer's theorem) that
$$\sum_{p\leq x}\frac{\log p}{p}=\log x+C+o(1)$$
for some constant $C$. From here, $(\ast)$ follows easily by partial summation.
A: In another answer here,
GH from MO shows PNT implies $\sum_{n \leq x} 1/p = \log x + C + o(1)$ for some constant $C$ (a hard step), and that implies $\sum_{p \leq x} 1/p = \log\log x + M + o(1/\log x)$ for some constant $M$ by partial summation. Here is a converse argument, from the estimate on $\sum_{p\leq x}1/p$ with error term $o(1/\log x)$ to PNT, thus establishing that there is an equivalence
$$
{\rm PNT} \Longleftrightarrow \sum_{p \leq x} \frac{1}{p} = \log\log x + M + o\left(\frac{1}{\log x}\right)
$$
for some constant $M$.
Set $A(x) = \sum_{p \leq x} 1/p$, so we assume $A(x) = \log\log x + M + o(1/\log x)$ for some $M$, as $x \to \infty$. Then
$$
\pi(x) = \sum_{p \leq x} 1 = \sum_{p \leq x} \frac{1}{p}p = \sum_{n \leq x} a_n n
$$
where $a_n = 1/p$ when $n = p$ is prime and $a_n = 0$ otherwise. Then $A(x) = \sum_{n \leq x} a_n$, so by partial summation
$$
\sum_{n \leq x} a_n n = A(x)x - \int_2^x A(y)\,dy
$$
for $x \geq 2$. Since $A(x) = \log\log x + M + o(1/\log x)$,
$$
A(x)x = x\log\log x + Mx + o\left(\frac{x}{\log x}\right)
$$
and
$$
\int_2^x A(y)\,dy = \int_2^x \log\log y\,dy + M(x-2) + \int_2^x o\left(\frac{1}{\log y}\right)\,dy.
$$
Using integration by parts,
$$
\int_2^x \log\log y\,dy = x\log\log x - 2\log\log 2 - \int_2^x\frac{dy}{\log y}.
$$
Putting these formulas into the expression for $\pi(x)$ as $\sum_{n \leq x} a_n n$, we get
$$
\pi(x) = \int_2^x \frac{dy}{\log y} + 2M + 2\log\log 2 + o\left(\frac{x}{\log x}\right) - \int_2^x o\left(\frac{1}{\log y}\right)\,dy.
$$
The constant $2M + 2\log\log 2$ can be absorbed into the $o(x/\log x)$ term.
Lastly, since $\int_2^x o(1/\log y)\,dy = o(\int_2^x dy/\log y)$ as $x \to \infty$, we have
$$
\pi(x) = \int_2^x \frac{dy}{\log y} + o\left(\frac{x}{\log x}\right), 
$$
which is PNT.
PS. The Wikipedia page you mention no longer refers to an equivalence between PNT and the estimate on $\sum_{p \leq x} 1/p$ with error term $o(1/\log x)$, but only a one-way implication from PNT to the estimate on $\sum_{p \leq x} 1/p$: the View History tab of that page shows an edit was made to it on Jan. 10, 2023 by RStanley31 (gee, who could that be... :)) with the comment "Result of equivalence is not found within the literature." Perhaps that page should be reverted to the original wording with a link to this page.
