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:
I didn't pay enough attention to the second case "final element and copy". I agree that this is more surprising.

In both the first and second case, the arrays are not the same object (``numbers2[1:] is not numbers2[1:]``). So they are treated the same by ``euclidean_distances``.

However, in the first case, both views point to the same place in the memory
``np.may_share_memory(numbers2[1:], numbers2[1:])``.

This is what makes the matrix matrix multiplication ``X @ Y.T`` behave differently. This operation is done by the underlying BLAS library (OpenBLAS, MKL, ...) which must have a slightly different behavior when the pointers point to the same address. The relative difference between the results of ``X @ Y.T`` is ~2e-16, which is expected if BLAS does not perform all operations in the exact same order depending on whether the pointers are the same or not. So in the end the explanation is still rounding errors but in a different place :)

Labels: Numerical Stability