Re: [PATCH libquadmath/PR68686]

2018-11-05 Thread Joseph Myers 
On Mon, 5 Nov 2018, Jakub Jelinek wrote:

> On Fri, Nov 02, 2018 at 11:43:03PM +, Joseph Myers wrote:
> > Here's an updated version of the patch that also updates most of the 
> > previously omitted libquadmath/math/ files that are based on glibc sources 
> > (not fmaq.c or rem_pio2q.c), including *gamma*.  It adds exp2q and 
> > issignalingq as new public interfaces, given how they are used in the 
> > current glibc versions of some of the functions already present in 
> > libquadmath, but doesn't add any other new functions from glibc.
> 
> LGTM (with a suitable ChangeLog).

Committed with the following logs.

Update most of libquadmath/math/ from glibc, automate update (PR 
libquadmath/68686).

libquadmath sources are mostly based on glibc sources at present, but
derived from them by a manual editing / substitution process and with
subsequent manual merges.  The manual effort involved in merges means
they are sometimes incomplete and long-delayed.

Since libquadmath was first created, glibc's support for this format
has undergone significant changes so that it can also be used in glibc
to provide *f128 functions for the _Float128 type from TS 18661-3.
This makes it significantly easier to use it for libquadmath in a more
automated fashion, since glibc has a float128_private.h header that
redefines many identifiers as macros as needed for building *f128
functions.

Simply using float128_private.h directly in libquadmath, with
unmodified glibc sources except for changing function names in that
one header to be *q instead of *f128, would be tricky, given its
dependence on lots of other glibc-internal headers (whereas
libquadmath supports non-glibc systems), and also given how some libm
functions in glibc are built from type-generic templates using a
further set of macros rather than from separate function
implementations for each type.

So instead this patch adds a script update-quadmath.py to convert
glibc sources into libquadmath ones, and the script reads
float128_private.h to identify many of the substitutions it should
make.  quadmath-imp.h is updated with various new internal
definitions, taken from glibc as needed; this is the main place
expected to need updating manually when subsequent merges from glibc
are done using the script.  No attempt is made to make the script
output match the details of existing formatting, although the
differences are of a size that makes a rough comparison (ignoring
whitespace) possible.

Two new public interfaces are added to libquadmath, exp2q and
issignalingq, at a new QUADMATH_1.2 symbol version, since those
interfaces are used internally by some of the glibc sources being
merged into libquadmath; although there is a new symbol version, no
change however is made to the libtool version in the libtool-version
file.  Although there are various other interfaces now in glibc libm
but not in libquadmath, this patch does nothing to add such interfaces
(although adding many of them would in fact be easy to do, given the
script).

One internal file (not providing any public interfaces),
math/isinf_nsq.c, is removed, as no longer used by anything in
libquadmath after the merge.

