BUG: FMA contractions/optimisations lead to inconsistent numerical results across different architectures

Describe the issue:

While generation of random numbers using NumPy Generators, we ran into the following problem:

Due to instruction selection differences across compilers, there are discrepancies within NumPy outputs across different architectures. (In the order of 1000s of ULPs on 32-bit architectures as well as linux-aarch64 and linux-ppc64le platforms). This is suspected to be happening due to Fused-Multiply-Add (FMA) optimisations at compiler level.

Linked comment: https://github.com/numba/numba/pull/8038#issuecomment-1165571368

cc @stuartarchibald

Reproduce the code example:

import numpy as np

rng_a = np.random.default_rng(1)
rng_b = np.random.default_rng(1)

scale = np.float64(3.)
loc = np.float64(1.5)
size = (100,)

# The following will undergo fma optimisations on fma enabled NumPy builds/architectures
ra = rng_a.normal(loc=loc, scale=scale, size=size)

# The following won't undergo fma optimizations
rb = loc + scale * rng_b.standard_normal(size=size, dtype=np.float64)

first = True

for (x, y) in zip(ra.flat, rb.flat):
    print('---')
    print(x)
    print(y)
    if first:
        first = False
        # check that the first result is the expected value
        np.testing.assert_array_max_ulp(x, 2.536752576194358, 1)
    np.testing.assert_array_max_ulp(x, y, 1)

Error message:

N.A.

Runtime information:

N.A.

Context for the issue:

No response

Start with NumPy Generators and run the provided comparison between rng.normal and loc + scale * standard_normal on the affected architectures and builds. Read the linked Numba comment about FMA optimisations and trace the normal-generation entry point; done means the reproduced cross-architecture numerical discrepancy is addressed and the example's ULP comparisons pass.