Let $d(a,b) = 1 - \frac{2\gcd(a,b)^3}{ab(a+b)}$ be a metric on natural numbers without $0$.

The metric space $X = \{x_0,x_1,\cdots,x_n\},n>2$ is isometric embeddable in $\mathbb{R}^n$ if and only if the matrix: $$M(x_0,x_1,\cdots,x_n) = (1/2 (d(x_0,x_i)^2+d(x_0,x_j)^2-d(x_i,x_j)^2))_{1 \le i,j \le n}$$ is positive semidefinite.

So my question is:

Is the matrix above for $d$ as above positive semidefinite for all choices of $x_i \in \mathbb{N}$? (Maybe it is possible to prove this using quadratic forms and then transform it to $\sum_{i} a_{ii} y_i^2$ showing then that $a_{ii}\ge 0$?

If it is so, then this would one allow to do euclidean geometry of natural numbers. For instance for three (pairwise distinct) points / natural numbers we would have:

1. a triangle

2. law of sines

3. law of cosines

4. All other theorems concerning triangles

Then in the limit three consecutive numbers / primes would build an equilateral triangle of side length $1$. Hence one could imagine primes ("in the limit") as an infinite dimensional simplex, which would be a funny thing to think about.

Thanks for your help.

Related question: https://math.stackexchange.com/questions/3385102/is-this-metric-matrix-positive-semidefinite

See Theorem 2.4 in https://books.google.de/books?id=7_DuCAAAQBAJ&printsec=frontcover&hl=de&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false for isometrically embedding of $(\mathbb{N},d)$ in a Hilbert space.

Edit: Here is some Sage code in case one wants to check this numerically for some examples:

def dABC(a,b):
    """ABC"""
    return 1- 2*gcd(a,b)**3/(a*b*(a+b))

def MM(xx,d=dABC):
    N = len(xx)
    return matrix([[1/2*(d(xx[0],xx[i])**2+d(xx[0],xx[j])**2-d(xx[i],xx[j])**2) for i in range(1,N)] for j in range(1,N)])

def skp(a,b,d=dABC):
    return 1/2*(d(a,1)**2+d(b,1)**2-d(a,b)**2)

def schur(M):
    from scipy.linalg import schur
    import numpy as np
    M_np = np.matrix(M,dtype='float64')
    A,B = schur(M_np,output="complex")
    return (matrix(np.asmatrix(A)),matrix(np.asmatrix(B)))

def createEmbedding(rr):
    M = MM(rr)
    n = len(rr)+1
    A,B = schur(M)
    E = diagonal_matrix([sqrt(x) for x in A.diagonal()])
    X = B*E
    ee = [ matrix([[i==j] for i in range(1,n-1)],ring=QQ) for j in range(1,n-1)]
    #print ee
    xx = [ X.transpose()*ee[i] for i in range(n-2)]
    return xx

N = 20
for i in primes(N):
    for j in primes(i+1,N):
        for k in primes(j+1,N):
            a = dABC(i,j)
            b = dABC(j,k)
            c = dABC(k,i)
            s = 1/2*(a+b+c)
            area = sqrt(s*(s-a)*(s-b)*(s-c)).n()
            alpha = pi.n()-arccos((skp(j,k)-skp(j,i)-skp(k,k)+skp(k,i))/(b*c))
            beta = pi.n()-arccos((skp(j,i)-skp(k,j)-skp(i,i)+skp(i,k))/(a*c))
            gamma = pi.n()-arccos((skp(j,k)-skp(k,i)-skp(j,j)+skp(j,i))/(b*a))
            print i,j,k,"area:",area, "sum:",(alpha+gamma+beta).n(),pi.n()
            print i,j,k,"sine law:",a/sin(alpha).n(),b/sin(beta).n(),c/sin(gamma).n()
            print i,j,k,"lengths:", a.n(),b.n(),c.n()
            print i,j,k,"cosine law: c", c**2.0,(a**2+b**2-2*a*b*cos(gamma)).n(),cos(gamma).n()
            print i,j,k,"cosine law: b", b**2.0,(c**2+a**2-2*c*a*cos(beta)).n(),cos(beta).n()
            print i,j,k,"cosine law: a", a**2.0,(c**2+b**2-2*c*b*cos(alpha)).n(),cos(alpha).n()
for n in range(2,101):
    print n, MM(range(1,n)).is_positive_definite()

Edit: Just out of curiosity: For $(a,b,c)=(1,2,2k+1)$, so $c \ge 3$ is odd, we get using the sum of angles in a triangle:

$$\alpha + \beta + \gamma = \pi$$

the following curious identity. For each odd $c \ge 3$ we have:

$$\operatorname{acos}(\frac{4 \, c^{5} + 28 \, c^{4} + 62 \, c^{3} + 2 \, c^{2} - 153 \, c - 135}{12 \, {\left(c + 2\right)}^{3} {\left(c + 1\right)} c} ) +$$ $$ \operatorname{acos}(\frac{14 \, c^{5} + 98 \, c^{4} + 226 \, c^{3} + 142 \, c^{2} - 135 \, c - 153}{18 \, {\left(c^{2} + 2 \, c - 1\right)} {\left(c + 2\right)}^{2} {\left(c + 1\right)}}) + $$ $$\operatorname{acos}(\frac{4 \, c^{6} + 24 \, c^{5} + 70 \, c^{4} + 156 \, c^{3} + 187 \, c^{2} - 18 \, c - 135}{12 \, {\left(c^{2} + 2 \, c - 1\right)} {\left(c + 2\right)} {\left(c + 1\right)}^{2} c}) = \pi$$