Re: [PATCH]libgfortran: Add additional long double math support for hppa64-hpux*
Jerry DeLisle <[email protected]> Thu, 23 Jul 2026 18:13:55 -0700
| Newsgroups | gmane.comp.gcc.patches,gmane.comp.gcc.fortran |
|---|---|
| Message-ID | <[email protected]> |
On 7/23/26 7:50 AM, John David Anglin wrote: > The attached patch completes the C99 long double support needed > for libgfortran on hppa*-hpux*. It is mostly a mechanical addition. > > The change fixes PR libfortran/126305. > > I added files to support ceill, floorl, frexpl, nextafterl, scalbnl > and truncl. These are implemented by bit twiddling and they are > now needed to successfully build libgfortran. I revised these files > to GCC coding standards to make them more readable. I only implemented > the 16-byte long double format. > > Other functions are implemented by calling routines in libquadmath. > These also have backup implementations using double support that > is limited in accuracy, etc. > > With libquadmath, the only fails are PR100914.f90, default_format_2.f90 > and default_format_denormal_2.f90. The fail of PR100914.f90 is due > to a missing complex.h header. With libquadmath support, norm2_3.f90 > also fails on hppa64-hp-hpux11.11. > > Note several files are copyright by Sun Microsystems. These are > ceill_16.c, floorl_16.c, frexpl_16.c, nextafterl_16.c, and > scalbnl_16.c. Variants of these files exist in glibc and > libquadmath (ceilq.c, floorq.c, frexpq.c, nextafterq.c and scalbnq.c. > > I believe this is consistent with the GCC Runtime Library Exception > but this needs to be confirmed by someone more knowledgeable about > copyright issues. > > Okay? This regressions tests OK for me on x86_64 so it does not break things here. Okay for mainline from my perspective Jerry > > Dave > --- > > libgfortran: Add additional long double math support for hppa64-hpux* > > On hppa-hpux, we lack support for C99 math routines but the long > double type conforms to the standard 16-byte IEEE (IEC 60559) > specification. This allows us to use the quad routines from > libquadmath. > > 2026-07-22 John David Anglin <[email protected]> > > libgfortran/ChangeLog: > > PR libfortran/126305 > * c99_protos.h (acosl, acoshl, asinl, asinhl,atan2l, atanl, > atanhl, ceill, cosl, coshl, expl, hypotl, logl, sinl, sinhl, > sqrtl, tanl, tanhl, truncl, nextafterl, powl, erfl, erfcl): > Add declarations. > * configure.ac (powl, erfl): Add GCC_CHECK_MATH_FUNCs. > * configure: Regenerate. > * config.h.in: Regenerate. > * intrinsics/c99_functions.c: Include math_imp.h. Move code > to manipulate the floating-point type to math_imp.h. > (erfl, erfcl, acosl, acoshl, asinl, asinhl, atan2l, atanl, > atanhl, ceill, cosl, coshl, expl, floorl, fmodl, hypotl, > logl, sinl, sinhl, sqrtl, tanl, tanhl, truncl, nextafterl, > powl, log10l): Implement. > (frexpl): Move code to intrinsics/frexpl_16.c. > (scalbnl): Move code to intrinsics/scalbnl_16.c. > * intrinsics/ceill_16.c: New. > * intrinsics/floorl_16.c: New. > * intrinsics/frexpl_16.c: New. > * intrinsics/math_imp.h: New. > * intrinsics/nextafterl_16.c: New. > * intrinsics/scalbnl_16.c: New. > * intrinsics/truncl_16.c: New. > * libgfortran.h (USE_LIBQUADLIB): Define on hppa-hpux. > > diff --git a/libgfortran/c99_protos.h b/libgfortran/c99_protos.h > index 952b251c31a..a233991b60d 100644 > --- a/libgfortran/c99_protos.h > +++ b/libgfortran/c99_protos.h > @@ -37,41 +37,84 @@ extern size_t strnlen(const char *, size_t); > extern float acosf(float); > #endif > > +/* On HPUX, some long double functions are mapped to functions in > + libquadmath, e.g., acosl(x) maps to acosq((__float128)x). */ > + > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ACOS)) && !defined(HAVE_ACOSL) > +#define HAVE_ACOSL 1 > +extern long double acosl(long double); > +#endif > + > #if HAVE_ACOSH && !HAVE_ACOSHF > #define HAVE_ACOSHF 1 > extern float acoshf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ACOS)) && !defined(HAVE_ACOSHL) > +#define HAVE_ACOSHL 1 > +extern long double acoshl(long double); > +#endif > + > #ifndef HAVE_ASINF > #define HAVE_ASINF 1 > extern float asinf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ASIN)) && !defined(HAVE_ASINL) > +#define HAVE_ASINL 1 > +extern long double asinl(long double); > +#endif > + > #if HAVE_ASINH && !HAVE_ASINHF > #define HAVE_ASINHF 1 > extern float asinhf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ASINH)) && !defined(HAVE_ASINHL) > +#define HAVE_ASINHL 1 > +extern long double asinhl(long double); > +#endif > + > #ifndef HAVE_ATAN2F > #define HAVE_ATAN2F 1 > extern float atan2f(float, float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATAN2)) && !defined(HAVE_ATAN2L) > +#define HAVE_ATAN2L 1 > +extern long double atan2l(long double, long double); > +#endif > + > #ifndef HAVE_ATANF > #define HAVE_ATANF 1 > extern float atanf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATAN)) && !defined(HAVE_ATANL) > +#define HAVE_ATANL 1 > +extern long double atanl(long double); > +#endif > + > #if HAVE_ATANH && !HAVE_ATANHF > #define HAVE_ATANHF 1 > extern float atanhf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATANH)) && !defined(HAVE_ATANHL) > +#define HAVE_ATANHL 1 > +extern long double atanhl (long double); > +#endif > + > #ifndef HAVE_CEILF > #define HAVE_CEILF 1 > extern float ceilf(float); > #endif > > +#if !defined(HAVE_CEILL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > +#define HAVE_CEILL 1 > +extern long double ceill(long double); > +#endif > + > #ifndef HAVE_COPYSIGNF > #define HAVE_COPYSIGNF 1 > extern float copysignf(float, float); > @@ -92,16 +135,31 @@ extern long double copysignl(long double, long double); > extern float cosf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_COS)) && !defined(HAVE_COSL) > +#define HAVE_COSL 1 > +extern long double cosl (long double); > +#endif > + > #ifndef HAVE_COSHF > #define HAVE_COSHF 1 > extern float coshf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_COSH)) && !defined(HAVE_COSHL) > +#define HAVE_COSHL 1 > +extern long double coshl(long double); > +#endif > + > #ifndef HAVE_EXPF > #define HAVE_EXPF 1 > extern float expf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_EXP)) && !defined(HAVE_EXPL) > +#define HAVE_EXPL 1 > +extern long double expl(long double); > +#endif > + > #ifndef HAVE_FABSF > #define HAVE_FABSF 1 > extern float fabsf(float); > @@ -142,7 +200,7 @@ extern long double fmodl (long double x, long double y); > extern float frexpf(float, int *); > #endif > > -#ifndef HAVE_FREXPL > +#if !defined(HAVE_FREXPL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > #define HAVE_FREXPL 1 > extern long double frexpl(long double, int *); > #endif > @@ -152,11 +210,21 @@ extern long double frexpl(long double, int *); > extern float hypotf(float, float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_HYPOT)) && !defined(HAVE_HYPOTL) > +#define HAVE_HYPOTL 1 > +extern long double hypotl(long double, long double); > +#endif > + > #ifndef HAVE_LOGF > #define HAVE_LOGF 1 > extern float logf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_LOG)) && !defined(HAVE_LOGL) > +#define HAVE_LOGL 1 > +extern long double logl(long double); > +#endif > + > #ifndef HAVE_LOG10F > #define HAVE_LOG10F 1 > extern float log10f(float); > @@ -172,7 +240,7 @@ extern double scalbn(double, int); > extern float scalbnf(float, int); > #endif > > -#ifndef HAVE_SCALBNL > +#if !defined(HAVE_SCALBNL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > #define HAVE_SCALBNL 1 > extern long double scalbnl(long double, int); > #endif > @@ -182,26 +250,51 @@ extern long double scalbnl(long double, int); > extern float sinf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_SIN)) && !defined(HAVE_SINL) > +#define HAVE_SINL 1 > +extern long double sinl(long double); > +#endif > + > #ifndef HAVE_SINHF > #define HAVE_SINHF 1 > extern float sinhf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_SINH)) && !defined(HAVE_SINHL) > +#define HAVE_SINHL 1 > +extern long double sinhl(long double); > +#endif > + > #ifndef HAVE_SQRTF > #define HAVE_SQRTF 1 > extern float sqrtf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_SQRT)) && !defined(HAVE_SQRTL) > +#define HAVE_SQRTL 1 > +extern long double sqrtl(long double); > +#endif > + > #ifndef HAVE_TANF > #define HAVE_TANF 1 > extern float tanf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_TAN)) && !defined(HAVE_TANL) > +#define HAVE_TANL 1 > +extern long double tanl(long double); > +#endif > + > #ifndef HAVE_TANHF > #define HAVE_TANHF 1 > extern float tanhf(float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_TANH)) && !defined(HAVE_TANHL) > +#define HAVE_TANHL 1 > +extern long double tanhl(long double); > +#endif > + > #ifndef HAVE_TRUNC > #define HAVE_TRUNC 1 > extern double trunc(double); > @@ -212,16 +305,31 @@ extern double trunc(double); > extern float truncf(float); > #endif > > +#if !defined(HAVE_TRUNCL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > +#define HAVE_TRUNCL 1 > +extern long double truncl(long double); > +#endif > + > #ifndef HAVE_NEXTAFTERF > #define HAVE_NEXTAFTERF 1 > extern float nextafterf(float, float); > #endif > > +#if !defined(HAVE_NEXTAFTERL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > +#define HAVE_NEXTAFTERL 1 > +extern long double nextafterl(long double, long double); > +#endif > + > #ifndef HAVE_POWF > #define HAVE_POWF 1 > extern float powf(float, float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_POW)) && !defined(HAVE_POWL) > +#define HAVE_POWL 1 > +extern long double powl(long double, long double); > +#endif > + > #ifndef HAVE_ROUND > #define HAVE_ROUND 1 > extern double round(double); > @@ -310,11 +418,21 @@ extern float ynf (int, float); > extern float erff (float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ERF)) && !defined(HAVE_ERFL) > +#define HAVE_ERFL 1 > +long double erfl(long double); > +#endif > + > #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF) > #define HAVE_ERFCF 1 > extern float erfcf (float); > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ERFC)) && !defined(HAVE_ERFCL) > +#define HAVE_ERFL 1 > +long double erfcl(long double); > +#endif > + > > > /* log10l is needed on all platforms for decimal I/O */ > diff --git a/libgfortran/configure.ac b/libgfortran/configure.ac > index 9acb323469e..c5f84a031cd 100644 > --- a/libgfortran/configure.ac > +++ b/libgfortran/configure.ac > @@ -470,6 +470,7 @@ GCC_CHECK_MATH_FUNC([nextafter]) > GCC_CHECK_MATH_FUNC([nextafterl]) > GCC_CHECK_MATH_FUNC([powf]) > GCC_CHECK_MATH_FUNC([pow]) > +GCC_CHECK_MATH_FUNC([powl]) > GCC_CHECK_MATH_FUNC([cpowf]) > GCC_CHECK_MATH_FUNC([cpow]) > GCC_CHECK_MATH_FUNC([cpowl]) > @@ -520,6 +521,7 @@ GCC_CHECK_MATH_FUNC([trunc]) > GCC_CHECK_MATH_FUNC([truncl]) > GCC_CHECK_MATH_FUNC([erff]) > GCC_CHECK_MATH_FUNC([erf]) > +GCC_CHECK_MATH_FUNC([erfl]) > GCC_CHECK_MATH_FUNC([erfcf]) > GCC_CHECK_MATH_FUNC([erfc]) > GCC_CHECK_MATH_FUNC([erfcl]) > diff --git a/libgfortran/intrinsics/c99_functions.c b/libgfortran/intrinsics/c99_functions.c > index 28acadb5f9f..ab3920702b9 100644 > --- a/libgfortran/intrinsics/c99_functions.c > +++ b/libgfortran/intrinsics/c99_functions.c > @@ -26,6 +26,7 @@ see the files COPYING3 and COPYING.RUNTIME respectively. If not, see > > #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW > #include "libgfortran.h" > +#include "math_imp.h" > > /* On a C99 system "I" (with I*I = -1) should be defined in complex.h; > if not, we define a fallback version here. */ > @@ -48,138 +49,6 @@ see the files COPYING3 and COPYING.RUNTIME respectively. If not, see > /* Prototypes are included to silence -Wstrict-prototypes > -Wmissing-prototypes. */ > > -/* Main union type we use to manipulate the floating-point type. */ > -typedef union > -{ > - long double value; > - > - struct > -#ifdef __MINGW32__ > - /* On mingw targets the ms-bitfields option is active by default. > - Therefore enforce gnu-bitfield style. */ > - __attribute__ ((gcc_struct)) > -#endif > - { > -#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > - unsigned negative:1; > - unsigned exponent:15; > - unsigned mantissa0:16; > - unsigned mantissa1:32; > - unsigned mantissa2:32; > - unsigned mantissa3:32; > -#else > - unsigned mantissa3:32; > - unsigned mantissa2:32; > - unsigned mantissa1:32; > - unsigned mantissa0:16; > - unsigned exponent:15; > - unsigned negative:1; > -#endif > - } ieee; > - > - struct > - { > -#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > - uint64_t high; > - uint64_t low; > -#else > - uint64_t low; > - uint64_t high; > -#endif > - } words64; > - > - struct > - { > -#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > - uint32_t w0; > - uint32_t w1; > - uint32_t w2; > - uint32_t w3; > -#else > - uint32_t w3; > - uint32_t w2; > - uint32_t w1; > - uint32_t w0; > -#endif > - } words32; > - > - struct > -#ifdef __MINGW32__ > - /* Make sure we are using gnu-style bitfield handling. */ > - __attribute__ ((gcc_struct)) > -#endif > - { > -#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > - unsigned negative:1; > - unsigned exponent:15; > - unsigned quiet_nan:1; > - unsigned mantissa0:15; > - unsigned mantissa1:32; > - unsigned mantissa2:32; > - unsigned mantissa3:32; > -#else > - unsigned mantissa3:32; > - unsigned mantissa2:32; > - unsigned mantissa1:32; > - unsigned mantissa0:15; > - unsigned quiet_nan:1; > - unsigned exponent:15; > - unsigned negative:1; > -#endif > - } ieee_nan; > - > -} ieee754_long_double; > - > - > -/* Get two 64 bit ints from a long double. */ > -#define GET_LDOUBLE_WORDS64(ix0,ix1,d) \ > -do { \ > - ieee754_long_double u; \ > - u.value = (d); \ > - (ix0) = u.words64.high; \ > - (ix1) = u.words64.low; \ > -} while (0) > - > -/* Set a long double from two 64 bit ints. */ > -#define SET_LDOUBLE_WORDS64(d,ix0,ix1) \ > -do { \ > - ieee754_long_double u; \ > - u.words64.high = (ix0); \ > - u.words64.low = (ix1); \ > - (d) = u.value; \ > -} while (0) > - > -/* Get the more significant 64 bits of a long double mantissa. */ > -#define GET_LDOUBLE_MSW64(v,d) \ > -do { \ > - ieee754_long_double u; \ > - u.value = (d); \ > - (v) = u.words64.high; \ > -} while (0) > - > -/* Set the more significant 64 bits of a long double mantissa from an int. */ > -#define SET_LDOUBLE_MSW64(d,v) \ > -do { \ > - ieee754_long_double u; \ > - u.value = (d); \ > - u.words64.high = (v); \ > - (d) = u.value; \ > -} while (0) > - > -/* Get the least significant 64 bits of a long double mantissa. */ > -#define GET_LDOUBLE_LSW64(v,d) \ > -do { \ > - ieee754_long_double u; \ > - u.value = (d); \ > - (v) = u.words64.low; \ > -} while (0) > - > -static const long double > -two114 = 2.0769187434139310514121985316880384E+34L, /* 0x4071000000000000, 0 */ > -twom114 = 4.8148248609680896326399448564623183E-35L, /* 0x3F8D000000000000, 0 */ > -huge = 1.0E+4900L, > -tiny = 1.0E-4900L; > - > /* Wrapper for systems without strnlen function. */ > > #ifndef HAVE_STRNLEN > @@ -276,6 +145,25 @@ erff (float x) > } > #endif > > +/* On HPUX, some long double functions are mapped to functions in > + libquadmath, e.g., erfl(x) maps to erfq((__float128)x). */ > + > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ERF)) && !defined(HAVE_ERFL) > +#define HAVE_ERFL 1 > +long double erfl (long double); > + > +long double > +erfl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) erfq ((__float128) x); > +#else > + return (long double) erf ((double) x); > +#endif > +} > +#endif > + > + > #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF) > #define HAVE_ERFCF 1 > float erfcf (float); > @@ -287,6 +175,21 @@ erfcf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ERFC)) && !defined(HAVE_ERFCL) > +#define HAVE_ERFCL 1 > +long double erfcl (long double); > + > +long double > +erfcl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) erfcq ((__float128) x); > +#else > + return (long double) erfc ((double) x); > +#endif > +} > +#endif > + > > #ifndef HAVE_ACOSF > #define HAVE_ACOSF 1 > @@ -299,6 +202,21 @@ acosf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ACOS)) && !defined(HAVE_ACOSL) > +#define HAVE_ACOSL 1 > +long double acosl (long double); > + > +long double > +acosl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) acosq ((__float128) x); > +#else > + return (long double) acos ((double) x); > +#endif > +} > +#endif > + > #if HAVE_ACOSH && !HAVE_ACOSHF > float acoshf (float x); > > @@ -309,6 +227,21 @@ acoshf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ACOS)) && !defined(HAVE_ACOSHL) > +#define HAVE_ACOSHL 1 > +long double acoshl (long double); > + > +long double > +acoshl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) acoshq ((__float128) x); > +#else > + return (long double) acosh ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_ASINF > #define HAVE_ASINF 1 > float asinf (float x); > @@ -320,6 +253,21 @@ asinf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ASIN)) && !defined(HAVE_ASINL) > +#define HAVE_ASINL 1 > +long double asinl(long double); > + > +long double > +asinl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) asinq ((__float128) x); > +#else > + return (long double) asin ((double) x); > +#endif > +} > +#endif > + > #if HAVE_ASINH && !HAVE_ASINHF > float asinhf (float x); > > @@ -330,6 +278,21 @@ asinhf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ASINH)) && !defined(HAVE_ASINHL) > +#define HAVE_ASINHL 1 > +long double asinhl(long double); > + > +long double > +asinhl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) asinhq ((__float128) x); > +#else > + return (long double) asinh ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_ATAN2F > #define HAVE_ATAN2F 1 > float atan2f (float y, float x); > @@ -341,6 +304,21 @@ atan2f (float y, float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATAN2)) && !defined(HAVE_ATAN2L) > +#define HAVE_ATAN2L 1 > +long double atan2l(long double, long double); > + > +long double > +atan2l (long double x, long double y) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) atan2q ((__float128) x, (__float128) y); > +#else > + return (long double) atan2 ((double) x, (double) y); > +#endif > +} > +#endif > + > #ifndef HAVE_ATANF > #define HAVE_ATANF 1 > float atanf (float x); > @@ -352,6 +330,21 @@ atanf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATAN)) && !defined(HAVE_ATANL) > +#define HAVE_ATANL 1 > +long double atanl (long double); > + > +long double > +atanl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) atanq ((__float128) x); > +#else > + return (long double) atan ((double) x); > +#endif > +} > +#endif > + > #if HAVE_ATANH && !HAVE_ATANHF > float atanhf (float x); > > @@ -362,6 +355,21 @@ atanhf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATANH)) && !defined(HAVE_ATANHL) > +#define HAVE_ATANHL 1 > +long double atanhl (long double); > + > +long double > +atanhl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) atanhq ((__float128) x); > +#else > + return (long double) atanh ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_CEILF > #define HAVE_CEILF 1 > float ceilf (float x); > @@ -373,6 +381,13 @@ ceilf (float x) > } > #endif > > +#if !defined(HAVE_CEILL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > +#define HAVE_CEILL 1 > +long double ceill (long double); > +#include "ceill_16.c" > +#endif > + > + > #if !defined(HAVE_COPYSIGN) && defined(HAVE_INLINE_BUILTIN_COPYSIGN) > #define HAVE_COPYSIGN 1 > double copysign (double x, double y); > @@ -417,6 +432,21 @@ cosf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_COS)) && !defined(HAVE_COSL) > +#define HAVE_COSL 1 > +long double cosl (long double); > + > +long double > +cosl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) cosq ((__float128) x); > +#else > + return (long double) cos ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_COSHF > #define HAVE_COSHF 1 > float coshf (float x); > @@ -428,6 +458,21 @@ coshf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_COSH)) && !defined(HAVE_COSHL) > +#define HAVE_COSHL 1 > +long double coshl (long double); > + > +long double > +coshl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) coshq ((__float128) x); > +#else > + return (long double) cosh ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_EXPF > #define HAVE_EXPF 1 > float expf (float x); > @@ -439,6 +484,21 @@ expf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_EXP)) && !defined(HAVE_EXPL) > +#define HAVE_EXPL 1 > +long double expl (long double); > + > +long double > +expl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) expq ((__float128) x); > +#else > + return (long double) exp ((double) x); > +#endif > +} > +#endif > + > #if !defined(HAVE_FABS) && defined(HAVE_INLINE_BUILTIN_FABS) > #define HAVE_FABS 1 > double fabs (double x); > @@ -483,6 +543,37 @@ floorf (float x) > } > #endif > > +#if !defined(HAVE_FLOORL) > +#define HAVE_FLOORL 1 > +long double floorl (long double); > + > +#if (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > +#include "floorl_16.c" > +#else > +long double > +floorl (long double x); > +{ > + /* Zero, possibly signed. */ > + if (x == 0) > + return x; > + > + /* Large magnitude. */ > + if (x > DBL_MAX || x < (-DBL_MAX)) > + return x; > + > + /* Small positive values. */ > + if (x >= 0 && x < DBL_MIN) > + return 0; > + > + /* Small negative values. */ > + if (x < 0 && x > (-DBL_MIN)) > + return -1; > + > + return floor (x); > +} > +#endif > +#endif > + > #ifndef HAVE_FMODF > #define HAVE_FMODF 1 > float fmodf (float x, float y); > @@ -494,6 +585,26 @@ fmodf (float x, float y) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_FLOORL)) && !defined(HAVE_FMODL) > +#define HAVE_FMODL 1 > +long double fmodl (long double, long double); > + > +long double > +fmodl (long double x, long double y) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) fmodq ((__float128) x, (__float128) y); > +#else > + if (y == 0.0L) > + return 0.0L; > + > + /* Need to check that the result has the same sign as x and magnitude > + less than the magnitude of y. */ > + return x - floorl (x / y) * y; > +#endif > +} > +#endif > + > #ifndef HAVE_FREXPF > #define HAVE_FREXPF 1 > float frexpf (float x, int *exp); > @@ -505,55 +616,10 @@ frexpf (float x, int *exp) > } > #endif > > -#if !defined(HAVE_FREXPL) > +#if !defined(HAVE_FREXPL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > #define HAVE_FREXPL 1 > -long double frexpl (long double x, int *eptr); > - > -/* s_frexpl.c -- long double version of s_frexp.c. > - * Conversion to IEEE quad long double by Jakub Jelinek, [email protected]. > - */ > - > -/* > - * ==================================================== > - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. > - * > - * Developed at SunPro, a Sun Microsystems, Inc. business. > - * Permission to use, copy, modify, and distribute this > - * software is freely granted, provided that this notice > - * is preserved. > - * ==================================================== > - */ > - > -long double > -frexpl (long double x, int *eptr) > -{ > - if (sizeof (long double) == 16) > - { > - uint64_t hx, lx, ix; > - > - GET_LDOUBLE_WORDS64(hx,lx,x); > - ix = 0x7fffffffffffffffULL&hx; > - *eptr = 0; > - if(ix >= 0x7fff000000000000ULL || ((ix|lx) == 0)) > - /* 0,inf,nan */ > - return x + x; > - if (ix < 0x0001000000000000ULL) > - { > - /* subnormal */ > - x *= two114; > - GET_LDOUBLE_MSW64(hx,x); > - ix = hx & 0x7fffffffffffffffULL; > - *eptr = -114; > - } > - *eptr += (ix>>48) - 16382; > - hx = (hx & 0x8000ffffffffffffULL) | 0x3ffe000000000000ULL; > - SET_LDOUBLE_MSW64(x,hx); > - return x; > - } > - else > - /* Intel 80 bit */ > - abort(); > -} > +long double frexpl (long double, int *); > +#include "frexpl_16.c" > #endif > > #ifndef HAVE_HYPOTF > @@ -567,6 +633,21 @@ hypotf (float x, float y) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_HYPOT)) && !defined(HAVE_HYPOTL) > +#define HAVE_HYPOTL 1 > +long double hypotl (long double, long double); > + > +long double > +hypotl (long double x, long double y) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) hypotq ((__float128) x, (__float128) y); > +#else > + return (long double) hypot ((double) x, (double) y); > +#endif > +} > +#endif > + > #ifndef HAVE_LOGF > #define HAVE_LOGF 1 > float logf (float x); > @@ -578,6 +659,21 @@ logf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_LOG)) && !defined(HAVE_LOGL) > +#define HAVE_LOGL 1 > +long double logl (long double); > + > +long double > +logl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) logq ((__float128) x); > +#else > + return (long double) log ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_LOG10F > #define HAVE_LOG10F 1 > float log10f (float x); > @@ -615,74 +711,10 @@ scalbnf (float x, int y) > } > #endif > > -#if !defined(HAVE_SCALBNL) > +#if !defined(HAVE_SCALBNL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > #define HAVE_SCALBNL 1 > -long double scalbnl (long double x, int n); > - > -/* s_scalbnl.c -- long double version of s_scalbn.c. > - * Conversion to IEEE quad long double by Jakub Jelinek, [email protected]. > - */ > - > -/* @(#)s_scalbn.c 5.1 93/09/24 */ > -/* > - * ==================================================== > - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. > - * > - * Developed at SunPro, a Sun Microsystems, Inc. business. > - * Permission to use, copy, modify, and distribute this > - * software is freely granted, provided that this notice > - * is preserved. > - * ==================================================== > - */ > - > -long double > -scalbnl (long double x, int n) > -{ > - if (sizeof (long double) == 16) > - { > - int64_t k, hx, lx; > - > - GET_LDOUBLE_WORDS64(hx,lx,x); > - > - /* extract exponent */ > - k = (hx >> 48) & 0x7fff; > - if (k == 0) > - { > - /* 0 or subnormal x */ > - if ((lx | (hx & 0x7fffffffffffffffULL)) == 0) > - return x; /* +-0 */ > - x *= two114; > - GET_LDOUBLE_MSW64(hx,x); > - k = ((hx >> 48) & 0x7fff) - 114; > - } > - if (k == 0x7fff) > - /* NaN or Inf */ > - return x+x; > - if (n < -50000) > - /*underflow*/ > - return tiny * copysignl (tiny, x); > - if (n > 50000 || k + n > 0x7ffe) > - /* overflow */ > - return huge * copysignl (huge, x); > - /* Now k and n are bounded we know that k = k+n does not overflow. */ > - k = k + n; > - if (k > 0) > - { > - /* normal result */ > - SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); > - return x; > - } > - if (k <= -114) > - /*underflow*/ > - return tiny * copysignl (tiny, x); > - k += 114; /* subnormal result */ > - SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); > - return x * twom114; > - } > - else > - /* Intel 80 bit */ > - abort(); > -} > +long double scalbnl (long double, int); > +#include "scalbnl_16.c" > #endif > > #ifndef HAVE_SINF > @@ -696,6 +728,21 @@ sinf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_SIN)) && !defined(HAVE_SINL) > +#define HAVE_SINL 1 > +long double sinl(long double); > + > +long double > +sinl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) sinq ((__float128) x); > +#else > + return (long double) sin ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_SINHF > #define HAVE_SINHF 1 > float sinhf (float x); > @@ -707,6 +754,21 @@ sinhf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_SINH)) && !defined(HAVE_SINHL) > +#define HAVE_SINHL 1 > +long double sinhl(long double); > + > +long double > +sinhl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) sinhq ((__float128) x); > +#else > + return (long double) sinh ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_SQRTF > #define HAVE_SQRTF 1 > float sqrtf (float x); > @@ -718,6 +780,21 @@ sqrtf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_SQRT)) && !defined(HAVE_SQRTL) > +#define HAVE_SQRTL 1 > +long double sqrtl(long double); > + > +long double > +sqrtl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) sqrtq ((__float128) x); > +#else > + return (long double) sqrt ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_TANF > #define HAVE_TANF 1 > float tanf (float x); > @@ -729,6 +806,21 @@ tanf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_TAN)) && !defined(HAVE_TANL) > +#define HAVE_TANL 1 > +long double tanl(long double); > + > +long double > +tanl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) tanq ((__float128) x); > +#else > + return (long double) tan ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_TANHF > #define HAVE_TANHF 1 > float tanhf (float x); > @@ -740,6 +832,21 @@ tanhf (float x) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_TANH)) && !defined(HAVE_TANHL) > +#define HAVE_TANHL 1 > +long double tanhl(long double); > + > +long double > +tanhl (long double x) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) tanhq ((__float128) x); > +#else > + return (long double) tanh ((double) x); > +#endif > +} > +#endif > + > #ifndef HAVE_TRUNC > #define HAVE_TRUNC 1 > double trunc (double x); > @@ -768,6 +875,12 @@ truncf (float x) > } > #endif > > +#if !defined(HAVE_TRUNCL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > +#define HAVE_TRUNCL 1 > +long double truncl(long double); > +#include "truncl_16.c" > +#endif > + > #ifndef HAVE_NEXTAFTERF > #define HAVE_NEXTAFTERF 1 > /* This is a portable implementation of nextafterf that is intended to be > @@ -833,6 +946,11 @@ nextafterf (float x, float y) > } > #endif > > +#if !defined(HAVE_NEXTAFTERL) && (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > +#define HAVE_NEXTAFTERL 1 > +long double nextafterl (long double, long double); > +#include "nextafterl_16.c" > +#endif > > #ifndef HAVE_POWF > #define HAVE_POWF 1 > @@ -845,6 +963,21 @@ powf (float x, float y) > } > #endif > > +#if (defined(USE_LIBQUADLIB) || defined(HAVE_POW)) && !defined(HAVE_POWL) > +#define HAVE_POWL 1 > +long double powl (long double, long double); > + > +long double > +powl (long double x, long double y) > +{ > +#if defined(USE_LIBQUADLIB) > + return (long double) powq ((__float128) x, (__float128) y); > +#else > + return (long double) pow ((double) x, (double) y); > +#endif > +} > +#endif > + > > #ifndef HAVE_ROUND > #define HAVE_ROUND 1 > @@ -1041,6 +1174,9 @@ long double log10l (long double x); > long double > log10l (long double x) > { > +#if defined(USE_LIBQUADLIB) > + return (long double) log10q ((__float128) x); > +#else > #if LDBL_MAX_EXP > DBL_MAX_EXP > if (x > DBL_MAX) > { > @@ -1070,51 +1206,7 @@ log10l (long double x) > } > #endif > return log10 (x); > -} > #endif > - > - > -#ifndef HAVE_FLOORL > -#define HAVE_FLOORL 1 > -long double floorl (long double x); > - > -long double > -floorl (long double x) > -{ > - /* Zero, possibly signed. */ > - if (x == 0) > - return x; > - > - /* Large magnitude. */ > - if (x > DBL_MAX || x < (-DBL_MAX)) > - return x; > - > - /* Small positive values. */ > - if (x >= 0 && x < DBL_MIN) > - return 0; > - > - /* Small negative values. */ > - if (x < 0 && x > (-DBL_MIN)) > - return -1; > - > - return floor (x); > -} > -#endif > - > - > -#ifndef HAVE_FMODL > -#define HAVE_FMODL 1 > -long double fmodl (long double x, long double y); > - > -long double > -fmodl (long double x, long double y) > -{ > - if (y == 0.0L) > - return 0.0L; > - > - /* Need to check that the result has the same sign as x and magnitude > - less than the magnitude of y. */ > - return x - floorl (x / y) * y; > } > #endif > > diff --git a/libgfortran/intrinsics/ceill_16.c b/libgfortran/intrinsics/ceill_16.c > new file mode 100644 > index 00000000000..40dd1fa09a6 > --- /dev/null > +++ b/libgfortran/intrinsics/ceill_16.c > @@ -0,0 +1,91 @@ > +/* s_ceill.c -- long double version of s_ceil.c. > + * Conversion to IEEE quad long double by Jakub Jelinek, [email protected]. > + */ > + > +/* > + * ==================================================== > + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. > + * > + * Developed at SunPro, a Sun Microsystems, Inc. business. > + * Permission to use, copy, modify, and distribute this > + * software is freely granted, provided that this notice > + * is preserved. > + * ==================================================== > + */ > + > +/* > + * ceill(x) > + * Return x rounded toward -inf to integral value > + * Method: > + * Bit twiddling. > + */ > + > +long double > +ceill (long double x) > +{ > + int64_t i0, i1, j0; > + uint64_t i, j; > + > + GET_LDOUBLE_WORDS64 (i0 ,i1 ,x); > + j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; > + if (j0 < 48) > + { > + if (j0 < 0) > + { > + /* return 0*sign(x) if |x|<1 */ > + if (i0 < 0) > + { > + i0 = 0x8000000000000000ULL; > + i1 = 0; > + } > + else if ((i0|i1) != 0) > + { > + i0 = 0x3fff000000000000ULL; > + i1 = 0; > + } > + } > + else > + { > + i = (0x0000ffffffffffffULL) >> j0; > + if (((i0 & i) | i1) == 0) > + /* x is integral */ > + return x; > + if (i0 > 0) > + i0 += (0x0001000000000000LL) >> j0; > + i0 &= (~i); > + i1 = 0; > + } > + } > + else if (j0 > 111) > + { > + if (j0 == 0x4000) > + /* inf or NaN */ > + return x+x; > + else > + /* x is integral */ > + return x; > + } > + else > + { > + i = -1ULL >> (j0 - 48); > + if ((i1 & i) == 0) > + /* x is integral */ > + return x; > + if (i0 > 0) > + { > + if (j0 == 48) > + i0 += 1; > + else > + { > + j = i1 + (1LL << (112 - j0)); > + if ((int64_t) j < i1) > + /* got a carry */ > + i0 += 1; > + i1 = j; > + } > + } > + i1 &= (~i); > + } > + SET_LDOUBLE_WORDS64 (x, i0, i1); > + return x; > +} > diff --git a/libgfortran/intrinsics/floorl_16.c b/libgfortran/intrinsics/floorl_16.c > new file mode 100644 > index 00000000000..d07fee1dc07 > --- /dev/null > +++ b/libgfortran/intrinsics/floorl_16.c > @@ -0,0 +1,88 @@ > +/* s_floorl.c -- long double version of s_floor.c. > + * Conversion to IEEE quad long double by Jakub Jelinek, [email protected]. > + */ > + > +/* > + * ==================================================== > + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. > + * > + * Developed at SunPro, a Sun Microsystems, Inc. business. > + * Permission to use, copy, modify, and distribute this > + * software is freely granted, provided that this notice > + * is preserved. > + * ==================================================== > + */ > + > +/* > + * floorl(x) > + * Return x rounded toward -inf to integral value > + * Method: > + * Bit twiddling. > + */ > + > +long double > +floorl (long double x) > +{ > + int64_t i0, i1, j0; > + uint64_t i, j; > + > + GET_LDOUBLE_WORDS64 (i0, i1, x); > + j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; > + if (j0 < 48) > + { > + if (j0 < 0) > + { > + /* return 0*sign(x) if |x|<1 */ > + if (i0 >= 0) > + i0 = i1 = 0; > + else if (((i0 & 0x7fffffffffffffffLL) | i1) != 0) > + { > + i0 = 0xbfff000000000000ULL; > + i1 = 0; > + } > + } > + else > + { > + i = (0x0000ffffffffffffULL) >> j0; > + if (((i0 & i) | i1) == 0) > + /* x is integral */ > + return x; > + if (i0 < 0) > + i0 += (0x0001000000000000LL) >> j0; > + i0 &= (~i); > + i1 = 0; > + } > + } > + else if (j0 > 111) > + { > + if (j0 == 0x4000) > + /* inf or NaN */ > + return x + x; > + else > + /* x is integral */ > + return x; > + } > + else > + { > + i = -1ULL >> (j0 - 48); > + if ((i1 & i) == 0) > + /* x is integral */ > + return x; > + if (i0 < 0) > + { > + if (j0 == 48) > + i0 += 1; > + else > + { > + j = i1 + (1LL << (112 - j0)); > + if ((int64_t) j < i1) > + /* got a carry */ > + i0 += 1; > + i1 = j; > + } > + } > + i1 &= (~i); > + } > + SET_LDOUBLE_WORDS64 (x, i0, i1); > + return x; > +} > diff --git a/libgfortran/intrinsics/frexpl_16.c b/libgfortran/intrinsics/frexpl_16.c > new file mode 100644 > index 00000000000..7f6bcc97f6b > --- /dev/null > +++ b/libgfortran/intrinsics/frexpl_16.c > @@ -0,0 +1,39 @@ > +/* s_frexpl.c -- long double version of s_frexp.c. > + * Conversion to IEEE quad long double by Jakub Jelinek, [email protected]. > + */ > + > +/* > + * ==================================================== > + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. > + * > + * Developed at SunPro, a Sun Microsystems, Inc. business. > + * Permission to use, copy, modify, and distribute this > + * software is freely granted, provided that this notice > + * is preserved. > + * ==================================================== > + */ > + > +long double > +frexpl (long double x, int *eptr) > +{ > + uint64_t hx, lx, ix; > + > + GET_LDOUBLE_WORDS64 (hx, lx, x); > + ix = 0x7fffffffffffffffULL&hx; > + *eptr = 0; > + if (ix >= 0x7fff000000000000ULL || ((ix | lx) == 0)) > + /* 0,inf,nan */ > + return x + x; > + if (ix < 0x0001000000000000ULL) > + { > + /* subnormal */ > + x *= two114; > + GET_LDOUBLE_MSW64 (hx, x); > + ix = hx & 0x7fffffffffffffffULL; > + *eptr = -114; > + } > + *eptr += (ix >> 48) - 16382; > + hx = (hx & 0x8000ffffffffffffULL) | 0x3ffe000000000000ULL; > + SET_LDOUBLE_MSW64 (x, hx); > + return x; > +} > diff --git a/libgfortran/intrinsics/math_imp.h b/libgfortran/intrinsics/math_imp.h > new file mode 100644 > index 00000000000..fe55648dfde > --- /dev/null > +++ b/libgfortran/intrinsics/math_imp.h > @@ -0,0 +1,170 @@ > +/* Header file for bit manipulation of long double type > + Copyright (C) 2004-2026 Free Software Foundation, Inc. > + > +This file is part of the GNU Fortran 95 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/>. */ > + > +#ifndef LIBGFORTRAN_MATH_IMP_H > +#define LIBGFORTRAN_MATH_IMP_H > + > +#if (__SIZEOF_LONG_DOUBLE__ == 16) && (__LDBL_IS_IEC_60559__ > 0) > + > +/* Main union type we use to manipulate the floating-point type. */ > +typedef union > +{ > + long double value; > + > + struct > +#ifdef __MINGW32__ > + /* On mingw targets the ms-bitfields option is active by default. > + Therefore enforce gnu-bitfield style. */ > + __attribute__ ((gcc_struct)) > +#endif > + { > +#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > + unsigned negative:1; > + unsigned exponent:15; > + unsigned mantissa0:16; > + unsigned mantissa1:32; > + unsigned mantissa2:32; > + unsigned mantissa3:32; > +#else > + unsigned mantissa3:32; > + unsigned mantissa2:32; > + unsigned mantissa1:32; > + unsigned mantissa0:16; > + unsigned exponent:15; > + unsigned negative:1; > +#endif > + } ieee; > + > + struct > + { > +#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > + uint64_t high; > + uint64_t low; > +#else > + uint64_t low; > + uint64_t high; > +#endif > + } words64; > + > + struct > + { > +#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > + uint32_t w0; > + uint32_t w1; > + uint32_t w2; > + uint32_t w3; > +#else > + uint32_t w3; > + uint32_t w2; > + uint32_t w1; > + uint32_t w0; > +#endif > + } words32; > + > + struct > +#ifdef __MINGW32__ > + /* Make sure we are using gnu-style bitfield handling. */ > + __attribute__ ((gcc_struct)) > +#endif > + { > +#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ > + unsigned negative:1; > + unsigned exponent:15; > + unsigned quiet_nan:1; > + unsigned mantissa0:15; > + unsigned mantissa1:32; > + unsigned mantissa2:32; > + unsigned mantissa3:32; > +#else > + unsigned mantissa3:32; > + unsigned mantissa2:32; > + unsigned mantissa1:32; > + unsigned mantissa0:15; > + unsigned quiet_nan:1; > + unsigned exponent:15; > + unsigned negative:1; > +#endif > + } ieee_nan; > + > +} ieee754_long_double; > + > + > +/* Get two 64 bit ints from a long double. */ > +#define GET_LDOUBLE_WORDS64(ix0,ix1,d) \ > +do { \ > + ieee754_long_double u; \ > + u.value = (d); \ > + (ix0) = u.words64.high; \ > + (ix1) = u.words64.low; \ > +} while (0) > + > +/* Set a long double from two 64 bit ints. */ > +#define SET_LDOUBLE_WORDS64(d,ix0,ix1) \ > +do { \ > + ieee754_long_double u; \ > + u.words64.high = (ix0); \ > + u.words64.low = (ix1); \ > + (d) = u.value; \ > +} while (0) > + > +/* Get the more significant 64 bits of a long double mantissa. */ > +#define GET_LDOUBLE_MSW64(v,d) \ > +do { \ > + ieee754_long_double u; \ > + u.value = (d); \ > + (v) = u.words64.high; \ > +} while (0) > + > +/* Set the more significant 64 bits of a long double mantissa from an int. */ > +#define SET_LDOUBLE_MSW64(d,v) \ > +do { \ > + ieee754_long_double u; \ > + u.value = (d); \ > + u.words64.high = (v); \ > + (d) = u.value; \ > +} while (0) > + > +/* Get the least significant 64 bits of a long double mantissa. */ > +#define GET_LDOUBLE_LSW64(v,d) \ > +do { \ > + ieee754_long_double u; \ > + u.value = (d); \ > + (v) = u.words64.low; \ > +} while (0) > + > +static const long double __attribute__ ((unused)) > +two114 = 2.0769187434139310514121985316880384E+34L, /* 0x4071000000000000, 0 */ > +twom114 = 4.8148248609680896326399448564623183E-35L, /* 0x3F8D000000000000, 0 */ > +huge = 1.0E+4900L, > +tiny = 1.0E-4900L; > +#endif > + > +#ifndef math_opt_barrier > +# define math_opt_barrier(x) \ > +({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; }) > +# define math_force_eval(x) \ > +({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); }) > +#endif > + > +#endif /* LIBGFORTRAN_MATH_IMP_H */ > diff --git a/libgfortran/intrinsics/nextafterl_16.c b/libgfortran/intrinsics/nextafterl_16.c > new file mode 100644 > index 00000000000..8334da5e3e7 > --- /dev/null > +++ b/libgfortran/intrinsics/nextafterl_16.c > @@ -0,0 +1,108 @@ > +/* s_nextafterl.c -- long double version of s_nextafter.c. > + * Conversion to IEEE quad long double by Jakub Jelinek, [email protected]. > + */ > + > +/* > + * ==================================================== > + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. > + * > + * Developed at SunPro, a Sun Microsystems, Inc. business. > + * Permission to use, copy, modify, and distribute this > + * software is freely granted, provided that this notice > + * is preserved. > + * ==================================================== > + */ > + > +/* IEEE functions > + * nextafterl(x,y) > + * return the next machine floating-point number of x in the > + * direction toward y. > + * Special cases: > + */ > + > +long double > +nextafterl (long double x, long double y) > +{ > + int64_t hx, hy, ix, iy; > + uint64_t lx, ly; > + > + GET_LDOUBLE_WORDS64 (hx, lx, x); > + GET_LDOUBLE_WORDS64 (hy, ly, y); > + ix = hx & 0x7fffffffffffffffLL; /* |x| */ > + iy = hy & 0x7fffffffffffffffLL; /* |y| */ > + > + if (((ix >= 0x7fff000000000000LL) && ((ix - 0x7fff000000000000LL) | lx) != 0) > + || ((iy >= 0x7fff000000000000LL) && ((iy - 0x7fff000000000000LL) | ly) != 0)) > + /* x or y is nan */ > + return x + y; > + if (x == y) > + /* x=y, return y */ > + return y; > + if ((ix | lx) == 0) > + { > + /* x == 0 */ > + long double u; > + > + /* return +-minsubnormal */ > + SET_LDOUBLE_WORDS64 (x, hy & 0x8000000000000000ULL, 1); > + u = math_opt_barrier (x); > + u = u * u; > + /* raise underflow flag */ > + math_force_eval (u); > + return x; > + } > + if (hx >= 0) > + { > + /* x > 0 */ > + if (hx > hy || ((hx == hy) && (lx > ly))) > + { > + /* x > y, x -= ulp */ > + if (lx == 0) > + hx--; > + lx--; > + } > + else > + { > + /* x < y, x += ulp */ > + lx++; > + if (lx == 0) > + hx++; > + } > + } > + else > + { > + /* x < 0 */ > + if (hy >= 0 || hx > hy || ((hx == hy) && (lx > ly))) > + { > + /* x < y, x -= ulp */ > + if (lx == 0) > + hx--; > + lx--; > + } > + else > + { > + /* x > y, x += ulp */ > + lx++; > + if (lx == 0) > + hx++; > + } > + } > + hy = hx & 0x7fff000000000000LL; > + if (hy == 0x7fff000000000000LL) > + { > + /* overflow */ > + long double u = x + x; > + math_force_eval (u); > + errno = ERANGE; > + } > + if (hy == 0) > + { > + /* underflow */ > + long double u = x * x; > + /* raise underflow flag */ > + math_force_eval (u); > + errno = ERANGE; > + } > + SET_LDOUBLE_WORDS64 (x, hx, lx); > + return x; > +} > diff --git a/libgfortran/intrinsics/scalbnl_16.c b/libgfortran/intrinsics/scalbnl_16.c > new file mode 100644 > index 00000000000..717f4d6a148 > --- /dev/null > +++ b/libgfortran/intrinsics/scalbnl_16.c > @@ -0,0 +1,61 @@ > +/* s_scalbnl.c -- long double version of s_scalbn.c. > + * Conversion to IEEE quad long double by Jakub Jelinek, [email protected]. > + */ > + > +/* @(#)s_scalbn.c 5.1 93/09/24 */ > +/* > + * ==================================================== > + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. > + * > + * Developed at SunPro, a Sun Microsystems, Inc. business. > + * Permission to use, copy, modify, and distribute this > + * software is freely granted, provided that this notice > + * is preserved. > + * ==================================================== > + */ > + > +long double > +scalbnl (long double x, int n) > +{ > + int64_t k, hx, lx; > + > + GET_LDOUBLE_WORDS64 (hx, lx, x); > + > + /* extract exponent */ > + k = (hx >> 48) & 0x7fff; > + if (k == 0) > + { > + /* 0 or subnormal x */ > + if ((lx | (hx & 0x7fffffffffffffffULL)) == 0) > + return x; /* +-0 */ > + x *= two114; > + GET_LDOUBLE_MSW64 (hx, x); > + k = ((hx >> 48) & 0x7fff) - 114; > + } > + if (k == 0x7fff) > + /* NaN or Inf */ > + return x + x; > + if (n < -50000) > + /*underflow*/ > + return tiny * copysignl (tiny, x); > + if (n > 50000 || k + n > 0x7ffe) > + /* overflow */ > + return huge * copysignl (huge, x); > + > + /* Now k and n are bounded we know that k = k + n does not overflow. */ > + k = k + n; > + if (k > 0) > + { > + /* normal result */ > + SET_LDOUBLE_MSW64 (x, (hx&0x8000ffffffffffffULL) | (k << 48)); > + return x; > + } > + if (k <= -114) > + /* underflow */ > + return tiny * copysignl (tiny, x); > + > + /* subnormal result */ > + k += 114; > + SET_LDOUBLE_MSW64 (x, (hx & 0x8000ffffffffffffULL) | (k << 48)); > + return x * twom114; > +} > diff --git a/libgfortran/intrinsics/truncl_16.c b/libgfortran/intrinsics/truncl_16.c > new file mode 100644 > index 00000000000..5632f69bd31 > --- /dev/null > +++ b/libgfortran/intrinsics/truncl_16.c > @@ -0,0 +1,50 @@ > +/* Truncate argument to nearest integral value not larger than the argument. > + Copyright (C) 1997-2018 Free Software Foundation, Inc. > + This file is part of the GNU C Library. > + Contributed by Ulrich Drepper <[email protected]>, 1997 and > + Jakub Jelinek <[email protected]>, 1999. > + > + The GNU C Library 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. > + > + The GNU C Library 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 the GNU C Library; if not, see > + <http://www.gnu.org/licenses/>. */ > + > +long double > +truncl (long double x) > +{ > + int32_t j0; > + uint64_t i0, i1, sx; > + > + GET_LDOUBLE_WORDS64 (i0, i1, x); > + sx = i0 & 0x8000000000000000ULL; > + j0 = ((i0 >> 48) & 0x7fff) - 0x3fff; > + if (j0 < 48) > + { > + if (j0 < 0) > + /* The magnitude of the number is < 1 so the result is +-0. */ > + SET_LDOUBLE_WORDS64 (x, sx, 0); > + else > + SET_LDOUBLE_WORDS64 (x, i0 & ~(0x0000ffffffffffffLL >> j0), 0); > + } > + else if (j0 > 111) > + { > + if (j0 == 0x4000) > + /* x is inf or NaN. */ > + return x + x; > + } > + else > + { > + SET_LDOUBLE_WORDS64 (x, i0, i1 & ~(0xffffffffffffffffULL >> (j0 - 48))); > + } > + > + return x; > +} > diff --git a/libgfortran/libgfortran.h b/libgfortran/libgfortran.h > index 8381cbd7cd4..3613b5d68c0 100644 > --- a/libgfortran/libgfortran.h > +++ b/libgfortran/libgfortran.h > @@ -59,6 +59,14 @@ see the files COPYING3 and COPYING.RUNTIME respectively. If not, see > #include <string.h> > #endif > > +/* Use libquadlib math routines. HP-UX on PA-RISC uses the 16-byte > +IEEE format for long double but doesn't implement any of the standard > +mathmetical routines. However, since the long double and __float128 > +types are identical, we can use the routines in libquadmath. */ > +#if defined(__hpux__) && defined(__hppa__) > +#define USE_LIBQUADLIB > +#endif > + > #if HAVE_COMPLEX_H > /* Must appear before math.h on VMS systems. */ > # include <complex.h>