No, this is just literally a property of chaos. Changes of 1e-16 will give O(1) differences in the trajectory. The difference of using muladd vs x*y + z is then enough to, over time after a Lyopunov time, give O(1) differences. You can only get exactly the same result out of chaotic systems if you have exactly, floating point exact, the same implementation.
Take a look at https://diffeqcallbacks.sciml.ai/stable/uncertainty_quantification/ to explore the uncertainty of such systems.
Everything that approximates a chaotic attractor is incorrect. You only have the shadowing lemma: the fact that there’s a trajectory with an epsilon different initial condition on the attractor that you happen to have pulled, but it can be O(1) different from the real trajectory. Numerical approximations to chaotic systems are only trustworthy in a probabilistic sense.
Here’s a nice post that shows some of the details: