Improve precision of 32-bit `gammaln`?

Description

In [2]: from jax.scipy.special import gammaln as gammaln_jax
In [5]: x = jax.numpy.asarray(2.00001)
In [7]: gammaln_jax(x)
Out[7]: Array(5.722046e-06, dtype=float32, weak_type=True)
...
In [1]: from scipy.special import gammaln
In [4]: x = np.asarray(2.00001)
In [5]: gammaln(x)
Out[5]: 4.227875597648359e-06

https://www.wolframalpha.com/input?i=ln%28Gamma%282.00001%29%29

Found in https://github.com/scipy/scipy/pull/20085#issuecomment-2117810225.

System info (python version, jaxlib version, accelerator, etc.)

jax:    0.4.27
jaxlib: 0.4.23.dev20240502
numpy:  1.26.4
python: 3.12.3 | packaged by conda-forge | (main, Apr 15 2024, 18:35:20) [Clang 16.0.6 ]
jax.devices (1 total, 1 local): [CpuDevice(id=0)]
process_count: 1
platform: uname_result(system='Darwin', node='Lucass-MacBook-Air-4.local', release='23.4.0', version='Darwin Kernel Version 23.4.0: Fri Mar 15 00:12:41 PDT 2024; root:xnu-10063.101.17~1/RELEASE_ARM64_T8103', machine='arm64')