various fixes detected with -Wdouble-promotion

Changes with respect to v1:
- added comment in e_j1f.c to explain the use of float is enough
Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>

math/w_jnf_compat.c

@@ -52,7 +52,7 @@ __ynf (int n, float x)
 	  feraiseexcept (FE_INVALID);
 	  return __kernel_standard_f (n, x, 113);
 	}
-      else if (x == 0.0)
+      else if (x == 0.0f)
 	{
 	  /* d = -one/(x-x) */
 	  feraiseexcept (FE_DIVBYZERO);

math/w_log2f_compat.c

@@ -29,7 +29,7 @@ __log2f_compat (float x)
 {
   if (__builtin_expect (islessequal (x, 0.0f), 0) && _LIB_VERSION != _IEEE_)
     {
-      if (x == 0.0)
+      if (x == 0.0f)
 	{
 	  feraiseexcept (FE_DIVBYZERO);
 	  return __kernel_standard_f (x, x, 148); /* log2(0) */

math/w_tgammaf_compat.c

@@ -26,7 +26,7 @@ __tgammaf(float x)
 	float y = __ieee754_gammaf_r(x, NULL);
 
 	if(__glibc_unlikely (!isfinite (y) || y == 0)
-	   && (isfinite (x) || (isinf (x) && x < 0.0))
+	   && (isfinite (x) || (isinf (x) && x < 0.0f))
 	   && _LIB_VERSION != _IEEE_) {
 	  if (x == (float)0.0)
 	    /* tgammaf pole */

sysdeps/ieee754/dbl-64/s_fmaf.c

@@ -45,7 +45,7 @@ __fmaf (float x, float y, float z)
 
   /* Ensure correct sign of an exact zero result by performing the
      addition in the original rounding mode in that case.  */
-  if (temp == -z)
+  if (temp == (double) -z)
     return (float) temp + z;
 
   union ieee754_double u;

sysdeps/ieee754/flt-32/e_acosf.c

@@ -112,9 +112,9 @@ __ieee754_acosf (float x)
   if (ax < (0x7eu << 24))
     {
       if (t == 0x328885a3u)
-	return 0x1.921fb6p+0f + 0x1p-25;
+	return 0x1.921fb6p+0f + 0x1p-25f;
       if (t == 0x39826222u)
-	return 0x1.920f6ap+0f + 0x1p-25;
+	return 0x1.920f6ap+0f + 0x1p-25f;
       double x2 = xs * xs;
       r = (pi2 - xs) - (xs * x2) * poly12 (x2, C0);
     }

sysdeps/ieee754/flt-32/e_atan2f.c

@@ -73,7 +73,7 @@ cr_atan2f_tiny (float y, float x)
   static const double c = -0x1.5555555555555p-2; /* -1/3 rounded to nearest */
   double zz = z * z;
   double cz = c * z;
-  e = e / x + cz * zz;
+  e = e / dx + cz * zz;
   uint64_t t = asuint64 (z);
   if ((t & UINT64_C(0xfffffff)) == 0) /* boundary case */
     {
@@ -139,7 +139,7 @@ __ieee754_atan2f (float y, float x)
 	  if (ux >> 31)
 	    return pi * sgn[uy >> 31];
 	  else
-	    return 0.0f * sgn[uy >> 31];
+	    return 0.0 * sgn[uy >> 31];
 	}
       if (yinf)
 	return pi2 * sgn[uy >> 31];
@@ -155,7 +155,7 @@ __ieee754_atan2f (float y, float x)
 	    return off[i];
 	}
       if (!(ux >> 31))
-	return 0.0f * sgn[uy >> 31];
+	return 0.0 * sgn[uy >> 31];
     }
   uint32_t gt = ay > ax;
   uint32_t i = (uy >> 31) * 4 + (ux >> 31) * 2 + gt;

sysdeps/ieee754/flt-32/e_coshf.c

@@ -77,7 +77,7 @@ __ieee754_coshf (float x)
   double rm = sm * (te - h * to);
   double r = rp + rm;
   float ub = r;
-  double lb = r - 1.45e-10 * r;
+  float lb = r - 1.45e-10 * r;
   if (__glibc_unlikely (ub != lb))
     {
       const double iln2h = 0x1.7154765p+5;

sysdeps/ieee754/flt-32/e_exp10f.c

@@ -163,7 +163,7 @@ __exp10f (float x)
       if (x < -0x2d.278d4p0f)
         return __math_uflowf (0);
 #if WANT_ERRNO_UFLOW
-      if (x < -0x2c.da7cfp0)
+      if (x < -0x2c.da7cfp0f)
         return __math_may_uflowf (0);
 #endif
       /* the smallest value such that 10^x >= 2^-126 (normal range)

sysdeps/ieee754/flt-32/e_j0f.c

@@ -203,7 +203,7 @@ j0f_asympt (float x)
     return 0xf.47039p-28f;
   double y = 1.0 / (double) x;
   double y2 = y * y;
-  double beta0 = 1.0f + y2 * (-0x1p-4f + 0x1.a8p-4 * y2);
+  double beta0 = 1.0 + y2 * (-0x1p-4 + 0x1.a8p-4 * y2);
   double alpha0 = y * (0x2p-4 - 0x1.0aaaaap-4 * y2);
   double h;
   int n;
@@ -484,7 +484,7 @@ y0f_asympt (float x)
     return 0x1.a48974p-40f;
   double y = 1.0 / (double) x;
   double y2 = y * y;
-  double beta0 = 1.0f + y2 * (-0x1p-4f + 0x1.a8p-4 * y2);
+  double beta0 = 1.0 + y2 * (-0x1p-4 + 0x1.a8p-4 * y2);
   double alpha0 = y * (0x2p-4 - 0x1.0aaaaap-4 * y2);
   double h;
   int n;
@@ -529,7 +529,7 @@ y0f_near_root (float x, float z)
   /* For degree 0 use a degree-5 polynomial, where the coefficients of
      degree 4 and 5 are hard-coded.  */
   float p6 = (index > 0) ? p[6]
-    : p[6] + y * (-0x3.a46c9p-4 + y * 0x3.735478p-4);
+    : p[6] + y * (-0x3.a46c9p-4f + y * 0x3.735478p-4f);
   float res = p[3] + y * (p[4] + y * (p[5] + y * p6));
   return res;
 }

sysdeps/ieee754/flt-32/e_j1f.c

@@ -207,7 +207,7 @@ j1f_asympt (float x)
     }
   double y = 1.0 / (double) x;
   double y2 = y * y;
-  double beta1 = 1.0f + y2 * (0x3p-4 - 0x3.18p-4 * y2);
+  double beta1 = 1.0 + y2 * (0x3p-4 - 0x3.18p-4 * y2);
   double alpha1;
   alpha1 = y * (-0x6p-4 + y2 * (0x2.ap-4 - 0x5.ef33333333334p-4 * y2));
   double h;
@@ -494,7 +494,7 @@ y1f_asympt (float x)
   float cst = 0xc.c422ap-4; /* sqrt(2/pi) rounded to nearest  */
   double y = 1.0 / (double) x;
   double y2 = y * y;
-  double beta1 = 1.0f + y2 * (0x3p-4 - 0x3.18p-4 * y2);
+  double beta1 = 1.0 + y2 * (0x3p-4 - 0x3.18p-4 * y2);
   double alpha1;
   alpha1 = y * (-0x6p-4 + y2 * (0x2.ap-4 - 0x5.ef33333333334p-4 * y2));
   double h;
@@ -537,10 +537,12 @@ y1f_near_root (float x, float z)
     return z;
   float xmid = p[1];
   float y = x - xmid, p6;
+  /* Using float rather than double in the constants below is enough to get
+     the desired accuracy (at most 9 ulps).  */
   if (index == 0)
-    p6 = p[6] + y * (-0x1.28043p-8 + y * 0x2.50e83p-8);
+    p6 = p[6] + y * (-0x1.28043p-8f + y * 0x2.50e83p-8f);
   else if (index == 1)
-    p6 = p[6] + y * -0xf.ff6b8p-12;
+    p6 = p[6] + y * -0xf.ff6b8p-12f;
   else
     p6 = p[6];
   return p[3] + y * (p[4] + y * (p[5] + y * p6));

sysdeps/ieee754/flt-32/e_jnf.c

@@ -61,7 +61,8 @@ __ieee754_jnf(int n, float x)
 	    b = __ieee754_j1f(x);
 	    for(i=1;i<n;i++){
 		temp = b;
-		b = b*((double)(i+i)/x) - a; /* avoid underflow */
+                /* the following computation should be done in double precision */
+		b = (double)b*((double)(i+i)/(double)x) - (double)a; /* avoid underflow */
 		a = temp;
 	    }
 	} else {
@@ -219,7 +220,8 @@ __ieee754_ynf(int n, float x)
 	GET_FLOAT_WORD(ib,b);
 	for(i=1;i<n&&ib!=0xff800000;i++){
 	    temp = b;
-	    b = ((double)(i+i)/x)*b - a;
+            /* the following computation should be done in double precision */
+	    b = ((double)(i+i)/(double)x)*(double)b - (double)a;
 	    GET_FLOAT_WORD(ib,b);
 	    a = temp;
 	}

sysdeps/ieee754/flt-32/e_lgammaf_r.c

@@ -181,6 +181,7 @@ __ieee754_lgammaf_r (float x, int *signgamp)
     *signgamp = 1 - ((((int) fx) & 1) << 1);
 
   double z = ax, f;
+  double s = x;
   if (__glibc_unlikely (ax < 0x1.52p-1f))
     {
       static const double rn[] =
@@ -198,7 +199,6 @@ __ieee754_lgammaf_r (float x, int *signgamp)
 	  -0x1.7dd25af0b83d4p+0, -0x1.36bf1880125fcp+0,
 	  -0x1.1379fc8023d9cp+0, -0x1.03712e41525d2p+0
 	};
-      double s = x;
       f = (c0 * s) * as_r8 (s, rn) / as_r8 (s, rd) - as_ln (z);
     }
   else
@@ -208,7 +208,7 @@ __ieee754_lgammaf_r (float x, int *signgamp)
 	  if (__glibc_unlikely (x > 0x1.895f1cp+121f))
 	    return math_narrow_eval (0x1p127f * 0x1p127f);
 	  /* |x|>=2**23, must be -integer */
-	  if (__glibc_unlikely (x < 0.0f && ax > 0x1p+23))
+	  if (__glibc_unlikely (x < 0.0f && ax > 0x1p+23f))
 	    return ax / 0.0f;
 	  double lz = as_ln (z);
 	  f = (z - 0.5) * (lz - 1) + 0x1.acfe390c97d69p-2;
@@ -286,7 +286,7 @@ __ieee754_lgammaf_r (float x, int *signgamp)
 	    {
 	      if (__glibc_unlikely (t < 0x40301b93u && t > 0x402f95c2u))
 		{
-		  double h = (x + 0x1.5fb410a1bd901p+1)
+		  double h = (s + 0x1.5fb410a1bd901p+1)
 		    - 0x1.a19a96d2e6f85p-54;
 		  double h2 = h * h;
 		  double h4 = h2 * h2;
@@ -302,7 +302,7 @@ __ieee754_lgammaf_r (float x, int *signgamp)
 		}
 	      else if (__glibc_unlikely (t > 0x401ceccbu && t < 0x401d95cau))
 		{
-		  double h = (x + 0x1.3a7fc9600f86cp+1)
+		  double h = (s + 0x1.3a7fc9600f86cp+1)
 		    + 0x1.55f64f98af8dp-55;
 		  double h2 = h * h;
 		  double h4 = h2 * h2;
@@ -319,7 +319,7 @@ __ieee754_lgammaf_r (float x, int *signgamp)
 		}
 	      else if (__glibc_unlikely (t > 0x40492009u && t < 0x404940efu))
 		{
-		  double h = (x + 0x1.9260dbc9e59afp+1)
+		  double h = (s + 0x1.9260dbc9e59afp+1)
 		    + 0x1.f717cd335a7b3p-53;
 		  double h2 = h * h;
 		  double h4 = h2 * h2;

sysdeps/ieee754/flt-32/e_powf.c

@@ -203,7 +203,7 @@ __powf (float x, float y)
 	}
     }
   double_t logx = log2_inline (ix);
-  double_t ylogx = y * logx; /* Note: cannot overflow, y is single prec.  */
+  double_t ylogx = (double) y * logx; /* Note: cannot overflow, y is single prec.  */
   if (__glibc_unlikely ((asuint64 (ylogx) >> 47 & 0xffff)
 			>= asuint64 (126.0 * POWF_SCALE) >> 47))
     {

sysdeps/ieee754/flt-32/e_remainderf.c

@@ -46,7 +46,7 @@ __ieee754_remainderf(float x, float p)
 	    x -= p;
 	  /* Make sure x is not -0. This can occur only when x = p
 	     and rounding direction is towards negative infinity. */
-	  else if (x == 0.0)
+	  else if (x == 0.0f)
 	    x = 0.0;
 	}
     }
