Inconsistent results from metrics.pairwise.euclidean_distances

### Describe the bug

Calculations on same set of floating point values return different results depending on object passed in.

### Steps/Code to Reproduce
```python
from sklearn.metrics import pairwise_distances
from sklearn.metrics.pairwise import euclidean_distances
import numpy as np

numbers2 = np.array([[0, 0, 0],
                     [4765960.395804689, -2726208.244406798, 3236212.738643452]])


def test(desc, arr1, arr2, idx):
    if arr1 is arr2:
        match="Object"
    elif (arr1 == arr2).all():
        match="Numeric"
    else:
        match="No"
        
    dist = euclidean_distances(arr1, arr2, squared=True)

    print("%s\t%s,%s\t" % (desc, match, str(dist.dtype)), dist[tuple(idx)])
        
print("Squared  distance directly")
test("final element only:\t", numbers2[1:], numbers2[1:], [0,0])
test("final element and copy:\t", numbers2[1:], numbers2[1:].copy(), [0,0])
print("\nOther cases:")
test("final element and copy:\t", numbers2[1:].copy(), numbers2[1:].copy(), [0,0])
test("same array:\t", numbers2, numbers2, [1,1])
test("array and copy:\t", numbers2, numbers2.copy(), [1,1])
test("array and slice:\t", numbers2[:], numbers2, [1,1])
test("array and final element:\t", numbers2[1:], numbers2, [0,1])
test("array and final w/ delta:\t", numbers2[1:], numbers2+1e-5, [0,1])
test("copy with delta:\t", numbers2, numbers2.copy() + 1e-5, [1,1])
test("both delta:\t", numbers2 + 1e-5, numbers2 + 1e-5, [1,1])

print("sklearn euclidean formuala") # https://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.euclidean_distances.html
def simple_dist(x, y):
    return np.sqrt(np.dot(x, x) - 2 * np.dot(x, y) + np.dot(y, y))
print("calculation:\t", simple_dist(numbers2[1], numbers2[1]))
```

### Expected Results

One of two results is expected:
1. All distances printed by `test` are 0
2. The distance for "same array" is 0, all other distances calculated are an identical value very close to 0




### Actual Results
```
Squared  distance directly
final element only:		Numeric,float64	 0.015625
final element and copy:		Numeric,float64	 0.0

Other cases:
final element and copy:		Numeric,float64	 0.0
same array:		Object,float64	 0.0
array and copy:		Numeric,float64	 0.015625
array and slice:		Numeric,float64	 0.015625
array and final element:		No,float64	 0.0
array and final w/ delta:		No,float64	 0.0
copy with delta:		No,float64	 0.0
both delta:		Numeric,float64	 0.0
sklearn euclidean formuala
calculation:	 0.0
```

### Versions

```shell
>>> import sklearn; sklearn.show_versions()

System:
    python: 3.9.16 (main, Dec  8 2022, 00:00:00)  [GCC 11.3.1 20221121 (Red Hat 11.3.1-4)]
executable: .../env/bin/python3
   machine: Linux-5.14.0-284.11.1.el9_2.x86_64-x86_64-with-glibc2.34

Python dependencies:
      sklearn: 1.2.2
          pip: 23.1.2
   setuptools: 53.0.0
        numpy: 1.19.5
        scipy: 1.5.4
       Cython: None
       pandas: 1.1.5
   matplotlib: 3.6.2
       joblib: 1.2.0
threadpoolctl: 2.2.0

Built with OpenMP: True

threadpoolctl info:
       filepath: .../env/lib/python3.9/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
         prefix: libgomp
       user_api: openmp
   internal_api: openmp
        version: None
    num_threads: 16

       filepath: .../env/lib/python3.9/site-packages/numpy.libs/libopenblasp-r0-09e95953.3.13.so
         prefix: libopenblas
       user_api: blas
   internal_api: openblas
        version: 0.3.13
    num_threads: 16
threading_layer: pthreads
   architecture: Prescott

       filepath: .../env/lib/python3.9/site-packages/scipy.libs/libopenblasp-r0-085ca80a.3.9.so
         prefix: libopenblas
       user_api: blas
   internal_api: openblas
        version: 0.3.9
    num_threads: 16
threading_layer: pthreads
   architecture: Prescott
```

COMMENT:
Thanks for the report @sa2329. There are 2 things here:

- The distance between numbers2[1:] and numbers2[:1] is not always 0. It's inherent to floating point arithmetic. We can only expect a result up to some precision. Working with float64 input, any operation is guaranteed up to a relative difference of ~2e-16. Repeating many operations can lead to higher relative difference. Here we are computing `d(x,x)²` and a good point for comparison is `||x||²`.
  ``||x||² = 40619662776008`` and `d(x,x) = 0.015625`, which represents a relative error of ~3e-16 which is coherent with my explanation above.

- The distance is not the same when both arrays are the same object compared to when they are not the same object. This is because we treat the 2 cases differently. Indeed, they don't have the same meaning. When both arrays are the same object, we know that the diagonal of the distance matrix should mathematically be 0 and so we enforce that. When they are not the same object we can't have this guarantee and so we report the computed distance which is subject to numerical precision issues.

Labels: Numerical Stability