Jakub Jelinek
133d0d422e
As mentioned by Joseph in PR105101, glibc 2.26 or later has on x86
(both -m32/-m64), powerpc64le, ia64 and mips support for
*f128 math/complex APIs plus strtof128 and strfromf128, and these APIs allow
us to avoid libquadmath for Fortran purposes on these architectures,
replace *q math/complex APIs, strtof128 instead of strtoflt128 and,
while strfromf128 unfortunately isn't a perfect replacement to
quadmath_snprintf, it can be made to work.
The advantage of this is that when configured against such glibcs
(2.26 is now almost 5 years old), we can avoid linking against an extra shared
library and the math support in glibc is maintained better than libquadmath.
We need both a compiler change (so that for glibc 2.26+ it uses *f128 APIs
instead of *q) and library change.
The above mentioned problem with strfromf128 is that the strfrom* functions
are severely restricted versions of snprintf. In libgfortran, we handle
!isfinite differently and just use snprintf/quadmath_snprintf for
%+-#.*{L,Q}{f,e} printing.
strfrom* doesn't allow +, -, # modifiers and it only supports .34 or
similar precision, not .* . The L/Q etc. letters are omitted.
The + is there to force + sign at the start if it is positive.
Workaround in the patch is to add the + at the start manually for
!signbit (val).
The - (left alignment instead of right) I don't understand why we need it,
when minimum field width isn't specified (for strfrom* can't be specified),
no padding is ever added anywhere I believe.
The # is to force adding . - workaround is to search for first . or e or '\0'
character, if it is '\0', just append ., if it is e, insert . before e and
memmove the rest (which is just a few bytes, e, +/- and at most a few digits)
one byte later.
The .* case is handled by creating the format string for strfrom* by
snprintf into a temporary buffer.
As requested, this patch also switches from using __float128 type in
libgfortran to _Float128 which is equivalent on all arches that support
__float128.
The change is done in a backwards compatible change, when GCC is configured
against glibc 2.26 or newer, libgfortran.so.5 itself doesn't link against
-lquadmath nor uses any libquadmath APIs, libgfortran.a doesn't use any
libquadmath APIs either. User programs and libraries when being linked
by gfortran driver are linked against -lgfortran and -lquadmath, but
the latter only in the --as-needed linker mode, which means it needs
to be around during linking and will be linked in if there are any
calls to math/complex functions with real(kind=16) or complex(kind=16)
in compilation units compiled by older versions of gcc, but if either
user code doesn't call those math/complex functions for the largest
supported kind, or the code is recompiled by gcc with this change in,
libquadmath won't be linked in.
2022-06-28 Jakub Jelinek <jakub@redhat.com>
gcc/fortran/
* gfortran.h (gfc_real_info): Add use_iec_60559 bitfield.
* trans-types.h (gfc_real16_use_iec_60559): Declare.
* trans-types.cc (gfc_real16_use_iec_60559): Define.
(gfc_init_kinds): When building powerpc64le-linux libgfortran
on glibc 2.26 to 2.31, set gfc_real16_use_iec_60559 and
use_iec_60559.
(gfc_build_real_type): Set gfc_real16_use_iec_60559 and use_iec_60559
on glibc 2.26 or later.
* trans-intrinsic.cc (gfc_build_intrinsic_lib_fndecls): Adjust
comment. Handle gfc_real16_use_iec_60559.
(gfc_get_intrinsic_lib_fndecl): Handle use_iec_60559.
libgfortran/
* configure.ac: Check for strtof128 and strfromf128.
Check for math and complex *f128 functions. Set
have_iec_60559_libc_support to yes if *f128 support is around, for
--enable-libquadmath-support default to "default" rather than yes if
have_iec_60559_libc_support is yes.
* acinclude.m4 (LIBGFOR_CHECK_FLOAT128): Test
_Float128/_Complex _Float128 rather than __float128 and
_Complex float __attribute__((mode(TC))). If libquadmath support
is defaulted and have_iec_60559_libc_support is yes, define and subst
USE_IEC_60559. Remove unused LIBGFOR_BUILD_QUAD conditional.
* Makefile.am (kinds.h): Pass @USE_IEC_60559@ as an extra
mk-kinds-h.sh argument.
* mk-kinds-h.sh: Accept 4th use_iec_60559 argument. Use
_Float128/_Complex _Float128 types instead of __float128 and
_Complex float __attribute__((mode(TC))), and if use_iec_60559 is yes,
use f128 suffix instead of q and define GFC_REAL_16_USE_IEC_60559.
* kinds-override.h: Use _Float128/_Complex _Float128 types instead of
__float128 and _Complex float __attribute__((mode(TC))), if
USE_IEC_60559 is defined, use f128 suffixes instead of q and
define GFC_REAL_17_USE_IEC_60559.
* libgfortran.h: Don't include quadmath_weak.h if USE_IEC_60559 is
defined.
(GFC_REAL_16_INFINITY, GFC_REAL_16_QUIET_NAN): Define
for GFC_REAL_16_USE_IEC_60559 differently.
* caf/single.c (convert_type): Use _Float128/_Complex _Float128
instead of __float128 and _Complex float __attribute__((mode(TC))).
For HAVE_GFC_REAL_10 when HAVE_GFC_REAL_16 isn't defined use
_Complex long double instead of long double.
* ieee/issignaling_fallback.h (ieee854_float128_shape_type): Use
_Float128 instead of __float128.
(__issignalingf128): Change argument type to _Float128.
(issignaling): Use _Float128 instead of __float128 in _Generic.
* intrinsics/cshift0.c (cshift0): Use _Float128 instead of __float128
in a comment. Fix a comment typo, logn double -> long double.
* intrinsics/erfc_scaled.c (_THRESH, _M_2_SQRTPI, _INF, _ERFC, _EXP):
Use different definitions if GFC_REAL_16_USE_IEC_60559.
(_THRESH, _M_2_SQRTPI): Use GFC_REAL_17_LITERAL macro.
(_ERFC, _EXP): Use different definitions if GFC_REAL_17_USE_IEC_60559.
* intrinsics/spread_generic.c (spread, spread_scalar): Use _Float128
instead of __float128 in a comment. Fix a comment typo,
logn double -> long double.
* intrinsics/trigd.c (ENABLE_SIND, ENABLE_COSD, ENABLE_TAND): Handle
GFC_REAL_16_USE_IEC_60559.
* intrinsics/pack_generic.c (pack): Use _Float128 instead of
__float128 in a comment. Fix a comment typo, logn double ->
long double.
* intrinsics/unpack_generic.c (unpack1, unpack0): Likewise.
* runtime/in_pack_generic.c (internal_pack): Likewise.
* runtime/in_unpack_generic.c (internal_unpack): Likewise.
* io/read.c (convert_real, convert_infnan): Handle
GFC_REAL_16_USE_IEC_60559 and GFC_REAL_17_USE_IEC_60559.
* io/transfer128.c (tmp1, tmp2): Don't define if libquadmath
isn't needed.
* io/write_float.def (gfor_strfromf128): New function.
(DTOA2Q, FDTOA2Q): Define differently if
GFC_REAL_16_USE_IEC_60559 or GFC_REAL_17_USE_IEC_60559.
* m4/mtype.m4: Use different suffix if GFC_REAL_16_USE_IEC_60559
or GFC_REAL_17_USE_IEC_60559.
* config.h.in: Regenerated.
* configure: Regenerated.
* Makefile.in: Regenerated.
* generated/bessel_r16.c: Regenerated.
* generated/bessel_r17.c: Regenerated.
* generated/norm2_r16.c: Regenerated.
* generated/norm2_r17.c: Regenerated.
177 lines
4.6 KiB
C
177 lines
4.6 KiB
C
/* Implementation of the ERFC_SCALED intrinsic.
|
|
Copyright (C) 2008-2022 Free Software Foundation, Inc.
|
|
|
|
This file is part of the GNU Fortran runtime library (libgfortran).
|
|
|
|
Libgfortran is free software; you can redistribute it and/or
|
|
modify it under the terms of the GNU General Public
|
|
License as published by the Free Software Foundation; either
|
|
version 3 of the License, or (at your option) any later version.
|
|
|
|
Libgfortran 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 General Public License for more details.
|
|
|
|
Under Section 7 of GPL version 3, you are granted additional
|
|
permissions described in the GCC Runtime Library Exception, version
|
|
3.1, as published by the Free Software Foundation.
|
|
|
|
You should have received a copy of the GNU General Public License and
|
|
a copy of the GCC Runtime Library Exception along with this program;
|
|
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
<http://www.gnu.org/licenses/>. */
|
|
|
|
#include "libgfortran.h"
|
|
|
|
/* This implementation of ERFC_SCALED is based on the netlib algorithm
|
|
available at http://www.netlib.org/specfun/erf */
|
|
|
|
#ifdef HAVE_GFC_REAL_4
|
|
#undef KIND
|
|
#define KIND 4
|
|
#include "erfc_scaled_inc.c"
|
|
#endif
|
|
|
|
#ifdef HAVE_GFC_REAL_8
|
|
#undef KIND
|
|
#define KIND 8
|
|
#include "erfc_scaled_inc.c"
|
|
#endif
|
|
|
|
#ifdef HAVE_GFC_REAL_10
|
|
#undef KIND
|
|
#define KIND 10
|
|
#include "erfc_scaled_inc.c"
|
|
#endif
|
|
|
|
#ifdef HAVE_GFC_REAL_16
|
|
|
|
/* For quadruple-precision, netlib's implementation is
|
|
not accurate enough. We provide another one. */
|
|
|
|
#ifdef GFC_REAL_16_IS_FLOAT128
|
|
|
|
# ifdef GFC_REAL_16_USE_IEC_60559
|
|
# define _THRESH -106.566990228185312813205074546585730F128
|
|
# define _M_2_SQRTPI M_2_SQRTPIf128
|
|
# define _INF __builtin_inff128()
|
|
# define _ERFC(x) erfcf128(x)
|
|
# define _EXP(x) expf128(x)
|
|
# else
|
|
# define _THRESH -106.566990228185312813205074546585730Q
|
|
# define _M_2_SQRTPI M_2_SQRTPIq
|
|
# define _INF __builtin_infq()
|
|
# define _ERFC(x) erfcq(x)
|
|
# define _EXP(x) expq(x)
|
|
# endif
|
|
|
|
#else
|
|
|
|
# define _THRESH -106.566990228185312813205074546585730L
|
|
# ifndef M_2_SQRTPIl
|
|
# define M_2_SQRTPIl 1.128379167095512573896158903121545172L
|
|
# endif
|
|
# define _M_2_SQRTPI M_2_SQRTPIl
|
|
# define _INF __builtin_infl()
|
|
# ifdef HAVE_ERFCL
|
|
# define _ERFC(x) erfcl(x)
|
|
# endif
|
|
# ifdef HAVE_EXPL
|
|
# define _EXP(x) expl(x)
|
|
# endif
|
|
|
|
#endif
|
|
|
|
#define ERFC_SCALED(k) \
|
|
GFC_REAL_ ## k \
|
|
erfc_scaled_r ## k (GFC_REAL_ ## k x) \
|
|
{ \
|
|
if (x < _THRESH) \
|
|
{ \
|
|
return _INF; \
|
|
} \
|
|
if (x < 12) \
|
|
{ \
|
|
/* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2). \
|
|
This is not perfect, but much better than netlib. */ \
|
|
return _ERFC(x) * _EXP(x * x); \
|
|
} \
|
|
else \
|
|
{ \
|
|
/* Calculate ERFC_SCALED(x) using a power series in 1/x: \
|
|
ERFC_SCALED(x) = 1 / (x * sqrt(pi)) \
|
|
* (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1)) \
|
|
/ (2 * x**2)**n) \
|
|
*/ \
|
|
GFC_REAL_ ## k sum = 0, oldsum; \
|
|
GFC_REAL_ ## k inv2x2 = 1 / (2 * x * x); \
|
|
GFC_REAL_ ## k fac = 1; \
|
|
int n = 1; \
|
|
\
|
|
while (n < 200) \
|
|
{ \
|
|
fac *= - (2*n - 1) * inv2x2; \
|
|
oldsum = sum; \
|
|
sum += fac; \
|
|
\
|
|
if (sum == oldsum) \
|
|
break; \
|
|
\
|
|
n++; \
|
|
} \
|
|
\
|
|
return (1 + sum) / x * (_M_2_SQRTPI / 2); \
|
|
} \
|
|
}
|
|
|
|
#if defined(_ERFC) && defined(_EXP)
|
|
|
|
extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
|
|
export_proto(erfc_scaled_r16);
|
|
|
|
ERFC_SCALED(16)
|
|
|
|
#endif
|
|
|
|
#undef _THRESH
|
|
#undef _M_2_SQRTPI
|
|
#undef _INF
|
|
#undef _ERFC
|
|
#undef _EXP
|
|
|
|
#endif
|
|
|
|
#ifdef HAVE_GFC_REAL_17
|
|
|
|
/* For quadruple-precision, netlib's implementation is
|
|
not accurate enough. We provide another one. */
|
|
|
|
# define _THRESH GFC_REAL_17_LITERAL(-106.566990228185312813205074546585730)
|
|
# define _M_2_SQRTPI GFC_REAL_17_LITERAL(M_2_SQRTPI)
|
|
# define _INF __builtin_inff128()
|
|
# ifdef POWER_IEEE128
|
|
# define _ERFC(x) __erfcieee128(x)
|
|
# define _EXP(x) __expieee128(x)
|
|
# elif defined(GFC_REAL_17_USE_IEC_60559)
|
|
# define _ERFC(x) erfcf128(x)
|
|
# define _EXP(x) expf128(x)
|
|
# else
|
|
# define _ERFC(x) erfcq(x)
|
|
# define _EXP(x) expq(x)
|
|
# endif
|
|
|
|
extern GFC_REAL_17 erfc_scaled_r17 (GFC_REAL_17);
|
|
export_proto(erfc_scaled_r17);
|
|
|
|
ERFC_SCALED(17)
|
|
|
|
#undef _THRESH
|
|
#undef _M_2_SQRTPI
|
|
#undef _INF
|
|
#undef _ERFC
|
|
#undef _EXP
|
|
#undef ERFC_SCALED
|
|
|
|
#endif
|