https://gcc.gnu.org/g:7c89181457dafe48203c0ce31a429c3e887dd378
commit r17-2676-g7c89181457dafe48203c0ce31a429c3e887dd378 Author: John David Anglin <[email protected]> Date: Thu Jul 23 22:06:50 2026 -0400 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: --- libgfortran/c99_protos.h | 122 +++++- libgfortran/config.h.in | 6 + libgfortran/configure | 96 +++++ libgfortran/configure.ac | 2 + libgfortran/intrinsics/c99_functions.c | 674 +++++++++++++++++++-------------- libgfortran/intrinsics/ceill_16.c | 91 +++++ libgfortran/intrinsics/floorl_16.c | 88 +++++ libgfortran/intrinsics/frexpl_16.c | 39 ++ libgfortran/intrinsics/math_imp.h | 170 +++++++++ libgfortran/intrinsics/nextafterl_16.c | 108 ++++++ libgfortran/intrinsics/scalbnl_16.c | 61 +++ libgfortran/intrinsics/truncl_16.c | 50 +++ libgfortran/libgfortran.h | 8 + 13 files changed, 1222 insertions(+), 293 deletions(-) diff --git a/libgfortran/c99_protos.h b/libgfortran/c99_protos.h index 952b251c31ab..a233991b60dc 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/config.h.in b/libgfortran/config.h.in index 240a1762de32..1d4bbfbb4a9e 100644 --- a/libgfortran/config.h.in +++ b/libgfortran/config.h.in @@ -435,6 +435,9 @@ /* Define to 1 if you have the `erff128' function. */ #undef HAVE_ERFF128 +/* Define to 1 if you have the `erfl' function. */ +#undef HAVE_ERFL + /* Define to 1 if you have the `exp' function. */ #undef HAVE_EXP @@ -810,6 +813,9 @@ /* Define to 1 if you have the `powf128' function. */ #undef HAVE_POWF128 +/* Define to 1 if you have the `powl' function. */ +#undef HAVE_POWL + /* Have PTHREAD_PRIO_INHERIT. */ #undef HAVE_PTHREAD_PRIO_INHERIT diff --git a/libgfortran/configure b/libgfortran/configure index 7820988523e3..6d31fa11be2a 100755 --- a/libgfortran/configure +++ b/libgfortran/configure @@ -22860,6 +22860,54 @@ _ACEOF + { $as_echo "$as_me:${as_lineno-$LINENO}: checking for powl" >&5 +$as_echo_n "checking for powl... " >&6; } +if ${gcc_cv_math_func_powl+:} false; then : + $as_echo_n "(cached) " >&6 +else + if test x$gcc_no_link = xyes; then + as_fn_error $? "Link tests are not allowed after GCC_NO_EXECUTABLES." "$LINENO" 5 +fi +cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + +#ifdef HAVE_COMPLEX_H +#include <complex.h> +#endif +#ifdef HAVE_MATH_H +#include <math.h> +#endif + +int (*ptr)() = (int (*)())powl; + +int +main () +{ + return 0; +} + +_ACEOF +if ac_fn_c_try_link "$LINENO"; then : + gcc_cv_math_func_powl=yes +else + gcc_cv_math_func_powl=no +fi +rm -f core conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +fi +{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $gcc_cv_math_func_powl" >&5 +$as_echo "$gcc_cv_math_func_powl" >&6; } + if test $gcc_cv_math_func_powl = yes; then + +cat >>confdefs.h <<_ACEOF +#define HAVE_POWL 1 +_ACEOF + + fi + + + + { $as_echo "$as_me:${as_lineno-$LINENO}: checking for cpowf" >&5 $as_echo_n "checking for cpowf... " >&6; } if ${gcc_cv_math_func_cpowf+:} false; then : @@ -25260,6 +25308,54 @@ _ACEOF + { $as_echo "$as_me:${as_lineno-$LINENO}: checking for erfl" >&5 +$as_echo_n "checking for erfl... " >&6; } +if ${gcc_cv_math_func_erfl+:} false; then : + $as_echo_n "(cached) " >&6 +else + if test x$gcc_no_link = xyes; then + as_fn_error $? "Link tests are not allowed after GCC_NO_EXECUTABLES." "$LINENO" 5 +fi +cat confdefs.h - <<_ACEOF >conftest.$ac_ext +/* end confdefs.h. */ + +#ifdef HAVE_COMPLEX_H +#include <complex.h> +#endif +#ifdef HAVE_MATH_H +#include <math.h> +#endif + +int (*ptr)() = (int (*)())erfl; + +int +main () +{ + return 0; +} + +_ACEOF +if ac_fn_c_try_link "$LINENO"; then : + gcc_cv_math_func_erfl=yes +else + gcc_cv_math_func_erfl=no +fi +rm -f core conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +fi +{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $gcc_cv_math_func_erfl" >&5 +$as_echo "$gcc_cv_math_func_erfl" >&6; } + if test $gcc_cv_math_func_erfl = yes; then + +cat >>confdefs.h <<_ACEOF +#define HAVE_ERFL 1 +_ACEOF + + fi + + + + { $as_echo "$as_me:${as_lineno-$LINENO}: checking for erfcf" >&5 $as_echo_n "checking for erfcf... " >&6; } if ${gcc_cv_math_func_erfcf+:} false; then : diff --git a/libgfortran/configure.ac b/libgfortran/configure.ac index 9acb323469e7..c5f84a031cd3 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 28acadb5f9f0..ab3920702b95 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 000000000000..40dd1fa09a6b --- /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 000000000000..d07fee1dc077 --- /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 000000000000..7f6bcc97f6bf --- /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 000000000000..fe55648dfdec --- /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 000000000000..8334da5e3e71 --- /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 000000000000..717f4d6a148c --- /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 000000000000..5632f69bd316 --- /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 8381cbd7cd48..3613b5d68c08 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>
