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) Sine law
3) Cosine law
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
    
    X = range(1,1600)
    for i in X:
        i,j,k = nth_prime(i),nth_prime(i+1),nth_prime(i+2)
        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 = acos((skp(j,k)-skp(j,i)-skp(k,k)+skp(k,i))/(b*c))
        print area,a/sin(alpha).n(),a.n(),b.n(),c.n()

    for n in range(2,101):
        print n, MM(range(1,n)).is_positive_definite()