@@ -63,7 +63,7 @@ __ieee754_remainderf(float x, float p)
 	  x -= p;
 	  if (x >= p_half)
 	    x -= p;
-	  else if (x == 0.0)
+	  else if (x == 0.0f)
 	    x = 0.0;
 	}
     }

sysdeps/ieee754/flt-32/e_sinhf.c

@@ -91,7 +91,7 @@ __ieee754_sinhf (float x)
   double rm = sm * (te - h * to);
   double r = rp - rm;
   float ub = r;
-  double lb = r - 1.52e-10 * r;
+  float lb = r - 1.52e-10 * r;
   if (__glibc_unlikely (ub != lb))
     {
       const double iln2h = 0x1.7154765p+5;

sysdeps/ieee754/flt-32/s_asinpif.c

@@ -65,7 +65,7 @@ __asinpif (float x)
       c0 += c2 * z4;
       c4 += c6 * z4;
       c0 += c4 * (z4 * z4);
-      if (__glibc_unlikely (ax != 0.0 && ax <= 0x1.921fb4p-126f))
+      if (__glibc_unlikely (ax != 0.0f && ax <= 0x1.921fb4p-126f))
 	__set_errno (ERANGE);
       return z * c0;
     }

sysdeps/ieee754/flt-32/s_atanpif.c