Conditionals in individual source files on  availability or
features are moved into quadmath-imp.h (providing trivial macro
versions of the functions if real implementations aren't available),
to simplify the substitutions in individual source files.  Note
however that I haven't tested for any configurations lacking ,
so further changes could well be needed there.

Two files in libquadmath/math/ are based on glibc sources but not
updated in this patch: fmaq.c and rem_pio2q.c.  Both could be updated
after further changes to the script (and quadmath-imp.h as needed); in
the case of rem_pio2q.c, based on two separate glibc source files,
those separate files would naturally be split out into separate
libquadmath source files in the process (as done in this patch with
expq_table.h and tanq_kernel.c, where previously two glibc source
files had been merged into one libquadmath source file).  complex.c,
nanq.c and sqrtq.c are not based on glibc sources (though four of the
(trivial) functions in complex.c could readily be replaced by instead
using the four corresponding files from glibc, if desired).

libquadmath also has printf/ and strtod/ sources based on glibc, also
mostly not updated for a long time.  Again the script could no doubt
be made to generate those automatically, although that would be a
larger change (effectively some completely separate logic in the
script, not sharing much if anything with the existing code).

Bootstrapped with no regressions on x86_64-pc-linux-gnu.

2018-11-05  Joseph Myers  

PR libquadmath/68686
* Makefile.am: (libquadmath_la_SOURCES): Remove math/isinf_nsq.c.
Add math/exp2q.c math/issignalingq.c math/lgammaq_neg.c
math/lgammaq_product.c math/tanq_kernel.c math/tgammaq_product.c
math/casinhq_kernel.c.
* Makefile.in: Regen

Re: [PATCH libquadmath/PR68686]

2018-11-05 Thread Ed Smith-Rowland 

On 11/5/18 1:19 PM, Joseph Myers wrote:

On Sun, 4 Nov 2018, Ed Smith-Rowland wrote:


I looked in glibc.  Unfortunately, I see how they have the same mistake:
glibc/math/w_tgammal_compat.c:
     long double
     __tgammal(long double x)
     {
         int local_signgam;
         long double y = __ieee754_gammal_r(x,&local_signgam);
     ...
     return local_signgam < 0 ? - y : y;
     }
I'm very sure this is where tgammaq came from.
Ditto for glibc/math/w_tgamma_compat.c and glibc/math/w_tgammaf_compat.c.

No, that's not a mistake.  __ieee754_gammal_r returns +/- the gamma
function and sets the integer pointed to by the second argument to
indicate whether to negate the result.  (This isn't a particularly good
interface design for tgamma, as opposed to lgamma; unfortunately
__gammal_r_finite, with this interface, is a public ABI.)

Excellent, I missed the replacement of  expq (lgammaq (x)) with the 
re-entrant lgamma with the good sign.

Thank you.  I'll let someone else check this off ;-)
FWIW, I'd like to see C++ at least return narrow types like
template
  struct lgamma_t
  {
    _Tp log_abs_gamma; // O something less verbose.
    int sign; // Maybe try for same size as _Tp.
  };
C could do something like also. That's a discussion for another forum.

Re: [PATCH libquadmath/PR68686]

2018-11-05 Thread Joseph Myers 
On Sun, 4 Nov 2018, Ed Smith-Rowland wrote:

> I *do* think a couple tests should be added to test-signgam-*.c to test
> alternation of signs:

The main tests for results of libm functions are in auto-libm-test-in 
(from which auto-libm-test-out-* are generated by gen-auto-libm-tests.c) 
and libm-test-*.inc.  I believe these already cover signgam setting 
thoroughly.  test-signgam-*.c are specifically for ISO C namespace issues 
(an ISO C program should be able to define its own variable called signgam 
and not have lgamma affect it, but an XSI program must be able to use the 
signgam variable defined in libm and have it affected by lgamma).

-- 
Joseph S. Myers
jos...@codesourcery.com

Re: [PATCH libquadmath/PR68686]

2018-11-05 Thread Joseph Myers 
On Sun, 4 Nov 2018, Ed Smith-Rowland wrote:

> I looked in glibc.  Unfortunately, I see how they have the same mistake:
> glibc/math/w_tgammal_compat.c:
>     long double
>     __tgammal(long double x)
>     {
>         int local_signgam;
>         long double y = __ieee754_gammal_r(x,&local_signgam);
>     ...
>     return local_signgam < 0 ? - y : y;
>     }
> I'm very sure this is where tgammaq came from.
> Ditto for glibc/math/w_tgamma_compat.c and glibc/math/w_tgammaf_compat.c.

No, that's not a mistake.  __ieee754_gammal_r returns +/- the gamma 
function and sets the integer pointed to by the second argument to 
indicate whether to negate the result.  (This isn't a particularly good 
interface design for tgamma, as opposed to lgamma; unfortunately 
__gammal_r_finite, with this interface, is a public ABI.)

-- 
Joseph S. Myers
jos...@codesourcery.com

Re: [PATCH libquadmath/PR68686]

2018-11-05 Thread Jakub Jelinek 
On Fri, Nov 02, 2018 at 11:43:03PM +, Joseph Myers wrote:
> Here's an updated version of the patch that also updates most of the 
> previously omitted libquadmath/math/ files that are based on glibc sources 
> (not fmaq.c or rem_pio2q.c), including *gamma*.  It adds exp2q and 
> issignalingq as new public interfaces, given how they are used in the 
> current glibc versions of some of the functions already present in 
> libquadmath, but doesn't add any other new functions from glibc.

LGTM (with a suitable ChangeLog).

Jakub

Re: [PATCH libquadmath/PR68686]

2018-11-05 Thread Joseph Myers 
On Sat, 3 Nov 2018, Jeff Law wrote:

> Note that Joseph's follow-up doesn't touch on the gamma problem AFAICT,
> but instead touches on the larger issues around trying to keep the
> quadmath implementations between glibc and gcc more in sync.

The second version of my patch 
 does address 
the gamma problem.

-- 
Joseph S. Myers
jos...@codesourcery.com

Re: [PATCH libquadmath/PR68686]

2018-11-04 Thread Ed Smith-Rowland 

On 11/3/18 10:09 PM, Jeff Law wrote:

On 10/23/18 7:45 PM, Ed Smith-Rowland wrote:

Greetings,

This is an almost trivial patch to get the correct sign for tgammaq.

I don't have a testcase as I don't know where to put one.

OK?

Ed Smith-Rowland



tgammaq.CL

2018-10-24  Edward Smith-Rowland  <3dw...@verizon.net>

PR libquadmath/68686
* math/tgammaq.c: Correct sign for negative argument.

I don't have the relevant background to evaluate this for correctness.
Can you refer back to any kind of documentation which indicates what the
sign of the return value ought to be?

Alternately, if you can point to the relevant code in glibc that handles
the resultant sign, that'd be useful too.

Note that Joseph's follow-up doesn't touch on the gamma problem AFAICT,
but instead touches on the larger issues around trying to keep the
quadmath implementations between glibc and gcc more in sync.

Jeff


I've looked at glibc lgamma, in particular signgam and I think those DTRT:


I'm pretty sure the lgamma that write to global signgam and the 
lgamma_r(x, int *signgam) DTRT.


The various __lgamma_neg* DTRT:

__lgamma_negX (REALTYPE x, int *signgamp)
{
  /* Determine the half-integer region X lies in, handle exact
 integers and determine the sign of the result.  */
  int i = __floorl (-2 * x);
  if ((i & 1) == 0 && i == -2 * x)
    return 1.0L / 0.0L;
  long double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
  i -= 4;
  *signgamp = ((i & 2) == 0 ? -1 : 1);
...

I think the various e_lgammaX_r.c are good too:

  if (se & 0x8000)
    {
  if (x < -2.0L && x > -33.0L)
    return __lgamma_negl (x, signgamp);
  t = sin_pi (x);
  if (t == zero)
    return one / fabsl (t);    /* -integer */
  nadj = __ieee754_logl (pi / fabsl (t * x));
  if (t < zero)
    *signgamp = -1;
  x = -x;
    }

I *do* think a couple tests should be added to test-signgam-*.c to test 
alternation of signs:

  signgam = 123;                        \
  c = FUNC (b);                        \
  if (signgam == 1)                        \

    puts ("PASS: " #FUNC " (-1.5) setting signgam");    \

  else                            \
    {                            \
      puts ("FAIL: " #FUNC " (-1.5) setting signgam");    \
      result = 1;                        \
    }                            \

Add to test lgamma_negX code paths...
  signgam = 123;                        \
  c = FUNC (b);                        \
  if (signgam == -1)                    \
    puts ("PASS: " #FUNC " (-34.5) setting signgam");    \
  else                            \
    {                            \
      puts ("FAIL: " #FUNC " (-34.5) setting signgam"); \
      result = 1;                        \
    }                            \
  signgam = 123;                        \
  c = FUNC (b);                        \
  if (signgam == 1)                        \
    puts ("PASS: " #FUNC " (-35.5) setting signgam");    \
  else                            \
    {                            \
      puts ("FAIL: " #FUNC " (-35.5) setting signgam"); \
      result = 1;                        \
    }                            \


I've not dealt with glibc directly.  Do I need separate Copyright and 
all that?  Is it similar to gcc in terms of devel?

Re: [PATCH libquadmath/PR68686]

2018-11-04 Thread Ed Smith-Rowland 

On 11/3/18 10:09 PM, Jeff Law wrote:

On 10/23/18 7:45 PM, Ed Smith-Rowland wrote:

Greetings,

This is an almost trivial patch to get the correct sign for tgammaq.

I don't have a testcase as I don't know where to put one.

OK?

Ed Smith-Rowland



tgammaq.CL

2018-10-24  Edward Smith-Rowland  <3dw...@verizon.net>

PR libquadmath/68686
* math/tgammaq.c: Correct sign for negative argument.

I don't have the relevant background to evaluate this for correctness.
Can you refer back to any kind of documentation which indicates what the
sign of the return value ought to be?

Alternately, if you can point to the relevant code in glibc that handles
the resultant sign, that'd be useful too.

Note that Joseph's follow-up doesn't touch on the gamma problem AFAICT,
but instead touches on the larger issues around trying to keep the
quadmath implementations between glibc and gcc more in sync.

Jeff

Thank you for (re)considering this.The maths here can be read from a 
very good NIST website:


For the functional formulas referred to below - https://dlmf.nist.gov/5.5

For cool pictures - https://dlmf.nist.gov/5.3#i
(Aside: the reciprocal gamma function is entire - possessing no poles or 
other analytic complications anywhere in the complex plane.  A rgamma 
function might be a nice addition to the C family and to TS 18661-1 for 
that matter.)


For a table of extrema (Table 5.4.1) - https://dlmf.nist.gov/5.4#iii
These would be nice tests for accuracy as well as sign.

TL;DR:
==

Either follow the factorial-like recursion Gamma(x+1) = x Gamma(x) [DLMF 
5.5.1] backwards:

    Gamma(x-1) = Gamma(x) / (x-1)
Given that Gamma(x) is positive for x > 0 this will give alternating 
negative signs if you start with, say, x=1/2 and keep going.

A cooler formula is [DLMF 5.5.3]:
    Gamma(x) Gamma(1-x) = pi / sin(pi x)
Start with x = 3/2 and continue with higher odd half integers to see the 
sign alternation.


GLibC
=
I looked in glibc.  Unfortunately, I see how they have the same mistake:
glibc/math/w_tgammal_compat.c:
    long double
    __tgammal(long double x)
    {
        int local_signgam;
        long double y = __ieee754_gammal_r(x,&local_signgam);
    ...
    return local_signgam < 0 ? - y : y;
    }
I'm very sure this is where tgammaq came from.
Ditto for glibc/math/w_tgamma_compat.c and glibc/math/w_tgammaf_compat.c.

This fix will need to be done upstream.
Ed

Re: [PATCH libquadmath/PR68686]

2018-11-03 Thread Jeff Law 
On 10/23/18 7:45 PM, Ed Smith-Rowland wrote:
> Greetings,
> 
> This is an almost trivial patch to get the correct sign for tgammaq.
> 
> I don't have a testcase as I don't know where to put one.
> 
> OK?
> 
> Ed Smith-Rowland
> 
> 
> 
> tgammaq.CL
> 
> 2018-10-24  Edward Smith-Rowland  <3dw...@verizon.net>
> 
>   PR libquadmath/68686
>   * math/tgammaq.c: Correct sign for negative argument.
I don't have the relevant background to evaluate this for correctness.
Can you refer back to any kind of documentation which indicates what the
sign of the return value ought to be?

Alternately, if you can point to the relevant code in glibc that handles
the resultant sign, that'd be useful too.

Note that Joseph's follow-up doesn't touch on the gamma problem AFAICT,
but instead touches on the larger issues around trying to keep the
quadmath implementations between glibc and gcc more in sync.

Jeff

Re: [PATCH libquadmath/PR68686]

2018-11-02 Thread Joseph Myers 
On Fri, 2 Nov 2018, Joseph Myers wrote:

> I think it would be best to move to having a script to generate 
> libquadmath sources automatically from glibc sources by appropriate 
> substitutions, so that while you might need to update the script or 
> quadmath-imp.h as part of updating libquadmath from glibc, you don't need 
> to merge lots of changes manually.
> 
> Specifically, any comments on the patch below (quadmath-imp.h changes and 
> new script shown, 6000 lines of diffs from running the script not shown)?  
> It doesn't yet update the *gamma* sources, but could be extended to do so.  
> (It also doesn't do anything with the parts of libquadmath outside of 
> libquadmath/math/, but again could be extended for that.)  Specifically, 
> the following files in libquadmath/math/ aren't yet updated by the script 
> (a few of these, e.g. sqrtq.c, aren't actually based on glibc sources at 
> all, while others just need the script to gain new features, or additional 
> source files to be added to libquadmath): cacoshq.c cacosq.c casinhq.c 
> complex.c expq.c fmaq.c ilogbq.c isinf_nsq.c lgammaq.c nanq.c rem_pio2q.c 
> sqrtq.c tanq.c tgammaq.c x2y2m1q.c.

Here's an updated version of the patch that also updates most of the 
previously omitted libquadmath/math/ files that are based on glibc sources 
(not fmaq.c or rem_pio2q.c), including *gamma*.  It adds exp2q and 
issignalingq as new public interfaces, given how they are used in the 
current glibc versions of some of the functions already present in 
libquadmath, but doesn't add any other new functions from glibc.

Index: libquadmath/Makefile.am
===
--- libquadmath/Makefile.am (revision 265750)
+++ libquadmath/Makefile.am (working copy)
@@ -44,7 +44,7 @@
 libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include
 
 libquadmath_la_SOURCES = \
-  math/x2y2m1q.c math/isinf_nsq.c math/acoshq.c math/fmodq.c \
+  math/x2y2m1q.c math/acoshq.c math/fmodq.c \
   math/acosq.c math/frexpq.c \
   math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \
   math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \
@@ -58,6 +58,8 @@
   math/tanhq.c math/expq.c math/modfq.c math/tanq.c math/fabsq.c \
   math/nanq.c math/tgammaq.c math/finiteq.c math/nextafterq.c \
   math/truncq.c math/floorq.c math/powq.c math/fmaq.c math/logbq.c \
+  math/exp2q.c math/issignalingq.c math/lgammaq_neg.c math/lgammaq_product.c \
+  math/tanq_kernel.c math/tgammaq_product.c math/casinhq_kernel.c \
   math/cacoshq.c math/cacosq.c math/casinhq.c math/casinq.c \
   math/catanhq.c math/catanq.c math/cimagq.c math/conjq.c math/cprojq.c \
   math/crealq.c math/fdimq.c math/fmaxq.c math/fminq.c math/ilogbq.c \
Index: libquadmath/libquadmath.texi
===
--- libquadmath/libquadmath.texi(revision 265750)
+++ libquadmath/libquadmath.texi(working copy)
@@ -157,6 +157,7 @@
 @item @code{cosq}: cosine function
 @item @code{erfq}: error function
 @item @code{erfcq}: complementary error function
+@item @code{exp2q}: base 2 exponential function
 @item @code{expq}: exponential function
 @item @code{expm1q}: exponential minus 1 function
 @need 800
@@ -173,6 +174,7 @@
 @item @code{ilogbq}: get exponent of the value
 @item @code{isinfq}: check for infinity
 @item @code{isnanq}: check for not a number
+@item @code{issignalingq}: check for signaling not a number
 @item @code{j0q}: Bessel function of the first kind, first order
 @item @code{j1q}: Bessel function of the first kind, second order
 @item @code{jnq}: Bessel function of the first kind, @var{n}-th order
Index: libquadmath/quadmath-imp.h
===
--- libquadmath/quadmath-imp.h  (revision 265750)
+++ libquadmath/quadmath-imp.h  (working copy)
@@ -21,10 +21,16 @@
 #ifndef QUADMATH_IMP_H
 #define QUADMATH_IMP_H
 
+#include 
+#include 
+#include 
 #include 
 #include 
 #include "quadmath.h"
 #include "config.h"
+#ifdef HAVE_FENV_H
+# include 
+#endif
 
 
 /* Under IEEE 754, an architecture may determine tininess of
@@ -36,7 +42,11 @@
 
 #define TININESS_AFTER_ROUNDING   1
 
+#define HIGH_ORDER_BIT_IS_SET_FOR_SNAN 0
 
+#define FIX_FLT128_LONG_CONVERT_OVERFLOW 0
+#define FIX_FLT128_LLONG_CONVERT_OVERFLOW 0
+
 /* Prototypes for internal functions.  */
 extern int32_t __quadmath_rem_pio2q (__float128, __float128 *);
 extern void __quadmath_kernel_sincosq (__float128, __float128, __float128 *,
@@ -43,9 +53,24 @@
   __float128 *, int);
 extern __float128 __quadmath_kernel_sinq (__float128, __float128, int);
 extern __float128 __quadmath_kernel_cosq (__float128, __float128);
+extern __float128 __quadmath_kernel_tanq (__float128, __float128, int);
+extern __float128 __quadmath_gamma_productq (__float128, __float128, int,
+__float128 *);
+extern __float128 __quadmath_g

Re: [PATCH libquadmath/PR68686]

2018-11-02 Thread Joseph Myers 
On Wed, 24 Oct 2018, Jakub Jelinek wrote:

> On Tue, Oct 23, 2018 at 09:45:13PM -0400, Ed Smith-Rowland wrote:
> > Greetings,
> > 
> > This is an almost trivial patch to get the correct sign for tgammaq.
> 
> Doesn't look trivial to me.  What happens if x is a NaN?  Or if x is outside
> of the range of int?
> 
> Generally, libquadmath follows what glibc does, I don't see such a change in
> there.  So, the right fix is probably port all the upstream glibc *gamma*
> changes, in PR65757 as the 
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757#c18
> comment and my patch submission says, I've left those out because they were
> too large and I didn't have spare cycles for that.

I think it would be best to move to having a script to generate 
libquadmath sources automatically from glibc sources by appropriate 
substitutions, so that while you might need to update the script or 
quadmath-imp.h as part of updating libquadmath from glibc, you don't need 
to merge lots of changes manually.

Specifically, any comments on the patch below (quadmath-imp.h changes and 
new script shown, 6000 lines of diffs from running the script not shown)?  
It doesn't yet update the *gamma* sources, but could be extended to do so.  
(It also doesn't do anything with the parts of libquadmath outside of 
libquadmath/math/, but again could be extended for that.)  Specifically, 
the following files in libquadmath/math/ aren't yet updated by the script 
(a few of these, e.g. sqrtq.c, aren't actually based on glibc sources at 
all, while others just need the script to gain new features, or additional 
source files to be added to libquadmath): cacoshq.c cacosq.c casinhq.c 
complex.c expq.c fmaq.c ilogbq.c isinf_nsq.c lgammaq.c nanq.c rem_pio2q.c 
sqrtq.c tanq.c tgammaq.c x2y2m1q.c.

(The script does not aim to match formatting, comments etc. with the 
existing libquadmath code in all cases, but the results should be similar 
enough for comparison to be reasonable, especially if you ignore 
whitespace; an important question is whether it's losing deliberate local 
changes in the libquadmath code, that are either relevant for portability 
or affect the generated code.)

Index: libquadmath/quadmath-imp.h
===
--- libquadmath/quadmath-imp.h  (revision 265724)
+++ libquadmath/quadmath-imp.h  (working copy)
@@ -21,10 +21,15 @@
 #ifndef QUADMATH_IMP_H
 #define QUADMATH_IMP_H
 
+#include 
+#include 
 #include 
 #include 
 #include "quadmath.h"
 #include "config.h"
+#ifdef HAVE_FENV_H
+# include 
+#endif
 
 
 /* Under IEEE 754, an architecture may determine tininess of
@@ -36,6 +41,8 @@
 
 #define TININESS_AFTER_ROUNDING   1
 
+#define FIX_FLT128_LONG_CONVERT_OVERFLOW 0
+#define FIX_FLT128_LLONG_CONVERT_OVERFLOW 0
 
 /* Prototypes for internal functions.  */
 extern int32_t __quadmath_rem_pio2q (__float128, __float128 *);
@@ -227,4 +234,50 @@
 }  \
   while (0)
 
+/* Likewise, for both real and imaginary parts of a complex
+   result.  */
+#define math_check_force_underflow_complex(x)  \
+  do   \
+{  \
+  __typeof (x) force_underflow_complex_tmp = (x);  \
+  math_check_force_underflow (__real__ force_underflow_complex_tmp); \
+  math_check_force_underflow (__imag__ force_underflow_complex_tmp); \
+}  \
+  while (0)
+
+#ifndef HAVE_FENV_H
+# define feraiseexcept(arg) ((void) 0)
+typedef int fenv_t;
+# define feholdexcept(arg) ((void) 0)
+# define fesetround(arg) ((void) 0)
+# define feupdateenv(arg) ((void) (arg))
+# define fesetenv(arg) ((void) (arg))
+# define fetestexcept(arg) 0
+# define feclearexcept(arg) ((void) 0)
+#else
+# ifndef HAVE_FEHOLDEXCEPT
+#  define feholdexcept(arg) ((void) 0)
+# endif
+# ifndef HAVE_FESETROUND
+#  define fesetround(arg) ((void) 0)
+# endif
+# ifndef HAVE_FEUPDATEENV
+#  define feupdateenv(arg) ((void) (arg))
+# endif
+# ifndef HAVE_FESETENV
+#  define fesetenv(arg) ((void) (arg))
+# endif
+# ifndef HAVE_FETESTEXCEPT
+#  define fetestexcept(arg) 0
+# endif
 #endif
+
+#ifndef __glibc_likely
+# define __glibc_likely(cond)  __builtin_expect ((cond), 1)
+#endif
+
+#ifndef __glibc_unlikely
+# define __glibc_unlikely(cond)__builtin_expect ((cond), 0)
+#endif
+
+#endif
Index: libquadmath/update-quadmath.py
===
--- libquadmath/update-quadmath.py  (nonexistent)
+++ libquadmath/update-quadmath.py  (working copy)
@@ -0,0 +1,198 @@
+#!/usr/bin/python3
+# Update libquadmath code from glibc sources.
+# Copyright (C) 2018 Free Software Foundation, Inc.
+# This file is part of the libquadmath library.
+#
+# Libquadmath is free software; you can redistribute it and/or
+# modify it u

Re: [PATCH libquadmath/PR68686]

2018-10-24 Thread Jakub Jelinek 
On Tue, Oct 23, 2018 at 09:45:13PM -0400, Ed Smith-Rowland wrote:
> Greetings,
> 
> This is an almost trivial patch to get the correct sign for tgammaq.

Doesn't look trivial to me.  What happens if x is a NaN?  Or if x is outside
of the range of int?

Generally, libquadmath follows what glibc does, I don't see such a change in
there.  So, the right fix is probably port all the upstream glibc *gamma*
changes, in PR65757 as the 
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757#c18
comment and my patch submission says, I've left those out because they were
too large and I didn't have spare cycles for that.

> --- libquadmath/math/tgammaq.c(revision 265345)
> +++ libquadmath/math/tgammaq.c(working copy)
> @@ -47,7 +47,9 @@
>  /* x == -Inf.  According to ISO this is NaN.  */
>  return x - x;
>  
> -  /* XXX FIXME.  */
>res = expq (lgammaq (x));
> -  return signbitq (x) ? -res : res;
> +  if (x > 0.0Q || ((int)(-x) & 1) == 1)
> +return res;
> +  else
> +return -res;
>  }


Jakub

[PATCH libquadmath/PR68686]

2018-10-23 Thread Ed Smith-Rowland 

Greetings,

This is an almost trivial patch to get the correct sign for tgammaq.

I don't have a testcase as I don't know where to put one.

OK?

Ed Smith-Rowland


2018-10-24  Edward Smith-Rowland  <3dw...@verizon.net>

	PR libquadmath/68686
	* math/tgammaq.c: Correct sign for negative argument.
Index: libquadmath/math/tgammaq.c
===
--- libquadmath/math/tgammaq.c	(revision 265345)
+++ libquadmath/math/tgammaq.c	(working copy)
@@ -47,7 +47,9 @@
 /* x == -Inf.  According to ISO this is NaN.  */
 return x - x;
 
-  /* XXX FIXME.  */
   res = expq (lgammaq (x));
-  return signbitq (x) ? -res : res;
+  if (x > 0.0Q || ((int)(-x) & 1) == 1)
+return res;
+  else
+return -res;
 }

Re: [PATCH libquadmath, pr68686] tgammaq(x) is always negative for noninteger x < 0

2017-11-14 Thread Steve Kargl 
Ed,

gfortran uses libquadmath for support of its REAL(16)
intrinsic subprograms.  I checked the list of intrinsics
and tgamma is current not in the list.  I don't have
time to try Fortran's C interop feature to see if an
ordinary user can access __float128.  While your patch
is probably useful for GCC, I doubr that gfortran in
its current state will use tgammal.  You'll need either
Jakub or Joseph (jsm28) to give an OK.

-- 
steve

On Mon, Nov 13, 2017 at 01:27:42PM -0500, Ed Smith-Rowland wrote:
> Here is a patch for tammaq for negative argument pr68686.
> 
> I know about depending on ports from upstream but this was done recently 
> and this (tgammaq) was left out.
> 
> This patch is basically a one-liner.
> 
> I have test cases but libquadmath doesn't have a testsuite.
> 
> One test just shows alternating signs for -0.5Q, -1.5Q, -2.5Q, etc.
> 
> Another test verifies the Gamma reflection formula:
> 
>   tgamma(1-x) * tgamma(x) - pi / sinq(pi*x) is tiny.
> 
> Builds on x86_64-linux and tests correctly offline.
> 
> 
> OK?
> 
> 

> 2017-11-10  Edward Smith-Rowland  <3dw...@verizon.net>
> 
>   PR libquadmath/68686
>   * math/tgammaq.c: Correct sign for negative argument.

> diff --git a/libquadmath/math/tgammaq.c b/libquadmath/math/tgammaq.c
> index a07d583..3080094 100644
> --- a/libquadmath/math/tgammaq.c
> +++ b/libquadmath/math/tgammaq.c
> @@ -47,7 +47,9 @@ tgammaq (__float128 x)
>  /* x == -Inf.  According to ISO this is NaN.  */
>  return x - x;
>  
> -  /* XXX FIXME.  */
>res = expq (lgammaq (x));
> -  return signbitq (x) ? -res : res;
> +  if (x > 0.0Q || ((int)(-x) & 1) == 1)
> +return res;
> +  else
> +return -res;
>  }

> #include 
> #include 
> 
> void
> test_alt_signs()
> {
>   __float128 result = tgammaq(-1.5Q);
>   assert(result > 0.0Q);
> 
>   result = tgammaq(-2.5Q);
>   assert(result < 0.0Q);
> 
>   result = tgammaq(-3.5Q);
>   assert(result > 0.0Q);
> 
>   result = tgammaq(-4.5Q);
>   assert(result < 0.0Q);
> }
> 
> /*
>  * Return |\Gamma(x) \Gamma(1 - x) - \frac{\pi}{sin(\pi x)}|
>  */
> __float128
> abs_delta(__float128 x)
> {
>   return fabsq(tgammaq(x) * tgammaq(1.0Q - x) - M_PIq / sinq(M_PIq * x));
> }
> 
> /*
>  * Test reflection: \Gamma(x) \Gamma(1 - x) = \frac{\pi}{sin(\pi x)}
>  */
> void
> test_reflection()
> {
>   assert(abs_delta(+1.5Q) < 100 * FLT128_EPSILON);
>   assert(abs_delta(+0.5Q) < 100 * FLT128_EPSILON);
>   assert(abs_delta(-0.5Q) < 100 * FLT128_EPSILON);
>   assert(abs_delta(-1.5Q) < 100 * FLT128_EPSILON);
>   assert(abs_delta(-2.5Q) < 100 * FLT128_EPSILON);
> }
> 
> int
> main()
> {
>   test_alt_signs();
> 
>   test_reflection();
> 
>   return 0;
> }


-- 
Steve
20170425 https://www.youtube.com/watch?v=VWUpyCsUKR4
20161221 https://www.youtube.com/watch?v=IbCHE-hONow

[PATCH libquadmath, pr68686] tgammaq(x) is always negative for noninteger x < 0

2017-11-13 Thread Ed Smith-Rowland 

Here is a patch for tammaq for negative argument pr68686.

I know about depending on ports from upstream but this was done recently 
and this (tgammaq) was left out.


This patch is basically a one-liner.

I have test cases but libquadmath doesn't have a testsuite.

One test just shows alternating signs for -0.5Q, -1.5Q, -2.5Q, etc.

Another test verifies the Gamma reflection formula:

 tgamma(1-x) * tgamma(x) - pi / sinq(pi*x) is tiny.

Builds on x86_64-linux and tests correctly offline.


OK?


2017-11-10  Edward Smith-Rowland  <3dw...@verizon.net>

PR libquadmath/68686
* math/tgammaq.c: Correct sign for negative argument.
diff --git a/libquadmath/math/tgammaq.c b/libquadmath/math/tgammaq.c
index a07d583..3080094 100644
--- a/libquadmath/math/tgammaq.c
+++ b/libquadmath/math/tgammaq.c
@@ -47,7 +47,9 @@ tgammaq (__float128 x)
 /* x == -Inf.  According to ISO this is NaN.  */
 return x - x;
 
-  /* XXX FIXME.  */
   res = expq (lgammaq (x));
-  return signbitq (x) ? -res : res;
+  if (x > 0.0Q || ((int)(-x) & 1) == 1)
+return res;
+  else
+return -res;
 }
#include 
#include 

void
test_alt_signs()
{
  __float128 result = tgammaq(-1.5Q);
  assert(result > 0.0Q);

  result = tgammaq(-2.5Q);
  assert(result < 0.0Q);

  result = tgammaq(-3.5Q);
  assert(result > 0.0Q);

  result = tgammaq(-4.5Q);
  assert(result < 0.0Q);
}

/*
 * Return |\Gamma(x) \Gamma(1 - x) - \frac{\pi}{sin(\pi x)}|
 */
__float128
abs_delta(__float128 x)
{
  return fabsq(tgammaq(x) * tgammaq(1.0Q - x) - M_PIq / sinq(M_PIq * x));
}

/*
 * Test reflection: \Gamma(x) \Gamma(1 - x) = \frac{\pi}{sin(\pi x)}
 */
void
test_reflection()
{
  assert(abs_delta(+1.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(+0.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(-0.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(-1.5Q) < 100 * FLT128_EPSILON);
  assert(abs_delta(-2.5Q) < 100 * FLT128_EPSILON);
}

int
main()
{
  test_alt_signs();

  test_reflection();

  return 0;
}