[PATCH]libgfortran: Add additional long double math support for hppa64-hpux*
John David Anglin <[email protected]> Thu, 23 Jul 2026 10:50:18 -0400
| Newsgroups | gmane.comp.gcc.patches,gmane.comp.gcc.fortran |
|---|---|
| Message-ID | <[email protected]> |
--9ziKhqtsx3XySAFl Content-Type: text/plain; charset=us-ascii Content-Disposition: inline Content-Transfer-Encoding: quoted-printable 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? 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.=20 (scalbnl): Move code to intrinsics/scalbnl_16.c.=20 * 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 =20 +/* 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 =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ACOS)) && !defined(HAVE_ACOSH= L) +#define HAVE_ACOSHL 1 +extern long double acoshl(long double); +#endif + #ifndef HAVE_ASINF #define HAVE_ASINF 1 extern float asinf(float); #endif =20 +#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 =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ASINH)) && !defined(HAVE_ASIN= HL) +#define HAVE_ASINHL 1 +extern long double asinhl(long double); +#endif + #ifndef HAVE_ATAN2F #define HAVE_ATAN2F 1 extern float atan2f(float, float); #endif =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATAN2)) && !defined(HAVE_ATAN= 2L) +#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 =20 +#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 =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATANH)) && !defined(HAVE_ATAN= HL) +#define HAVE_ATANHL 1 +extern long double atanhl (long double); +#endif + #ifndef HAVE_CEILF #define HAVE_CEILF 1 extern float ceilf(float); #endif =20 +#if !defined(HAVE_CEILL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 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 =20 +#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 =20 +#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 =20 +#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 =20 -#ifndef HAVE_FREXPL +#if !defined(HAVE_FREXPL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 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 =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_HYPOT)) && !defined(HAVE_HYPO= TL) +#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 =20 +#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 =20 -#ifndef HAVE_SCALBNL +#if !defined(HAVE_SCALBNL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 16) && (__LDB= L_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 =20 +#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 =20 +#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 =20 +#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 =20 +#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 =20 +#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 =20 +#if !defined(HAVE_TRUNCL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 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 =20 +#if !defined(HAVE_NEXTAFTERL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 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 =20 +#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 =20 +#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 =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ERFC)) && !defined(HAVE_ERFCL) +#define HAVE_ERFL 1 +long double erfcl(long double); +#endif + =20 =20 /* 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/intrinsic= s/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 =20 #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW #include "libgfortran.h" +#include "math_imp.h" =20 /* On a C99 system "I" (with I*I =3D -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= =2E If not, see /* Prototypes are included to silence -Wstrict-prototypes -Wmissing-prototypes. */ =20 -/* 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__ =3D=3D __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__ =3D=3D __ORDER_BIG_ENDIAN__ - uint64_t high; - uint64_t low; -#else - uint64_t low; - uint64_t high; -#endif - } words64; - - struct - { -#if __BYTE_ORDER__ =3D=3D __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__ =3D=3D __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 =3D (d); \ - (ix0) =3D u.words64.high; \ - (ix1) =3D 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 =3D (ix0); \ - u.words64.low =3D (ix1); \ - (d) =3D 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 =3D (d); \ - (v) =3D 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 =3D (d); \ - u.words64.high =3D (v); \ - (d) =3D 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 =3D (d); \ - (v) =3D u.words64.low; \ -} while (0) - -static const long double -two114 =3D 2.0769187434139310514121985316880384E+34L, /* 0x407100000000000= 0, 0 */ -twom114 =3D 4.8148248609680896326399448564623183E-35L, /* 0x3F8D0000000000= 00, 0 */ -huge =3D 1.0E+4900L, -tiny =3D 1.0E-4900L; - /* Wrapper for systems without strnlen function. */ =20 #ifndef HAVE_STRNLEN @@ -276,6 +145,25 @@ erff (float x) } #endif =20 +/* 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 =20 +#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 + =20 #ifndef HAVE_ACOSF #define HAVE_ACOSF 1 @@ -299,6 +202,21 @@ acosf (float x) } #endif =20 +#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); =20 @@ -309,6 +227,21 @@ acoshf (float x) } #endif =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ACOS)) && !defined(HAVE_ACOSH= L) +#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 =20 +#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); =20 @@ -330,6 +278,21 @@ asinhf (float x) } #endif =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ASINH)) && !defined(HAVE_ASIN= HL) +#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 =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATAN2)) && !defined(HAVE_ATAN= 2L) +#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 =20 +#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); =20 @@ -362,6 +355,21 @@ atanhf (float x) } #endif =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_ATANH)) && !defined(HAVE_ATAN= HL) +#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 =20 +#if !defined(HAVE_CEILL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 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 =20 +#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 =20 +#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 =20 +#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 =20 +#if !defined(HAVE_FLOORL) +#define HAVE_FLOORL 1 +long double floorl (long double); + +#if (__SIZEOF_LONG_DOUBLE__ =3D=3D 16) && (__LDBL_IS_IEC_60559__ > 0) +#include "floorl_16.c" +#else +long double +floorl (long double x); +{ + /* Zero, possibly signed. */ + if (x =3D=3D 0) + return x; + + /* Large magnitude. */ + if (x > DBL_MAX || x < (-DBL_MAX)) + return x; + + /* Small positive values. */ + if (x >=3D 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 =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_FLOORL)) && !defined(HAVE_FMO= DL) +#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 =3D=3D 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 =20 -#if !defined(HAVE_FREXPL) +#if !defined(HAVE_FREXPL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 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]. - */ - -/* - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D - * 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. - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D - */ - -long double -frexpl (long double x, int *eptr) -{ - if (sizeof (long double) =3D=3D 16) - { - uint64_t hx, lx, ix; - - GET_LDOUBLE_WORDS64(hx,lx,x); - ix =3D 0x7fffffffffffffffULL&hx; - *eptr =3D 0; - if(ix >=3D 0x7fff000000000000ULL || ((ix|lx) =3D=3D 0)) - /* 0,inf,nan */ - return x + x; - if (ix < 0x0001000000000000ULL) - { - /* subnormal */ - x *=3D two114; - GET_LDOUBLE_MSW64(hx,x); - ix =3D hx & 0x7fffffffffffffffULL; - *eptr =3D -114; - } - *eptr +=3D (ix>>48) - 16382; - hx =3D (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 =20 #ifndef HAVE_HYPOTF @@ -567,6 +633,21 @@ hypotf (float x, float y) } #endif =20 +#if (defined(USE_LIBQUADLIB) || defined(HAVE_HYPOT)) && !defined(HAVE_HYPO= TL) +#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 =20 +#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 =20 -#if !defined(HAVE_SCALBNL) +#if !defined(HAVE_SCALBNL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 16) && (__LDB= L_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 */ -/* - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D - * 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. - * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D - */ - -long double -scalbnl (long double x, int n) -{ - if (sizeof (long double) =3D=3D 16) - { - int64_t k, hx, lx; - - GET_LDOUBLE_WORDS64(hx,lx,x); - - /* extract exponent */ - k =3D (hx >> 48) & 0x7fff; - if (k =3D=3D 0) - { - /* 0 or subnormal x */ - if ((lx | (hx & 0x7fffffffffffffffULL)) =3D=3D 0) - return x; /* +-0 */ - x *=3D two114; - GET_LDOUBLE_MSW64(hx,x); - k =3D ((hx >> 48) & 0x7fff) - 114; - } - if (k =3D=3D 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 =3D k+n does not overflow.= */ - k =3D k + n; - if (k > 0) - { - /* normal result */ - SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); - return x; - } - if (k <=3D -114) - /*underflow*/ - return tiny * copysignl (tiny, x); - k +=3D 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 =20 #ifndef HAVE_SINF @@ -696,6 +728,21 @@ sinf (float x) } #endif =20 +#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 =20 +#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 =20 +#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 =20 +#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 =20 +#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 =20 +#if !defined(HAVE_TRUNCL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 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 =20 +#if !defined(HAVE_NEXTAFTERL) && (__SIZEOF_LONG_DOUBLE__ =3D=3D 16) && (__= LDBL_IS_IEC_60559__ > 0) +#define HAVE_NEXTAFTERL 1 +long double nextafterl (long double, long double); +#include "nextafterl_16.c" +#endif =20 #ifndef HAVE_POWF #define HAVE_POWF 1 @@ -845,6 +963,21 @@ powf (float x, float y) } #endif =20 +#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 + =20 #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 =3D=3D 0) - return x; - - /* Large magnitude. */ - if (x > DBL_MAX || x < (-DBL_MAX)) - return x; - - /* Small positive values. */ - if (x >=3D 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 =3D=3D 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 =20 diff --git a/libgfortran/intrinsics/ceill_16.c b/libgfortran/intrinsics/cei= ll_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]. + */ + +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +/* + * 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 =3D ((i0 >> 48) & 0x7fff) - 0x3fff; + if (j0 < 48) + { + if (j0 < 0) + { + /* return 0*sign(x) if |x|<1 */ + if (i0 < 0) + { + i0 =3D 0x8000000000000000ULL; + i1 =3D 0; + } + else if ((i0|i1) !=3D 0) + { + i0 =3D 0x3fff000000000000ULL; + i1 =3D 0; + } + } + else + { + i =3D (0x0000ffffffffffffULL) >> j0; + if (((i0 & i) | i1) =3D=3D 0) + /* x is integral */ + return x; + if (i0 > 0) + i0 +=3D (0x0001000000000000LL) >> j0; + i0 &=3D (~i); + i1 =3D 0; + } + } + else if (j0 > 111) + { + if (j0 =3D=3D 0x4000) + /* inf or NaN */ + return x+x; + else + /* x is integral */ + return x; + } + else + { + i =3D -1ULL >> (j0 - 48); + if ((i1 & i) =3D=3D 0) + /* x is integral */ + return x; + if (i0 > 0) + { + if (j0 =3D=3D 48) + i0 +=3D 1; + else + { + j =3D i1 + (1LL << (112 - j0)); + if ((int64_t) j < i1) + /* got a carry */ + i0 +=3D 1; + i1 =3D j; + } + } + i1 &=3D (~i); + } + SET_LDOUBLE_WORDS64 (x, i0, i1); + return x; +} diff --git a/libgfortran/intrinsics/floorl_16.c b/libgfortran/intrinsics/fl= oorl_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]. + */ + +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +/* + * 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 =3D ((i0 >> 48) & 0x7fff) - 0x3fff; + if (j0 < 48) + { + if (j0 < 0) + { + /* return 0*sign(x) if |x|<1 */ + if (i0 >=3D 0) + i0 =3D i1 =3D 0; + else if (((i0 & 0x7fffffffffffffffLL) | i1) !=3D 0) + { + i0 =3D 0xbfff000000000000ULL; + i1 =3D 0; + } + } + else + { + i =3D (0x0000ffffffffffffULL) >> j0; + if (((i0 & i) | i1) =3D=3D 0) + /* x is integral */ + return x; + if (i0 < 0) + i0 +=3D (0x0001000000000000LL) >> j0; + i0 &=3D (~i); + i1 =3D 0; + } + } + else if (j0 > 111) + { + if (j0 =3D=3D 0x4000) + /* inf or NaN */ + return x + x; + else + /* x is integral */ + return x; + } + else + { + i =3D -1ULL >> (j0 - 48); + if ((i1 & i) =3D=3D 0) + /* x is integral */ + return x; + if (i0 < 0) + { + if (j0 =3D=3D 48) + i0 +=3D 1; + else + { + j =3D i1 + (1LL << (112 - j0)); + if ((int64_t) j < i1) + /* got a carry */ + i0 +=3D 1; + i1 =3D j; + } + } + i1 &=3D (~i); + } + SET_LDOUBLE_WORDS64 (x, i0, i1); + return x; +} diff --git a/libgfortran/intrinsics/frexpl_16.c b/libgfortran/intrinsics/fr= expl_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]. + */ + +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +long double +frexpl (long double x, int *eptr) +{ + uint64_t hx, lx, ix; + + GET_LDOUBLE_WORDS64 (hx, lx, x); + ix =3D 0x7fffffffffffffffULL&hx; + *eptr =3D 0; + if (ix >=3D 0x7fff000000000000ULL || ((ix | lx) =3D=3D 0)) + /* 0,inf,nan */ + return x + x; + if (ix < 0x0001000000000000ULL) + { + /* subnormal */ + x *=3D two114; + GET_LDOUBLE_MSW64 (hx, x); + ix =3D hx & 0x7fffffffffffffffULL; + *eptr =3D -114; + } + *eptr +=3D (ix >> 48) - 16382; + hx =3D (hx & 0x8000ffffffffffffULL) | 0x3ffe000000000000ULL; + SET_LDOUBLE_MSW64 (x, hx); + return x; +} diff --git a/libgfortran/intrinsics/math_imp.h b/libgfortran/intrinsics/mat= h_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__ =3D=3D 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__ =3D=3D __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__ =3D=3D __ORDER_BIG_ENDIAN__ + uint64_t high; + uint64_t low; +#else + uint64_t low; + uint64_t high; +#endif + } words64; + + struct + { +#if __BYTE_ORDER__ =3D=3D __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__ =3D=3D __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 =3D (d); \ + (ix0) =3D u.words64.high; \ + (ix1) =3D 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 =3D (ix0); \ + u.words64.low =3D (ix1); \ + (d) =3D 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 =3D (d); \ + (v) =3D 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 =3D (d); \ + u.words64.high =3D (v); \ + (d) =3D 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 =3D (d); \ + (v) =3D u.words64.low; \ +} while (0) + +static const long double __attribute__ ((unused)) +two114 =3D 2.0769187434139310514121985316880384E+34L, /* 0x407100000000000= 0, 0 */ +twom114 =3D 4.8148248609680896326399448564623183E-35L, /* 0x3F8D0000000000= 00, 0 */ +huge =3D 1.0E+4900L, +tiny =3D 1.0E-4900L; +#endif + +#ifndef math_opt_barrier +# define math_opt_barrier(x) \ +({ __typeof (x) __x =3D (x); __asm ("" : "+m" (__x)); __x; }) +# define math_force_eval(x) \ +({ __typeof (x) __x =3D (x); __asm __volatile__ ("" : : "m" (__x)); }) +#endif + +#endif /* LIBGFORTRAN_MATH_IMP_H */ diff --git a/libgfortran/intrinsics/nextafterl_16.c b/libgfortran/intrinsic= s/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]. + */ + +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +/* 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 =3D hx & 0x7fffffffffffffffLL; /* |x| */ + iy =3D hy & 0x7fffffffffffffffLL; /* |y| */ + + if (((ix >=3D 0x7fff000000000000LL) && ((ix - 0x7fff000000000000LL) | lx= ) !=3D 0) + || ((iy >=3D 0x7fff000000000000LL) && ((iy - 0x7fff000000000000LL) | ly)= !=3D 0)) + /* x or y is nan */ + return x + y; + if (x =3D=3D y) + /* x=3Dy, return y */ + return y; + if ((ix | lx) =3D=3D 0) + { + /* x =3D=3D 0 */ + long double u; + + /* return +-minsubnormal */ + SET_LDOUBLE_WORDS64 (x, hy & 0x8000000000000000ULL, 1); + u =3D math_opt_barrier (x); + u =3D u * u; + /* raise underflow flag */ + math_force_eval (u); + return x; + } + if (hx >=3D 0) + { + /* x > 0 */ + if (hx > hy || ((hx =3D=3D hy) && (lx > ly))) + { + /* x > y, x -=3D ulp */ + if (lx =3D=3D 0) + hx--; + lx--; + } + else + { + /* x < y, x +=3D ulp */ + lx++; + if (lx =3D=3D 0) + hx++; + } + } + else + { + /* x < 0 */ + if (hy >=3D 0 || hx > hy || ((hx =3D=3D hy) && (lx > ly))) + { + /* x < y, x -=3D ulp */ + if (lx =3D=3D 0) + hx--; + lx--; + } + else + { + /* x > y, x +=3D ulp */ + lx++; + if (lx =3D=3D 0) + hx++; + } + } + hy =3D hx & 0x7fff000000000000LL; + if (hy =3D=3D 0x7fff000000000000LL) + { + /* overflow */ + long double u =3D x + x; + math_force_eval (u); + errno =3D ERANGE; + } + if (hy =3D=3D 0) + { + /* underflow */ + long double u =3D x * x; + /* raise underflow flag */ + math_force_eval (u); + errno =3D ERANGE; + } + SET_LDOUBLE_WORDS64 (x, hx, lx); + return x; +} diff --git a/libgfortran/intrinsics/scalbnl_16.c b/libgfortran/intrinsics/s= calbnl_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 */ +/* + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + * 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. + * =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D + */ + +long double +scalbnl (long double x, int n) +{ + int64_t k, hx, lx; + + GET_LDOUBLE_WORDS64 (hx, lx, x); + + /* extract exponent */ + k =3D (hx >> 48) & 0x7fff; + if (k =3D=3D 0) + { + /* 0 or subnormal x */ + if ((lx | (hx & 0x7fffffffffffffffULL)) =3D=3D 0) + return x; /* +-0 */ + x *=3D two114; + GET_LDOUBLE_MSW64 (hx, x); + k =3D ((hx >> 48) & 0x7fff) - 114; + } + if (k =3D=3D 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 =3D k + n does not overflow. = */ + k =3D k + n; + if (k > 0) + { + /* normal result */ + SET_LDOUBLE_MSW64 (x, (hx&0x8000ffffffffffffULL) | (k << 48)); + return x; + } + if (k <=3D -114) + /* underflow */ + return tiny * copysignl (tiny, x); + + /* subnormal result */ + k +=3D 114; + SET_LDOUBLE_MSW64 (x, (hx & 0x8000ffffffffffffULL) | (k << 48)); + return x * twom114; +} diff --git a/libgfortran/intrinsics/truncl_16.c b/libgfortran/intrinsics/tr= uncl_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 argumen= t. + 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 =3D i0 & 0x8000000000000000ULL; + j0 =3D ((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 =3D=3D 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 =20 +/* 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> --9ziKhqtsx3XySAFl Content-Type: application/pgp-signature; name=signature.asc -----BEGIN PGP SIGNATURE----- iQIzBAABCAAdFiEEnRzl+6e9+DTrEhyEXb/Nrl8ZTfEFAmpiKiAACgkQXb/Nrl8Z TfHm7Q//eB1Tpqkp4GYEDKYT6WTBpwkmHADjmafz16Qe8aagOBVHQ89cQo07iYOc d/amiB16cnJkDGZA2Epj0UDgtTCHoSw5hquZx+SxtF8a5syDGFu7Xi3LWBkeAQig lVkMkVUtrRsMnN94XWyfW6G2a6opWkO9Rs1HolmG3XyAwzg92d5UrPZeup3p/KHH 2B3rt89ONmiRfrkEPyNkgKEx83Cn4dTi6ko0F35rTnqInuqT5CsiC7EzbxD9/Kjm 6k5tyZs2zmCIJ08qpBCUKgW7/BobY26bOntOUjjriHaRYKmavfiU4DyV0+4D+C0J Dfuj/sudOJ8EgBKHlxpzMoS6mmdGvTzuvw5zm8hDKaxnaxY84h1Huo+T7GEeqni1 7lKMe3BhTwmNttkQGZTKipExOd305lfuzznjVl2AZ0dibv7c5QsBfzCt5jnSJbc3 W9owmSQpKqA+KQPr881YuYgMpVrU4WLp0likqoVuHZzjLtLZk9qhvGX+jr3wu1mx w42T7Vm2/sUoh+2pOv9ZNwGxivk8yT4g/spHylcCO0xl7qsFyZ6MOw4TPSRjdfLq gnpGpFkymG0Jb2Fchz8WB0qQeC4civFSVyk9MzaSvTsnCupM3DARxuP3ssMQOBCo 48RSMEifCitPzDIhuyZNGkKrAKdqFocAKFFCOMfj7g4TDo5sCzE= =shYg -----END PGP SIGNATURE----- --9ziKhqtsx3XySAFl--