[Bug c/107370] New: long double precision is wrong in ARM 64

[Bug c/107370] New: long double precision is wrong in ARM 64

[jacob at jacob dot remcomp.fr]

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

            Bug ID: 107370
           Summary: long double precision is wrong in ARM 64
           Product: gcc
           Version: 10.2.1
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: c
          Assignee: unassigned at gcc dot gnu.org
          Reporter: jacob at jacob dot remcomp.fr
  Target Milestone: ---

This program
#include <stdio.h>
#include <math.h>
int main(void) { printf("%.20Lg \n",sqrtl(2.0L)); }
produces the output:
1.4142135623730951455 
The four last digits are WRONG. Correct output is:
1.414213562373095048801688724209698

The long double precision is identical to double precision. It should be approx
30 digits but it is only 16.


gcc -v
MACHINE ARM 64 (raspberry pi)
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib/gcc/aarch64-linux-gnu/10/lto-wrapper
Target: aarch64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Debian 10.2.1-6'
--with-bugurl=file:///usr/share/doc/gcc-10/README.Bugs
--enable-languages=c,ada,c++,go,d,fortran,objc,obj-c++,m2 --prefix=/usr
--with-gcc-major-version-only --program-suffix=-10
--program-prefix=aarch64-linux-gnu- --enable-shared --enable-linker-build-id
--libexecdir=/usr/lib --without-included-gettext --enable-threads=posix
--libdir=/usr/lib --enable-nls --enable-bootstrap --enable-clocale=gnu
--enable-libstdcxx-debug --enable-libstdcxx-time=yes
--with-default-libstdcxx-abi=new --enable-gnu-unique-object
--disable-libquadmath --disable-libquadmath-support --enable-plugin
--enable-default-pie --with-system-zlib --enable-libphobos-checking=release
--with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch
--enable-fix-cortex-a53-843419 --disable-werror --enable-checking=release
--build=aarch64-linux-gnu --host=aarch64-linux-gnu --target=aarch64-linux-gnu
--with-build-config=bootstrap-lto-lean --enable-link-mutex
Thread model: posix
Supported LTO compression algorithms: zlib zstd
gcc version 10.2.1 20210110 (Debian 10.2.1-6) 

uname -a
Linux raspi 5.15.61-v8+ #1579 SMP PREEMPT Fri Aug 26 11:16:44 BST 2022 aarch64
GNU/Linux

Thanks in advance for your attention.

[Bug c/107370] long double precision is wrong in ARM 64

[jacob at jacob dot remcomp.fr]

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

--- Comment #1 from jacob navia <jacob at jacob dot remcomp.fr> ---
This program produces correct results:
#include <stdio.h>
#include <math.h>
long double nsqrt(long double x)
{
        double m = sqrt(x);
        long double r = m;

        r = (r+x/r)/2.0L;

        r = (r+x/r)/2.0L;
        r = (r+x/r)/2.0L;
        r = (r+x/r)/2.0L;
        r = (r+x/r)/2.0L;
        return r;
}

int main(void) { printf("%.30Lg \n",nsqrt(2.0L)); }
produces:
1.41421356237309504880168872421 
All digits que OK.

jacob

[Bug target/107370] long double precision is wrong in ARM 64

[pinskia at gcc dot gnu.org]

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

--- Comment #2 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
I get the following output from ubuntu 1804's gcc and the trunk of gcc from
20220808:

ubuntu@ubuntu:~/src/upstream-gcc\# ~/upstream-gcc/bin/gcc t54.c -O2
ubuntu@ubuntu:~/src/upstream-gcc\# ./a.out
1.4142135623730950488

trunk gcc:
.LC1:
        .word   325511829
        .word   -922176773
        .word   -429395012
        .word   1073703433

ubuntu gcc:
.LC0:
        .word   325511829
        .word   3372790523
        .word   3865572284
        .word   1073703433

Both of these are the same, just unsigned vs signed output.

As far as I can tell, this is not a GCC issue. It might be a binutils issue or
a glibc issue with printf.

[Bug target/107370] long double precision is wrong in ARM 64

[jacob at jacob dot remcomp.fr]

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

--- Comment #3 from jacob navia <jacob at jacob dot remcomp.fr> ---
  1 trunk gcc:
  2 .LC1:
  3        .word   325511829  # 0x1366EA95 <<<--- SHOULD BE 325508205
  4        .word   -922176773 # 0xC908B2FB OK
  5        .word   -429395012 # 0xE667F3BC OK
  6        .word   1073703433 # 0x3FFF6A09 OK

This data is wrong, I repeat, the first number (line 3) should be 325508205 or 
0x1366DC6D.

It CAN'T be a printf issue because when you use my "nsqrt" function IT PRINTS
CORRECTLY. 

THEN

printf is able to print correctly when the input is correct. The sqrtl function
under ARM 64 bits produces only double precision, instead of quadruple
precision.

That's all I am saying. 

And the intention here (from my side) is to help you correct this problem, not
to minimize the quality of gcc, a great software!

jacob

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

[pinskia at gcc dot gnu.org]

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

