[Bug middle-end/107370] long double sqrtl constant folding is wrong

["jacob at jacob dot remcomp.fr" <gcc-bugzilla@gcc.gnu.org>]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107370

--- Comment #15 from jacob navia <jacob at jacob dot remcomp.fr> ---
> Le 24 oct. 2022 à 22:21, jakub at gcc dot gnu.org <gcc-bugzilla@gcc.gnu.org> a écrit :
> 
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107370
> 
> --- Comment #12 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
> I have to agree with Joseph (on x86_64-linux I get):
> #define _GNU_SOURCE
> #include <stdlib.h>
> #include <stdio.h>
> 
> int
> main ()
> {
>  volatile _Float128 x = __builtin_sqrtf128(2.0f128);
>  char buf[256];
>  strfromf128 (buf, 128, "%a", x);
>  strfromf128 (buf + 128, 128, "%.36f", x);
>  printf ("%s\t%s\n", buf, buf + 128);
>  x = x * x;
>  strfromf128 (buf, 128, "%a", x);
>  strfromf128 (buf + 128, 128, "%.36f", x);
>  printf ("%s\t%s\n", buf, buf + 128);
>  // The value you want to see
>  x = 0x1.6a09e667f3bcc908b2fb1366dc6dp+0f128;
>  strfromf128 (buf, 128, "%a", x);
>  strfromf128 (buf + 128, 128, "%.36f", x);
>  printf ("%s\t%s\n", buf, buf + 128);
>  x = x * x; places.

>  strfromf128 (buf, 128, "%a", x);
>  strfromf128 (buf + 128, 128, "%.36f", x);
>  printf ("%s\t%s\n", buf, buf + 128);
> }
> 
> 0x1.6a09e667f3bcc908b2fb1366ea95p+0     1.414213562373095048801688724209697984
> 0x1.ffffffffffffffffffffffffffffp+0     1.999999999999999999999999999999999807
> 0x1.6a09e667f3bcc908b2fb1366dc6dp+0     1.414213562373095048801688724209000027
> 0x1.ffffffffffffffffffffffffd7f4p+0     1.999999999999999999999999999998025537
> 
> Clearly the value GCC uses for sqrtf128(2.0f128) multiplied by itself is 1ulp
> smaller than 2.0,
> while the value you want to see is 10252ulps smaller than 2.0.
> 
> -- 
> You are receiving this mail because:
> You reported the bug.

Sorry but the program I used was compiled with gcc. It suffers from lack of
precision too.

In ANY case, in decimal now, the square root of 2 is
1.414 213 562 373 095 048 801 688 724 209 698 078 GP PARI calculator
1.414 213 562 373 095 145 547 462 185 874 gcc: only the first 15 digits are OK
                      ^^^^^ 
In ANOTHER example:
#include <stdio.h>
#include <math.h>

int main(void) { printf("%.30Lg \n",sqrtl(powl(2,225))); }

This produces:
7.34301663720716943338259962711e+33
But bc produces 
~/lcc/book/test % bc
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
scale=32
sqrt(2^225)
7343016637207168931428032607349397.48923080385562326502615381523522

To make things clearer I put now the results side by side
7.34301663720716943338259962711e+33 gcc
734301663720716 8931428032607349397.48923080385562326502615381523522 bc

As you can see, after the 15th digits all others are wrong!!!

Jacob