Does every large prime satisfy $\sum_{\substack{ab\equiv1\pmod p}}\frac{1}{\sqrt{ab}}\longrightarrow 5 ? $

My experimental observation suggests that the sum of the reciprocals of the square roots of the products of all multiplicative-inverse pairs modulo a prime $p$ tends to $5$ as $p \to \infty$.

More specifically, for each prime $p$, consider the pairs $(a,b)$ satisfying $1\le a,b\le p-1$ and $ab\equiv1\pmod p$. Then

$$ \sum_{\substack{1\le a,b\le p-1\\ ab\equiv1\pmod p}}\frac{1}{\sqrt{ab}} \longrightarrow 5. $$

Can this be proved?

Source code:

%%time
import math, random

P = {
    10: 29, 100: 541, 1_000: 7919, 10_000: 104729,
    100_000: 1299709, 1_000_000: 15485863,
    10_000_000: 179424673, 100_000_000: 2038074743,
    1_000_000_000: 22801763489
}

def S_exact(p):
    inv = [0] * p
    inv[1] = 1
    s = 1.0
    for a in range(2, p):
        inv[a] = p - (p//a) * inv[p % a] % p
        s += 1 / math.sqrt(a * inv[a])
    return s

def S_mc(p, N=300_000):
    lo, hi = math.sqrt(2), math.sqrt(p)
    total = 0.0
    for _ in range(N):
        a = int((lo + (hi-lo)*random.random())**2)
        b = pow(a, -1, p)
        q = ((math.sqrt(a+1)-math.sqrt(a)) +
             (math.sqrt(b+1)-math.sqrt(b))) / (2*(hi-lo))
        total += 1 / (math.sqrt(a*b) * q)
    return 1 + total/N

for n, p in P.items():
    s = S_exact(p) if p < 2_000_000 else S_mc(p)
    print(f"n={n:>10,}   p={p:>14,}   S(p)={s:.8f}   error={s-5:+.8f}")