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>