@@ -48,7 +48,7 @@ __atanpif (float x)
 	}
       /* Warning: 0x1.45f306p-2f / x underflows for |x| >= 0x1.45f306p+124 */
       if (fabsf (x) >= 0x1.45f306p+124f)
-	return f - copysign (0x1p-26f, x);
+	return f - copysignf (0x1p-26f, x);
       else
 	return f - 0x1.45f306p-2f / x;
     }
@@ -63,7 +63,7 @@ __atanpif (float x)
 	    __set_errno (ERANGE);
 	  return rsx;
 	}
-      return sx - (0x1.5555555555555p-2 * sx) * (x * x);
+      return sx - (0x1.5555555555555p-2 * sx) * (z * z);
     }
   uint32_t ax = t & (~0u >> 1);
   if (__glibc_unlikely (ax == 0x3fa267ddu))

sysdeps/ieee754/flt-32/s_erfcf.c

@@ -131,7 +131,7 @@ __erfcf (float xf)
 	  0x1.20dd750429b6dp+0, -0x1.812746b03610bp-2, 0x1.ce2f218831d2fp-4,
 	  -0x1.b82c609607dcbp-6, 0x1.553af09b8008ep-8
 	};
-      double f0 = xf
+      double f0 = (double) xf
 	    * (c[0] + x2 * (c[1] + x2 * (c[2] + x2 * (c[3] + x2 * (c[4])))));
       return 1.0 - f0;
     }

