Also: add basic tests to detect fpu miscompilation.
The program sets the fpu to 53 bit mode to get consistent results. It
turns out that doing this on musl breaks fmt_fp (-> printf) which uses
long double to format floats. Note this is an issue even on 64 bit musl,
although setting fpu to 53 bits is unnecessary on 64 bits because of the
default -mfpmath=sse.
The way Configure works: it tries different combinations of flags which
exercise different methods to change the fpu control word until one
works, meaning doubles are effectively 53 bits. The fix here is to try
first NOT touching the fpu control word. On x86_64 using the default
-mfpmath=sse this will succeed bypassing all fpucw modification.
For i686 using -mfpmath=387 (always if sse2 not available) the code will
fall back to changing the fpu control word as before. Since this is not
a problem for glibc, everything still works. I expect i686-musl to still
be broken, but we don't support that arch.
A simple way to test the bug on musl is to run
$ sympow -curve '[0,0,0,0,1]' -analrank
...
Analytic Rank is 0 : L-value 7.01091e-01
When the fpucw is changed, on musl this prints "7.01092e-01" instead.
The actual value computed to 128 bits with pari is:
$ echo 'elllseries(ellinit([0,0,0,0,1]),1)' | gp -q
0.70109105266272713058750953952514706773
so the glibc output is the correct one.
Again: what is broken is not computing but printing, as fmt_fp in musl
uses long double, which means messing with fpu control word breaks it.
This is important since sagemath parses sympow output.
This affects sagemath doctests as in
$ sage -t src/sage/lfunctions/sympow.py
...
(3 failures)
After this commit, the doctest in question passes on the 3 supported
archs.
Érico Nogueira
Omar Zeghouani
674Y3r
mobinmob
Duncaen
DaErich
Andrew J. Hesford
Leah Neukirchen
Đoàn Trần Công Danh
Gonzalo Tornaría
Michal Vasilek
meator
Clyhtsuriva
KawaiiAmber
Ramdziana F Y
Daniel Eyßer
Daniel Lewan
Ján Priner
Frank Steinborn
Ranme Fritzvilliams