Andrew Pinski <pinskia at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |NEW
           Keywords|                            |wrong-code
     Ever confirmed|0                           |1
          Component|target                      |middle-end
            Summary|long double precision is    |long double sqrtl constant
                   |wrong in ARM 64             |folding is wrong
   Last reconfirmed|                            |2022-10-24

--- Comment #4 from Andrew Pinski <pinskia at gcc dot gnu.org> ---

Hmm, this might be a bug in mpfr though.  GCC uses MPFR to do constant folding.


Without constant folding libm produces:
1.4142135623730950488

That is:
```
#include <stdio.h>
#include <math.h>
int main(void) { volatile long double a = 2.0L; printf("%.20Lg \n",sqrtl(a)); }
```

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

[jakub at gcc dot gnu.org]

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

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jakub at gcc dot gnu.org

--- Comment #5 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
On powerpc64le-linux with -O2 -nostdinc -mlong-double-128 -mabi=ieeelongdouble
int main() { __builtin_printf("%.20Lg \n",__builtin_sqrtl(2.0L)); }
we also emit on the trunk:
        .long   325511829
        .long   -922176773
        .long   -429395012
        .long   1073703433
And similarly for:
long double foo (long double x) { return __builtin_hypotl (x, x); }
with -O2 -nostdinc -mlong-double-128 -mabi=ieeelongdouble -ffast-math where we
optimize the hypotl to fabs(x) * M_PI2l (dconst_sqrt2) it is the same constant.
const REAL_VALUE_TYPE *
dconst_sqrt2_ptr (void)
{
  static REAL_VALUE_TYPE value;

  /* Initialize mathematical constants for constant folding builtins.
     These constants need to be given to at least 160 bits precision.  */
  if (value.cl == rvc_zero)
    {
      mpfr_t m;
      mpfr_init2 (m, SIGNIFICAND_BITS);
      mpfr_sqrt_ui (m, 2, MPFR_RNDN);
      real_from_mpfr (&value, m, NULL_TREE, MPFR_RNDN);
      mpfr_clear (m);
    }
  return &value;
}

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #6 from jacob navia <jacob at jacob dot remcomp.fr> ---
 mpfr_init2 (m, SIGNIFICAND_BITS);

How much is that SIGNIFICAND_BITS? I guess that is the problem.

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #7 from jacob navia <jacob at jacob dot remcomp.fr> ---
 mpfr_init2 (m, SIGNIFICAND_BITS);

How much is that SIGNIFICAND_BITS? I guess that is the problem.

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

[joseph at codesourcery dot com]

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

--- Comment #8 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
On Mon, 24 Oct 2022, jacob at jacob dot remcomp.fr via Gcc-bugs wrote:

> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107370
> 
> --- Comment #3 from jacob navia <jacob at jacob dot remcomp.fr> ---
>   1 trunk gcc:
>   2 .LC1:
>   3        .word   325511829  # 0x1366EA95 <<<--- SHOULD BE 325508205
>   4        .word   -922176773 # 0xC908B2FB OK
>   5        .word   -429395012 # 0xE667F3BC OK
>   6        .word   1073703433 # 0x3FFF6A09 OK
> 
> This data is wrong, I repeat, the first number (line 3) should be 325508205 or 
> 0x1366DC6D.

Why do you think that number is wrong?  If I compute the square root of 
2**225 using GMP (so not involving MPFR at all, just integer square root 
in GMP), I get 1366ea95 as the low 32 bits (and the next bit is a 0, so 
rounding toward 0 is correct in this case).

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #9 from jacob navia <jacob at jacob dot remcomp.fr> ---
Actually, this looks that a bug in all versions of gcc!
Using the "GP" calculator I obtain:
1.414213562373095048801688724209698078
Using the "bc" calculator

~ % 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=33
sqrt(2)
1.414213562373095048801688724209698

All agree: gcc is just wrong.

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #10 from jacob navia <jacob at jacob dot remcomp.fr> ---
Please let's use 30 digits after the decimal point not 20. Long double gives 32
digits with 128 bits.

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #11 from jacob navia <jacob at jacob dot remcomp.fr> ---
Please let's use 30 digits after the decimal point not 20. Long double gives 32
digits with 128 bits.

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

[jakub at gcc dot gnu.org]

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;
  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.

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

[jakub at gcc dot gnu.org]

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

--- Comment #13 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
glibc math.h has:
# define M_SQRT2f128    __f128 (1.414213562373095048801688724209698079) /*
sqrt(2) */

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

[jakub at gcc dot gnu.org]

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

