From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gcc-bugzilla@gcc.gnu.org>
Received: by sourceware.org (Postfix, from userid 48)
 id 7A267385840C; Wed, 30 Mar 2022 12:15:18 +0000 (GMT)
DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 7A267385840C
From: "jakub at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>
To: gcc-bugs@gcc.gnu.org
Subject: [Bug libquadmath/105101] incorrect rounding for sqrtq
Date: Wed, 30 Mar 2022 12:15:18 +0000
X-Bugzilla-Reason: CC
X-Bugzilla-Type: changed
X-Bugzilla-Watch-Reason: None
X-Bugzilla-Product: gcc
X-Bugzilla-Component: libquadmath
X-Bugzilla-Version: unknown
X-Bugzilla-Keywords: wrong-code
X-Bugzilla-Severity: normal
X-Bugzilla-Who: jakub at gcc dot gnu.org
X-Bugzilla-Status: NEW
X-Bugzilla-Resolution: 
X-Bugzilla-Priority: P3
X-Bugzilla-Assigned-To: unassigned at gcc dot gnu.org
X-Bugzilla-Target-Milestone: ---
X-Bugzilla-Flags: 
X-Bugzilla-Changed-Fields: cc
Message-ID: <bug-105101-4-kfMlbedA2e@http.gcc.gnu.org/bugzilla/>
In-Reply-To: <bug-105101-4@http.gcc.gnu.org/bugzilla/>
References: <bug-105101-4@http.gcc.gnu.org/bugzilla/>
Content-Type: text/plain; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable
X-Bugzilla-URL: http://gcc.gnu.org/bugzilla/
Auto-Submitted: auto-generated
MIME-Version: 1.0
X-BeenThere: gcc-bugs@gcc.gnu.org
X-Mailman-Version: 2.1.29
Precedence: list
List-Id: Gcc-bugs mailing list <gcc-bugs.gcc.gnu.org>
List-Unsubscribe: <https://gcc.gnu.org/mailman/options/gcc-bugs>,
 <mailto:gcc-bugs-request@gcc.gnu.org?subject=unsubscribe>
List-Archive: <https://gcc.gnu.org/pipermail/gcc-bugs/>
List-Post: <mailto:gcc-bugs@gcc.gnu.org>
List-Help: <mailto:gcc-bugs-request@gcc.gnu.org?subject=help>
List-Subscribe: <https://gcc.gnu.org/mailman/listinfo/gcc-bugs>,
 <mailto:gcc-bugs-request@gcc.gnu.org?subject=subscribe>
X-List-Received-Date: Wed, 30 Mar 2022 12:15:18 -0000

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

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

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

--- Comment #1 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
Trying
#include <math.h>
#include <stdio.h>

int
main ()
{
  volatile long double ld =3D 0x2.4e5658b5f6d2a1ec21e3829dd7f8p+0L;
  ld =3D sqrtl (ld);
  printf ("%.50La\n", ld);
  return 0;
}
on s390x-linux shows:
0x1.84bfedcba255aeaf7e99184e6f3b0000000000000000000000p+0
but that is implemented using a hw instruction.
#include <quadmath.h>
#include <stdio.h>

int
main ()
{
  volatile __float128 ld =3D 0x2.4e5658b5f6d2a1ec21e3829dd7f8p+0Q;
  ld =3D sqrtq (ld);
  char buf[70];
  quadmath_snprintf (buf, 69, "%.50Qa", ld);
  printf ("%s\n", buf);
  return 0;
}
on x86_64 indeed prints ...6f3a
The first testcase also prints ...6f3b on powerpc64le-linux, but again
that doesn't say much, because powerpc64le-linux __ieee754_sqrtf128
comes from sysdeps/powerpc/powerpc64/le/fpu/e_sqrtf128.c where in my case i=
t is
implemented using soft-fp (and when glibc is configured for power9 or later
using hw instruction).

Seems libquadmath sqrtq.c isn't based on any glibc code, instead it uses sq=
rt
or sqrtl for initial approximation and then
    /* Two Newton iterations.  */
    y -=3D 0.5q * (y - x / y);
    y -=3D 0.5q * (y - x / y);
(or for sqrtl one iteration).
I bet that doesn't and can't guarantee correct rounding.
So, do we need to switch to soft-fp based implementation for it (we have a =
copy
already in libgcc/soft-fp), or do we have other options?=