Make expm1 detection work correctly

[heplesser]

In many of our models, we use expm1() to compute propagator elements. Since expm1() historically was not guaranteed to be available in math.h or cmath, we let CMake check if expm1() was available in math.h and if not, we provided our own implementation.

It now turns out that this CMake mechanism was incorrectly implemented: The detection worked correctly, but HAVE_EXPM1 was never set in config.h due to an extra argument to the pertaining CMake function that got things messed up. This PR fixes this, so that HAVE_EXPM1 is now properly set.

One could wonder whether we need the check at all, given that according to CPPReference is part of the standard; but a Perplexity query gave a more ambiguous result.

Now the problem here is that this fix will change NEST results, because the expm1() function from C++ may return different results from our implementation, which has been used in all NEST experiments over at least the past decade due to the bug described and fixed here.

To test the consequences, I performed simulations with the brunel_alpha_nest.py example and a modified version using the precise iaf_psc_alpha_ps neuron. For the normal iaf_psc_alpha neuron, I did not find a single spike differed when recording from 500 spikes for 10s. But for the precise neurons, already the time of the first spikes differed by about 10 eps with "macroscopic" differences in spike patterns later in the simulation. The precise neurons are a much harder test since they evaluate expm1() essentially for every incoming spike and thus also probe a much wider range of argument values, especially values close to 0. The implementation from math.h, btw, gives about 25% faster network simulation times for the ps-case.

How do we proceed here?

[heplesser]

I have thought about this a bit more especially in light of

Blanco, W., Lopes, P. H., De S. Souza, A. A., & Mascagni, M. (2020).
Non-replicability circumstances in a neural network model with Hodgkin-Huxley-type neurons.
Journal of Computational Neuroscience, 48(3), 357–363. https://doi.org/10.1007/s10827-020-00748-3

They showed that the exp() function can return different results on different computers (IEEE754 only applies to simple arithmetic). So it is reasonable to assume that also expm1() and all other higher functions can return system-dependent results. Making a change that is within the range of this system-dependent variation should be acceptable.

But we might want to give the user a choice to force the old behavior with a compiler switch. I am just not sure about the name, maybe --Dwith-classic-expm1 or -Dwith-nest-expm1?

[diesmann]

Nice finding HEP. I think if expm1 is in the standard we should normally use this, as you say. Still, also as you say, having the opportunity to reproduce our old results in special situations is of high value. So a compiler switch is a good option. I think with-classic-expm1 is the better choice because with-nest-expm1 suggests a bit that this is extra good.

[heplesser]

Nice finding HEP. I think if expm1 is in the standard we should normally use this, as you say. Still, also as you say, having the opportunity to reproduce our old results in special situations is of high value. So a compiler switch is a good option. I think with-classic-expm1 is the better choice because with-nest-expm1 suggests a bit that this is extra good.

Good point, but I am not really happy with "classic" either—what does it really mean? Any other ideas?

[heplesser]

Explored some more: The differences between our handcrafted solution and std::expm1 are of the same order as differences in std::expm1 between Clang and GCC when compiling with -O3 (which may mean that arithmetic is not IEEE754-conformant).

[heplesser]

I did some more experiments (with some help from duck.ai/Haiku 4.5) as follows:

• Compute expm1() using both our code and std::expm1() for logarithmically spaced x with $10^{-25} < |x| < 10$ for $n=10^7$ values.
• Compute the difference between our and the standard result.
• Do this both with Apple Clang 21 and gcc 15, building with -O3

Results show the following:

• For about 24% of values tested, our implementation and std differ
• Relative differences never exceed $4 \epsilon$ in either direction
• The mean of the absolute relative differences are around $0.75\epsilon$
• Comparing results obtained with clang with results obtained with gcc, I find
  • About 8000 out of $10^7$ values of $x$ differ when executing the same function compiled with different compilers
  • Relative differences are between $-3\epsilon$ and $2.5\epsilon$
  • The mean of the absolute relative differences for the cases with non-zero differences is $1\epsilon$ for our code and $0.5\epsilon$ for std::expm1().
• Results are very similar when compiling while suppressing agressive math optimizations (clang: -O1 -fno-fast-math -ffp-contract=off -fexcess-precision=standard, gcc: -O1 -fno-fast-math -fexcess-precision=standard -frounding-math).

The difference between different compilers are thus of the same order as between the different implementations. All my data are obtained on Apple Silicon, Intel-ish CPUs with 80 bit floating point registers may add more variation in results for the same operations.

Therefore, I think that it is justified to simply correct the detection mechanism in CMake without adding an extra switch. If someone really wants to test things, they can modify config.h.

I will modernize the CMake mechanism a little because the function currently used there is deprecated. I will also switch from math.h to cmath.