sysdeps/ieee754/flt-32/s_log10p1f.c

@@ -155,7 +155,7 @@ __log10p1f (float x)
 	      0x1.63c62378fa3dbp-3, 0x1.fca4139a42374p-4
 	    };
 	  float ret = z * ((c[0] + z2 * c[1]) + z4 * (c[2] + z2 * c[3]));
-	  if (x != 0.0f && ret == 0.0)
+	  if (x != 0.0f && ret == 0.0f)
 	    __set_errno (ERANGE);
 	  return ret;
 	}

sysdeps/ieee754/flt-32/s_log1pf.c

@@ -154,21 +154,21 @@ __log1pf (float x)
 	  double tr = rl + Lh;
 	  if (__glibc_unlikely ((asuint64 (tr) & 0xfffffffll) == 0))
 	    {
-	      if (x == -0x1.247ab0p-6)
+	      if (x == -0x1.247ab0p-6f)
 		return -0x1.271f0ep-6f - 0x1p-31f;
-	      if (x == -0x1.3a415ep-5)
+	      if (x == -0x1.3a415ep-5f)
 		return -0x1.407112p-5f + 0x1p-30f;
-	      if (x == 0x1.fb035ap-2)
+	      if (x == 0x1.fb035ap-2f)
 		return 0x1.9bddc2p-2f + 0x1p-27f;
 	      tr += 64 * (rl + (Lh - tr));
 	    }
 	  else if (rl + (Lh - tr) == 0.0)
 	    {
-	      if (x == 0x1.b7fd86p-4)
+	      if (x == 0x1.b7fd86p-4f)
 		return 0x1.a1ece2p-4f + 0x1p-29f;
-	      if (x == -0x1.3a415ep-5)
+	      if (x == -0x1.3a415ep-5f)
 		return -0x1.407112p-5f + 0x1p-30f;
-	      if (x == 0x1.43c7e2p-6)
+	      if (x == 0x1.43c7e2p-6f)
 		return 0x1.409f80p-6f + 0x1p-31f;
 	    }
 	  ub = tr;

sysdeps/ieee754/flt-32/s_remquof.c

@@ -85,7 +85,7 @@ __remquof (float x, float y, int *quo)
     }
   else
     {
-      float y_half = 0.5 * y;
+      float y_half = 0.5f * y;
       if (x > y_half)
 	{
 	  x -= y;