The data type of input data for LinearRegression class will affect the results

### Describe the bug

Our team just used the class sklearn.linear_model.LinearRegression to do multi-linear regression. And, we found out that the same data, which means the values are identical for each element, with different data formats will return different coefficients. The array we use are all integer values but saved as float32 and float64, separately. However, the array in float64 as input get us the expected result, but float32 doesn't.

The code below is a simple reproducible example.

### Steps/Code to Reproduce

```python
In [1]: from sklearn.linear_model import LinearRegression
    ...: import numpy as np
    ...: rng = np.random.default_rng(0)
    ...: x_1_array = rng.integers(low=1000, high=10000, size=10_000_000)
    ...: x_2_array = rng.integers(low=1000, high=10000, size=10_000_000)
    ...: x_3_array = rng.integers(low=1000, high=10000, size=10_000_000)
    ...: y_array = 3 * x_1_array + x_2_array - 5 * x_3_array
    ...: intercept_data = rng.integers(low=-3, high=3)
    ...: print(f'True intercept {intercept_data}')
    ...: y_array += intercept_data
    ...: X_train = np.stack([x_1_array,x_2_array,x_3_array],axis=-1)
    ...: y_train = y_array
    ...: linreg = LinearRegression()
    ...: model = linreg.fit(X_train.astype(np.float64), y_train)
    ...: print(f'Float 64 Intercept: {linreg.intercept_}')
    ...: print(f'Float 64 Coefficient: {linreg.coef_}')
    ...: model = linreg.fit(X_train.astype(np.float32), y_train)
    ...: print(f'Float 32 Intercept: {linreg.intercept_}')
    ...: print(f'Float 32 Coefficient: {linreg.coef_}')
```

### Expected Results

![image](https://user-images.githubusercontent.com/51532833/168492124-d98f9663-4899-47ee-83b2-4952462c37a6.png)

### Actual Results

![image](https://user-images.githubusercontent.com/51532833/168492136-4e0f4100-9e48-4f58-aa2a-97e826f06c61.png)


### Versions

```shell
System:
    python: 3.7.9 (v3.7.9:13c94747c7, Aug 15 2020, 01:31:08)  [Clang 6.0 (clang-600.0.57)]
executable: /Users/lazy/PycharmProjects/projects/bin/python
   machine: Darwin-21.4.0-x86_64-i386-64bit

Python dependencies:
          pip: 21.3
   setuptools: 58.2.0
      sklearn: 1.0.2
        numpy: 1.21.4
        scipy: 1.7.2
       Cython: None
       pandas: 1.3.4
   matplotlib: 3.4.3
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True
```

COMMENT:
Hello everyone 

I’m a beginner in contributing, and I’ve been following this issue for a while.  
I totally agree that this is mainly a numerical precision limitation rather than a direct bug.

I was thinking of contributing in two small ways:

1. **Documentation improvement** — adding a short note under the *Ordinary Least Squares* section explaining that using `float32` inputs may lead to small differences in coefficients or intercepts compared to `float64`, with a link to this issue and a short mitigation tip (`float64` or `StandardScaler`).

2. **Runtime warning** — adding a `UserWarning` inside `LinearRegression.fit()` when `X.dtype == np.float32`, to gently alert users without changing model behavior.

I understand this won’t “fix” the underlying mathematical limitation, but it could help users better understand and avoid confusion when they see slightly different results.

If this sounds good, could you please assign this issue to me?  
I’d be happy to start working on it

Labels: Bug, module:linear_model, float32