[gcc r13-9604] libquadmath: Fix up THREEp96 constant in expq

[Jakub Jelinek]

https://gcc.gnu.org/g:b4988146684a9813aab174aedf82704b42f2a0db

commit r13-9604-gb4988146684a9813aab174aedf82704b42f2a0db
Author: Jakub Jelinek <jakub@redhat.com>
Date:   Wed Apr 9 22:01:30 2025 +0200

    libquadmath: Fix up THREEp96 constant in expq
    
    Here is a cherry-pick from glibc [BZ #32411] fix.
    
    As mentioned by the reporter in a pull request against gcc-mirror,
    the THREEp96 constant in e_expl.c is incorrect, it is actually 0x3.p+94f128
    rather than 0x3.p+96f128.
    
    The algorithm uses that to compute the t2 integer (tval2), by whose
    delta it adjusts the x+xl pair and then in the result uses the precomputed
    exp value for that entry.
    Using 0x3.p+94f128 rather than 0x3.p+96f128 results in tval2 sometimes
    being one smaller, sometimes one larger than the desired value, thus can mean
    the x+xl pair after adjustment will be larger in absolute value than it
    should be.
    
    DesWursters created a test program for this
    https://github.com/DesWurstes/comparefloats
    and his results were
    total: 1135000000 not_equal: 4322 earlier_score: 674 later_score: 3648
    I've modified this so with
    https://sourceware.org/bugzilla/show_bug.cgi?id=32411#c3
    so that it actually tests pseudo-random _Float128 values with range
    (-16384.,16384) with strong bias on values larger than 0.0002 in absolute
    value (so that tval1/tval2 aren't zero most of the time) and that gave
    total: 10000000000 not_equal: 29861 earlier_score: 4606 later_score: 25255
    So, in both cases, in most cases the change doesn't result in any differences,
    and in those rare cases where does, about 85% have smaller ulp than without
    the patch.
    Additionally I've tried
    https://sourceware.org/bugzilla/show_bug.cgi?id=32411#c4
    and in 2 billion iterations it didn't find any case where x+xl after the
    adjustments without this change would be smaller in absolute value compared
    to x+xl after the adjustments with this change.
    
    2025-04-09  Jakub Jelinek  <jakub@redhat.com>
    
            * math/expq.c (C): Fix up THREEp96 constant.
    
    (cherry picked from commit e081ced345c45581a4891361c08e50e07720239e)

Diff:
---
 libquadmath/math/expq.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/libquadmath/math/expq.c b/libquadmath/math/expq.c
index 1aaad91c2446..171598824f12 100644
--- a/libquadmath/math/expq.c
+++ b/libquadmath/math/expq.c
@@ -74,7 +74,7 @@ static const __float128 C[] = {
 
 /* 3x2^96 */
 #define THREEp96 C[2]
- 59421121885698253195157962752.0Q,
+ 237684487542793012780631851008.0Q,
 
 /* 3x2^103 */
 #define THREEp103 C[3]