| Message ID | bca16da8-46e7-a274-3724-6472aa6803bd@verizon.net |
|---|---|
| State | New |
| Headers | show |
| Series | [libquadmath/PR68686] | expand |
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
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 <errno.h> +#include <limits.h> #include <stdint.h> #include <stdlib.h> #include "quadmath.h" #include "config.h" +#ifdef HAVE_FENV_H +# include <fenv.h> +#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 under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# Libquadmath is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with libquadmath; if not, see +# <https://www.gnu.org/licenses/>. + +# Usage: update-quadmath.py glibc_srcdir quadmath_srcdir + +import argparse +import os.path +import re + + +def replace_in_file(repl_map, src, dest): + """Apply the replacements in repl_map to the file src, producing dest.""" + with open(src, 'r') as src_file: + text = src_file.read() + for re_src, re_repl in sorted(repl_map.items()): + text = re.sub(re_src, re_repl, text) + text = text.rstrip() + '\n' + with open(dest, 'w') as dest_file: + dest_file.write(text) + + +def update_sources(glibc_srcdir, quadmath_srcdir): + """Update libquadmath sources.""" + glibc_ldbl128 = os.path.join(glibc_srcdir, 'sysdeps/ieee754/ldbl-128') + glibc_math = os.path.join(glibc_srcdir, 'math') + quadmath_math = os.path.join(quadmath_srcdir, 'math') + float128_h = os.path.join(glibc_srcdir, + 'sysdeps/ieee754/float128/float128_private.h') + repl_map = {} + # Use float128_private.h to get an initial list of names to + # replace for libquadmath. + repl_names = {} + with open(float128_h, 'r') as header: + for line in header: + line = line.strip() + if not line.startswith('#define '): + continue + match = re.fullmatch('^#define[ \t]+([a-zA-Z0-9_]+)' + '[ \t]+([a-zA-Z0-9_]+)', line) + if not match: + continue + macro = match.group(1) + result = match.group(2) + result = result.replace('f128', 'q') + result = result.replace('__ieee754_', '') + if result not in ('__expq_table', '__sincosq_table', + '__builtin_signbit'): + result = result.replace('__', '') + result = result.replace('_do_not_use', '') + if result in ('rem_pio2q', 'kernel_sincosq', 'kernel_sinq', + 'kernel_cosq', 'kernel_tanq', 'x2y2m1q'): + # Internal function names, for which the above removal + # of leading '__' was inappropriate and a leading + # '__quadmath_' needs adding instead. + result = '__quadmath_' + result + if result == 'ieee854_float128_shape_type': + result = 'ieee854_float128' + if result == 'HUGE_VAL_F128': + result = 'HUGE_VALQ' + repl_names[macro] = result + # More such names that aren't simply defined as object-like macros + # in float128_private.h. + repl_names['_Float128'] = '__float128' + repl_names['parts32'] = 'words32' + for macro in ('GET_LDOUBLE_LSW64', 'GET_LDOUBLE_MSW64', + 'GET_LDOUBLE_WORDS64', 'SET_LDOUBLE_LSW64', + 'SET_LDOUBLE_MSW64', 'SET_LDOUBLE_WORDS64'): + repl_names[macro] = macro.replace('LDOUBLE', 'FLT128') + # The classication macros are replaced. + for macro in ('FP_NAN', 'FP_INFINITE', 'FP_ZERO', 'FP_SUBNORMAL', + 'FP_NORMAL'): + repl_names[macro] = 'QUAD' + macro + for macro in ('fpclassify', 'signbit', 'isnan', 'isinf'): + repl_names[macro] = macro + 'q' + repl_names['isfinite'] = 'finiteq' + # Map comparison macros to the __builtin forms. + for macro in ('isgreater', 'isgreaterequal', 'isless', 'islessequal', + 'islessgreater', 'isunordered'): + repl_names[macro] = '__builtin_' + macro + # Replace macros used in type-generic templates in glibc. + repl_names['FLOAT'] = '__float128' + repl_names['CFLOAT'] = '__complex128' + repl_names['M_NAN'] = 'nanq ("")' + repl_names['M_HUGE_VAL'] = 'HUGE_VALQ' + repl_names['INFINITY'] = '__builtin_inf ()' + for macro in ('MIN_EXP', 'MAX_EXP', 'MIN', 'MAX', 'MANT_DIG', 'EPSILON'): + repl_names['M_%s' % macro] = 'FLT128_%s' % macro + for macro in ('COPYSIGN', 'FABS', 'SINCOS', 'SCALBN', 'LOG1P', 'ATAN2', + 'COSH', 'EXP', 'HYPOT', 'LOG', 'SINH', 'SQRT'): + repl_names['M_%s' % macro] = macro.lower() + 'q' + # Each such name is replaced when it appears as a whole word. + for macro in repl_names: + repl_map[r'\b%s\b' % macro] = repl_names[macro] + # issignalingq not yet defined in libquadmath. + repl_map[r'\bissignaling *\(.*?\)'] = '0' + # Likewise SET_RESTORE_ROUNDL (should map to SET_RESTORE_ROUNDF128 + # above once supported). + repl_map[r'\bSET_RESTORE_ROUNDL *\(.*?\)'] = '(void) 0' + # Also replace the L macro for constants; likewise M_LIT and M_MLIT. + repl_map[r'\bL *\((.*?)\)'] = r'\1Q' + repl_map[r'\bM_LIT *\((.*?)\)'] = r'\1Q' + repl_map[r'\bM_MLIT *\((.*?)\)'] = r'\1q' + # M_DECL_FUNC and M_SUF need similar replacements. + repl_map[r'\bM_DECL_FUNC *\((?:__)?(.*?)\)'] = r'\1q' + repl_map[r'\bM_SUF *\((?:__)?(?:ieee754_)?(.*?)\)'] = r'\1q' + # Further adjustments are then needed for certain internal + # functions called via M_SUF. + repl_map[r'\bx2y2m1q\b'] = '__quadmath_x2y2m1q' + # Replace calls to __set_errno. + repl_map[r'\b__set_errno *\((.*?)\)'] = r'errno = \1' + # Different names used in union. + repl_map[r'\.d\b'] = '.value' + # Calls to alias and hidden_def macros are all eliminated. + for macro in ('strong_alias', 'weak_alias', 'libm_alias_ldouble', + 'declare_mgen_alias', 'libm_hidden_def', 'mathx_hidden_def'): + repl_map[r'\b%s *\(.*?\);?' % macro] = '' + # Replace all #includes with a single include of quadmath-imp.h. + repl_map['(\n+#include[^\n]*)+\n+'] = '\n\n#include "quadmath-imp.h"\n\n' + # Omitted from this list because code comes from more than one + # glibc source file: exp, rem_pio2, tan. Omitted because of a + # union not currently provided in libquadmath: fma. Omitted + # because of fallback macro definitions not currently handled + # here: ilogb. + ldbl_files = { + 'e_acoshl.c': 'acoshq.c', 'e_acosl.c': 'acosq.c', + 's_asinhl.c': 'asinhq.c', 'e_asinl.c': 'asinq.c', + 'e_atan2l.c': 'atan2q.c', 'e_atanhl.c': 'atanhq.c', + 's_atanl.c': 'atanq.c', 's_cbrtl.c': 'cbrtq.c', 's_ceill.c': 'ceilq.c', + 's_copysignl.c': 'copysignq.c', 'e_coshl.c': 'coshq.c', + 's_cosl.c': 'cosq.c', 'k_cosl.c': 'cosq_kernel.c', + 's_erfl.c': 'erfq.c', 's_expm1l.c': 'expm1q.c', 's_fabsl.c': 'fabsq.c', + 's_finitel.c': 'finiteq.c', 's_floorl.c': 'floorq.c', + 'e_fmodl.c': 'fmodq.c', 's_frexpl.c': 'frexpq.c', + 'e_hypotl.c': 'hypotq.c', 's_isinfl.c': 'isinfq.c', + 's_isnanl.c': 'isnanq.c', 'e_j0l.c': 'j0q.c', 'e_j1l.c': 'j1q.c', + 'e_jnl.c': 'jnq.c', 's_llrintl.c': 'llrintq.c', + 's_llroundl.c': 'llroundq.c', 'e_log10l.c': 'log10q.c', + 's_log1pl.c': 'log1pq.c', 'e_log2l.c': 'log2q.c', + 's_logbl.c': 'logbq.c', 'e_logl.c': 'logq.c', 's_lrintl.c': 'lrintq.c', + 's_lroundl.c': 'lroundq.c', 's_modfl.c': 'modfq.c', + 's_nearbyintl.c': 'nearbyintq.c', 's_nextafterl.c': 'nextafterq.c', + 'e_powl.c': 'powq.c', 'e_remainderl.c': 'remainderq.c', + 's_remquol.c': 'remquoq.c', 's_rintl.c': 'rintq.c', + 's_roundl.c': 'roundq.c', 's_scalblnl.c': 'scalblnq.c', + 's_scalbnl.c': 'scalbnq.c', 's_signbitl.c': 'signbitq.c', + 't_sincosl.c': 'sincos_table.c', 's_sincosl.c': 'sincosq.c', + 'k_sincosl.c': 'sincosq_kernel.c', 'e_sinhl.c': 'sinhq.c', + 's_sinl.c': 'sinq.c', 'k_sinl.c': 'sinq_kernel.c', + 's_tanhl.c': 'tanhq.c', 's_truncl.c': 'truncq.c' + } + # Omitted from this list because of dependency on __kernel_casinh, + # not currently in libquadmath: cacosh cacos casinh. + template_files = { + 's_casin_template.c': 'casinq.c', + 's_catanh_template.c': 'catanhq.c', 's_catan_template.c': 'catanq.c', + 's_ccosh_template.c': 'ccoshq.c', 's_cexp_template.c': 'cexpq.c', + 'cimag_template.c': 'cimagq.c', 's_clog10_template.c': 'clog10q.c', + 's_clog_template.c': 'clogq.c', 'conj_template.c': 'conjq.c', + 's_cproj_template.c': 'cprojq.c', 'creal_template.c': 'crealq.c', + 's_csinh_template.c': 'csinhq.c', 's_csin_template.c': 'csinq.c', + 's_csqrt_template.c': 'csqrtq.c', 's_ctanh_template.c': 'ctanhq.c', + 's_ctan_template.c': 'ctanq.c', 's_fdim_template.c': 'fdimq.c', + 's_fmax_template.c': 'fmaxq.c', 's_fmin_template.c': 'fminq.c', + 's_ldexp_template.c': 'ldexpq.c' + } + for src, dest in ldbl_files.items(): + replace_in_file(repl_map, os.path.join(glibc_ldbl128, src), + os.path.join(quadmath_math, dest)) + for src, dest in template_files.items(): + replace_in_file(repl_map, os.path.join(glibc_math, src), + os.path.join(quadmath_math, dest)) + +def main(): + parser = argparse.ArgumentParser(description='Update libquadmath code.') + parser.add_argument('glibc_srcdir', help='glibc source directory') + parser.add_argument('quadmath_srcdir', help='libquadmath source directory') + args = parser.parse_args() + update_sources(args.glibc_srcdir, args.quadmath_srcdir) + + +if __name__ == '__main__': + main() Property changes on: libquadmath/update-quadmath.py
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 <errno.h> +#include <limits.h> +#include <stdbool.h> #include <stdint.h> #include <stdlib.h> #include "quadmath.h" #include "config.h" +#ifdef HAVE_FENV_H +# include <fenv.h> +#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_gammaq_r (__float128, int *); +extern __float128 __quadmath_lgamma_negq (__float128, int *); +extern __float128 __quadmath_lgamma_productq (__float128, __float128, + __float128, int); +extern __float128 __quadmath_lgammaq_r (__float128, int *); extern __float128 __quadmath_x2y2m1q (__float128 x, __float128 y); -extern int __quadmath_isinf_nsq (__float128 x); +extern __complex128 __quadmath_kernel_casinhq (__complex128, int); +static inline void +mul_splitq (__float128 *hi, __float128 *lo, __float128 x, __float128 y) +{ + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = fmaq (x, y, -*hi); +} @@ -227,4 +252,86 @@ } \ 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 + +#if defined HAVE_FENV_H && defined HAVE_FESETROUND && defined HAVE_FEUPDATEENV +struct rm_ctx +{ + fenv_t env; + bool updated_status; +}; + +# define SET_RESTORE_ROUNDF128(RM) \ + struct rm_ctx ctx __attribute__((cleanup (libc_feresetround_ctx))); \ + libc_feholdsetround_ctx (&ctx, (RM)) + +static inline __attribute__ ((always_inline)) void +libc_feholdsetround_ctx (struct rm_ctx *ctx, int round) +{ + ctx->updated_status = false; + + /* Update rounding mode only if different. */ + if (__glibc_unlikely (round != fegetround ())) + { + ctx->updated_status = true; + fegetenv (&ctx->env); + fesetround (round); + } +} + +static inline __attribute__ ((always_inline)) void +libc_feresetround_ctx (struct rm_ctx *ctx) +{ + /* Restore the rounding mode if updated. */ + if (__glibc_unlikely (ctx->updated_status)) + feupdateenv (&ctx->env); +} +#else +# define SET_RESTORE_ROUNDF128(RM) ((void) 0) +#endif + +#endif Index: libquadmath/quadmath.h =================================================================== --- libquadmath/quadmath.h (revision 265750) +++ libquadmath/quadmath.h (working copy) @@ -58,6 +58,7 @@ extern __float128 cosq (__float128) __quadmath_throw; extern __float128 erfq (__float128) __quadmath_throw; extern __float128 erfcq (__float128) __quadmath_throw; +extern __float128 exp2q (__float128) __quadmath_throw; extern __float128 expq (__float128) __quadmath_throw; extern __float128 expm1q (__float128) __quadmath_throw; extern __float128 fabsq (__float128) __quadmath_throw; @@ -73,6 +74,7 @@ extern int isinfq (__float128) __quadmath_throw; extern int ilogbq (__float128) __quadmath_throw; extern int isnanq (__float128) __quadmath_throw; +extern int issignalingq (__float128) __quadmath_throw; extern __float128 j0q (__float128) __quadmath_throw; extern __float128 j1q (__float128) __quadmath_throw; extern __float128 jnq (int, __float128) __quadmath_throw; Index: libquadmath/quadmath.map =================================================================== --- libquadmath/quadmath.map (revision 265750) +++ libquadmath/quadmath.map (working copy) @@ -101,3 +101,9 @@ global: logbq; } QUADMATH_1.0; + +QUADMATH_1.2 { + global: + exp2q; + issignalingq; +} QUADMATH_1.1; Index: libquadmath/quadmath_weak.h =================================================================== --- libquadmath/quadmath_weak.h (revision 265750) +++ libquadmath/quadmath_weak.h (working copy) @@ -50,6 +50,7 @@ __qmath3 (cosq) __qmath3 (erfq) __qmath3 (erfcq) +__qmath3 (exp2q) __qmath3 (expq) __qmath3 (expm1q) __qmath3 (fabsq) @@ -65,6 +66,7 @@ __qmath3 (ilogbq) __qmath3 (isinfq) __qmath3 (isnanq) +__qmath3 (issignalingq) __qmath3 (j0q) __qmath3 (j1q) __qmath3 (jnq) Index: libquadmath/update-quadmath.py =================================================================== --- libquadmath/update-quadmath.py (nonexistent) +++ libquadmath/update-quadmath.py (working copy) @@ -0,0 +1,255 @@ +#!/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 under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# Libquadmath is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with libquadmath; if not, see +# <https://www.gnu.org/licenses/>. + +# Usage: update-quadmath.py glibc_srcdir quadmath_srcdir + +import argparse +from collections import defaultdict +import os.path +import re + + +def replace_in_file(repl_map, extra_map, src, dest): + """Apply the replacements in repl_map, then those in extra_map, to the + file src, producing dest.""" + with open(src, 'r') as src_file: + text = src_file.read() + for re_src, re_repl in sorted(repl_map.items()): + text = re.sub(re_src, re_repl, text) + for re_src, re_repl in sorted(extra_map.items()): + text = re.sub(re_src, re_repl, text) + text = text.rstrip() + '\n' + with open(dest, 'w') as dest_file: + dest_file.write(text) + + +def update_sources(glibc_srcdir, quadmath_srcdir): + """Update libquadmath sources.""" + glibc_ldbl128 = os.path.join(glibc_srcdir, 'sysdeps/ieee754/ldbl-128') + glibc_math = os.path.join(glibc_srcdir, 'math') + quadmath_math = os.path.join(quadmath_srcdir, 'math') + float128_h = os.path.join(glibc_srcdir, + 'sysdeps/ieee754/float128/float128_private.h') + repl_map = {} + # Use float128_private.h to get an initial list of names to + # replace for libquadmath. + repl_names = {} + with open(float128_h, 'r') as header: + for line in header: + line = line.strip() + if not line.startswith('#define '): + continue + match = re.fullmatch('^#define[ \t]+([a-zA-Z0-9_]+)' + '[ \t]+([a-zA-Z0-9_]+)', line) + if not match: + continue + macro = match.group(1) + result = match.group(2) + result = result.replace('f128', 'q') + result = result.replace('__ieee754_', '') + if result not in ('__expq_table', '__sincosq_table', + '__builtin_signbit'): + result = result.replace('__', '') + result = result.replace('_do_not_use', '') + if result in ('rem_pio2q', 'kernel_sincosq', 'kernel_sinq', + 'kernel_cosq', 'kernel_tanq', 'gammaq_r', + 'gamma_productq', 'lgamma_negq', 'lgamma_productq', + 'lgammaq_r', 'x2y2m1q'): + # Internal function names, for which the above removal + # of leading '__' was inappropriate and a leading + # '__quadmath_' needs adding instead. In the + # libquadmath context, lgammaq_r is an internal name. + result = '__quadmath_' + result + if result == 'ieee854_float128_shape_type': + result = 'ieee854_float128' + if result == 'HUGE_VAL_F128': + result = 'HUGE_VALQ' + repl_names[macro] = result + # More such names that aren't simply defined as object-like macros + # in float128_private.h. + repl_names['_Float128'] = '__float128' + repl_names['SET_RESTORE_ROUNDL'] = 'SET_RESTORE_ROUNDF128' + repl_names['parts32'] = 'words32' + for macro in ('GET_LDOUBLE_LSW64', 'GET_LDOUBLE_MSW64', + 'GET_LDOUBLE_WORDS64', 'SET_LDOUBLE_LSW64', + 'SET_LDOUBLE_MSW64', 'SET_LDOUBLE_WORDS64'): + repl_names[macro] = macro.replace('LDOUBLE', 'FLT128') + # The classication macros are replaced. + for macro in ('FP_NAN', 'FP_INFINITE', 'FP_ZERO', 'FP_SUBNORMAL', + 'FP_NORMAL'): + repl_names[macro] = 'QUAD' + macro + for macro in ('fpclassify', 'signbit', 'isnan', 'isinf', 'issignaling'): + repl_names[macro] = macro + 'q' + repl_names['isfinite'] = 'finiteq' + # Map comparison macros to the __builtin forms. + for macro in ('isgreater', 'isgreaterequal', 'isless', 'islessequal', + 'islessgreater', 'isunordered'): + repl_names[macro] = '__builtin_' + macro + # Replace macros used in type-generic templates in glibc. + repl_names['FLOAT'] = '__float128' + repl_names['CFLOAT'] = '__complex128' + repl_names['M_NAN'] = 'nanq ("")' + repl_names['M_HUGE_VAL'] = 'HUGE_VALQ' + repl_names['INFINITY'] = '__builtin_inf ()' + for macro in ('MIN_EXP', 'MAX_EXP', 'MIN', 'MAX', 'MANT_DIG', 'EPSILON'): + repl_names['M_%s' % macro] = 'FLT128_%s' % macro + for macro in ('COPYSIGN', 'FABS', 'SINCOS', 'SCALBN', 'LOG1P', 'ATAN2', + 'COSH', 'EXP', 'HYPOT', 'LOG', 'SINH', 'SQRT'): + repl_names['M_%s' % macro] = macro.lower() + 'q' + # Each such name is replaced when it appears as a whole word. + for macro in repl_names: + repl_map[r'\b%s\b' % macro] = repl_names[macro] + # Also replace the L macro for constants; likewise M_LIT and M_MLIT. + repl_map[r'\bL *\((.*?)\)'] = r'\1Q' + repl_map[r'\bM_LIT *\((.*?)\)'] = r'\1Q' + repl_map[r'\bM_MLIT *\((.*?)\)'] = r'\1q' + # M_DECL_FUNC and M_SUF need similar replacements. + repl_map[r'\bM_DECL_FUNC *\((?:__)?(?:ieee754_)?(.*?)\)'] = r'\1q' + repl_map[r'\bM_SUF *\((?:__)?(?:ieee754_)?(.*?)\)'] = r'\1q' + # Further adjustments are then needed for certain internal + # functions called via M_SUF. + repl_map[r'\bx2y2m1q\b'] = '__quadmath_x2y2m1q' + repl_map[r'\bkernel_casinhq\b'] = '__quadmath_kernel_casinhq' + # Replace calls to __set_errno. + repl_map[r'\b__set_errno *\((.*?)\)'] = r'errno = \1' + # Eliminate glibc diagnostic macros. + repl_map[r' *\bDIAG_PUSH_NEEDS_COMMENT;'] = '' + repl_map[r' *\bDIAG_IGNORE_NEEDS_COMMENT *\(.*?\);'] = '' + repl_map[r' *\bDIAG_POP_NEEDS_COMMENT;'] = '' + # Different names used in union. + repl_map[r'\.d\b'] = '.value' + repl_map[r'\bunion ieee854_float128\b'] = 'ieee854_float128' + # Calls to alias and hidden_def macros are all eliminated. + for macro in ('strong_alias', 'weak_alias', 'libm_alias_ldouble', + 'declare_mgen_alias', 'declare_mgen_finite_alias', + 'libm_hidden_def', 'mathx_hidden_def'): + repl_map[r'\b%s *\(.*?\);?' % macro] = '' + # Replace all #includes with a single include of quadmath-imp.h. + repl_map['(\n+#include[^\n]*)+\n+'] = '\n\n#include "quadmath-imp.h"\n\n' + # Omitted from this list because code comes from more than one + # glibc source file: rem_pio2. Omitted because of a union not + # currently provided in libquadmath: fma. + ldbl_files = { + 'e_acoshl.c': 'acoshq.c', 'e_acosl.c': 'acosq.c', + 's_asinhl.c': 'asinhq.c', 'e_asinl.c': 'asinq.c', + 'e_atan2l.c': 'atan2q.c', 'e_atanhl.c': 'atanhq.c', + 's_atanl.c': 'atanq.c', 's_cbrtl.c': 'cbrtq.c', 's_ceill.c': 'ceilq.c', + 's_copysignl.c': 'copysignq.c', 'e_coshl.c': 'coshq.c', + 's_cosl.c': 'cosq.c', 'k_cosl.c': 'cosq_kernel.c', + 's_erfl.c': 'erfq.c', 's_expm1l.c': 'expm1q.c', 'e_expl.c': 'expq.c', + 't_expl.h': 'expq_table.h', 's_fabsl.c': 'fabsq.c', + 's_finitel.c': 'finiteq.c', 's_floorl.c': 'floorq.c', + 'e_fmodl.c': 'fmodq.c', 's_frexpl.c': 'frexpq.c', + 'e_lgammal_r.c': 'lgammaq.c', 'lgamma_negl.c': 'lgammaq_neg.c', + 'lgamma_productl.c': 'lgammaq_product.c', 'e_hypotl.c': 'hypotq.c', + 'e_ilogbl.c': 'ilogbq.c', 's_isinfl.c': 'isinfq.c', + 's_isnanl.c': 'isnanq.c', 's_issignalingl.c': 'issignalingq.c', + 'e_j0l.c': 'j0q.c', 'e_j1l.c': 'j1q.c', 'e_jnl.c': 'jnq.c', + 's_llrintl.c': 'llrintq.c', 's_llroundl.c': 'llroundq.c', + 'e_log10l.c': 'log10q.c', 's_log1pl.c': 'log1pq.c', + 'e_log2l.c': 'log2q.c', 's_logbl.c': 'logbq.c', 'e_logl.c': 'logq.c', + 's_lrintl.c': 'lrintq.c', 's_lroundl.c': 'lroundq.c', + 's_modfl.c': 'modfq.c', 's_nearbyintl.c': 'nearbyintq.c', + 's_nextafterl.c': 'nextafterq.c', 'e_powl.c': 'powq.c', + 'e_remainderl.c': 'remainderq.c', 's_remquol.c': 'remquoq.c', + 's_rintl.c': 'rintq.c', 's_roundl.c': 'roundq.c', + 's_scalblnl.c': 'scalblnq.c', 's_scalbnl.c': 'scalbnq.c', + 's_signbitl.c': 'signbitq.c', 't_sincosl.c': 'sincos_table.c', + 's_sincosl.c': 'sincosq.c', 'k_sincosl.c': 'sincosq_kernel.c', + 'e_sinhl.c': 'sinhq.c', 's_sinl.c': 'sinq.c', + 'k_sinl.c': 'sinq_kernel.c', 's_tanhl.c': 'tanhq.c', + 's_tanl.c': 'tanq.c', 'k_tanl.c': 'tanq_kernel.c', + 'e_gammal_r.c': 'tgammaq.c', 'gamma_productl.c': 'tgammaq_product.c', + 's_truncl.c': 'truncq.c', 'x2y2m1l.c': 'x2y2m1q.c' + } + template_files = { + 's_cacosh_template.c': 'cacoshq.c', 's_cacos_template.c': 'cacosq.c', + 's_casinh_template.c': 'casinhq.c', + 'k_casinh_template.c': 'casinhq_kernel.c', + 's_casin_template.c': 'casinq.c', 's_catanh_template.c': 'catanhq.c', + 's_catan_template.c': 'catanq.c', 's_ccosh_template.c': 'ccoshq.c', + 's_cexp_template.c': 'cexpq.c', 'cimag_template.c': 'cimagq.c', + 's_clog10_template.c': 'clog10q.c', 's_clog_template.c': 'clogq.c', + 'conj_template.c': 'conjq.c', 's_cproj_template.c': 'cprojq.c', + 'creal_template.c': 'crealq.c', 's_csinh_template.c': 'csinhq.c', + 's_csin_template.c': 'csinq.c', 's_csqrt_template.c': 'csqrtq.c', + 's_ctanh_template.c': 'ctanhq.c', 's_ctan_template.c': 'ctanq.c', + 'e_exp2_template.c': 'exp2q.c', 's_fdim_template.c': 'fdimq.c', + 's_fmax_template.c': 'fmaxq.c', 's_fmin_template.c': 'fminq.c', + 's_ldexp_template.c': 'ldexpq.c' + } + # Some files have extra substitutions to apply. + extra_maps = defaultdict(dict) + extra_maps['expq.c'] = {r'#include "quadmath-imp\.h"\n': + '#include "quadmath-imp.h"\n' + '#include "expq_table.h"\n'} + extra_maps['ilogbq.c'] = {r'#include "quadmath-imp\.h"\n': + '#include <math.h>\n' + '#include "quadmath-imp.h"\n' + '#ifndef FP_ILOGB0\n' + '# define FP_ILOGB0 INT_MIN\n' + '#endif\n' + '#ifndef FP_ILOGBNAN\n' + '# define FP_ILOGBNAN INT_MAX\n' + '#endif\n', + r'return ([A-Z0-9_]+);': + r'{ errno = EDOM; feraiseexcept (FE_INVALID); ' + r'return \1; }'} + extra_maps['lgammaq.c'] = {r'#include "quadmath-imp\.h"\n': + '#include "quadmath-imp.h"\n' + '#ifdef HAVE_MATH_H_SIGNGAM\n' + '# include <math.h>\n' + '#endif\n' + '__float128\n' + 'lgammaq (__float128 x)\n' + '{\n' + '#ifndef HAVE_MATH_H_SIGNGAM\n' + ' int signgam;\n' + '#endif\n' + ' return __quadmath_lgammaq_r (x, &signgam);\n' + '}\n'} + extra_maps['tgammaq.c'] = {r'#include "quadmath-imp\.h"\n': + '#include "quadmath-imp.h"\n' + '__float128\n' + 'tgammaq (__float128 x)\n' + '{\n' + ' int sign;\n' + ' __float128 ret;\n' + ' ret = __quadmath_gammaq_r (x, &sign);\n' + ' return sign < 0 ? -ret : ret;\n' + '}\n'} + for src, dest in ldbl_files.items(): + replace_in_file(repl_map, extra_maps[dest], + os.path.join(glibc_ldbl128, src), + os.path.join(quadmath_math, dest)) + for src, dest in template_files.items(): + replace_in_file(repl_map, extra_maps[dest], + os.path.join(glibc_math, src), + os.path.join(quadmath_math, dest)) + +def main(): + parser = argparse.ArgumentParser(description='Update libquadmath code.') + parser.add_argument('glibc_srcdir', help='glibc source directory') + parser.add_argument('quadmath_srcdir', help='libquadmath source directory') + args = parser.parse_args() + update_sources(args.glibc_srcdir, args.quadmath_srcdir) + + +if __name__ == '__main__': + main() Property changes on: libquadmath/update-quadmath.py
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 <3dw4rd@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
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 <3dw4rd@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
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 <3dw4rd@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?
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 <https://gcc.gnu.org/ml/gcc-patches/2018-11/msg00136.html> does address the gamma problem.
On Fri, Nov 02, 2018 at 11:43:03PM +0000, 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
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.)
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).
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<typename _Tp> 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.
On Mon, 5 Nov 2018, Jakub Jelinek wrote: > On Fri, Nov 02, 2018 at 11:43:03PM +0000, 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 <fenv.h> 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 <fenv.h>, 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 <joseph@codesourcery.com> 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: Regenerate. * libquadmath.texi (exp2q, issignalingq): Document. * quadmath-imp.h: Include <errno.h>, <limits.h>, <stdbool.h> and <fenv.h>. (HIGH_ORDER_BIT_IS_SET_FOR_SNAN, FIX_FLT128_LONG_CONVERT_OVERFLOW) (FIX_FLT128_LLONG_CONVERT_OVERFLOW, __quadmath_kernel_tanq) (__quadmath_gamma_productq, __quadmath_gammaq_r) (__quadmath_lgamma_negq, __quadmath_lgamma_productq) (__quadmath_lgammaq_r, __quadmath_kernel_casinhq, mul_splitq) (math_check_force_underflow_complex, __glibc_likely) (__glibc_unlikely, struct rm_ctx, SET_RESTORE_ROUNDF128) (libc_feholdsetround_ctx, libc_feresetround_ctx): New. (feraiseexcept, fenv_t, feholdexcept, fesetround, feupdateenv) (fesetenv, fetestexcept, feclearexcept): Define if not supported through <fenv.h>. (__quadmath_isinf_nsq): Remove. * quadmath.h (exp2q, issignalingq): New. * quadmath.map (QUADMATH_1.2): New. * quadmath_weak.h (exp2q, issignalingq): New. * update-quadmath.py: New file. * math/isinf_nsq.c: Remove file. * math/casinhq_kernel.c, math/exp2q.c, math/expq_table.h, math/issignalingq.c, math/lgammaq_neg.c, math/lgammaq_product.c, math/tanq_kernel.c, math/tgammaq_product.c: New files. Generated from glibc sources with update-quadmath.py. * math/acoshq.c, math/acosq.c, math/asinhq.c, math/asinq.c, math/atan2q.c, math/atanhq.c, math/atanq.c, math/cacoshq.c, math/cacosq.c, math/casinhq.c, math/casinq.c, math/catanhq.c, math/catanq.c, math/cbrtq.c, math/ccoshq.c, math/ceilq.c, math/cexpq.c, math/cimagq.c, math/clog10q.c, math/clogq.c, math/conjq.c, math/copysignq.c, math/coshq.c, math/cosq.c, math/cosq_kernel.c, math/cprojq.c, math/crealq.c, math/csinhq.c, math/csinq.c, math/csqrtq.c, math/ctanhq.c, math/ctanq.c, math/erfq.c, math/expm1q.c, math/expq.c, math/fabsq.c, math/fdimq.c, math/finiteq.c, math/floorq.c, math/fmaxq.c, math/fminq.c, math/fmodq.c, math/frexpq.c, math/hypotq.c, math/ilogbq.c, math/isinfq.c, math/isnanq.c, math/j0q.c, math/j1q.c, math/jnq.c, math/ldexpq.c, math/lgammaq.c, math/llrintq.c, math/llroundq.c, math/log10q.c, math/log1pq.c, math/log2q.c, math/logbq.c, math/logq.c, math/lrintq.c, math/lroundq.c, math/modfq.c, math/nearbyintq.c, math/nextafterq.c, math/powq.c, math/remainderq.c, math/remquoq.c, math/rintq.c, math/roundq.c, math/scalblnq.c, math/scalbnq.c, math/signbitq.c, math/sincos_table.c, math/sincosq.c, math/sincosq_kernel.c, math/sinhq.c, math/sinq.c, math/sinq_kernel.c, math/tanhq.c, math/tanq.c, math/tgammaq.c, math/truncq.c, math/x2y2m1q.c: Regenerate from glibc sources with update-quadmath.py.
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; }