--- Comment #14 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
(In reply to jacob navia from comment #7)
>  mpfr_init2 (m, SIGNIFICAND_BITS);
> 
> How much is that SIGNIFICAND_BITS? I guess that is the problem.

160 for 32-bit hosts, 192 for 64-bit hosts.

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

[jacob at jacob dot remcomp.fr]

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

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #16 from jacob navia <jacob at jacob dot remcomp.fr> ---
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

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #17 from jacob navia <jacob at jacob dot remcomp.fr> ---
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


> Le 24 oct. 2022 à 22:50, gcc-bugzilla@gcc.gnu.org a écrit :
> 
> Attachments with a MIME type of "text/html" are not allowed on this installation. 
> 
> jacob navia wrote:
>> 
>> 
>>> 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

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #18 from jacob navia <jacob at jacob dot remcomp.fr> ---
I can’t reply you because you do not accept accented characters, just plain
ascii, and my mailer uses a French
Accentuated character to say « You wrote » and the whole message is translated
into mime/html.

I posted my reply in the bug mailer 
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107370
<https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107370>

Jacob

> Le 24 oct. 2022 à 22:50, gcc-bugzilla@gcc.gnu.org a écrit :
> 
> Attachments with a MIME type of "text/html" are not allowed on this installation. 
> 
> jacob navia wrote:
>> 
>> 
>>> 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

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #19 from jacob navia <jacob at jacob dot remcomp.fr> ---
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

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

[jakub at gcc dot gnu.org]

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

--- Comment #20 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
All I'm arguing is that everything I see from gcc seems to be correct.
30 digits precision is not enough for IEEE quad, 36 is enough.
If I try your sqrt (2^255) case on x86_64, I see:
#define _GNU_SOURCE
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int
main ()
{
  volatile _Float128 x = __builtin_sqrtf128(powf128(2, 225));
  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);
}

0x1.6a09e667f3bcc908b2fb1366ea95p+112  
7343016637207168931428032607349397.000000000000000000000000000000000000
0x1.ffffffffffffffffffffffffffffp+224  
53919893334301279589334030174039256154977430310253516431710891278336.000000000000000000000000000000000000
so again, x * x is 1ulp from 0x1.0p+225.

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #21 from jacob navia <jacob at jacob dot remcomp.fr> ---
When I post I come to a page

MIDAIR COLLISION DETECTED

and that provokes these repeated messages

Sorry but there is a bug in the bug reporting software...

Gosh!

:-)

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #22 from jacob navia <jacob at jacob dot remcomp.fr> ---
(In reply to Jakub Jelinek from comment #20)
> All I'm arguing is that everything I see from gcc seems to be correct.
> 30 digits precision is not enough for IEEE quad, 36 is enough.
> If I try your sqrt (2^255) case on x86_64, I see:
> #define _GNU_SOURCE
> #include <stdlib.h>
> #include <stdio.h>
> #include <math.h>
> 
> int
> main ()
> {
>   volatile _Float128 x = __builtin_sqrtf128(powf128(2, 225));
>   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);
> }
> 
> 0x1.6a09e667f3bcc908b2fb1366ea95p+112
> 7343016637207168931428032607349397.000000000000000000000000000000000000
> 0x1.ffffffffffffffffffffffffffffp+224
> 53919893334301279589334030174039256154977430310253516431710891278336.
> 000000000000000000000000000000000000
> so again, x * x is 1ulp from 0x1.0p+225.

OK, but you are NOT RUNNING IN AN ARM 64 do you? Because that's the machine I
am using...

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

[jakub at gcc dot gnu.org]

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

--- Comment #23 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
> OK, but you are NOT RUNNING IN AN ARM 64 do you? Because that's the machine
> I am using...

So that would point to glibc's printf then.
Try to use %a to see the hex values, those are much more easily verifiable (you
can e.g. look at the bytes in memory in a debugger etc.).

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

[jacob at jacob dot remcomp.fr]

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

--- Comment #24 from jacob navia <jacob at jacob dot remcomp.fr> ---
Sorry, you are right. This is a terrible mistake FROM MY PART!

I am using the raspberry pi through ssh and I had several shells open, using my
Macintosh ARM 64 as the machine where I started those shells.

But I got confused with the Macintosh shells and thought that I was using the
raspberry pi when I was using the Macintosh instead.

The macintosh compiler also is called "gcc" but actually is calling clang, that
doesn't support long doubles of 128 bits but only 64 bits, what explains the
loss of precision after the 15th bits!!!

The bug was in MY SIDE and I have taken your time for nothing. I am very sorry
about this situation. I do not know how to repair it.

My sincere excuses.

jacob

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

[pinskia at gcc dot gnu.org]

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

Andrew Pinski <pinskia at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
         Resolution|---                         |INVALID
             Status|NEW                         |RESOLVED

--- Comment #25 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
Note the darwin (mac OS X) aarch64 (arm64) ABI says long double is 64bits. It
would be implemented that way even in GCC for Mac OS X ARM64.
clang/llvm for aarch64 other than Mac OS X does support 128bit long double
though as required by the ABI too.

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

[pinskia at gcc dot gnu.org]

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

--- Comment #26 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
(In reply to Andrew Pinski from comment #25)
> Note the darwin (mac OS X) aarch64 (arm64) ABI says long double is 64bits.
> It would be implemented that way even in GCC for Mac OS X ARM64.
> clang/llvm for aarch64 other than Mac OS X does support 128bit long double
> though as required by the ABI too.

PR 107315 records the adding of support of 128bit float for Darwin (as
_Float128).