Update copyright dates with scripts/update-copyrights.

[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128ibm / e_expl.c

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c
index 1b994cd7a92645325d63851db59d6f40b6a514a2..69faad2afdac39cfce6aa2867cf0c0a144181ba7 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/e_expl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_expl.c
@@ -1,5 +1,5 @@
 /* Quad-precision floating point e^x.
-   Copyright (C) 1999-2013 Free Software Foundation, Inc.
+   Copyright (C) 1999-2019 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
    Partly based on double-precision code
@@ -65,7 +65,10 @@
 #include <fenv.h>
 #include <inttypes.h>
 #include <math_private.h>
-#include <sysdeps/ieee754/ldbl-128/t_expl.h>
+#include <fenv_private.h>
+
+
+#include "t_expl.h"
 
 static const long double C[] = {
 /* Smallest integer x for which e^x overflows.  */
@@ -131,81 +134,80 @@ static const long double C[] = {
  1.98412698413981650382436541785404286E-04L,
 };
 
+/* Avoid local PLT entry use from (int) roundl (...) being converted
+   to a call to lroundl in the case of 32-bit long and roundl not
+   inlined.  */
+long int lroundl (long double) asm ("__lroundl");
+
 long double
 __ieee754_expl (long double x)
 {
+  long double result, x22;
+  union ibm_extended_long_double ex2_u, scale_u;
+  int unsafe;
+
   /* Check for usual case.  */
   if (isless (x, himark) && isgreater (x, lomark))
     {
-      int tval1, tval2, unsafe, n_i, exponent2;
-      long double x22, n, result, xl;
-      union ibm_extended_long_double ex2_u, scale_u;
-      fenv_t oldenv;
-
-      feholdexcept (&oldenv);
-#ifdef FE_TONEAREST
-      fesetround (FE_TONEAREST);
-#endif
+      int tval1, tval2, n_i, exponent2;
+      long double n, xl;
 
-      n = __roundl (x*M_1_LN2);
+      SET_RESTORE_ROUND (FE_TONEAREST);
+
+      n = roundl (x*M_1_LN2);
       x = x-n*M_LN2_0;
       xl = n*M_LN2_1;
 
-      tval1 = __roundl (x*TWO8);
+      tval1 = roundl (x*TWO8);
       x -= __expl_table[T_EXPL_ARG1+2*tval1];
       xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
 
-      tval2 = __roundl (x*TWO15);
+      tval2 = roundl (x*TWO15);
       x -= __expl_table[T_EXPL_ARG2+2*tval2];
       xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
 
       x = x + xl;
 
       /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]).  */
-      ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
-               * __expl_table[T_EXPL_RES2 + tval2];
+      ex2_u.ld = (__expl_table[T_EXPL_RES1 + tval1]
+                 * __expl_table[T_EXPL_RES2 + tval2]);
       n_i = (int)n;
       /* 'unsafe' is 1 iff n_1 != 0.  */
       unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1;
-      ex2_u.ieee.exponent += n_i >> unsafe;
+      ex2_u.d[0].ieee.exponent += n_i >> unsafe;
       /* Fortunately, there are no subnormal lowpart doubles in
         __expl_table, only normal values and zeros.
         But after scaling it can be subnormal.  */
-      exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe);
-      if (ex2_u.ieee.exponent2 == 0)
-       /* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */;
+      exponent2 = ex2_u.d[1].ieee.exponent + (n_i >> unsafe);
+      if (ex2_u.d[1].ieee.exponent == 0)
+       /* assert ((ex2_u.d[1].ieee.mantissa0|ex2_u.d[1].ieee.mantissa1) == 0) */;
       else if (exponent2 > 0)
-       ex2_u.ieee.exponent2 = exponent2;
+       ex2_u.d[1].ieee.exponent = exponent2;
       else if (exponent2 <= -54)
        {
-         ex2_u.ieee.exponent2 = 0;
-         ex2_u.ieee.mantissa2 = 0;
-         ex2_u.ieee.mantissa3 = 0;
+         ex2_u.d[1].ieee.exponent = 0;
+         ex2_u.d[1].ieee.mantissa0 = 0;
+         ex2_u.d[1].ieee.mantissa1 = 0;
        }
       else
        {
          static const double
            two54 = 1.80143985094819840000e+16, /* 4350000000000000 */
            twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */
-         ex2_u.dd[1] *= two54;
-         ex2_u.ieee.exponent2 += n_i >> unsafe;
-         ex2_u.dd[1] *= twom54;
+         ex2_u.d[1].d *= two54;
+         ex2_u.d[1].ieee.exponent += n_i >> unsafe;
+         ex2_u.d[1].d *= twom54;
        }
 
       /* Compute scale = 2^n_1.  */
-      scale_u.d = 1.0L;
-      scale_u.ieee.exponent += n_i - (n_i >> unsafe);
+      scale_u.ld = 1.0L;
+      scale_u.d[0].ieee.exponent += n_i - (n_i >> unsafe);
 
       /* Approximate e^x2 - 1, using a seventh-degree polynomial,
         with maximum error in [-2^-16-2^-53,2^-16+2^-53]
         less than 4.8e-39.  */
       x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
 
-      /* Return result.  */
-      fesetenv (&oldenv);
-
-      result = x22 * ex2_u.d + ex2_u.d;
-
       /* Now we can test whether the result is ultimate or if we are unsure.
         In the later case we should probably call a mpn based routine to give
         the ultimate result.
@@ -235,15 +237,11 @@ __ieee754_expl (long double x)
            return __ieee754_expl_proc2 (origx);
          }
        */
-      if (!unsafe)
-       return result;
-      else
-       return result * scale_u.d;
     }
   /* Exceptional cases:  */
   else if (isless (x, himark))
     {
-      if (__isinfl (x))
+      if (isinf (x))
        /* e^-inf == 0, with no error.  */
        return 0;
       else
@@ -253,5 +251,10 @@ __ieee754_expl (long double x)
   else
     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
     return TWO1023*x;
+
+  result = x22 * ex2_u.ld + ex2_u.ld;
+  if (!unsafe)
+    return result;
+  return result * scale_u.ld;
 }
 strong_alias (__ieee754_expl, __expl_finite)