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:
If this problem is impacting many users, we might consider investing efforts in implementing alternative metric (e.g. "stable_euclidean") that would not rely on the GEMM trick. However to keep good performance, this would probably require some kind of manual SIMD tooling, either based on the numpy SIMD infrastructure if it is made available as a library or on XSIMD (see #26010).

Labels: Numerical Stability