Re: [PATCH] libquadmath: add quad support for trig-pi functions
On Thu, Jul 03, 2025 at 02:43:43PM +0200, Michael Matz wrote:
> Hello,
>
> On Thu, 3 Jul 2025, Yuao Ma wrote:
>
> > This patch adds the required function for Fortran trigonometric functions to
> > work with glibc versions prior to 2.26. It's based on glibc source commit
> > 632d895f3e5d98162f77b9c3c1da4ec19968b671.
> >
> > I've built it successfully on my end. Documentation is also included.
> >
> > Please take a look when you have a moment.
>
> +__float128
> +cospiq (__float128 x)
> +{
> ...
> + if (__builtin_islessequal (x, 0.25Q))
> +return cosq (M_PIq * x);
>
> Isn't the whole raison d'etre for the trig-pi functions that the internal
> argument reduction against multiples of pi becomes trivial and hence (a)
> performant, and (b) doesn't introduce rounding artifacts? Expressing the
> trig-pi functions in terms of their counterparts completely defeats this
> purpose. The other way around would be more sensible for the cases where
> it works, but the above doesn't seem very attractive.
glibc has
FLOAT
M_DECL_FUNC (__cospi) (FLOAT x)
{
if (isless (M_FABS (x), M_EPSILON))
return M_LIT (1.0);
if (__glibc_unlikely (isinf (x)))
__set_errno (EDOM);
x = M_FABS (x - M_LIT (2.0) * M_SUF (round) (M_LIT (0.5) * x));
if (islessequal (x, M_LIT (0.25)))
return M_SUF (__cos) (M_MLIT (M_PI) * x);
else if (x == M_LIT (0.5))
return M_LIT (0.0);
else if (islessequal (x, M_LIT (0.75)))
return M_SUF (__sin) (M_MLIT (M_PI) * (M_LIT (0.5) - x));
else
return -M_SUF (__cos) (M_MLIT (M_PI) * (M_LIT (1.0) - x));
}
for this case, so if it is incorrect, has too bad precision
or there are better ways to do it, then perhaps it should be changed
on the glibc side first and then ported to libquadmath.
Jakub
Re: [PATCH] libquadmath: add quad support for trig-pi functions
On Thu, 3 Jul 2025, Jakub Jelinek wrote: > > Isn't the whole raison d'etre for the trig-pi functions that the internal > > argument reduction against multiples of pi becomes trivial and hence (a) > > performant, and (b) doesn't introduce rounding artifacts? Expressing the > > trig-pi functions in terms of their counterparts completely defeats this > > purpose. The other way around would be more sensible for the cases where > > it works, but the above doesn't seem very attractive. > x = M_FABS (x - M_LIT (2.0) * M_SUF (round) (M_LIT (0.5) * x)); In particular, this is what trivial range reduction looks like: no need to do multiple-precision multiplication with the relevant bits of a multiple-precision value of 1/pi, just round to the nearest integer (typically a single instruction). -- Joseph S. Myers josmy...@redhat.com
Re: [PATCH] libquadmath: add quad support for trig-pi functions
Hi all,
On 7/3/2025 9:21 PM, Jakub Jelinek wrote:
On Thu, Jul 03, 2025 at 02:43:43PM +0200, Michael Matz wrote:
Hello,
On Thu, 3 Jul 2025, Yuao Ma wrote:
This patch adds the required function for Fortran trigonometric functions to
work with glibc versions prior to 2.26. It's based on glibc source commit
632d895f3e5d98162f77b9c3c1da4ec19968b671.
I've built it successfully on my end. Documentation is also included.
Please take a look when you have a moment.
+__float128
+cospiq (__float128 x)
+{
...
+ if (__builtin_islessequal (x, 0.25Q))
+return cosq (M_PIq * x);
Isn't the whole raison d'etre for the trig-pi functions that the internal
argument reduction against multiples of pi becomes trivial and hence (a)
performant, and (b) doesn't introduce rounding artifacts? Expressing the
trig-pi functions in terms of their counterparts completely defeats this
purpose. The other way around would be more sensible for the cases where
it works, but the above doesn't seem very attractive.
glibc has
FLOAT
M_DECL_FUNC (__cospi) (FLOAT x)
{
if (isless (M_FABS (x), M_EPSILON))
return M_LIT (1.0);
if (__glibc_unlikely (isinf (x)))
__set_errno (EDOM);
x = M_FABS (x - M_LIT (2.0) * M_SUF (round) (M_LIT (0.5) * x));
if (islessequal (x, M_LIT (0.25)))
return M_SUF (__cos) (M_MLIT (M_PI) * x);
else if (x == M_LIT (0.5))
return M_LIT (0.0);
else if (islessequal (x, M_LIT (0.75)))
return M_SUF (__sin) (M_MLIT (M_PI) * (M_LIT (0.5) - x));
else
return -M_SUF (__cos) (M_MLIT (M_PI) * (M_LIT (1.0) - x));
}
for this case, so if it is incorrect, has too bad precision
or there are better ways to do it, then perhaps it should be changed
on the glibc side first and then ported to libquadmath.
That's exactly what I want to say. I only touched the generation script;
the other main changes are borrowed from glibc. Therefore, it should be
at least as precise as the glibc generic implementation. If the current
implementation isn't sufficient, we can update glibc first and then
mirror the changes back. This is only used for systems with glibc
versions prior to 2.26. Modern systems will use the implementation
directly from glibc, which should avoid this potential issue.
Regarding math function implementation, I think LLVM-libc provides a
great example. While it may not be as performant as its glibc
equivalent, its structure is very clear. Even for a newcomer like me,
who isn't very familiar with computational mathematics, I can grasp the
basic ideas of implementation and testing, such as nan/inf handling,
polynomial curve fitting, exception values, and range reduction. I'd
love to learn from and contribute to the glibc generic implementation as
Joseph suggested, but it seems it will take more time to get familiar
with the glibc source code. Hopefully, I'll be able to do that one day.
Regarding the patch itself, the previous version was missing the
declaration in the header file. The newly attached version fixes this.
Is this patch good to go, or does it need further amendment? Looking
forward to your further review comments.
Thanks,
Yuao
From 719886689fd21555c9ba1bd63a44111408c0f01d Mon Sep 17 00:00:00 2001
From: Yuao Ma
Date: Thu, 3 Jul 2025 23:01:38 +0800
Subject: [PATCH] libquadmath: add quad support for trig-pi functions
This function is required for Fortran trigonometric functions with glibc <2.26.
Use glibc commit 632d895f3e5d98162f77b9c3c1da4ec19968b671.
libquadmath/ChangeLog:
* Makefile.am: Add sources to makefile.
* Makefile.in: Regen makefile.
* libquadmath.texi: Add doc for trig-pi funcs.
* quadmath.h (acospiq, asinpiq, atanpiq, atan2piq, cospiq, sinpiq,
tanpiq): New.
* update-quadmath.py: Update generation script.
* math/acospiq.c: New file.
* math/asinpiq.c: New file.
* math/atan2piq.c: New file.
* math/atanpiq.c: New file.
* math/cospiq.c: New file.
* math/sinpiq.c: New file.
* math/tanpiq.c: New file.
Signed-off-by: Yuao Ma
---
libquadmath/Makefile.am| 2 ++
libquadmath/Makefile.in| 26 --
libquadmath/libquadmath.texi | 7
libquadmath/math/acospiq.c | 33 ++
libquadmath/math/asinpiq.c | 40 ++
libquadmath/math/atan2piq.c| 36
libquadmath/math/atanpiq.c | 35 +++
libquadmath/math/cospiq.c | 37
libquadmath/math/sinpiq.c | 44
libquadmath/math/tanpiq.c | 62 ++
libquadmath/quadmath.h | 7
libquadmath/update-quadmath.py | 48 ++
12 files changed, 359 insertions(+), 18 deletions(-)
create mode 100644 libquadmath/math/acospiq.c
create mode 100644 libquadmath/math/asinpiq.c
create mode 100644 libquadmath/math/atan2piq.c
create mode 100644 libquadmath/math/atanpiq.c
create mode 1
Re: [PATCH] libquadmath: add quad support for trig-pi functions
On Thu, Jul 03, 2025 at 02:43:43PM +0200, Michael Matz wrote:
> Hello,
>
> On Thu, 3 Jul 2025, Yuao Ma wrote:
>
> > This patch adds the required function for Fortran trigonometric functions to
> > work with glibc versions prior to 2.26. It's based on glibc source commit
> > 632d895f3e5d98162f77b9c3c1da4ec19968b671.
> >
> > I've built it successfully on my end. Documentation is also included.
> >
> > Please take a look when you have a moment.
>
> +__float128
> +cospiq (__float128 x)
> +{
> ...
> + if (__builtin_islessequal (x, 0.25Q))
> +return cosq (M_PIq * x);
>
> Isn't the whole raison d'etre for the trig-pi functions that the internal
> argument reduction against multiples of pi becomes trivial and hence (a)
> performant, and (b) doesn't introduce rounding artifacts? Expressing the
> trig-pi functions in terms of their counterparts completely defeats this
> purpose. The other way around would be more sensible for the cases where
> it works, but the above doesn't seem very attractive.
>
It's more than just easier range reduction. There are special
cases that give exact results. For example, if x is integer,
then sinpi(x) = +-0, exactly.
It's also telling when looking at reported accuracy
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
--
Steve
[PATCH] fortran: Add the preliminary code of MOVE_ALLOC arguments
From: Mikael Morin
Regression-tested on aarch64-unknown-linux-gnu.
OK for master?
-- >8 --
Add the preliminary code produced for the evaluation of the FROM and TO
arguments of the MOVE_ALLOC intrinsic before using their values.
Before this change, the preliminary code was ignored and dropped,
limiting the validity of the implementation of MOVE_ALLOC to simple
cases without preliminary code.
This change also adds the cleanup code of the same arguments. It
doesn't make any difference on the testcase though. Because of the
limited set of arguments that are allowed (variables or components
without subreference), it is possible that the cleanup code is actually
guaranteed to be empty. At least adding the cleanup code makes the
array case consistent with the scalar case.
gcc/fortran/ChangeLog:
* trans-intrinsic.cc (conv_intrinsic_move_alloc): Add pre and
post code for the FROM and TO arguments.
gcc/testsuite/ChangeLog:
* gfortran.dg/move_alloc_20.f03: New test.
---
gcc/fortran/trans-intrinsic.cc | 5 +
gcc/testsuite/gfortran.dg/move_alloc_20.f03 | 151
2 files changed, 156 insertions(+)
create mode 100644 gcc/testsuite/gfortran.dg/move_alloc_20.f03
diff --git a/gcc/fortran/trans-intrinsic.cc b/gcc/fortran/trans-intrinsic.cc
index f1bfd3eee51..be984271d6a 100644
--- a/gcc/fortran/trans-intrinsic.cc
+++ b/gcc/fortran/trans-intrinsic.cc
@@ -13101,6 +13101,8 @@ conv_intrinsic_move_alloc (gfc_code *code)
}
gfc_conv_expr_descriptor (&to_se, to_expr);
gfc_conv_expr_descriptor (&from_se, from_expr);
+ gfc_add_block_to_block (&block, &to_se.pre);
+ gfc_add_block_to_block (&block, &from_se.pre);
/* For coarrays, call SYNC ALL if TO is already deallocated as MOVE_ALLOC
is an image control "statement", cf. IR F08/0040 in 12-006A. */
@@ -13174,6 +13176,9 @@ conv_intrinsic_move_alloc (gfc_code *code)
if (fin_label)
gfc_add_expr_to_block (&block, build1_v (LABEL_EXPR, fin_label));
+ gfc_add_block_to_block (&block, &to_se.post);
+ gfc_add_block_to_block (&block, &from_se.post);
+
return gfc_finish_block (&block);
}
diff --git a/gcc/testsuite/gfortran.dg/move_alloc_20.f03
b/gcc/testsuite/gfortran.dg/move_alloc_20.f03
new file mode 100644
index 000..20403c30028
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/move_alloc_20.f03
@@ -0,0 +1,151 @@
+! { dg-do run }
+!
+! Check the presence of the pre and post code of the FROM and TO arguments
+! of the MOVE_ALLOC intrinsic subroutine.
+
+module m
+ implicit none
+ type :: t
+integer, allocatable :: a(:)
+ end type
+end module
+
+module pre
+ use m
+ implicit none
+ private
+ public :: check_pre
+
+contains
+
+ subroutine check_pre
+integer, parameter :: n = 5
+type(t) :: x(n)
+integer, allocatable :: tmp(:)
+integer :: array(4) = [ -1, 0, 1, 2 ]
+integer :: i
+
+if (allocated(tmp)) error stop 1
+
+tmp = [17]
+
+if (.not. allocated(tmp)) error stop 11
+if (any(shape(tmp) /= [1])) error stop 12
+if (any(tmp /= [17])) error stop 13
+do i=1,n
+ if (allocated(x(i)%a)) error stop 14
+end do
+
+! Check that the index of X is properly computed for the evaluation of TO.
+call move_alloc(tmp, x(sum(array))%a)
+
+do i=1,n
+ if (i == 2) cycle
+ if (allocated(x(i)%a)) error stop 21
+end do
+if (.not. allocated(x(2)%a)) error stop 22
+if (any(shape(x(2)%a) /= [1])) error stop 23
+if (any(x(2)%a /= [17])) error stop 24
+if (allocated(tmp)) error stop 25
+
+! Check that the index of X is properly computed for the evaluation of
FROM.
+call move_alloc(x(sum(array))%a, tmp)
+
+if (.not. allocated(tmp)) error stop 31
+if (any(shape(tmp) /= [1])) error stop 32
+if (any(tmp /= [17])) error stop 33
+do i=1,n
+ if (allocated(x(i)%a)) error stop 34
+end do
+ end subroutine
+
+end module
+
+module post
+ use m
+ implicit none
+ private
+ public :: check_post
+ integer, parameter :: n = 5
+ type(t), target :: x(n)
+ type :: u
+integer :: a
+ contains
+final :: finalize
+ end type
+ integer :: finalization_count = 0
+
+contains
+
+ function idx(arg)
+type(u) :: arg
+integer :: idx
+idx = mod(arg%a, n)
+ end function
+
+ subroutine check_post
+type(u) :: y
+integer, allocatable :: tmp(:)
+integer, target :: array(4) = [ -1, 0, 1, 2 ]
+integer :: i
+
+y%a = 12
+
+if (allocated(tmp)) error stop 1
+
+tmp = [37]
+
+if (.not. allocated(tmp)) error stop 11
+if (any(shape(tmp) /= [1])) error stop 12
+if (any(tmp /= [37])) error stop 13
+if (finalization_count /= 0) error stop 14
+do i=1,n
+ if (allocated(x(i)%a)) error stop 15
+end do
+
+! Check that the cleanup code for the evaluation of TO is properly
+! executed after MOVE_ALLOC: the result of GET_U should be finalized.
+call move_alloc(tmp, x(idx(get_u(y)))%a)
+
+do i=1,n
+ if (i == 2) cycle
+
Re: [PATCH] fortran: Add the preliminary code of MOVE_ALLOC arguments
On Thu, Jul 03, 2025 at 10:12:52PM +0200, Mikael Morin wrote: > From: Mikael Morin > > Regression-tested on aarch64-unknown-linux-gnu. > OK for master? > Yes. Almost looks obvious once someone finds and fixes the issue. Thanks for the patch. -- Steve
Re: [PATCH] libquadmath: add quad support for trig-pi functions
On Thu, 3 Jul 2025, Michael Matz wrote: > Yes. And then the above is multiplied by PI, passed to cos/sin and that > one then tries to figure out the multiple of PI (i.e. the 'x' above) again > via range reduction (not a _terribly_ slow one anymore in a good > implementation, because of the limited input range, but still). That's > backwards. People are encouraged to contribute format-specific implementations of these functions to glibc (for example, from CORE-MATH), after adding inputs to the glibc benchmarks so it can be verified if new implementations are a performance improvement. Indeed, that's already been done for various functions for binary32 format. The point of type-generic implementations is not to be optimal for performance or accuracy but to get the API supported across all glibc configurations within a reasonable time (including formats such as IBM long double that probably no-one cares enough about to write specific implementations for now), so that further improvements can be made incrementally. -- Joseph S. Myers josmy...@redhat.com
Re: [PATCH] libquadmath: add quad support for trig-pi functions
On Jul 03 2025, Michael Matz wrote: > Yes. And then the above is multiplied by PI, passed to cos/sin and that > one then tries to figure out the multiple of PI (i.e. the 'x' above) again > via range reduction (not a _terribly_ slow one anymore in a good > implementation, because of the limited input range, but still). That "range reduction" consists of a simple compare against the point where range reduction is required, which is cheap. -- Andreas Schwab, sch...@linux-m68k.org GPG Key fingerprint = 7578 EB47 D4E5 4D69 2510 2552 DF73 E780 A9DA AEC1 "And now for something completely different."
Re: [PATCH] libquadmath: add quad support for trig-pi functions
Hello, On Thu, 3 Jul 2025, Joseph Myers wrote: > > > Isn't the whole raison d'etre for the trig-pi functions that the internal > > > argument reduction against multiples of pi becomes trivial and hence (a) > > > performant, and (b) doesn't introduce rounding artifacts? Expressing the > > > trig-pi functions in terms of their counterparts completely defeats this > > > purpose. The other way around would be more sensible for the cases where > > > it works, but the above doesn't seem very attractive. > > > x = M_FABS (x - M_LIT (2.0) * M_SUF (round) (M_LIT (0.5) * x)); > > In particular, this is what trivial range reduction looks like: no need to > do multiple-precision multiplication with the relevant bits of a > multiple-precision value of 1/pi, just round to the nearest integer > (typically a single instruction). Yes. And then the above is multiplied by PI, passed to cos/sin and that one then tries to figure out the multiple of PI (i.e. the 'x' above) again via range reduction (not a _terribly_ slow one anymore in a good implementation, because of the limited input range, but still). That's backwards. Ciao, Michael.
Re: [PATCH] libquadmath: add quad support for trig-pi functions
Hello,
On Thu, 3 Jul 2025, Yuao Ma wrote:
> This patch adds the required function for Fortran trigonometric functions to
> work with glibc versions prior to 2.26. It's based on glibc source commit
> 632d895f3e5d98162f77b9c3c1da4ec19968b671.
>
> I've built it successfully on my end. Documentation is also included.
>
> Please take a look when you have a moment.
+__float128
+cospiq (__float128 x)
+{
...
+ if (__builtin_islessequal (x, 0.25Q))
+return cosq (M_PIq * x);
Isn't the whole raison d'etre for the trig-pi functions that the internal
argument reduction against multiples of pi becomes trivial and hence (a)
performant, and (b) doesn't introduce rounding artifacts? Expressing the
trig-pi functions in terms of their counterparts completely defeats this
purpose. The other way around would be more sensible for the cases where
it works, but the above doesn't seem very attractive.
Ciao,
Michael.
Re: [Fortran, Patch, PR120843, v3] Fix reject valid, because of inconformable coranks
Hi Jerry, thanks for the review and the ok. Committed as gcc-16-1967-g15413e05eb9. And special thanks for the kind words in the private mail you send me. It's very much appreciated that you even applied a translator to translate it to German. Thank you very much. I have set myself a reminder to backport this patch to gcc-15 is about a week. Regards, Andre On Wed, 2 Jul 2025 09:40:59 -0700 Jerry D wrote: > On 7/2/25 3:14 AM, Andre Vehreschild wrote: > > Hi all, > > > > I successfully created a big mess with the previous patch. First of all by > > applying an outdated one and secondly by adding the conformance checks for > > coranks in a3f1cdd8ed46f9816b31ab162ae4dac547d34ebc. Checking the standard > > even using AI (haha) to figure if coranks of an expression have > > restrictions on them, failed. I found nothing. AI fantasized about > > restrictions that did not exist. Therefore the current approach is to > > remove the conformance check and just use the computed coranks in > > expressions to prevent recomputaion whenever they needed. > > > > Jerry, Harald: Sorry for all the bother and all my mistakes. I am really > > sorry to have wasted your time. > > > > The patch has been regtested fine on x86_64-pc-linux-gnu / F41. Ok for > > mainline and later backport to gcc-15? > > > > Regards, > > Andre > > > --- snip --- > > With this fixer patch, I can successfully compile Toon's test case. > > The patch also regression tests here OK. > > OK to push. > > Jerry -- Andre Vehreschild * Email: vehre ad gmx dot de
