Trigonometry / Euclidean Geometry for natural numbers? 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()

Second 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$$
Third edit:
I think the main property which distinguishes $d$ for example from the Jaccard or other metrics is the proven property ( https://mathoverflow.net/a/342921/6671) :
For all $a \neq b, a\neq c$ we have:
$$d(a,b)+d(a,c) > 1$$
I have tested other metrics with this property and they also seem to embedd in Euclidean Space. On the other hand metrics who do not have this property do not seem to embedd. So I think this is the point to be taken into consideration.
If someone has an idea how to exploit this property that would be very nice!
 A: It can be done for the metric 
$$d(a,b)^2 = 1 -  \frac{(a,b)}{\sqrt{ab}},$$
and other similar ones like $d(a,b)^2 = 1 -  \frac{(a,b)^2}{ab}$, with some twists in the construction.
Suppose we want to embed $1,2,..., n$ in $\mathbb{R}^n$. We will first embed these in $\mathbb{R}^m$, where $m = lcm(1,2,...,n)$.
For each natural number $k\in\{1,...,n\}$ map it to the vector $v_k \in \mathbb{R}^m$ whose $i$-th entry is equal to $\sqrt{k}$ if $i$ is a multiple of $k$ and $0$ otherwise. Noticing that the vectors $v_a$ and $v_b$ are only both non-zero at the entries multiple of $[a,b] = lcm(a,b)$ we get:
$$\|v_a-v_b\|^2 = (\frac{m}{a}-\frac{m}{[a,b]})a+(\frac{m}{b}-\frac{m}{[a,b]})b+\frac{m}{[a,b]}(\sqrt{a}-\sqrt{b})^2$$
$$=2m(1 - \frac{\sqrt{ab}}{[a,b]}) = 2m(1-\frac{(a,b)}{\sqrt{ab}}).$$
This means that after normalization by $2m$ we get the desired embedding. For an embedding in $\mathbb{R}^n$ take the induced embedding in the subspace $span(v_1,...v_n)$.

Another nice embedding straight to a Hilbert space follows from the identity for any natural numbers $a,b$
$$\int_0^1 \psi(at)\psi(bt) dt = \frac{1}{12} \frac{(a,b)^2}{ab}.$$
Where $\psi(t) = t - \lfloor t \rfloor  - \frac{1}{2}$ is the sawtooth function. Therefore for any natural numbers $a,b$
$$\|\psi(at) - \psi(bt)\|_{L^2}^2 = \frac{1}{6}(1-\frac{(a,b)^2}{ab}).$$
So the embedding $\mathbb{N} \hookrightarrow L^2([0,1])$ taking $n \mapsto \psi(nt)$ (also normalizing by $\frac{1}{6}$) preserves this metric!
From the point of view of Fourier series this construction is similar to the previous one, noticing that $\psi(nt)$ only has non-zero Fourier coefficients at the entries divisible by $n$.
A: Given that the set of integers has fractal dimension -1, I would not be surprized that such trigonometry is possible, it would be trigonometry on a manifold of negative dimension. Particularly, the set of prime numbers would play the role of the -2-sphere.
A: I just add references for a mentioned identity by previous answer that is a known identity by Franel ([1]), and Landau ([2]). 
It is an identity related in the study of an equivalent form to the Riemann hypothesis.
References:
[1] J. Franel, Les suites de Farey et le problème des nombres premiers, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse (1924), pp 198-201. 
[2] Edmund Landau, Bemerkungen zu der vorstehenden Abhandlung von Herrn Franel, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, (1924), Mathematisch-Physikalische Klasse, pp 202–206.
A: The answer given are already very good. I just wanted to point, that there are also an infinte family of metrics who embed in euclidean space.
First notice that for finite sets $X,Y$ contained in a larger finite set $Z$, the symmetric difference metric $d(X,Y) = \sqrt{|X|+|Y|-2|X \cap Y|}$ can be embedded in euclidean space by listing the elements of $Z$ in an ordered way and the vector $\phi(X)$ is a binary vector with the $i$-th entry $=1$ if $z_i \in X$ and $0$ otherwise.
Then $|\phi(X)| = |X|$ and $|\phi(X)-\phi(Y)|^2 = d(X,Y)$, which shows the embedding.
By considering the sets $X_a = \{ a/k | 1 \le k \le a \}$ and noticing that
$|X_a \cap X_b| = \gcd(a,b)$ we get the metric on natural numbers:
$d(a,b) = \sqrt{|X_a|+|X_b|-2|X_a\cap X_b|} = \sqrt{a+b-2 \gcd(a,b)}$ which can be embedded as was shown in Euclidean space.
On the other hand if we consider sets $X_a$ such that $|X_a| = \sigma_k(a)$ where $k \ge 0$ and $\sigma_k(a) = \sum_{d|a}d^k$, which are not difficult to construct, and such that $|X_a \cap X_b| = \sigma(\gcd(a,b))$, we get an infinte sequences of metrics, which can be embedded in Euclidean space:
$$d_{\sigma,k}(a,b) = \sqrt{\sigma_k(a)+\sigma_k(b)-2\sigma_k(\gcd(a,b))}$$.
For $k=0$ and $\tau(a) = \sigma_0(a)$ we observe,that primes $p$ have norm equal to one:
$$|p|:= d_{\sigma,0}(1,p) = \sqrt{1+2-2\cdot 1}=1$$
Hence under this metric all primes are on the unit sphere.
Especially for $k=1$ and $p,q,r$ three pairwise distinct primes, we get invoking the law of cosines and $d_1(p,q)^2 = p+q$ the following nice formula:
$$\pi = \operatorname{acos}(\frac{r}{\sqrt{(p+r)(q+r)}})+\operatorname{acos}(\frac{q}{\sqrt{(p+q)(q+r)}})+\operatorname{acos}(\frac{p}{\sqrt{(p+r)(p+q)}})$$
A: By the answer given by @DenisSerre on this question:
The abc-conjecture as an inequality for inner-products?
the metric:
$$d_L(a,b) = \sqrt{1-\frac{2\gcd(a,b)}{a+b}}$$
can be embedded in Euclidean Space and plays a role in the abc-conjecture.
A: Let $e_d$ be the $d$-th standard-basis vector in the Hilbert space $H=l_2(\mathbb{R})$.
Let $h(n)  = J_2(n)$ be the second Jordan totient function.
Define:
$$\phi(n) = \frac{1}{n} \sum_{d|n}\sqrt{h(d)} e_d$$.
Then we have:
$$ \left < \phi(a),\phi(b) \right > = \frac{1}{ab} \sum_{d|a,d|b} \sqrt{h(d)}^2 = $$
$$=\frac{1}{ab} \sum_{d|\gcd(a,b)} \sqrt{h(d)}^2 $$
$$=\frac{1}{ab} \sum_{d|\gcd(a,b)} h(d) =\frac{\gcd(a,b)^2}{ab}=:k(a,b)$$
So this is another embedding of the natural numbers in the Hilbert space of series.
Related question:
A geometric approach to the odd perfect number problem?
