Test results depend on `auto_diff` (thus `sympy`) implementation

When dealing with the formatting of `auto_diff`, I relinted the python code as well and generated new `auto_diff_*.f90` using an up to date python environment.
This fails 3 testcases, see [here](https://testhub.mesastar.org/formatting%2Fauto_diff/commits/ff34670).

given that only the generated `auto_diff_*.f90` files changed, the fails must come from different implementations of the functions as generated by `sympy`. I used `sympy-1.11.1` from `conda-forge`, and would like to know what version built the current `auto_diff_*.f90` files.

subtle changes are _eg_ how `sympy` simplifies expressions:
```
q0 = powm1(sqrt(x%val - 1))*powm1(sqrt(x%val + 1))
```
in my branch (formatting/auto_diff) vs.
```
q0 = powm1(sqrt(pow2(x%val) - 1))
```
in `main`.
These are obviously analytically identical, but numerically not necessarily.

I have tested the suite on a separate branch changing only `auto_diff_real_2var_order3` (used by `skye_eos`) [here](https://testhub.mesastar.org/auto_diff_tests/commits/4c1b79a) and found the same discrepant results as on the formatting branch.

In general it is worrying that our code depends on the implementation of `sympy`, although we do provide the fortran code as-is.
I don't think there is an obvious 'fix' to this issue other than freezing the `sympy` version we use.

Start with the generated auto_diff_*.f90 files and the discrepant test reports linked from the formatting/auto_diff and auto_diff_tests branches. Compare auto_diff_real_2var_order3, including its use by skye_eos, and determine which SymPy version generated the current files. Done means the source of the numerical differences and a reproducible dependency or generation decision are documented.

Written by the indexing model from the issue text.